基于MODIS影像的亚象元积雪覆盖率提取方法研究Word格式文档下载.docx

基于MODIS影像的亚象元积雪覆盖率提取方法研究Word格式文档下载.docx

《基于MODIS影像的亚象元积雪覆盖率提取方法研究Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《基于MODIS影像的亚象元积雪覆盖率提取方法研究Word格式文档下载.docx（10页珍藏版）》请在冰豆网上搜索。

数据预处理

数据预处理包括波段反射率计算以及投影转换TM由于ETM+数据已经过投影转换,故只需计算其各波段的反射率值.以ENVI413为例,使用CalibrationUtilites中的LandsatTM功能计算其反射率MODIS对于MODIS1B数据,首先使用MODIStools的bow2tie工具对图像去“蝴蝶结”效应,接着使用MODIS

1B影像自带的控制点信息对图像进行几何校正.然后根据公式

（1）计算各波段的反射率值:

ρ=reflectance-scales×

（DN2reflectance-offsets）,

（1）

式中,ρ为反射率;

DN为像元值;

reflectance-scales为各波段缩放系数;

reflectance-offsets为偏置系数.得到波段反射率以后,从MODIS图像中裁剪出包含ETM+数据范围的区域,并将其转换为UTM投影.在转换过程中为保持波段反射率不变,故使用最邻近插值法.为了与LandsatETM+数据进行对比分析,还需要对裁剪后的MODIS影像进行空间配准.通过MODIS1、4、3波段合成影像和ETM+4、3、2合成影像的对比分析选择控制点,并对MODIS影像进行配准.得到相应的MODIS图像后,分别计算雪盖指数（NDSI,（ρ4-ρ6）P（ρ4+ρ6））和植被指数（NDVI,（ρ2-ρ1）P（ρ2+ρ1））.

雪盖率计算

对于ETM+数据,首先利用SNOWMAP方法[4]进行分类.SNOWMAP算法的核心内容是NDSI,首先设定NDSI的阈值为0140,当像元NDSI≥0140,且ρ4>

0111时判定为雪像元.其中,ρ4为LandsatETM+4波段的反射率.生成的二值积雪分类图中,1代表积雪像元,0代表非积雪像元,像元大小为30m;

然后,

选择3个研究区,分别从分类图像和MODIS影像中提取对应区域.为了计算MODIS雪盖率,首先在Matlab中分别导入MODIS图像和ETM+积雪分类图;

然后统计每个MODIS像元对应的区域中所包含的ETM+分类像元总数及其中值为1的积雪像元个数,按公式

（2）计算MODIS像元雪盖率（fraction）.

其中,ns为MODIS像元内的ETM+二值积雪分类图中积雪像元的个数;

n为MODIS像元内ETM+二值

积雪分类图的像元个数.

误差分析

（1）位于中间部分的

NDSI值相对较稀疏,从而导致模型拟合结果受到影响;

（2）从

空间分布来看,这些数据主要集中在高积雪区与非积雪区的

过渡地带,而这些地区地形起伏较大,所以山体的阴影也会对

模型反演的精度产生影响;

（3）在过渡地区地表覆盖复杂,仅

仅使用NDVI并不能完全反映地表的影响;

（4）MODIS数据和

ETM+在空间分辨率上的较大差异,也会导致雪盖率的计算产

生误差,所以这些区域的雪盖率值波动较大。

-------------------------------MODIS亚像元积雪覆盖率提取方法1.pdf

针对目前应用广泛的MODIS传感

器数据,充分利用了雪盖指数在积雪监测中的重要性,并在考虑了地表覆盖的情况下,建立了像元雪

盖率与雪盖指数、植被指数之间的线性关系模型,并利用ETM+数据对模型估算的雪盖率进行了验

证.结果表明,该方法能有效地提取亚像元尺度的积雪信息

.

1引言

􀀁

􀀁

积雪是地表覆盖的重要组成部分,积雪表面的高反射特性使其成为决定地球辐射平衡的一个关键因子,也使其成为区域与全球气候变化研究的重要内容之一.同时,积雪的空间分布还是山区与季节性积雪区水文、气象模型的重要输入因子[1].因此,研究积雪分布及其特性,对于水资源利用、灾害分析、大气环流分析以及气候演变有着重要的意义

[2

近年来,遥感资料已经成为积雪监测的主要信息来源,其中最主要的有AVHRR（AdvancedVer􀀁

yHighResolutionRadiometer）、MODIS（Moder􀀁

ateResolutionImagingSpectroradiometer）、TM/ETM+（ThematicMapper/EnhancedThematicMapper）等.针对这些数据源,多种方法已经用于积雪识别与雪盖制图,包括:

目视判读,多光谱图像分类,反射率特性计算法,决策树等[3-4].这些方法的一个共同特点就是,将遥感图像分为两个部分,即积雪区与非积雪区.也就是说,在进行积雪像元识别时,都假设所有的像元是纯像元,一个像元对应一种地物类型.因此,当遥感图像的空间分辨率较低时,混合像元的问题就比较突出,造成分类结果的精度不够理想[5].针对混合像元问题,一些学者提出了光谱线性分解模型,通过选取恰当的地表覆盖类型组合,就能在一定程度上提高最终结果的精度[6]

光谱线性分解模型在精度上有所提高,但是它通常使用最小二乘法来估算每个像元中的组分比例,当遥感数据量较大时,就会耗费大量的计算时间,并且这种方法的自动化程度也很低,很难满足一些业务上的需要.针对这样的问题,本文引入了两个基本输入参数（雪盖指数与植被指数）,提出了一种简化的算法来计算单个像元中积雪所占的比例（即像元雪盖率,snowfraction）.采用研究区内的MODIS与ETM+数据,对算法中的参数进行了计算,最后对算法的结果进行了比较和分析

2研究区域与数据准备

该区域是高原主要的积雪区之一,且地形比较复杂,植被覆盖类型多样,主要有高寒草甸、高山灌丛、裸岩、积雪等.从一定程度上来讲,该区域具有良好的代表性,对其研究将有助于分析遥感数据在整个青藏高原积雪信息提取中的实用性以及精度.

本文收集了4景ETM+影像数据,相应的参数如表1（数据已经过投影转换,为UTM投影,空间分辨率为30m）.同时,选择了包含以上区域两个时期内的MODIS1B500m分辨率影像数据（本文使用了MODIS1、2、4、6波段数据,其中1、2波段分辨率为250m,4、6波段分辨率为500m,为了统一精度,这里使用了500m分辨率数据）.

3􀀁

数据预处理

数据预处理包括波段反射率计算以及投影转换.对于本文所使用的ETM+数据,因其已经过投影转换,故只需要计算其各波段的反射率值,使用下面的公式将其像元灰度值转换为反射率值!

式中:

a、b为定标常数;

L为辐射亮度;

d为日地距离;

E0为大气顶层太阳辐照度;

!

太阳天顶角.对于MODIS1B数据,首先使用MODIStools的bow-tie工具对图像的边缘重叠效应进行改正,接着使用MODIS1BGeoreference命令对图像进行

了几何校正（以ENVI软件为例）.然后,使用下面的公式将MODIS各波段的像元灰度值转换为反射率值[7]:

􀀂

为反射率;

DN为像元值;

各波段比例系

数（reflectance_scales）与偏置系数（reflectance_off􀀁

sets）见表3（只列出了文中使用的几个波段参数

值）.

在得到波段反射率以后,从MODIS图像中裁减出包含ETM+数据范围的区域,并将其转换为UTM投影.对应MODIS图像,使用下面的公式计算其雪盖指数（NDSI）与植被指数（NDVI）.

4􀀁

研究方法

本文以ETM+数据获取的雪盖率作为#真值∃,

建立起MODIS数据雪盖率与NDSI、NDVI之间的线性关系.其中,雪盖率#真值∃计算需要经过以下几个步骤:

首先,通过监督分类的方法对ETM+数据进行积雪信息的提取,并将积雪像元的值设定为1;

统计各MODIS像元所包含的ETM+像元

总数以及像元值为1的像元数（积雪像元）;

使用式（6）来计算MODIS像元雪盖率（fraction）:

从图2的散点图可以看出,雪盖率存在于两个主要的区间,即高值区（>

0.8）和低值区（<

0.1）,其中高值区表明了积雪是一种大面积比较连续的地表覆盖,而低值区的存在除了研究区内积雪较少的原因以外,还在一定程度上受计算误差的影响（其中包括MODIS与ETM+数据之间的坐标配准误

差）.针对NDSI与雪盖率之间的关系,Bartonet

al.[8]和Salomonsonetal.[9]先后给出了曲线

（fraction=0.18+0.37∀NDSI+0.26∀

（NDSI）2）和直线（fraction=0.06+1.21∀NDSI）

拟合模型.其中,直线模型有两种建立方法,分别为fraction=a+b∀NDSI（算法A）或者NDSI=a+b∀fraction（算法B）,图2的曲线和直线分别表示了这3种算法,较细的直线表示算法A,较粗的直线表示算法B.从图2还可以看出,对于曲线模型,其在低值区的拟合效果不佳,将高估该区域的积雪值;

同样,对于算法A,其在高值区低估了积雪值,在低值区却高估了积雪值.对于算法B,则能得到比较好的结果,但是它也低估了高值区的积雪系数.造成以上误差的主要原因,分析认为:

对于高值区,在积雪覆盖超过一定值后,影响NDSI

的主要因素不再是积雪的比例,而是积雪的性质,其中积雪颗粒的影响较为严重（积雪颗粒大小影响了积雪在红外光谱区间内的反射率,较细的积雪具有较高的反射率,这样就会得到较低的NDSI值）;

对于积雪的低值区,由于地表覆盖的复杂性,使得NDSI与积雪比例之间的相关性降低.通过以上的分析,本文认为只使用单因子来进行雪盖率的反演还存在一定的不足,因此,可以结合NDVI值来进行雪盖率的反演.在考虑了地表覆盖对雪盖率的影响以后,使用下面的方法来计算雪盖率:

使用区域1和2得到的积雪信息作为雪盖率的#真值∃,通过线性回归的方法,本文对上式中的3个参数进行了拟合计算为了不增加计算的复杂性,最后的参数只是对上面的各区域拟合参数进行了平均（对于Salomon􀀁

son提出的算法,将算法B作变形,使用NDSI值

来表示雪盖率,并将其与算法A的参数进行平均

5􀀁

算法结果验证及其精度分析

图4为MODIS日积雪图与本文算法计算的雪盖率图.其中,MODIS积雪图使用像元值200来表示积雪像元,它只是在像元层次上给出了区域的积雪分布情况,而使用雪盖率的方法可以得到像元内积雪所占的比例,它能清楚地表现积雪分布的强

弱情况,尤其是在积雪边缘区域,雪盖率更能反应积雪覆盖渐变的细节信息.􀀁

针对图4所示的区域,使用本文算法计算了雪盖率,并对结果进行了统计和分析,相关的误差分布见图5

从图5可以看出,研究区域内雪盖率估算误差

分布的主要区间集中在[-0.25,0.1].研究区积雪覆盖率为0.598,由算法计算的像元雪盖率统计得到的积雪覆盖率为0.509,则区域内积雪覆盖率的相对误差为14.8%.从算法误差与NDSI之间的散点图（图6）可以看出,对于雪盖率的估计,明显地存在一个低值区,也就是在NDSI>

0.6时,算

法低估了像元的雪盖率.造成这种结果的原因可能

是由于山谷内的地表覆盖情况复杂,地形变化明显,尤其是植被（森林）覆盖率较高,同时还受一定的山体阴影影响,这样采用相同参数的模型必然会得到差异较大的结果.通过分析,可以得出以下规律:

1）误差分布与积雪状况关系密切.对于几乎完全覆盖的积雪区域,也就是NDSI的高值区,所有的算法总是低估了积雪的比例;

同时积雪物理性质的差异还会使算法的通用性降低,对于不同的区域,算法中各参数的拟合值稍微有所不同;

2）误差的大小与地表覆盖有一定的关系.不同的地表覆盖情况影响着算法中各参数的值,而且地表覆盖的复杂情况也将降低NDSI与雪盖率之间的相关性.

6􀀁

结论与讨论

本文在研究了雪盖率与雪盖指数、植被指数之间的关系后,提出了一种简单的方法来估算MO􀀁

DIS像元雪盖率,并对估算结果进行了评价和分析.本文算法在积雪覆盖的两种极端情况都很难给出准确的估算值,即对于完全积雪区估值偏低,而对于稀疏积雪区,算法结果误差波动较大.造成这种结果的主要原因可能是:

在完全积雪区,当ND􀀁

SI>

0.8时,它的变化主要是由积雪物理特性引起的,这时NDSI值与雪盖率之间的相关性减弱,导致雪盖率估算误差增大;

对于稀疏的积雪区域,误

差主要是由于地表覆盖的多样性引起的.从本文的结果来看,只是单纯地使用NDVI作为这种多样性的一个指标还不够充分.不过,NDSI与雪盖率之间的相关性还是比较明显,在以后的工作中还可以考虑使用分段函数的方法来计算雪盖率,以提高

估算精度.总的来讲,使用本文提出的雪盖率估算方法还是能够得到较高精度的结果,这对于大面积的积雪实时监测具有一定的实用性.