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

基于先验知识进行混合象元信息分解* 973项目资助万华伟北京师范大学遥感与GIS研究中心，资源与环境科学系，环境遥感与数字城市北京市重点实验室北京，100875摘要：混合象元一直是影响遥感图像精度的一个重要因素。本文通过使用地物先验知识，利用高分辨率图像的分类图做限制条件，计算低分辨率的图像的象元地物的百分比，然后通过线性模型做混合象元分解，然后又通过分类图提供的地物信息作为先验知识，使用约束性的最小二乘法反演，从而得到比较精确的地物的值。我们利用了MODIS43B3反照率产品作为低分辨率图像，ASTER的可见光和近红外波段作为高分辨率图像，进行了实验，然后与地表的实测相对比，得到比较令人满意的结果。关键词：混合象元，线性分解，先验知识，遥感1 引言在我们获得的遥感影象上，基本上都存在着混合象元。根据估计，在我们使用的MODIS的ALBEDO产品（1*1 km2）上，90%以上为混合象元。这样，我们在使用这些数据时，如果不加考虑就直接使用，无疑会大大地降低我们研究的精度，而且引起信息量的丢失。因此对混合象元的分解是提高遥感信息的精度的一个重要方式1。目前的混合象元处理有很多方法，当前混合象元有两个基本方法：1）模糊分类方法；2）确定一个象素的光谱响应和其中端元的光谱与其所占面积比的关系1。对于混合象元的分解，Boris Zhukov 等提出了约束和非约束的MMT算法，提出利用高分辨率图像和低分辨率数据作融合2，但是他们在反演时都没有充分地利用前面由分类图所增加的先验知识。李小文等指出，当地表的先验知识存在时，不加利用是可惜的3。在本文中，我们在对混合象元进行分解时，充分利用先验知识，在使用高分辨率作为限制条件的同时，利用引入先验知识的最小二乘分解，加入了根据分类图得到的“边界”。并且利用MODIS的反照率产品（分辨率为1km）做的实验。使用了同地区ASTER的图像前三波段（分辨率为15m）作为限制条件。通过进行混合象元的分解得到植被的反照率，将得到的数据与地表实测值以及如果直接使用混合象元的值作的对比以及误差分析，得出结论，经过引入先验知识的混合象元分解的误差比直接使用混合象元值低了很多，将误差从26.58%降低到16.89%。2 算法描述21 图像的配准算法的流程图及图像的表示如Fig1低分辨率图像（含研究需要的信息）1km，记LI高分辨率图像（含地物详细的空间分布信息）15m，记HI1515聚合窗聚合中间尺度图像225m（由高分辨率图像生成），记MHI找控制点，配准中间尺度图像225m（由低分辨率图像生成），记MLI具有相同的坐标系Fig1：算法的流程图及图像的表示我们首先将低分辨率图像和高分辨率图像配准，在这里，由于两幅图像的分辨率相差太多（1km和15m），所以我们选取一中间尺度图像，我们先将高分辨率图像进行向下的尺度转换，我们直接使用1515的聚合窗来进行重新采样。这样形成一幅225m分辨率的图像，我们再用它与1km的进行配准，这样配准的精度会比原来直接配准的精度提高许多。这样就生成了同一坐标系下的同为225m分辨率的两幅图像。22 找出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。23 找出HI和MHI之间的对应关系因为MHI由HI直接聚合而成，因此，它们之间的对应关系非常简单，MHI上的1个象元对应HI上1515225个象元，这些象元是规则分布的，我们取MHI上任意的象元(i,j)，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上的一个象元。24 找出LI和HI之间的对应关系并求出其面积百分比由于MLI与MHI已经配准到同一坐标系下，也就是说通过以上两步找到的LI与MLI和HI与MHI之间的关系，已经找出了，LI与HI之间的关系。LI上的一个象元，对应着MI上的NUM225个象元。我们在上面找到了LI上一个象元与HI上的NUM225个象元以后，由于，HI图像为分类图像，其中，每个象素所属的地物类别都已知，这样，我们也就知道了LI上这一象元的各类地物的象素数，然后除以总的象元数，便得到了该象元的每类地物的百分比。，其中，classi代表i类的面积百分比，代表class=i的象元数。25 使用最小二乘模型进行亚象元分解目前大部分反演采用的方法是最小二乘方法，其代价函数常用的形式为：其中，M为测量总数，为观测数据，为模型参数，为由参数得到的模型结果4。理论上说，当观测数据足够多时，搜索能产生使s最小的一组参数值，也就是参数的最佳估计值。然而在实际情况的，由于传统的最小二乘方法并不考虑其参数的物理边界，有可能出现无物理意义的值，因此有发展了约束化的最小二乘方法，我们在这里，利用地物分类后的所属类别所在的物理含义的边界，作为先验知识加以约束。将代价函数修改为：，其中，y，x为观测值和参量值；M，N为样本数和参量数；为地物的先验值8。这样，我们得到的目标的值就会比纯粹的使用混合象元的值精确的多。3 实验数据收集与处理31 实验数据描述在本文中，我们使用的低分辨率数据是MODIS的反照率产品MOD43B3，MODIS（中分辨率成像光谱仪）是美国宇航局研制大型空间遥感仪器，它在36个相互配准的光谱波段、以中等分辨率水平（0.25Km1Km）、每12天观测地球表面一次，获取陆地和海洋温度、初级生产率、陆地表面覆盖、云、汽溶胶、水汽和火情等目标的图像，其17波段均为研究陆地和云性质的5。它的标准的陆地产品，如Table1所示：陆地MOD09MOD11MOD12MOD13MOD14MOD15MOD17MOD43陆地表面反射比陆地表面温度陆地覆盖植被指数火情叶面指数和部分光合作用辐射净初级生产率/光合作用BRDF/反照率Table1MODIS的标准陆地产品我们在这里的实验数据，低分辨率数据为MODIS的反照率产品，mod43b3，其分辨率为1km，它是由16天的不同角度的数据反演得到的，这里，我们通过网上订购2001年第113天的一景图像。它的波段分布Table2所示，我们采用其0.35.0um宽波段的黑半球反照率。我们使用的低分辨率约束图像Aster使用2001年8月23日成象的中国河北栾城地区的一景图像。我们的地面实测数据是栾城生态试验田（，），为0.35um的宽波段实际反照率。我们看植被的反射率曲线Fig2，可以知道在0.32.5um之间，其能量大约占到0.35um的99以上，在大多数的辐射能量中，对于植被的反射率，在2500nm之后是可以忽略不计的6。因此我们可以使用他们两个直接对比。波段序号位置和宽度 （nm）12620670841876345674594795455651230125016281652210521558910300-700 700-5000300-5000Table2MOD43B3波段分布Fig1 典型植被光谱32 数据处理过程321 对MODIS数据进行处理MODIS的存储格式为HDF格式，它采用分层次的串块型数据格式-HDF（Hierarchical Data Format）。这种数据格式更适用于大数据量的快速传输、存储和提取7。但是目前大多数的通用遥感处理软件还不支持这种数据格式，我们先使用MODIS TOOL将其按波段转出。322 对aster图像进行分类对aster图像进行分类，得到15m分辨率的分类图。分类的精度对于我们混合象元分解影响比较大，所以我们应该尽量使分类达到一个高的精度。在这里，我们先通过目视确定地物有哪些类，然后确定类数，使用监督分类。我们的ASTER图像，为2001年8月23日成像，这一时期，当地的地表覆盖主要有农作物，还有城镇和水体。在这里我们主要是对这三类地物进行亚象元分解，以消除城镇和水体的影响。323 图像配准MOD43B3为1km尺度，而aster为15m尺度，直接配准误差会很大，所以我们先对aster进行向上的尺度转换，进行15*15的聚合，分辨率变为225m。然后从MODIS和aster图像上找到30个控制点，以聚合后的ASTER图像为基图，以MODIS图像为要纠正的图，进行图像的纠正，使得两幅图像在同一坐标系下。配准的结果如Fig2。MODIS-225m分辨率ASTER-225m分辨率（标准假彩色）Fig2.配准的结果324 混合象元分解通过上面提到的搜索算法，记录modis1km和modis225m的对应关系。即一对多的关系，一个象元对应多个象元，我们记作NUM。然后又把modis225和aster15m的关系找到，这样就找到了，mod43b3产品的一个象元所对应的各个类的部分。我们通过计算，classi代表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.551667381630.5937250.406275161990.5456210.454379391700.6386930.361307171950.6772220.322778401700.6335800.366420181910.6856300.314370411710.5665190.433481191910.4933330.506667421560.6274510.372549201880.6377780.362222431580.6810460.318954211930.6922220.307778441850.6319440.368056221850.5961110.403889451680.6544440.345556231820.4483330.551667324 基于先验知识的最小二乘反演在这里，为了反演的更加精确，我们把只包括植被和城镇两类的象元选择出来，根据公式使用代价函数，其中，y,x为观测值和残量值；M,N为样本数和参量数；为参量的先验值；我们在这里先不考虑样本之间的相关性，我们假定它们之间相互独立，对于结果同等重要。我们先设定其硬边界，植被的反照率在宽波段0.35um的先验值为100200之间。这样我们得到的结果，植被的反照率为0.1732。4 结果与讨论我们通过以上处理得到反演结果如Table2，我们又选取了植被占90%以上的象元的值取均值，通过这三种数据的比较，我们计算它们的误差，这里，定义误差Error=（计算得到的值-实测值）/实测值地面实测数据经纬度定位单点值Error最小二乘反演的植被反照率Error0.20840.15326.58%0.173216.89%Table2.反演结果比较这样，我们反演得到的值的误差为16.89%，而如果直接使用混合象元会造成26.58%的误差，可见我们经过反演使误差减少了很多，不过在这里我们分析之所以误差这么大的原因，大概是因为a） 配准的误差，配准的精度虽然已经很高，但是可能仍然会引起误差。b） 分类图的误差，我们使用的8月份的图像进行非监督分类，而且只是分成了水体，城镇和植被三类，虽然经过实地考察，8月份和3月份的地表覆盖大致相同，但仍有可能形成较大误差，采取一种更好的分类方法值得我们进一步研究。c） 地面实测的为0.34um的宽波段实际反照率，而modis我们取的其0.32.5um宽波段的黑半球反照率，虽然说植被在0.32.5um占了0.34um的99以上的能量，但是可能还是存在一定的误差。不过总的来说，由于我们加入了先验知识，使得我们的混合象元分解更为精确，得到的值也比直接使用混合象元更为可靠。参考文献1 Wu JianpingAn Attempt to Observe Paths of Particle of Wind Flowing over Buildings Utilizing Simplified Aerial Photogrammetry1992，ACRS。2 Zhukov, B., Oertel, D., Lanzl, F., Reinhackel, G., Unmixing-Based Multisensor Multiresolution Image Fusion. IEEE Transactions on Geoscience & Remote Sensing, 37(3), 1212-26, May 1999。3 李小文, 王锦地, 胡宝新, 等先验知识在遥感反演中的作用. 中国科学, D辑, 1998, 28(1): 6772。4 阎广建，朱重光，王锦地等遥感反演中约束最优化方法的拓展遥感学报，2002，Vol.6 No.2:81-87。5 /modis/modis-intro.htm.中分辨率成像光谱仪（MODIS）简介。6 Christophe FRANCOIS，Catherine OTTLE, Albert OLIOSO, Laurent PREVOT,Nadine BRUGUIER, Yannick DUCROSConversion of 4001100 nm vegetation albedo measurements into total shortwave broadband albedo using a canopy radiative transfer modelAgronomie 22 (2002) 611618. INRA, EDP Sciences, 20027 刘闯,葛成辉美国对地观测系统(EOS)中分辨率成像光谱仪(MODIS)遥感数据的特点与应用，遥感信息，2000-3：45-48。8 唐世浩，朱启疆，闫广建. 遥感地表参量反演的理论与方法北京师范大学学报（自然科学版）, 2001，Vol.37 No.2：266-273。Practice of Quantitative Remote Sensing Model Library Based on COM TechniqueAbstract：With the development of remote sensing, new models are available continuously. In order to extend the practicability of the model library, the authors introduce component object model (COM) technique. COM is a software architecture that allows the components made by different software vendors to be combined into a variety of applications. Remote sensing model library is composed of three parts, common objects, model objects and accessorial objects. Common objects include input/output procedure, solar angle calculating procedure in bi-directional reflectance, and metadata about all the models. Model objects include the models contributing to quantitative remote sensing applications, which comprise system models, simulant models and application models. Accessorial objects include prior knowledge, measurement data and image data. In the article, the executable project is validated with an instance in th