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

积雪是地球表面最为活跃的自然要素之一,其特征(如积雪面积、雪深、积雪反照率等) 是全球能量平衡模型中的重要输入参数 中分辨率成像光谱仪(MODIS) 是EOS 系统中主要的传感器之一,它具有较高的空间分 辨率、光谱分辨率和时间分辨率、虽然MODIS 全球每日500m 空间分辨率雪产品较以往雪盖制图精度有了一定的提高,但是混合像元的问题仍比较突出. 针对混合像元问题,一些学者提出了光谱线性分解模型,通过选取恰当的地表覆盖类型组合,就能在一定程度上提高最终结果的精度[8 ] . 但是这种模型运算量较大, 不适合积雪覆盖率的业务化制图. Salomoson[9 ] 等提出基于归一化雪盖指数(normalized difference snow index ,NDSI) 的MODIS 像元雪盖率(FRA) 反演简化方法、中国学者曹云刚[10 ] 考虑地表覆盖,建立了像元雪盖率与雪盖指数、植被指数之间的线性关系模型. 然而精确反演雪盖率依旧非常困难,其中一个重要原因在于积雪本身的物理性质也会对NDSI 产生影响;另一方面,地表覆盖的复杂性也会降低NDSI 与积雪比例间的相关性. 若对数据不加区分统一建模,则必然会对模型的精度产生 影响. 因此,本文尝试将雪盖率按NDSI 分类,然后对不同类别分别进行建模、选择研究区域需要在一定程度上,研究区具有很好的代表性 本文使用MODIS 1B 500m 分辨率数据反演雪盖率 (MODIS02HKM) (图1) ,其中1、2 波段为MODIS 250m 原 始影像经重采样到500m 后得到. 如表1 选择同一区域、 同一天ETM + 数据作为检验数据(数据已经过投影转 换,为UTM投影,空间分辨率为30m) . 数据预处理 数据预处理包括波段反射率计算以及投影转换TM由于ETM + 数据已经过投影转换,故只需计算其各波段的反射率值. 以ENVI413 为例,使用Calibration Utilites 中的Landsat TM功能计算其反射率MODIS对于MODIS 1B 数据,首先使用MODIS tools 的bow2tie 工具对图像去“蝴蝶结”效应,接着使用MODIS 1B 影像自带的控制点信息对图像进行几何校正. 然后根据公式(1) 计算各波段的反射率值: ρ = reflectance- scales ×(DN2reflectance- offsets) , (1) 式中,ρ为反射率;DN 为像元值;reflectance - scales 为各波段缩放系数;reflectance - offsets 为偏置系数.得到波段反射率以后,从MODIS 图像中裁剪出包含ETM+ 数据范围的区域,并将其转换为UTM 投影. 在转换过程中为保持波段反射率不变,故使用最邻近插值法. 为了与Landsat ETM + 数据进行对比分析,还需要对裁剪后的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 为Landsat ETM+ 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+ 在空间分辨率上的较大差异,也会导