基于先验知识进行混合象元信息分解

基于先验知识进行混合象元信息分解*万华伟北京师范大学遥感与GIS研究中心，资源与环境科学系，环境遥感与数字城市北京市重点实验室北京，100875摘要：混合象元一直是影响遥感图像精度的一个重要因素。本文通过使用地物先验知识，利用高分辨率图像的分类图做限制条件，计算低分辨率的图像的象元地物的百分比，然后通过线性模型做混合象元分解，然后又通过分类图提供的地物信息作为先验知识，使用约束性的最小二乘法反演，从而得到比较精确的地物的值。我们利用了MODIS43B3反照率产品作为低分辨率图像，ASTER的可见光和近红外波段作为高分辨率图像，进行了实验，然后与地表的实测相对比，得到比较令人满意的结果。关键词：混合象元，线性分解，先验知识，遥感1引言在我们获得的遥感影象上，基本上都存在着混合象元。根据估计，在我们使用的MODIS的ALBEDO产品（1*1km2）上，90%以上为混合象元。这样，我们在使用这些数据时，如果不加考虑就直接使用，无疑会大大地降低我们研究的精度，而且引起信息量的丢失。因此对混合象元的分解是提高遥感信息的精度的一个重要方式[1]。目前的混合象元处理有很多方法，当前混合象元有两个基本方法：1）模糊分类方法；2）确定一个象素的光谱响应和其中端元的光谱与其所占面积比的关系[1]。对于混合象元的分解，BorisZhukov等提出了约束和非约束的MMT算法，提出利用高分辨率图像和低分辨率数据作融合[2]，但是他们在反演时都没有充分地利用前面由分类图所增加的先验知识。李小文等指出，当地表的先验知识存在时，不加利用是可惜的[3]。在本文中，我们在对混合象元进行分解时，充分利用先验知识，在使用高分辨率作为限制条件的同时，利用引入先验知识的最小二乘分解，加入了根据分类图得到的“边界”。并且利用MODIS的反照率产品（分辨率为1km）做的实验。使用了同地区ASTER的图像前三波段（分辨率为15m）作为限制条件。通过进行混合象元的分解得到植被的反照率，将得到的数据与地表实测值以及如果直接使用混合象元的值作的对比以及误差分析，得出结论，经过引入先验知识的混合象元分解的误差比直接使用混合象元值低了很多，将误差从26.58%降低到16.89%。2算法描述2．1图像的配准算法的流程图及图像的表示如Fig1*973项目资助Fig1：算法的流程图及图像的表示我们首先将低分辨率图像和高分辨率图像配准，在这里，由于两幅图像的分辨率相差太多（1km和15m），所以我们选取一中间尺度图像，我们先将高分辨率图像进行向下的尺度转换，我们直接使用15×15的聚合窗来进行重新采样。这样形成一幅225m分辨率的图像，我们再用它与1km的进行配准，这样配准的精度会比原来直接配准的精度提高许多。这样就生成了同一坐标系下的同为225m分辨率的两幅图像。2．2找出LI和MLI之间的对应关系首先找到LI与MLI的对应关系，MLI是通过最近邻插值到MHI的坐标系上的、由LI生成的图像。最近邻法是最简单的一种插值方法，这种插值的过程中，DN是是不发生变化的，只与最邻近的象元有关，和其他的象元均无关系。采用搜索算法来找到它们之间的关系，此算法的依据是插值的过程中，虽然由LI上一个象元生成的MLI的象元是不规则的，但是是连续的。1．建立一个和MLI一样大小的mark文件；2．从（1，1）开始，设（1，1）为基本点，记录其value，然后搜索其8个邻域，如果值为value，则记录其行列值，将该象元加入该基本点的一组，并将mark值设为真，再继续搜索它的8个邻域象元。如果值与value不等，则停止。3．行列值增加，记录第一个mark不为真的象元为基本点，重复步骤2。直至所有的象元的mark值都为真。这样，便可以找到一系列的由LI上一个象元生成的MLI上的多个象元，我们记做NUM。2．3找出HI和MHI之间的对应关系因为MHI由HI直接聚合而成，因此，它们之间的对应关系非常简单，MHI上的1个象元对应HI上15×15＝225个象元，这些象元是规则分布的，我们取MHI上任意的象元(i,j)，低分辨率图像（含研究需要的信息）1km，记LI高分辨率图像（含地物详细的空间分布信息）15m，记HI15×15聚合窗聚合中间尺度图像225m（由高分辨率图像生成），记MHI找控制点，配准中间尺度图像225m（由低分辨率图像生成），记MLI具有相同的坐标系i,j为它们的行列号，则对应HI上的象元的四个边界点为：（i,j）：左上角：（（i-1）*15+1,(j-1)*15+1）右上角：（（i-1）*15+1,(j-1)*15+15）左下角：（（i-1）*15+15,(j-1)*15+1）右下角：（（i-1）*15+15,(j-1)*15+15）这四个边界点之间的225个象元，便对应着MHI上的一个象元。2．4找出LI和HI之间的对应关系并求出其面积百分比由于MLI与MHI已经配准到同一坐标系下，也就是说通过以上两步找到的LI与MLI和HI与MHI之间的关系，已经找出了，LI与HI之间的关系。LI上的一个象元，对应着MI上的NUM×225个象元。我们在上面找到了LI上一个象元与HI上的NUM×225个象元以后，由于，HI图像为分类图像，其中，每个象素所属的地物类别都已知，这样，我们也就知道了LI上这一象元的各类地物的象素数，然后除以总的象元数，便得到了该象元的每类地物的百分比。delta(class=i)class[i]=NUM*225，其中，class[i]代表i类的面积百分比，delta(class=i)代表class=i的象元数。2．5使用最小二乘模型进行亚象元分解目前大部分反演采用的方法是最小二乘方法，其代价函数常用的形式为：21(())Miiiisyfx其中，M为测量总数，iy为观测数据，ix为模型参数，()iifx为由参数ix得到的模型结果[4]。理论上说，当观测数据足够多时，搜索能产生使s最小的一组参数值，也就是参数的最佳估计值。然而在实际情况的，由于传统的最小二乘方法并不考虑其参数的物理边界，有可能出现无物理意义的值，因此有发展了约束化的最小二乘方法，我们在这里，利用地物分类后的所属类别所在的物理含义的边界，作为先验知识加以约束。将代价函数修改为：2211(())()MNiiijjiisyfxxx，其中，y，x为观测值和参量值；M，N为样本数和参量数；xj为地物的先验值[8]。这样，我们得到的目标的值就会比纯粹的使用混合象元的值精确的多。3实验数据收集与处理3．1实验数据描述在本文中，我们使用的低分辨率数据是MODIS的反照率产品MOD43B3，MODIS（中分辨率成像光谱仪）是美国宇航局研制大型空间遥感仪器，它在36个相互配准的光谱波段、以中等分辨率水平（0.25Km～1Km）、每1～2天观测地球表面一次，获取陆地和海洋温度、初级生产率、陆地表面覆盖、云、汽溶胶、水汽和火情等目标的图像，其1－7波段均为研究陆地和云性质的[5]。它的标准的陆地产品，如Table1所示：陆地MOD09MOD11MOD12MOD13MOD14MOD15MOD17MOD43陆地表面反射比陆地表面温度陆地覆盖植被指数火情叶面指数和部分光合作用辐射净初级生产率/光合作用BRDF/反照率Table1．MODIS的标准陆地产品我们在这里的实验数据，低分辨率数据为MODIS的反照率产品，mod43b3，其分辨率为1km，它是由16天的不同角度的数据反演得到的，这里，我们通过网上订购2001年第113天的一景图像。它的波段分布Table2所示，我们采用其0.3～5.0um宽波段的黑半球反照率。我们使用的低分辨率约束图像Aster使用2001年8月23日成象的中国河北栾城地区的一景图像。我们的地面实测数据是栾城生态试验田（11440。''，''3750。），为0.3～5um的宽波段实际反照率。我们看植被的反射率曲线Fig2，可以知道在0.3～2.5um之间，其能量大约占到0.3～5um的99％以上，在大多数的辐射能量中，对于植被的反射率，在2500nm之后是可以忽略不计的[6]。因此我们可以使用他们两个直接对比。波段序号位置和宽度（nm）12620～670841～87634567459～479545～5651230～12501628～16522105～21558910300-700700-5000300-5000Table2MOD43B3波段分布Fig1典型植被光谱3．2数据处理过程3．2．1对MODIS数据进行处理MODIS的存储格式为HDF格式，它采用分层次的串块型数据格式-HDF（HierarchicalDataFormat）。这种数据格式更适用于大数据量的快速传输、存储和提取[7]。但是目前大多数的通用遥感处理软件还不支持这种数据格式，我们先使用MODISTOOL将其按波段转出。3．2．2对aster图像进行分类对aster图像进行分类，得到15m分辨率的分类图。分类的精度对于我们混合象元分解影响比较大，所以我们应该尽量使分类达到一个高的精度。在这里，我们先通过目视确定地物有哪些类，然后确定类数，使用监督分类。我们的ASTER图像，为2001年8月23日成像，这一时期，当地的地表覆盖主要有农作物，还有城镇和水体。在这里我们主要是对这三类地物进行亚象元分解，以消除城镇和水体的影响。3．2．3图像配准MOD43B3为1km尺度，而aster为15m尺度，直接配准误差会很大，所以我们先对aster进行向上的尺度转换，进行15*15的聚合，分辨率变为225m。然后从MODIS和aster图像上找到30个控制点，以聚合后的ASTER图像为基图，以MODIS图像为要纠正的图，进行图像的纠正，使得两幅图像在同一坐标系下。配准的结果如Fig2。MODIS-225m分辨率ASTER-225m分辨率（标准假彩色）Fig2.配准的结果3．2．4混合象元分解通过上面提到的搜索算法，记录modis1km和modis225m的对应关系。即一对多的关系，一个象元对应多个象元，我们记作NUM。然后又把modis225和aster15m的关系找到，这样就找到了，mod43b3产品的一个象元所对应的各个类的部分。我们通过计算delta(class=i)class[i]=NUM*225，class[i]代表i类的面积百分比我们便完成了对MODIS象元的分解。可以得到每个象元的各类所占的面积百分比。在这里我们选择植被和城镇的信息尽量分布均匀的象元。选取如下45个混合象元进行反演：序号图像值植被百分比城市百分比序号图像值植被百分比城市百分比11730.6337780.366222241910.3277780.67222221740.6088890.391111251940.3008330.69916731920.4367410.563259261770.2405560.75944441620.5858330.414167271870.1730560.82694451670.6948150.305185281920.4557040.54429661720.6941670.305833292040.1525000.84750071850.6989630.301037301900.6151110.38488981730.6133330.386667312070.6497220.35027891690.4636110.536389321970.5759480.424052101850.6882960.311704331910.5366670.463333112030.6538890.346111341800.6125930.387407121930.5431110.456889351840.6063700.393630131860.6686110.331389361930.5913890.408611141900.6894120.310588371810.4936110.506389151940.4483330.551667381