基于MODIS影像的亚象元积雪覆盖率提取方法研究

积雪是地球表面最为活跃的自然要素之一,其特征(如积雪面积、雪深、积雪反照率等)是全球能量平衡模型中的重要输入参数中分辨率成像光谱仪(MODIS)是EOS系统中主要的传感器之一， 它具有较高的空间分辨率、光谱分辨率和时间分辨率、虽然MODIS全球每日500m空间分辨率雪产品较以往雪盖制图精度有了一定的提高 ,但是混合像元的问题仍比较突出.针对混合像元问题,一些学者提出了光谱线性分解模型通过选取恰当的地表覆盖类型组合,就能在一定程度上提高最终结果的精度[8].但是这种模型运算量较大,不适合积雪盖率的业务化制图.Salomoson[9]等提出基于归一化雪盖指数(normalizeddifferencesnowindex,NDSI)的MODIS像元雪盖率(FRA)反演简化方法、中国学者曹云刚[10]考虑地表覆盖,建立了像元雪盖率与雪盖指数、植被指数之间的线性关系模型.然而精确反演雪盖率依旧非常困难其中一个重要原因在于积雪本身的物理性质也会对NDSI产生影响另一方面地表覆盖的复杂性也会降低NDSI与积雪比例间的相关性.若对数据不加区分统一建模,则必然会对模型的精度产生影响.因此,本文尝试将雪盖率按NDSI分类,然后对不同类别分别进行建模、选择研究区域需要在一定程度上,研究区具有很好的代表性本文使用MODIS1B500m分辨率数据反演雪盖率(1),1、2波段为MODIS250m原500m后得到.1同一天ETM+数据作为检验数据(数据已经过投影转换UTM投影,30m).表1“luhat/f：™+和MODIS影像数据信息校感器类型LandsaL.ETM+

成像时间2001/10?09

穫盖范围那，4197'38.4779N

便用波段沁in 2,4t5EOS'MDDK 2001/10/09数据预处理

100-4S13—103-3857E35-49715514N994468 4067E

500m 1,2,4?5.6数据预处理包括波段反射率计算以及投影转换 TM由于ETM+数据已经过投影转换故只需计算其各波段的反射率值 .以ENVI413为例使用CalibratenUtilites中的LandsatTM功能计MODISMODIS1B数据,MODIStoolsbow2tie工具对图像去“蝴蝶结”效应，接着使用MODIS1B影像自带的控制点信息对图像进行几何校正(1)计算各波段的反射率值:

.然后根据公式p=reflectanee-scalesx(DN2reflectanee-offsets),(1)式中,p射率;DN为像元值;refleetanee-scaleseflectanee-offsets为偏置系数.得到波段反射率以后，从MODIS图ETM+数据范围的区域，并将其转换为UTM投影.转换过程中为保持波段反射率不变，故使用最邻近插值法.为了与LandsatETM+数据进行对比分析，还需要对裁剪后的MODIS影像进行空间配准.MODIS1、4、3ETM4、3、2合成影像的对比分析选择控制点MODIS影像进行配准得到相应的MODIS图像后，分别计算雪盖指数(NDSI,(p4-p6)P(p4+p6))和植被指数(NDVI,(p2-p1)P(p2+p1)).雪盖率计算ETM+数据，SNOWMAP方法［4］进行分类.SNOWMAPNDSI,NDSI的阈值为0140NDSI0140,p4>0111时判定为雪像元,p4为LandsatETM+4波段的反射率,1代表积雪像元,0代表非积雪像元，30m;然后3个研究区分MODISMODIS率MatlabMODISETM+积雪分类图；MODIS像元对应的区域中所包含的ETM+1的积雪像元个数MODIS像元雪盖率(fraction).fraction=丛?"其中，nsMODISETM+二值积雪分类图中积雪像元的个数；nMODISETM类图的像元个数.误差分析(1)位于中间部分的NDSI值相对较稀疏，从而导致模型拟合结果受到影响；(2)从空间分布来看这些数据主要集中在高积雪区与非积雪区的过渡地带而这些地区地形起伏较大所以山体的阴影也会对模型反演的精度产生影响;(3)在过渡地区地表覆盖复杂仅仅使用NDVI并不能完全反映地表的影响;(4)MODIS数据和ETM+差异也会导致雪盖率的计算产生误差所以这些区域的雪盖率值波动较大。------------------------------MODIS♘像元积雪覆盖率提取方法1.pdf针对目前应用广泛的MODIS传感器数据,充分利用了雪盖指数在积雪监测中的重要性,并在考虑了地表覆盖的情况下,建立了像元雪盖率与雪盖指数、植被指数之间的线性关系模型,并利用ETM+据对模型估算的雪盖率进行了验证.结果表明,该方法能有效地提取♘像元尺度的积雪信息1 引言积雪是地表覆盖的重要组成部分,积雪表面的高反射特性使其成为决定地球辐射平衡的一个关键因子,也使其成为区域与全球气候变化研究的重要内容之一.同时,积雪的空间分布还是山区与季节性积雪区水文、气象模型的重要输入因子[1]因此,研究积雪分布及其特性,对于水资源利用、灾害分析、大气环流分析以及气候演变有着重要的意义[2近年来,遥感资料已经成为积雪监测的主要信息来源,其中最主要的有AVHRR(AdvancedVeryHighResolutionRadiometer)、MODIS(ModerateResolutionImagingSpectroradiometer)、TM/ETM+(ThematicMapper/EnhancedThematicMapper)等.针对这些数据源,多种方法已经用于积雪识别与雪盖制图,包括:目视判读,多光谱图像分类,反射率特性计算法,决策树等[3-4].这些方法的一个共同特点就是,将遥感图像分为两个部分,即积雪区与非积雪区.也就是说,在进行积雪像元识别时,都假设所有的像元是纯像元,一个像元对应一种地物类型.因此,当遥感图像的空间分辨率较低时,混合像元的问题就比较突出,造成分类结果的精度不够理想[5].针对混合像元问题,一些学者提出了光谱线性分解模型,通过选取恰当的地表覆盖类型组合,就能在一定程度上提高最终结果的精度[6]光谱线性分解模型在精度上有所提高,但是它通常使用最小二乘法来研究区域与数据准2估算每个像元中的组分比例,当遥感数据量较大时,就会耗费大量的计算时间,并且这种方法的自动化程度也很低,很难满足一些业务上的需要.针对这样的问题,本文引入了两个基本输入参数（雪盖指数与植被指数）,提出了一种简化的算法来计算单个像元中积雪所占的比例（即像元雪盖率,snowfraction）.采用研究区内的MODIS与ETM+数据,对算法中的参数进行了计算,最后对算法的结果进行了比较和分析 研究区域与数据准2域是高原主要的积雪区之一,且地形比较复杂,植被覆盖类型多样,4ETM+影像数据1（UTM投影30m）.同时选择了包含以上区域两个时期内的MODIS1B500m（MODIS1、2、、61、2250m,4、6波段分辨率为500m,500m）.1相关ETM+数据信息Table1InformalionoftheETM+■dala黑域黑域成像日期轨道号太阳方位角12<1)]1021pl36rO3&15],23M7太阳高度角43.U44歹2(1)11021P136rO39150.225744.178^31023149.714(f44.70H2斗2(1)11U23pl34rO4O15],7295142.435ff2KTVI+各波段耳值Table2ValtHLofdifferentwavebantfeETM+波段号/<i】/(*1叶*匸]yfw*315514]OU5225.77K2,073 数据预处理数据预处理包括波段反射率计算以及投影转换.对于本文所使用的ETM+数据，因其已经过投影转换，故只需要计算其各波段的反射率值，使用下面的公式将其像元灰度值转换为反射率(I)£(>xoosO式中:a、b为定标常数;L为辐射亮度；d为日地距离；E0为大气顶层太阳辐照度；！太阳天顶角.MODIS1B数据，MODIStoolsbow-tie工具对图像的边缘重叠效应进行改正,接着使用MODIS1BGeoreferenee命令对图像进行了几何校正（ENVI软件为例）.然后，使用下面的公式将MODIS各波段的像元灰度值转换为反射率值[7]:P=reflectan x（Z）A—refleclance_offsets）式中: 为反射率；DN为像元值各波段比例系数（refleetance_scales） 与偏置系数（refleetance_offsets）见表3（只列出了文中使用的几个波段参数表3MODIS 据部分波段定标参数值）值）.60+000()3238{)波段号比例系数倾直系数1(1(XW51370.00X)31970()4{)在得到波段反射率以后，MODISETM+数据范围的区域UTM投影.MODIS图像，使用下面的公式计算其雪盖指数（NDSI）与植被指数（NDVIAOS/二NDVI=式中:PMODISi波段反射率.4 研究方法

(51本文以ETM+数据获取的雪盖率作为#真值？，(51MODISNDSI、NDVI之间的线性关系.，雪盖率真值？计算需要经过以下几个步骤首先通过监督分类的ETM+数据进行积雪信息的提取，MODISETM+1的像元数（积雪像元）；使用式（6）来计算MODIS像元雪盖率raction）:fraction= 积？;像元数像元总数

(6/2的散点图可以看出雪盖率存在于两个主要的区间，即高值区（0.8）和低值区（0.1）,其中高值区表明了积雪是一种大面积比较连续的地表覆盖,而低值区的存在除了研究区内积雪较少的原因以外，还在一定程度上受计算误差的影响（其中包括MODIS与ETM+数据之间的坐标配准误差）.针对NDSI与雪盖率之间的关系，Bartonetal.[8]Salomonsonetal.[9]（fraction=0.18+0.37 ?NDSI+0.26 ?（NDSI）2）和直线（fraction=0.06+1.21 ?NDSI）其中fraction=a+bDSI（A）NDSI=a+b?fraction（B）,2的3种算法，AB.2其在低值区的同样A,其在高值B,则能得到造成以上误差的影响NDSI其中积雪颗（NDSI值）;对于积雪的低值区，由于地表覆盖的复杂性NDSI与积雪比例之间的相关性降低.通过以上的分析本文认为只使用单因子来进行雪盖率的反演还存在一定的不足因此NDVI值来进的方法来计算雪盖率:fraction=</+/>xNDS1+rxADVI(7)12得到的积雪信息作为雪盖率的#真值？，通过线性回归的方法，本文对上式中的3个参数进行了拟合计算为了不增加计算的复杂性最后的参数只是对上面的各区域拟合参数进行了平均 (对Salomonson提出的算法，B作变形，NDSI示雪盖率A的参数进行平均X5 使用各区域算法平均后得到的最终结果Table5Th<^(ruitioiiKforiliffc^rrntalgoritLms最后拟合结果Baiton fradinti-0 24+it63xADS/+0.31+SalcmonAon

fraetinn=0.1K+0T

97xAI)Si帛文算法f(X+(K1x\-Q«A5 算法结果验证及其精度分析4MODIS日积雪图与本文算法计算的雪盖率图其中，MODIS积雪图使用像元值200来表示积雪像元，它只是在像元层次上给出了区域的积雪分布情况而使用雪盖率的方法可以得到像元积雪所占的比例，它能清楚地表现积雪分布的强弱情况，尤其是在积雪边缘区域雪盖率更能反应积雪覆盖渐变的细节信息. 针对图4所示的区域，使用本文算法计算了雪盖率并对结果进行了统计和分析，5从图可以看出，研究区域内雪盖率估算误差分布的主要区间集中在[-0.25,0.1]. 研究区积雪覆盖率为0.598,由算法计算的像元雪盖率统计得到的积雪覆盖率为 0.509,14.8%.NDSI（6）NDSI0.6率.造成这种结果的原因可能(h)无驾积雪 0 0.5 1.04MODISH积雪图(a)与雪盖率图(b)Pig.4