分层变端元混合像元分解的新疆北部积雪分量制图研究

刘艳，杨耘，李杨

（1．中国气象局乌鲁木齐沙漠气象研究所，新疆 乌鲁木齐830002；2．长安大学地质工程与测绘学院，陕西 西安710054）

具有不同波谱属性的物质出现在同一个像元内时会出现波谱混合现象，这类像元便称为混合像元［1］。因传感器空间分辨率的限制和被探测目标的多样性，遥感影像中往往会产生混合像元［2］。理论上，形成混合像元的原因主要有以下几方面：1）单一成分物质的光谱、几何结构及其在像元中的分布；2）大气传输过程中的混合效应；3）传感器本身的混合效应［3］。其中，1）是线性效应，指混合像元中各端元间相互独立互不影响时，混合像元光谱是该像元内各端元光谱的线性叠加［2］；2）和3）为非线性效应，它是由于端元（纯像元）间的散射传输路径和遥感仪器的混合效应所引入的光谱非线性叠加［2］。大气校正可以对2）进行修正，而传感器本身的混合效应可以通过仪器校准、定标加以部分克服。根据混合像元光谱的产生机理，研究人员建立了许多分解模型，主要包括线性、非线性分解模型、模糊监督和神经网络等。线性分解模型因物理含义明确、建模简单，得到了广泛应用［4－9］。

新疆北疆地区冬季多云天气较多，再加上地形和植被的影响，导致该地区积雪分布不均，这使得高精度的积雪遥感监测难以实现。从多光谱、高分辨率影像（ETM＋、QuickBird等）中虽能获取中、高空间分辨率的积雪覆盖数据，但这类影像具有如下特点：商业遥感影像价格昂贵且影像覆盖范围不大，时间分辨率有限。而新疆地域广阔，冬天地面积雪随时间变化频繁。因此，高空间分辨率影像难以在气象业务服务中得到广泛应用。相对于高空间分辨率影像，中、低空间分辨率影像（MODIS、FY－3）可全球免费获取，时间分辨率高，易获取大范围晴空条件下的影像数据，便于在业务中使用。但是，因传感器空间分辨率的限制，MODIS等中低分辨率影像存在大量的混合像元，需要对混合像元进行分解，才能提高制图精度。有学者利用中等分辨率成像光谱仪MODIS积雪数据对新疆或天山南、北地区进行积雪、冰川和融雪径流研究［10－13］，但提取积雪面积以像元为单位，考虑混合像元问题的研究甚少。陈晓娜等［14］以MOD02HKM为基础，通过线性分解模型对新疆天山中段MODIS影像进行分解，从中提取积雪面积；延昊和张国平［15］从甚高分辨率辐射计（advanced very high resolution radiometer，AVHRR）数据反演积雪盖度；Painter等［16］利用机载可见光／红外成像光谱仪（airborne visible infrared imaging spectrometer，AVIRIS）影像结合DISORT模型并考虑了亚像元对雪粒径反演的影响，发展了MEMSCAG（MODIS snow－cover area and gain size）模型同步反演雪粒径，取得了较好效果。通过遥感分析积雪覆盖、草场植被和成片森林等非立体空间目标时，由于端元间散射效应弱，线性效应在混合光谱的成因中起绝对主导作用［2］。但是，线性分解模型采用固定端元对整个影像进行解混，在实际应用中一个像元可能由2、3或更多个端元组成，采用固定端元的解混方法会产生端元分解过剩的问题，特别是新疆北疆这种土地覆盖和积雪分布状况。

多端元混合像元分解法 （multiple endmember spectral mixture analysis，MESMA）是传统混合像元分解法（spectral mixture analysis，SMA）的一个扩展，它的基本思想是基于像元为每类地物选取多条光谱，并以此生成多个端元组合（每个端元组合由不同地物的光谱组成），然后对每个像元搜索最小二乘误差最小的端元组合，进而反演出每个像元的端元分量［17］。目前，MESMA被广泛应用于高光谱遥感影像城市土地覆盖分类中［17－22］。但是，MESMA仍存在如下问题：用于解混的端元模型中若端元数过少，导致模型难以达到满意的解译结果；如果模型中包含的端元个数过多，测量光谱和模拟光谱间的轻微偏离通常分配给模型所用端元，实际上这样的端元并不存在［21］。在此背景下，针对新疆北疆地区特有的地理及气候条件，本文利用MODIS数据，提出了新的分层变端元混合像元分解（hierarchical dynamic endmember spectral mixture analysis，DESMA）方法，进行积雪分量的遥感制图研究。

1 材料与方法

1．1 研究区和数据

本文所用遥感影像是2012年1月6日获取的、覆盖新疆北疆地区的MODIS 1B数据，空间分辨率为500m，包含7个波段（MODIS 1～7）。影像中有一定量的积雪覆盖，适合积雪分量遥感制图研究（图1）。数据进行了FLAASH（fast line－of－sight atmospheric analysis of spectral hypercube）大气校正、辐射校正和几何纠正等预处理。同时，收集同期研究区HJ－1BCCD数据和MODIS 8d合成最大积雪覆盖产品（MOD10A2）。

图1 研究区域位置Fig．1 Location of the study area in Northern Xinjiang

1．2 分层变端元混合像元分解方法基本思想

该方法的基本思想是根据地类特征对影像进行逐级分类，并逐级逐类地解混，将前一级解混所得A类地物的空间分布作为下一级其子类分解时的空间限制，采用既定模型仅在A区域内对影像进行解混（图2）。在每一级分解时均假定像元在某一波段的反射率是以构成该像元的端元（纯像元）反射率及其所占像元面积比例（分量）为权重的线性组合。其中，线性混合像元分解模型的表达式［23］为：

冬季研究区内耕地、草场和水体等均匀地表会被积雪完全覆盖，云杉、荒漠植被和积雪混合分布，部分区域存在裸露地表。针对地物这种分布特征，第1级时将影像分为积雪－植被和土壤两大类。然后仅对积雪－植被区采用积雪和植被两端元模型进行解混，再消除山体阴影的影响。最后在植被区实施植被和积雪两端元模型分解。综合上述3级分解结果，最终获取积雪分量空间分布。因滤去了第1级中与积雪无关地类，端元数减少，降低了地物错分性。

图2 分层变端元混合像元分解法流程图Fig．2 Workflow of hierarchical dynamic endmember spectral mixture analysis（DESMA）

1．3 分层变端元混合像元分解制图

1．3．1 建立端元库 研究区主要地物有云杉、荒漠植被、耕地、裸土、草场、水体等，冬季大部分地物被积雪覆盖。因此，针对研究区地物特征可建立植被－积雪－土壤的端元模型。可选端元库包含影像端元和参考端元两种。影像端元由研究区MODIS影像获取。结合野外雪深GPS记录点，在ENVI图3环境中实现MODIS R（6）G（2）B（1）假彩色合成积雪（蓝色）、植被（墨绿色）、裸土（红色）的感兴趣区域选择，每类地物对应1～3个像元的矩形区域，由ENVI／VIPER TOOL［21］实现影像端元库建立，包含91条光谱，其中积雪光谱19条、植被光谱49条、土壤光谱23条。参考端元源自全波段地物光谱仪测量积雪、荒漠植被、裸土光谱数据。测量光谱经MODIS传感器响应函数（spectral response function，SRF）［22］转换后利用 ENVI／Spectral Library Builder建立参考端元库，包含13条积雪、2条植被和4条土壤光谱。

图3 MODIS R（6）G（2）B（1）假彩色合成图Fig．3 False color composite image of MODIS of the study［Red（6）Green（2）Blue（1）］

1．3．2 初选端元 利用EAR（Endmember Average RMSE）计算A端元光谱集内某一条光谱用A端元内其他光谱来分解产生的残差，EAR越低，表明这条光谱的代表性越好，很高则证明这条光谱可能是离群点、没有代表性［22］。因此，计算可选端元库中各端元的EAR，挑选EAR较小的光谱。结果选出14条积雪光谱，6条植被光谱和8条土壤光谱。EAR计算公式：

式中，i为端元光谱；j为参考光谱；n为参考光谱总数；RMSE为均方根误差。

1．3．3 优选端元 EAR初选后的端元光谱虽能很好地模拟其所在类的其他端元光谱，能否很好地对影像中所属地类进行分解，需进一步实验分析［22］。以积雪端元为例，首先，对积雪端元库中各条光谱的EAR进行升序排列，然后，采用EAR最低的端元1对整个影像进行2端元分解，得到该端元分解影像像元数，其后引入端元2，利用端元1和端元2对影像解混，依次类推，得到各端元光谱组合分解像元数（表1），判断各组合中每条光谱分解像元数是否大于44022（即分解像元数占整个影像像元总数的20%）。

1．3．4 逐级逐类解混和制图 为了解释地物因卫星传感器角度、地形起伏和其他一些阴影，所有模型都包含阴影端元，用来表征端元亮度变化，而阴影端元并非地物组成成分。因此，需要对阴影进行归一化以消除阴影端元，即用每个像元的非阴影端元分量除以全部非阴影端元分量的总和，最后得到每个像元实际组成端元的分量。

2 结果与分析

2．1 典型端元库

大气校正后MODIS通道1存在部分像元反射率大于1的反射率超饱和现象和通道5存在坏数据的情况。因此，本文选用MODIS 2、3、4、6和7通道反射光谱数据构建端元光谱库。依据上述端元库建立原则，最终选取6条积雪（表1中1、7、9、10、13和14）、5条土壤和4条植被光谱建立了适合研究区的积雪－植被－土壤端元库（图4a，b，c）。

表1 积雪端元优选原则Table 1 Criteria used to select snow endmembers for DESMA library

图4 优选积雪（a），土壤（b），植被（c）端元的光谱Fig．4 Spectra contained in optimal endmember library：snow（a），soil（b），vegetation（c）

2．2 分级分解结果

由于仅收集到研究区同期HJ－1BCCD数据，缺少可用的近红外数据IRS，仅利用可见光波段采用非监督分类和手动校正方法，提取积雪覆盖，同时，利用夏季HJ－1BCCD获取植被分布图（图5），将其与分级结果进行比对。

2．2．1 第1级分解结果 第1级分解时，2端元或3端元分解模型对分量最大和最小值不做限制。首先，将包含6条积雪、5条土壤和4条植被光谱的端元输入第1光谱库中，分别进行积雪－阴影、土壤－阴影和植被－阴影的2端元分解，端元分解次数为15；将包含积雪和植被的端元光谱输入第1光谱库，土壤端元光谱输入第2光谱库，进行积雪－土壤－阴影和植被－土壤－阴影的3端元分解，端元分解次数为50。然后，比较2端元和3端元分解模型的RMSE，如果RMSE变化大于0．1，选择3端元分解结果，否则选择2端元。比较发现，2端元模型RMSE整体较小，表明选择2端元即可以完成分解，分解像元数达100%。合并2端元模型生成的积雪－阴影和植被－阴影类，生成积雪－植被区类型图（图6），吉木乃、黑山头和福海附近存在裸露地表，阿勒泰和富蕴山区部分山坡无积雪覆盖，与HJ－1BCCD分类图（图7）整体分类结果一致。但是，由于地形起伏和阴影影响，阴坡积雪区域呈现出明显的非积雪特性，因此，仅利用HJ－1BCCD可见光通道信息进行分类，这部分区域被分成了非积雪区，造成了漏分。

图5 研究区HJ－1BCCD假彩色合成图Fig．5 False color composite image of HJ－1BCCD

图6 第1级：2－端元和3－端元模型解混分类图Fig．6 Merged result of 1st layer derived from 2－endmember and 3－endmember model of DESMA

2．2．2 第2级分解结果 第1级分解后得到积雪－植被的空间分布，将其对整个MODIS影像进行掩膜运算，得到积雪－植被区的MODIS数据，将包含积雪和植被的端元光谱输入第1光谱库中，进行积雪－阴影和植被－阴影的2端元分解，模型分解像元最小和最大分量控制在（－0．05，1．05），端元分解次数为10（图8），吉木乃和黑山头附近存在荒漠植被，哈巴河、阿勒泰和富蕴中山带有云杉分布和山体阴影，前山带靠近古尔班通古特沙漠边缘区域有荒漠植被，与 HJ－1BCCD R（4）G（3）B（2）假彩色合成图（红色区域为植被）（图5）分类结果一致。发现在积雪和植被混合区出现漏分现象。

图7 研究区HJ－1BCCD积雪覆盖图Fig．7 Snow－cover map of HJ－1BCCD of the study area

图8 第2级：掩膜后积雪－植被区MODIS影像2端元分解分类图Fig．8 Classification of 2nd layer derived from 2－endmember model of DESMA for MODIS image in snow and vegetation area

2．2．3 第3级分解结果 第2级分解后得到植被区主要植被类型有荒漠植被和天山云杉两类，考虑山体阴影与云杉易错分，引入山体阴影端元，对掩膜后的植被区MODIS数据进行植被－阴影和山体阴影－阴影2端元模型，模型分解像元最小和最大分量控制在（－0．05，1．05）且满足RMSE≥0．0025，将山体阴影从植被区中去除（图9）。

图9 修正了部分山体阴影的第2级分类图Fig．9 Classification result of level 2excluded part of the mountain shadows

修正了部分山体阴影的植被区，还会存在部分积雪，将积雪光谱输入第1光谱库，植被光谱输入第2光谱库，对植被区MODIS数据进行积雪－植被－阴影的3端元分解，端元分解次数为24。

2．3 积雪分量制图与精度验证

叠加第2级积雪分量图和第3级山体阴影修正后植被区积雪－植被－阴影3端元分解归一化后所得积雪分量，得到研究区积雪分量图。2012年1月在北疆地区进行了雪深、积雪覆盖度的野外观测（图10和表3），辅助记录了样点周围1km场景。将积雪分量与29个实测点进行比对，其中25个点正确，分量准确率为87%（表3）。

同时，收集同期MODIS 8d合成最大积雪覆盖产品（MOD10A2）与研究区积雪分量图进行比对分析（图11），结果发现，29个实测点上MOD10A2均监测为积雪（200表征像元被积雪覆盖），但是MOD10A2无雪区面积（白色区域）比本文提出的分层变端元混合像元分解生成的积雪分量图中无雪区域面积大很多（白色区域），MOD10A2中无雪区比例为8%，本文计算积雪分量图中无雪区域占3%。冬季研究区内湖泊被积雪覆盖，MOD10A2未对湖泊和湖冰进行积雪识别，占研究区面积比例为1%。同时，分层变端元混合像元分解法考虑区域地貌特征和积雪分布状况，以积雪为核心地类，各分级分类时，与积雪相关的地类利用2端元或3端元模型进行解混，为下一级积雪分量的定量分析提供相对纯净的数据环境，在山区云杉、山体阴影和积雪以及荒漠植被、积雪交错分布的地貌环境下，与MOD10A2仅从像元尺度上提供积雪和非积雪空间信息相比，较大程度地提高了积雪识别率，真实反映了积雪的空间分布特征。

3 结论与讨论

针对研究区土地覆盖类型，提出了基于最小EAR和少端元分解模型的MESMA方法来实现以积雪为核心的地类分级分类方法。该方法通过逐层删除非积雪地类，最终挑选出积雪，其适用条件宽泛，积雪覆盖区均能适用。实验结果表明，变端元解混时选择较少的端元个数更有利于提高积雪制图精度。相对于MESMA［17－22］，本文采用的最小EAR和选择2－端元分解模型构建典型端元库的方法，可使积雪制图精度高达87%，即选择较少的端元个数更有利于提高。

表3 模型计算值与实测值比较Table 3 Comparison between the result from our model and field data

图10 研究区积雪分量与积雪盖度测量点位置分布图Fig．10 Snow fraction map generated from DESMA of the study area and distribution of the measurement point of snow cover

图11 研究区MODIS积雪识别产品（MOD10A2）Fig．11 MODIS snow products of MOD10A2of the study

同时，本文方法具有很好的扩展性，可推广应用到我国风云三号中分辨率光谱成像仪（FY－3／MERSI）数据，有利于提升中低分辨率影像积雪覆盖精度，为气象业务工作中的融雪径流模拟、雪灾监测等提供数据和技术支持。在上述应用中，MODIS和FY－3／MERSI数据中积雪占混合像元比例多少会直接影响到积雪、植被和土壤等端元光谱的提取，进而影响到积雪分量求解结果。因此，当使用上述中分辨率遥感影像数据时，该方法适用于新疆北疆地区冬季积雪累积期的混合像元分解，这是由于该区域覆盖了积雪、植被、裸土和戈壁等地貌，地貌类型相对简单。随着遥感影像空间分辨率的提高，如采用ETM＋等影像时，积雪占混合像元的比例对于积雪等端元光谱的提取精度影响会随之降低，该方法适用性也会有所提高。

［1］邓书斌，陈秋锦．基于MTMF的混合像元分解方法研究［A］．第十七届中国遥感大会摘要集［C］．南京：中国遥感应用协会，2010：69－73．

［2］潘明忠，元洪兴，肖功海，等．基于高精度端元的混合像元线性分解模型研究［J］．红外与毫米波学报，2010，（5）：357－361．

［3］阎国倩，赵云升，宁艳玲，等．混合像元光谱特征影响因子分析［J］．光谱学与光谱分析，2009，29（12）：3358－3361．

［4］Chabrillat S，Pinet P C，Ceuleneer G，etal．Ronda peridotite massif：Methodology for its geological mapping and lithological discrimination from airborne hyper spectral data［J］．International Journal of Remote Sensing，2000，21：2363－2388．

［5］李慧，陈健飞，余明．线性光谱混合模型的ASTER影像植被应用分析［J］．地球信息科学，2005，7（1）：103－106．

［6］Small C．Estimation of urban vegetation abundance by spectral mixture analysis［J］．International Journal of Remote Sensing，2001，22（7）：1305－1334．

［7］张熙川，赵英时．应用线性光谱混合模型快速评价土地退化的方法研究［J］．中国科学院研究生院学报，1999，16（2）：169－176．

［8］万军，蔡运龙．应用线性光谱分离技术研究喀斯特地区土地覆被变化——以贵州省关岭县为例［J］．地理研究，2003，22（4）：440－443．

［9］邹蒲，王云鹏，王志石，等．基于ETM＋图像的混合像元线性分解方法在澳门植被信息提取中的应用及效果评价［J］．华南师范大学学报（自然科学版），2007，（7）：131－136．

［10］王玮，黄晓东，吕志邦，等．基于MODIS和AMER－E资料的青藏高原牧区雪被制图研究［J］．草业学报，2013，22（4）：227－238．

［11］冯琦胜，张学通，梁天刚．基于MOD10A1和AMSR－E的北疆牧区积雪动态监测研究［J］．草业学报，2009，18（1）：125－133．

［12］裴欢，房世峰，覃志豪，等．基于遥感的新疆北疆积雪盖度及雪深监测［J］．自然灾害学报，2008，17（5）：52－57．

［13］李宝林，张一驰，周成虎．天山开都河流域雪盖消融曲线研究［J］．资源科学，2004，26（6）：23－29．

［14］陈晓娜，包安明，张红利，等．基于混合像元分解的MODIS积雪面积信息提取及其精度评价——以天山中段为例［J］．应用气象学报，2004，15（6）：665－671．

［15］延昊，张国平．混合像元分解法提取积雪盖度［J］．应用气象学报，2004，15（6）：655－661．

［16］Painter T H，Dozier J，Robiters D，etal．Retrieval of subpixel snow－covered area and grain size from imaging spectrometer data［J］．Remote Sensing of Environment，2003，25（1）：64－77．

［17］Roberts D A，Gardner M，Church R，etal．Mapping chaparral in the Santa Monica mountains using multiple endmember spectral mixture models［J］．Remote Sensing of Environment，1998，65（3）：267－279．

［18］Roberts D，Halligan K，Dennison P．VIPER Tools User Manual Version 1．5［EB／OL］．（2007－12－11）．http：／／www．vipertools．org．

［19］Dennison P E，Halligan K Q，Roberts D A．A comparison of error metrics and constraints for multiple endmember spectral mixture analysis and spectral angle mapper［J］．Remote Sensing of Environment，2004，93（3）：359－367．

［20］Song C．Spectral mixture analysis for subpixel vegetation fractions in the urban environment：How to incorporate endmember variability［J］．Remote Sensing of Environment，2005，95（2）：248－263．

［21］马孟莉，朱艳，李文龙，等．基于分层多端元混合像元分解的水稻面积信息提取［J］．农业工程学报，2012，8（2）：154－159．

［22］Frankea J，Robertsb D，Halliganc K，etal．Hierarchical multiple endmember spectral mixture analysis（MESMA）of hyperspectral imagery for urban environments［J］．Remote Sensing of Environment，2009，93（3）：1712－1719．

［23］饶萍．EOS－MODIS像元组分分解中端元的选择与改进［A］．《测绘通报》测绘科学前沿技术论坛摘要集［M］．北京：测绘出版社，2008：1－10．