基于混合像元分解的天山典型地区冰雪变化监测

金 鑫，柯长青

(南京大学地理信息科学系，南京 210093)

基于混合像元分解的天山典型地区冰雪变化监测

金 鑫，柯长青

(南京大学地理信息科学系，南京 210093)

针对中低分辨率遥感图像中存在大量混合像元，而传统的图像分类方法存在只能将某个像元归到某一类中，不能正确反映混合像元实际情况的问题。以新疆天山典型冰川覆盖区为例，根据TM/ETM+遥感图像的光谱特征，结合天山地区地表覆盖特点，在线性混合像元分解方法基础上，设计一种符合冰川地区特点的“冰雪－植被－裸露山体－阴影”端元组分模型。通过选择合适的端元并将其反射率值代入改进后的且满足约束条件的线性混合像元分解模型，得到各端元组分丰度图，进而精确提取出冰雪信息并计算其面积。1989年TM和2000年ETM+遥感图像冰雪信息提取结果表明，运用线性混合像元分解模型能很好地监测实验区的冰雪覆盖变化情况。

冰雪覆盖;混合像元分解;线性混合模型;端元选择;天山地区

0 引言

冰川对地表径流及地下水有着直接的影响，而且可揭示大范围多尺度气候波动的响应特征［1］。冰雪覆盖变化监测具有重要的理论和现实意义。天山作为中国主要的冰川区，是我国西北部干旱区最重要的淡水资源［2］。乌鲁木齐河发源于天山天格尔峰1号冰川，它是乌鲁木齐市的主要水源。我国已有天山1号冰川近50 a的观测数据［3］，但采用的传统手段主要为实地野外考察以及冰川观测站的定点观测。这种方法具有周期长、工作量大、费用高、人为误差影响显著等弱点，尤其在自然气候条件较为恶劣的高山冰川地区，实地考察观测极为困难。

随着空间技术和传感器技术的发展，卫星遥感手段已经成为监测地表覆盖变化最主要的途径之一。然而，以传统的目视解译结合监督分类为主的方法，耗费时间长且分类结果因人而异;大多数遥感图像分类算法都是基于纯净像元的分类理念，其分类结果往往难以准确反映地表覆盖的真实情况。

混合像元分解方法正是针对遥感图像中大量存在的混合像元［4－5］，将每一混合像元进行分解并求得端元组分(即覆盖类型)占像元的丰度(即百分含量)，这样由混合像元的归属而产生的错分、误分问题也就迎刃而解［6－7］。其中，考虑混合像元的反射率、端元的光谱特征和丰度之间的响应关系时，根据其他地面特性和影像特征影响的不同，形成了各种各样基于混合像元分解方法的光谱混合模型，包括线性光谱混合模型，基于决策树、人工神经网络或边界提取的光谱混合模型等［8－10］。然而，光谱混合模型并不完全决定混合像元分解结果的精度，因为混合像元分解模型的求解是根据像元的反射率、端元的光谱特征，通过最小二乘法求解端元丰度来实现的，其结果的精度在很大程度上取决于端元组分的选择［11］。结合冰川研究区的具体情况选择地表覆盖模型及端元组分，这对于精确获取冰雪覆盖信息具有重要意义。

本文在线性混合像元分解模型的基础上，结合天山地区地表覆盖特征，设计出符合冰雪地区特点的端元组分选择方法，利用线性混合像元分解模型对该区1989年和2000年获取的遥感图像中的冰雪信息进行精确提取并计算面积，分析冰雪覆盖的变化特征，并对冰雪信息提取误差进行探讨。

1 研究方法

1.1 技术路线

基于端元的线性混合像元分解的主要流程包括:遥感图像预处理、最小噪声分离(minimum noise fraction，MNF)变换、2维及3维散点图查看、端元组分选择与确定、像元纯度指数(pixel purity index，PPI)计算、端元组分光谱反射率计算、线性混合像元分解及其误差分析等(图1)。其中，端元组分的选择及地表覆盖类型的确定由人机交互判读完成，其余部分则利用ENVI等软件自动完成。

图1 基于端元的线性混合像元分解技术流程Fig.1 Flow chart of mixed pixel decomposition using linear model based on end－member

1.2 线性光谱混合模型

线性混合模型是在混合像元分解提取地物信息中最为常用的模型。它利用线性关系表达遥感系统中一个像元内各地物的类型、比例与地物的光谱响应，其基本假设是组成混合像元的几种不同地物的光谱以线性方式组合成混合像元的光谱，即假定混合像元内各成分光谱之间是独立的［12－14］。

线性混合模型的数学物理原理是假设每个到达传感器的光子只与一种地表覆盖类型有关。如果传感器采集到n个波段的数据，地表覆盖类型数量为c，则线性模型定义为式中:向量 X(n×1)为多光谱的像元观测值;向量f(c×1)为地表覆盖类型 i的端元组分比例;M(n×c)代表各波段上地表覆盖类型i的端元光谱值。端元光谱值指对纯像元观测得到的，例如只包括单一的地表覆盖类型的像元。

为了使这个线性混合模型符合客观规律，需要满足2个约束条件:

①像元各端元组分比例总和为1，即

②各端元比例不能为负数，即

式中j=1，…，c。如果所有端元组分光谱之间都是线性无关的，那么当c=n+1时，由式(1)(2)确定的方程组总有唯一解;当c更大时，也可以求得无限个解。然而实际上并不能得到精确解，原因是存在传感器的噪声。因此将式(1)改为［15－16］

式中:r为噪声，假定噪声r服从多元正态分布，均值为0，以Nf为协方差矩阵(依赖于像元的组成成分而定)。各端元组分对Nf的贡献率与其所占比例fi有关，因此也依赖于相应的地表覆盖类型的协方差矩阵Ni，其中

应用线性混合模型的关键是尽可能减少每个像元的误差，使均方根误差RMS为最小，即

式中:P为像元总数;Rmnj为j波段m行n列像元的误差值。

解算该模型最常用的方法是最小二乘法，即用

求解f。线性光谱混合模型建立在像元内相同地物都有相同的光谱特征以及植被指数具有线性可加性的基础上，并忽略了地表多次反射的情况，模型简单，物理含义明确，实际运算简便快速，对解决像元内的混合现象有较好的效果。

为了满足式(2)(3)约束条件，本文利用IDL编程对ENVI软件的最小二乘法求解线性混合像元分解模型进行改进，以得到合理的f值从而有利于冰雪信息的提取。

2 实验

2.1 遥感数据

以新疆天山1号冰川所在地为实验区。遥感数据包括1989年9月3日获取的TM和2000年9月17日获取的ETM+图像数据。图像质量较好，没有云的覆盖，经过了辐射校正与几何纠正处理。结合本研究的实际需要，选取TM和ETM+图像的第1，2，3，4，5，7 波段作为参加像元分解的波段，并对 2景图像进行几何配准，以便做变化分析。

2.2 端元选择

端元组分的选取直接影响线性混合像元分解模型的分类精度。本文采用影像端元选择方法，直接从待分类的图像中选择端元，然后不断修改、调整。这样，一方面可以省去投入人力和物力通过野外波谱测量获得影像波谱库(尤其是在条件恶劣的冰川地区);另一方面，从图像本身获取的端元也与图像数据具有相同的度量尺度。采用散点图和PPI端元选择方法进行影像端元选择。先对图像数据进行MNF变换，以减少噪声干扰和数据维数;然后对图像中的像元进行PPI迭代运算，找出图像中相对纯净的像元。

2.2.1 地表覆盖模型

依据实验区自然地理特征及遥感图像的成像特点，该区常年被冰雪覆盖［17］;在9月份有高山草甸、苔草、禾本科等植被生长，也有部分裸露的山体岩石;另外高山地区在阳光照射下产生了不少阴影区。由于提取目标主要是冰雪覆盖信息，而且为便于线性混合像元分解模型的求解，暂不考虑植被类别内部细类之间的光谱特征差异，因此将地表覆盖模型确定为“冰雪－植被－裸露山体－阴影”端元模型。

2.2.2 MNF 变换

MNF变换用以判定图像数据内在的维数(即波段数)，分离数据中的噪声，减少随后处理中的计算量。原6个波段的TM图像经MNF变换后，可以得到6个分量的MNF图像及其特征值图像，同时可以求得各分量占图像总信息量的百分比(表1)。

表1 MNF变换后各分量特征值及所占百分比Tab.1 Eigenvalue and percentage of each band transformed by MNF

从表1看出，经MNF变换后前3个分量共占了图像信息量的91.01%，集中了图像的主要信息;而最后2个分量总共只占图像信息量的3.86%，已表现出明显的噪声特征。同时，经MNF变换后的各分量之间没有相关性，消除了噪声。因此，采用前3个分量选择端元组分［18－19］。通过制作MNF前3个分量的两两波段的散点图(图2)，结合TM图像分析，可以在散点图上大致判别出各个端元组分在特征空间中的位置和分布情况。

图2 各端元组分在MNF1，2组合散点图中的位置Fig.2 Position of each end －member in the scattering chart of the MNF1，2

2.2.3 PPI计算

根据散点图可初步确定端元组分，但图2中像元多且分布复杂，人工勾画端元会遇到很大困难，因此需通过选择纯净像元来缩小端元选择的范围。PPI计算可用于在多光谱和高光谱图像中寻找波谱最“纯”的像元。PPI图像中像元值表示在某些映射下该像元被认为是纯净像元的次数，这些次数同时也显示出接近每个像元的数据集的局部凸率以及每个像元与数据凸面的接近程度。像元值较高(影像较亮)的像元接近N维数据集群凸出的犄角，因而比像元值较低的像元相对更纯净一些;相应地，像元值为零的像元肯定不是纯净像元［18－19］。通过给出适当的阈值可以提取出相对纯净的像元(图3)。

图3 像元纯度指数图像Fig.3 Image of the pixel purity index

PPI计算后，可选作端元的纯净像元的数量大大减少。统计结果显示，PPI＞500的像元占像元总数的比例仅为7.50%;因此，本文选取PPI=500作为阈值，这样可以在减少数据复杂度的同时选出相对纯净的像元作为端元。

2.2.4 端元组分确定与反射率计算

根据计算结果选择MNF变换后分量中的纯净像元，以PPI＞500的像元作为端元，通过ENVI软件中的N维可视化器来完成(图4)。

图4 N维可视化器选取的端元Fig.4 End －members selected by N dimensional visualizer

对各端元在各个分量上的光谱反射率计算其均值，可以得到各端元光谱反射率曲线。

2.3 线性混合像元分解

获得端元组分的类型及其光谱反射率信息后，即可代入选用的线性混合像元分解模型中，用带约束的最小二乘法求解。在使用最小二乘法时需满足所有端元组分的丰度大于0且总和等于1。本文使用ENVI中的IDL编程实现满足该条件的模型求解。最后得到的结果是各端元组分的丰度图，并以一系列灰阶图像的形式出现。各组分丰度图的像元值在(0，1］之间，表示在某个像元内某组分所占面积的百分比。

2.4 误差分析

线性混合像元分解误差主要来自3个方面:①端元像元的选择。由于端元像元来自图像本身，同时又采用人工方法选择端元像元，得到的端元光谱特征可能有误差。②阈值选择的误差。选择纯净像元时使用了人工输入的PPI阈值参数，有人为误差。③图像本身的误差。遥感图像在获取过程中存在辐射畸变或几何畸变。

式(5)计算出的均方根误差图像中负值和较大值对应着较大的误差，象征着错误的端元。结果表明94.54%的像元值误差为0＜RMS＜5，说明误差在允许范围之内［20］。

3 结果分析

3.1 冰雪覆盖信息提取

对1989年获取的TM和2000年获取的ETM+图像做如上处理后，可得到2个时期内各端元组分(冰雪、植被、裸露山体、阴影)的丰度图，从而可以提取冰雪覆盖丰度(图5)和统计冰雪面积(表2)。

图5 1989年和2000年冰雪覆盖丰度图Fig.5 Images of the snow cover abundance in 1989 and 2000

表2 1989年和2000年冰雪覆盖统计Tab.2 Statistics of the snow cover in 1989 and 2000

3.2 冰雪覆盖变化分析

通过将2000年与1989年冰雪覆盖图像作相减运算，得到冰雪覆盖丰度变化图(图6)。图中红色表示1989—2000年间冰雪减少的像元，绿色表示冰雪增加的像元。冰雪覆盖丰度变化统计见表3。

表3 冰雪覆盖丰度变化统计表Tab.3 Statistics of the snow cover abundance change

图6 冰雪覆盖丰度变化图Fig.6 Image of the snow cover abundance change

从图6、表2和表3中可以看出:

1)1989—2000年间有84.06%的像元失去了冰雪覆盖，占了全部像元的绝大多数;有34.19%的像元失去了90%以上的冰雪丰度，即失去了这些像元面积90%以上的冰雪覆盖;实验区冰雪覆盖的百分比从37.68%减少到28.46%，变化幅度达9.22%。产生这种现象的原因有:①气候变化是主要原因。在全球变暖的大背景条件下，天山冰川地区气候也有局部区域的响应现象，使得在1989—2000年间气温不断升高;气温升高导致冰川融化消退，冰雪覆盖量减少，这与天山1号冰川站长期实地观测的结果相符合［1］。②冰川消融是海拔较低地区的冰雪覆盖减少的重要原因，因为冰雪的正常融化往往是先从雪线以下的地区(即气温高于零度的地区)开始。从图6可知，冰雪覆盖减少区域出现在海拔较低的地区。

2)有15.94%的像元有了新的冰雪覆盖。产生这种现象的原因可能有:①冰川的移动。天山冰川每年移动5m，局部区域可能达数十米［3］，这使得各像元上冰雪面积各有增减，当有冰川移动进入某像元区域时其冰雪丰度就有增加。②山区局部区域的特殊地形所致。某些特殊地形如山谷沟涧，利于冰雪的堆积且不易流逝(如在山体的背阴面，冰雪堆积后不易融化，因此会有冰雪的汇积)。如图6右下方的绿色区域，表示该区域冰雪量有所增加;而结合TM假彩色合成图像可看出，该区域正是处于山谷位置且在背阴面，这可能使冰雪更易累积。③2个时相遥感图像的成像时间不同。天山地区于8月26日进入冰期［17］，而1989年图像成像于9月3日，2000年图像成像于9月17日，后者比前者晚了14 d;因此随着雪季的逼近，冰雪覆盖量会逐日增加也可能是原因之一。

3)此外，有56.01%的像元冰雪覆盖变化百分比在0.1之内，说明实验区内有半数以上区域的冰雪覆盖量没有太大变化。这表明天山冰川地区气候变化不是非常剧烈，冰雪总量基本上相对稳定。

4 结论

本文以天山典型冰川区为实验区，通过线性混合像元分解实现了对冰川地区冰雪覆盖变化的遥感监测。主要结论如下:

1)利用混合像元分解技术和线性模型成功提取了冰川地区的冰雪覆盖信息，证明该方法用于冰川地区冰雪覆盖变化遥感监测可以获得较好的效果。

2)提出适合于冰川地区的“冰雪－植被－裸露山体－阴影”端元模型，改善了端元组分的选取，为冰川地区提供了较好的混合像元分解模型。

3)进行了1989—2000年间天山典型冰川区内冰雪覆盖变化监测并探讨了冰雪覆盖变化的成因。

［1］上官冬辉，刘时银，丁永建，等.玉龙喀什河源区32年来冰川变化遥感监测［J］.地理学报，2004，59(6):855 －862.Shangguan D H，Liu S Y，Ding Y J，et al.Glacier Changes at the Head of Yurungkaz River in the West Kunlun Mountains in the Past 32 Years［J］.Acta Geographica Sinica，2004，59(6):855 －862(in Chinese with English Abstract).

［2］李忠勤，沈永平，王飞腾，等.冰川消融对气候变化的响应——以乌鲁木齐河源1号冰川为例［J］.冰川冻土，2007，29(3):334－342.Li Z Q，Shen Y P，Wang F T，et al.Response of Glacier Melting to Climate Change——Take Rümqi Glacier No.1 as an Example［J］.Journal of Glaciology and Geocryology，2007，29(3):334 －342(in Chinese with English Abstract).

［3］焦克勤，井哲帆，韩添丁，等.42a来天山乌鲁木齐河源1号冰川变化及趋势预测［J］.冰川冻土，2004，26(3):254 －260.Jiao K Q，Jing Z F，Han T D，et al.Variation of the Glacier No.1 at the Headwaters of the Rümqi River in the Tianshan Mountains during the Past 42 Years and Its Trend Prediction［J］.Journal of Glaciology and Geocryology，2004，26(3):254 － 260(in Chinese with English Abstract).

［4］Gebbinck M S.Decomposition of Mixed Pixels in Remote Sensing Images to Improve the Area Estimation of Agricultural Fields［D］.Nijmegen:Katholieke Universiteit Nijmegen，1998.

［5］吕长春，王忠武，钱少猛.混合像元分解模型综述［J］.遥感信息，2003(3):55－60.Lv C C，Wang Z W，Qian S M.A Review of Pixel Unmixing Models［J］.Remote Sensing Information，2003(3):55 － 60(in Chinese with English Abstract).

［6］惠巍巍，衣德萍，廖彩霞，等.混合像元分解研究综述［J］.林业科技情报，2007，39(1):2 －3.Hui W W，Yi D P，Liao C X，et al.The Study of Decomposing Mix Element［J］.Forestry Science and Technology Information，2007，39(1):2－3(in Chinese with English Abstract).

［7］邹 蒲.遥感混合像元分解技术在土地覆盖与生态环境变化中应用研究［D］.广州:中国科学院研究生院(广州地球化学研究所)，2007:31－46.Zhou P.Unpixel Technique of Remote Sensing Image and Its Applications in Land Cover and Ecological Environmental Changes:A Case Study in Huadu，Guangzhou［D］.Guangzhou:Graduate School of Chinese Academy of Sciences(Guangzhou Institute of Geochemistry)，2007:31 －46(in Chinese with English Abstract).

［8］岳文泽，吴次芳.基于混合光谱分解的城市不透水面分布估算［J］.遥感学报，2007，11(6):914 －922.Yue W Z，Wu C F.Urban Impervious Surface Distribution Estimation by Spectral Mixture Analysis［J］.Journal of Remote Sensing，2007，11(6):914 －922(in Chinese with English Abstract).

［9］李 霞，王 飞，徐德斌，等.基于混合像元分解提取大豆种植面积的应用探讨［J］.农业工程学报，2008，24(1):213 －217.Li X，Wang F，Xu D B，et al.Application Research on the Method for Extracting Soybean Covered Areas Based on the Pixel Unmixing［J］.Transactions of the Chinese Society of Agricultural Engineering，2008，24(1):213 －217(in Chinese with English Abstract).

［10］Plaza A，Martinez P，Perez R，et al.A New Approach to Mixed Pixel Classification of Hyperspectral Imagery Based on Extended Morphological Profiles［J］.Pattern Recognition，2004，37(2004):1097－1116.

［11］李 素，李文正，周建军，等.遥感影像混合像元分解中的端元选择方法综述［J］.地理与地理信息科学，2007，23(5):35－42.Li S，Li W Z，Zhou J J，et al.A Review on Endmember Selection Methods in the Course of Mixed pixel Decomposition of Remote Sensing Images［J］.Geography and Geo － Information Science，2007，23(5):35 －42(in Chinese with English Abstract).

［12］Quarmby N A，Townshend J R，Settle J J，et al.Linear Mixture Modeling Applied to AVHRR Data for Crop Area Estimation［J］.International Journal of Remote Sensing，1992，13(3):415 －425.

［13］Settle J J，Drake N A.Linear Mixing and the Estimation of Ground Cover Proportions［J］.International Journal of Remote Sensing，1993，14(6):1159 －1177.

［14］Kerdiles H，Grondona M O.NOAA － AVHRR NDVI Decomposition and Subpixel Classification Using Linear Mixing in the Argentinean Pampa［J］.International Journal of Remote Sensing，1995，16(7):1303－1325.

［15］Klein Gebbinck M S，Schouten T E，Settle J J.Decomposition of Mixed Pixels Based on a Physical Linear Mixture Model［J］.Image and Signal Processing for Remote Sensing II，1995，2579:104 －115.

［16］Cross A M，Settle J J，Drark N A，et al.Subpixel Measurement of Tropical Forest Cover Using AVHRR Data［J］.International Journal of Remote Sensing，1991，12(5):1119 －1129.

［17］胡汝骥.中国天山自然地理［M］.北京:中国环境科学出版社，2004:155－180.Hu R J.China’s Tianshan Mountain Natural Geography［M］.Beijing:China Environmental Science Press，2004:155 － 180(in Chinese).

［18］李 素，李文正，周建军，等.ETM+影像亚像元级城市土地覆盖组分丰度提取——以南京市为例［J］.地理与地理信息科学，2008，24(2):17 －24.Li S，Li W Z，Zhou J J，et al.Urban Land Cover Composition Abundance Extraction of Sub－pixel Level Based on ETM+Image:A Case Study of Nanjing［J］.Geography and Geo － Information Science，2008，24(2):17 －24(in Chinese with English Abstract).

［19］张志新，邓孺孺，李 灏，等.基于混合像元分解的南方地区植被覆盖度遥感监测——以广州市为例［J］.国土资源遥感，2011(3):88－94.Zhang Z X，Deng R R，Li H，et al.Remote Sensing Monitoring of Vegetation Coverage in Southern China Based on Pixel Unmixing:A Case Study of Guangzhou City［J］.Remote Sensing for Land and Resources，2011(3):88 －94(in Chinese with English Abstract).

［20］郑有飞，范旻昊，张雪芬，等.基于MODIS遥感数据的混合像元分解技术研究和应用［J］.南京气象学院学报，2008，31(2):145－150.Zheng Y F，Fan W H，Zhang X F，et al.Pixel Unmixing Technology of MODIS Remote Sensing Data［J］.Journal of Nanjing Institute of Meteorology，2008，31(2):145 － 150(in Chinese with English Abstract).

Monitoring of Snow Cover Changes in Tianshan Mountains Based on Mixed Pixel Decomposition

JIN Xin，KE Chang－qing
(Department of Geographical Information Sciences，Nanjing University，Nanjing 210093，China)

Mixed pixels are abundant in medium － low resolution images，but the traditional methods for image classification could only assign pixels to one class，with the ignorance of the mixed pixels.To tackle this problem，the authors selected the typical area in Tianshan Glacier of Xinjiang as an experimental area.Based on the theory of mixed pixel decomposition and the principle of the linear model and taking into account the spectral characteristics of TM/ETM+image as well as the land cover characteristics of Tianshan area，the authors developed an endmember composition model suitable for the glacier area，i.e.，Snow － Vegetation － Rock － Shade model.After the appropriate end－members were selected and the reflectance values were substituted into the improved linear mixed pixel decomposition model，which satisfied the constraints，the abundance image of individual end － member was calculated and the snow cover information was easily and precisely extracted.The extraction results of snow cover in 1989 and 2000 demonstrate that the mixed pixel decomposition and the linear model could be used to monitor the snow cover changes in the glacier area.

snow cover;mixed pixel decomposition;linear mixture model;end－member selection;Tianshan Mountains

TP 79

A

1001－070X(2012)04－0146－06

2011－12－30;

2012－03－12

国家自然科学基金项目(编号:40971044)、国家科技支撑计划项目(编号:2012BAH28B02)、教育部新世纪优秀人才支持计划项目(编号:NCET－08－0276)和江苏高校优势学科建设工程项目共同资助。

10.6046/gtzyyg.2012.04.24

金 鑫(1989－)，男，硕士研究生，主要研究方向为遥感与地理信息系统。E－mail:jiangxi_jinxin@163.com。

柯长青(1969－)，男，教授，主要研究方向为遥感、GIS及其应用。E－mail:kecq@nju.edu.cn。

(责任编辑:邢 宇)