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

1、基于先验知识进行混合象元信息分解 *万华伟北 京 师 范 大 学 遥 感 与 GIS 研 究 中 心 ，资 源 与 环 境 科 学 系 ，环 境 遥 感 与 数 字 城 市 北 京 市 重 点 实验 室北京， 100875摘要：混合象元一直是影响遥感图像精度的一个重要因素。本文通过使用地物先验知识，利用高分辨率图像的分类图做限制条件，计算低分辨率的图像的象元地物的百分比，然后通过线性模型做混合象元分解，然后又通过分类图提供的地物信息作为先验知识，使用约束性的最小二乘法反演，从而得到比较精确的地物的值。我们利用了 MODIS43B3 反照率产品作为低分辨率图像，ASTER 的可见光和近红外波段作

2、为高分辨率图像，进行了实验，然后与地表的实测相对比，得到比较令人满意的结果。关键词：混合象元，线性分解，先验知识，遥感1 引言在我们获得的遥感影象上，基本上都存在着混合象元。根据估计，在我们使用的MODIS 的 ALBEDO 产品（1*1 km2）上，90% 以上为混合象元。这样，我们在使用这些数据时，如果不加考虑就直接使用，无疑会大大地降低我们研究的精度，而且引起信息量的丢失。因此对混合象元的分解是提高遥感信息的精度的一个重要方式 1。目前的混合象元处理有很多方法，当前混合象元有两个基本方法：1）模糊分类方法；2）确定一个象素的光谱响应和其中端元的光谱与其所占面积比的关系 1。对于混合象元的

3、分解，Boris Zhukov 等提出了约束和非约束的 MMT 算法，提出利用高分辨率图像和低分辨率数据作融合 2，但是他们在反演时都没有充分地利用前面由分类图所增加的先验知识。李小文等指出，当地表的先验知识存在时，不加利用是可惜的 3。在本文中，我们在对混合象元进行分解时，充分利用先验知识，在使用高分辨率作为限制条件的同时，利用引入先验知识的最小二乘分解，加入了根据分类图得到的“边界” 。并且利用 MODIS 的反照率产品（分辨率为 1km）做的实验。使用了同地区 ASTER 的图像前三波段（分辨率为 15m）作为限制条件。通过进行混合象元的分解得到植被的反照率，将得到的数据与地表实测值以及

4、如果直接使用混合象元的值作的对比以及误差分析，得出结论，经过引入先验知识的混合象元分解的误差比直接使用混合象元值低了很多，将误差从 26.58%降低到 16.89%。2 算法描述21 图像的配准* 973 项目资助算法的流程图及图像的表示如 Fig1低分辨率图像（含研究需要的信息）1km，记 LI高分辨率图像（含地物详细的空间分布信息）15m，记HI1515 聚合窗聚合中间尺度图像225m（由高分辨率图像生成） ，记 MHI找控制点，配准中间尺度图像225m（由低分辨率图像生成） ，记 MLI具有相同的坐标系Fig1：算法的流程图及图像的表示我们首先将低分辨率图像和高分辨率图像配准，在这里，由

5、于两幅图像的分辨率相差太多（1km 和 15m） ，所以我们选取一中间尺度图像，我们先将高分辨率图像进行向下的尺度转换，我们直接使用 1515 的聚合窗来进行重新采样。这样形成一幅 225m 分辨率的图像，我们再用它与 1km 的进行配准，这样配准的精度会比原来直接配准的精度提高许多。这样就生成了同一坐标系下的同为 225m 分辨率的两幅图像。22 找出 LI 和 MLI 之间的对应关系首先找到 LI 与 MLI 的对应关系，MLI 是通过最近邻插值到 MHI 的坐标系上的、由 LI生成的图像。最近邻法是最简单的一种插值方法，这种插值的过程中，DN 是是不发生变化的，只与最邻近的象元有关，和其

6、他的象元均无关系。采用搜索算法来找到它们之间的关系，此算法的依据是插值的过程中，虽然由 LI 上一个象元生成的 MLI 的象元是不规则的，但是是连续的。1建立一个和 MLI 一样大小的 mark 文件；2从（1，1）开始，设（1，1）为基本点，记录其 value，然后搜索其 8 个邻域，如果值为 value，则记录其行列值，将该象元加入该基本点的一组，并将 mark 值设为真，再继续搜索它的 8 个邻域象元。如果值与 value 不等，则停止。3行列值增加，记录第一个 mark 不为真的象元为基本点，重复步骤 2。直至所有的象元的 mark 值都为真。这样，便可以找到一系列的由 LI 上一个象

7、元生成的 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）这四个边界点之间的

8、225 个象元，便对应着 MHI 上的一个象元。24 找出 LI 和 HI 之间的对应关系并求出其面积百分比由于 MLI 与 MHI 已经配准到同一坐标系下，也就是说通过以上两步找到的 LI 与 MLI和 HI 与 MHI 之间的关系，已经找出了，LI 与 HI 之间的关系。LI 上的一个象元，对应着MI 上的 NUM225 个象元。我们在上面找到了 LI 上一个象元与 HI 上的 NUM225 个象元以后，由于，HI 图像为分类图像，其中，每个象素所属的地物类别都已知，这样，我们也就知道了 LI 上这一象元的各类地物的象素数，然后除以总的象元数，便得到了该象元的每类地物的百分比。，其中，cl

9、assi代表 i 类的面积百分比， 代表delta(cs=i)clsiNUM*25 delta(cs=i)class=i 的象元数。25 使用最小二乘模型进行亚象元分解目前大部分反演采用的方法是最小二乘方法，其代价函数常用的形式为： 21()Miiisyfx其中，M 为测量总数， 为观测数据， 为模型参数， 为由参数 得到的模型结果i ix()ifxix4。理论上说，当观测数据足够多时，搜索能产生使 s 最小的一组参数值，也就是参数的最佳估计值。然而在实际情况的，由于传统的最小二乘方法并不考虑其参数的物理边界，有可能出现无物理意义的值，因此有发展了约束化的最小二乘方法，我们在这里，利用地物分类

10、后的所属类别所在的物理含义的边界，作为先验知识加以约束。将代价函数修改为：，2211()()MNiijji isyfxx其中，y，x 为观测值和参量值；M，N 为样本数和参量数；为地物的先验值 8。这样，我们得到的目标的值就会比纯粹的使用混合象元的值精确的多。3 实验数据收集与处理31 实验数据描述在本文中，我们使用的低分辨率数据是 MODIS 的反照率产品MOD43B3，MODIS（中分辨率成像光谱仪）是美国宇航局研制大型空间遥感仪器，它在36 个相互配准的光谱波段、以中等分辨率水平（0.25Km1Km） 、每 12 天观测地球表面一次，获取陆地和海洋温度、初级生产率、陆地表面覆盖、云、汽溶

11、胶、水汽和火情等目标的图像，其 17 波段均为研究陆地和云性质的 5。它的标准的陆地产品，如 Table1 所示：陆地MOD09MOD11MOD12MOD13MOD14MOD15MOD17MOD43陆地表面反射比陆地表面温度陆地覆盖植被指数火情叶面指数和部分光合作用辐射净初级生产率/光合作用BRDF/反照率Table1MODIS 的标准陆地产品我们在这里的实验数据，低分辨率数据为MODIS的反照率产品，mod43b3，其分辨率为1km，它是由16天的不同角度的数据反演得到的，这里，我们通过网上订购2001年第113天的一景图像。它的波段分布Table2 所示，我们采用其0.3 5.0um宽波段

12、的黑半球反照率。我们使用的低分辨率约束图像Aster使用2001年8月23日成象的中国河北栾城地区的一景图像。我们的地面实测数据是栾城生态试验田（ ， ） ，为0.35um的宽波段实140。 375。际反照率。我们看植被的反射率曲线Fig2，可以知道在0.3 2.5um 之间，其能量大约占到0.35um的99以上，在大多数的辐射能量中，对于植被的反射率，在2500nm之后是可以忽略不计的 6。因此我们可以使用他们两个直接对比。波段序号 位置和宽度 （nm）12620670841876345674594795455651230125016281652210521558910300-700 700

13、-5000300-5000Table2 MOD43B3 波段分布 Fig1 典型植被光谱32 数据处理过程321 对 MODIS 数据进行处理MODIS 的存储格式为 HDF 格式，它采用分层次的串块型数据格式 -HDF（Hierarchical Data Format） 。这种数据格式更适用于大数据量的快速传输、存储和提取 7。但是目前大多数的通用遥感处理软件还不支持这种数据格式，我们先使用 MODIS TOOL 将其按波段转出。322 对 aster 图像进行分类对 aster 图像进行分类，得到 15m 分辨率的分类图。分类的精度对于我们混合象元分解影响比较大，所以我们应该尽量使分类达到

14、一个高的精度。在这里，我们先通过目视确定地物有哪些类，然后确定类数，使用监督分类。我们的 ASTER 图像，为 2001 年 8 月 23日成像，这一时期，当地的地表覆盖主要有农作物，还有城镇和水体。在这里我们主要是对这三类地物进行亚象元分解，以消除城镇和水体的影响。323 图像配准MOD43B3 为 1km 尺度，而 aster 为 15m 尺度，直接配准误差会很大，所以我们先对aster 进行向上的尺度转换，进行 15*15 的聚合，分辨率变为 225m。然后从 MODIS 和 aster图像上找到 30 个控制点，以聚合后的 ASTER 图像为基图，以 MODIS 图像为要纠正的图，进行

15、图像的纠正，使得两幅图像在同一坐标系下。配准的结果如 Fig2。MODIS-225m 分辨率 ASTER-225m 分辨率（标准假彩色）Fig2.配准的结果324 混合象元分解通过上面提到的搜索算法，记录 modis1km 和 modis225m 的对应关系。即一对多的关系，一个象元对应多个象元，我们记作 NUM。然后又把 modis225 和 aster15m 的关系找到，这样就找到了，mod43b3 产品的一个象元所对应的各个类的部分。我们通过计算，classi代表 i 类的面积百分比delta(cs=i)clsiNUM*25我们便完成了对 MODIS 象元的分解。可以得到每个象元的各类所

16、占的面积百分比。在这里我们选择植被和城镇的信息尽量分布均匀的象元。选取如下 45 个混合象元进行反演：序号 图像值 植被百分比城市百分比序号 图像值 植被百分比城市百分比1 173 0.633778 0.366222 24 191 0.327778 0.6722222 174 0.608889 0.391111 25 194 0.300833 0.6991673 192 0.436741 0.563259 26 177 0.240556 0.7594444 162 0.585833 0.414167 27 187 0.173056 0.8269445 167 0.694815 0.305185

17、 28 192 0.455704 0.5442966 172 0.694167 0.305833 29 204 0.152500 0.8475007 185 0.698963 0.301037 30 190 0.615111 0.3848898 173 0.613333 0.386667 31 207 0.649722 0.3502789 169 0.463611 0.536389 32 197 0.575948 0.42405210 185 0.688296 0.311704 33 191 0.536667 0.46333311 203 0.653889 0.346111 34 180 0.

18、612593 0.38740712 193 0.543111 0.456889 35 184 0.606370 0.39363013 186 0.668611 0.331389 36 193 0.591389 0.40861114 190 0.689412 0.310588 37 181 0.493611 0.50638915 194 0.448333 0.551667 38 163 0.593725 0.40627516 199 0.545621 0.454379 39 170 0.638693 0.36130717 195 0.677222 0.322778 40 170 0.633580

19、 0.36642018 191 0.685630 0.314370 41 171 0.566519 0.43348119 191 0.493333 0.506667 42 156 0.627451 0.37254920 188 0.637778 0.362222 43 158 0.681046 0.31895421 193 0.692222 0.307778 44 185 0.631944 0.36805622 185 0.596111 0.403889 45 168 0.654444 0.34555623 182 0.448333 0.551667324 基于先验知识的最小二乘反演在这里，为

20、了反演的更加精确，我们把只包括植被和城镇两类的象元选择出来，根据公式 iveg, city,AlbedolAlbedoivegictyPrcntPrn使用代价函数，其中，y,x 为观测值和残量值；M,N 为样本2211()()MNiijji isyfxx数和参量数； 为参量的先验值；我们在这里先不考虑样本之间的相关性，我们假定j它们之间相互独立，对于结果同等重要。我们先设定其硬边界，植被的反照率在宽波段0.35um 的先验值为 100200 之间。这样我们得到的结果，植被的反照率为 0.1732。4 结果与讨论我们通过以上处理得到反演结果如 Table2，我们又选取了植被占 90%以上的象元的

21、值取均值，通过这三种数据的比较，我们计算它们的误差，这里，定义误差Error=（计算得到的值-实测值）/实测值Table2.反演结果比较这样，我们反演得到的值的误差为 16.89%，而如果直接使用混合象元会造成 26.58%的误差，可见我们经过反演使误差减少了很多，不过在这里我们分析之所以误差这么大的原因，大概是因为a） 配准的误差，配准的精度虽然已经很高，但是可能仍然会引起误差。b） 分类图的误差，我们使用的 8 月份的图像进行非监督分类，而且只是分成了水体，城镇和植被三类，虽然经过实地考察，8 月份和 3 月份的地表覆盖大致相同，但仍有可能形成较大误差，采取一种更好的分类方法值得我们进一步

22、研究。c） 地面实测的为 0.34um 的宽波段实际反照率，而 modis 我们取的其 0.32.5um宽波段的黑半球反照率，虽然说植被在 0.32.5um 占了 0.34um 的 99以上的能量，但是可能还是存在一定的误差。不过总的来说，由于我们加入了先验知识，使得我们的混合象元分解更为精确，得到的值也比直接使用混合象元更为可靠。地面实测数据 经纬度定位单点值Error 最小二乘反演的植被反照率Error0.2084 0.153 26.58% 0.1732 16.89%参考文献1 Wu JianpingAn Attempt to Observe Paths of Particle of Wi

23、nd 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 李小文, 王锦地, 胡宝新, 等先验知识在遥感反演中的作用. 中国科学

24、, D 辑, 1998, 28(1): 6772。4 阎广建，朱重光，王锦地等遥感反演中约束最优化方法的拓展遥感学报，2002，Vol.6 No.2:81-87。5 http:/ ）简介。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 cano

25、py radiative transfer model Agronomie 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 TechniqueAbstr

26、act：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

27、 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 mode

28、ls. 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

29、an instance in the end. The spectrum of typical land surface observed in nadir viewing direction is simulated in pixel scale. In this process, crop model, PROSPECT model and SAIL model are used to calculate the spectrum character of the pixel. Crop model is used to simulate leaf area index, and PROS

30、PECT model is to simulate reflectance and transmission of leaves. The final result is calculated with SAIL model. The simulated spectrum of winter wheat observed in nadir viewing direction is at 11 oclock Apr.15, 2001.In the article, each model and each common object are designed to be components. The model library based on COM technique adapts to the progress and is propitious to be expanded and modified.Key Words: Quantitative Remote Sensing, Model Library, COM Technique