祁连冰沟流域浅雪层光谱特征分析与遥感算法反演
梁慧, 黄晓东, 王云龙, 高金龙, 马晓芳, 梁天刚
草地农业生态系统国家重点实验室 兰州大学草地农业科技学院,甘肃 兰州 730020
第一作者:梁慧(1991-),女,甘肃环县人,在读硕士生,研究方向为草地遥感与地理信息系统E-mail:liangh16@lzu.edu.cn
摘要

祁连山冰沟流域地形复杂,积雪深度较浅且破碎化严重,针对MODIS标准积雪面积比例产品在该地区监测精度较差的问题,本研究基于冰沟流域浅雪光谱特征分析及结合野外实测经验,探索浅雪的光谱特征对MODIS浅雪面积比例提取精度的影响;然后通过线性回归法、线性混合像元分解法及BP神经网络模型3种方法分别构建了针对研究区的MODIS积雪制图算法,并利用同时相的Landsat 8 OLI二值积雪数据作为真值对上述3种制图方法进行精度验证。结果表明,1)浅雪的光谱反射率对基于NDSI阈值法的MODIS浅雪提取精度几乎没影响,MODIS提取浅雪精度差的主要原因为该地区复杂的地形而导致的积雪分布破碎化,即混合像元的大量存在;2)利用BP神经网络模型反演积雪面积比例的最佳输入参数组合为(ρ17)+NDSI+DEM;3)线性混合像元分解模型在该研究区的积雪面积比例提取精度较低,BP神经网络模型精度最好;4)在地形复杂区域,多因素模型(BP神经网络模型)相对于单因素模型(一元线性回归模型)具有较好的积雪面积比例提取精度和稳定性,可以为研究区MODIS积雪面积比例的反演提供一种理想的方法。

关键词: 光谱反射率; 积雪面积比例; MODIS; Landsat 8 OLI; 积雪制图; 线性回归; BP神经网络
中图分类号:S127 文献标志码:A 文章编号:1001-0629(2017)07-1353-12 doi: 10.11829/j.issn.1001-0629.2017-0166
Analysis of thin snow spectral characteristic and retrieval algorithm construction of the fractional snow cover in Qilian Binggou Basin
Liang Hui, Huang Xiao-dong, Wang Yun-long, Gao Jin-long, Ma Xiao-fang, Liang Tian-gang
State Key Laboratory of Grassland and Agro-ecosystems, College of Pastoral Agriculture Science and Technology, Lanzhou University, Lanzhou 730020, China
Corresponding author: Huang Xiao-dong E-mail: huangxd@lzu.edu.cn
Abstract

The terrain in Qilian Binggou Basin is relatively complex, the snow depth is generally thin and snow distribution is fragmentized. Aiming at the problem of poor monitoring accuracy of MODIS snow products in this area, This study explored the influence of spectral characteristic of thin snow on the retrieval accuracy of MODIS data based on thin snow spectral characteristic analysis and combined with field survey experience; Then three MODIS fractional snow cover retrieval models are constructed through linear regression, linear mixed pixel unmixing and artificial neural network, and the snow map retrieved from Landsat 8 OLI image is taken as the ground truth to validate the three models’ accuracy respectively. The results show that: 1) The spectral reflectance of thin snow almost has no effect on MODIS snow retrieval accuracy based on the NDSI threshold method in the area. The poor accuracy of MODIS retrieval from thin snow is mainly due to the fragmentation of snow distribution caused by the complex terrain in the area, that is, the existence of a quantity of mixed pixels. 2) The best input parameters combination of BP artificial neural network model for MODIS fractional snow cover retrieval is(ρ17)+NDSI+DEM. 3) The linear mixed pixel unmixing model has the lowest accuracy and the BP artificial neural network model has the best accuracy in terms of snow cover extraction in the study area. 4) Multi-factor model (BP artificial neural network) has better accuracy and stability of snow coverage extraction compared with the single factor model (unary linear regression model) in complex terrain, it is an ideal method for the retrieval of fractional snow cover in the study area.

Keyword: spectral reflectance; fractional snow cover; MODIS; Landsat 8 OLI; Snow Mapping; linear regression model; BP artificial neural network

积雪是地表最活跃的自然因素之一, 其特征如积雪分布、积雪面积、雪深等是全球能量平衡、气候、水文以及生态模型中的重要输入参数[1]。山区积雪不仅是气候、水文、生态环境等的重要影响因素, 而且积雪形成的径流也是中下游地区重要的淡水资源, 可以及时补给河流与地下水。祁连山区地形复杂, 积雪深度较浅, 每年大量的积雪融水是河西走廊绿洲农业灌溉所需水资源的重要保障, 是河西走廊经济建设的基础[2]。研究表明, 在河西走廊农业、工业和生活用水有50%来自于祁连山区的融雪[3]; 季节性积雪变化不仅对青藏高原地区高寒草地植物的多样性发育与维持、群落结构与功能产生直接影响[4, 5], 而且也将对地球所有生态系统都产生潜在影响[6]。因此, 及时了解祁连山积雪资源变化状况, 准确监测该地区积雪变化, 制作高精度的雪盖图, 完善积雪观测系统, 对河西走廊干旱地区显得尤为重要, 而且对水文、气候、能量平衡以及建立准确和长期的积雪产品数据库都有重要意义。目前, 卫星遥感是进行大尺度积雪实时监测与信息提取的唯一技术手段[7]。自20世纪60年代初期, 卫星遥感技术以其宏观、快速、周期性、多尺度、多谱段、多时相等优势[8], 在积雪动态监测中发挥着重要作用[9]。在积雪面积比例反演方面得到广泛应用的光学传感器主要有Landsat TM/ETM+/OLI、NOAA/AVHRR、EOS/MODIS等, 其中新一代中分辨率辐射扫描仪MODIS(moderate resolution imaging spectroradiometer)因共享性和时间序列的完整性而成为大区域积雪监测研究所广泛使用的数据源[10, 11, 12]。但MODIS数据最高空间分辨率为250 m, 对复杂地形下的破碎化山区积雪监测来说, 数据中混合像元现象普遍存在[13]

目前, 国内外的遥感积雪制图主要以传统的积雪二值分类图为主, 即影像上的所有像元都被当成纯净像元, 有雪或无雪。积雪像元识别的方法主要有亮度阈值法[14]、监督分类法以及归一化雪被指数NDSI(normalized difference of snow index)法[15], 在这3种方法中, NDSI是最常用到的方法。基于Landsat TM影像的SNOMAP方法[16]提取的积雪图精度很高, 原因在于该方法的核心是NDSI。SNOMAP方法作为MODIS全球积雪产品的标准算法, 将NDSI的阈值设置为40%, 提取的积雪面积比例精度较高。但验证区在北美和欧洲, 其在混合像元等问题上仍然存在限制, 对祁连山区小于4 cm的浅雪层几乎没有识别能力, 造成MODIS积雪产品对祁连山区积雪提取精度很差, 很难满足复杂地形下积雪面积监测精度日益增高的需求, 很难为水文、气候等模型提供更加准确的输入因子。所以, 找出该地区浅雪提取精度差的主要原因、提高其整体积雪识别精度已显得尤为必要。

鉴于此, 本研究以祁连山冰沟流域为研究区, 针对NDSI阈值设置的原理, 分析浅雪的光谱特征, 重点探索:1)浅雪光谱特征对MODIS浅雪面积比例提取精度的影响, 从而为积雪面积比例模型的构建提供一定的依据; 2)在混合像元存在的事实下, 线性回归模型、线性混合像元分解模型以及BP(error back propagtion)神经网络模型对研究区积雪面积比例提取的精度和稳定性; 3)BP神经网络模型在研究区积雪面积比例反演方面的适用性, 旨在发展面向该研究区的积雪制图算法, 解决目前国内外流行的积雪产品在该地区积雪监测精度较差的问题, 建立一套高精度的适合祁连山区等复杂地形下的积雪制图算法, 为准确评价祁连山乃至整个青藏高原积雪变化奠定基础, 同时为水文、气候模型提供更准确的输入因子。

1 研究区概况

研究区系祁连山冰沟流域, 地处祁连山中段北坡黑河流域上游东支的一个较大二级支流。流域经纬度范围为100° 12'-100° 18' E, 38° 01'-38° 04' N, 平均海拔3 920 m, 总面积为30.48 km2。流域的自然景观具有明显垂直地带性, 森林带的海拔为2 800~3 300 m, 灌丛草甸带为3 400~3 700 m, 海拔4 000 m以上多为无植被的高山荒漠带, 属我国西部高山高原多年冻土区[17]。流域属于大陆性气候, 降水量丰富, 年均降水量达774 mm, 但年度分配不均匀, 主要集中在夏季(6月-8月)。流域内积雪属于季节性积雪, 深度平均约0.5 m, 最深0.8~1 m[18]。研究区地形复杂, 山区坡度较大[2]图1为研究区的Landsat 8陆地成像仪(Operational Land Imager, OLI)的假彩色合成影像图。在研究区选择两个区域, 区域1用来构建MODIS积雪面积比例监测模型, 区域2用于优化BP神经网络模型的输入参数, 整个研究区用来验证3种MODIS积雪面积比例模型的精度。

图1 祁连山冰沟流域Landsat 8 OLI假彩色合成影像Fig. 1 The Landsat 8 OLI false color composite image in Qilian Binggou Basin

2 数据与方法
2.1 研究数据

2.1.1 光谱数据 浅雪光谱数据的测量时间为2017年1月4日和5日。所用仪器为美国ASD公司生产的便携式地物光谱仪, 其光谱测量范围为350-2 500 nm, 光谱分辨率可达1 nm。在研究区选取18个积雪覆盖比例从0到1变化的典型样点, 样点观测范围设置直径为50 cm的圆, 下垫面类型均为岩甸。经过白板校正, 在距离雪面50 cm的垂直高度上测量积雪光谱, 同时测量了对应样点的雪深。为了降低随机误差, 光谱及雪深测量在每个样点处均进行了3次。期间用数码相机在对应样点处进行拍照, 并记录相应的编号。

与其它地表覆盖物相比, 积雪具有在可见光波段较高和在短波红外波段较低的反射率特性[19]。基于NDSI的MODIS积雪自动检测算法, 充分利用了这种光谱特点, 使用第4波段MODIS 4(0.545 0-0.565 0 μ m)和第6波段MODIS 6(1.628 0-1.652 0 μ m)的反射率(R4, R6)计算NDSI, 采用一套分组决策测试方法来检测积雪[20]。NDSI是较高反射率和较低短波红外反射率波段的一种组合指数:

NDSI= (R4-R6)(R4+R6)(1)

式中:R4表示积雪在可见光波段(MODIS第4波段)的反射率; R6表示积雪在短波红外波段(MODIS第6波段)的反射率。

一般而言, 积雪具有比地表其它地物类型高的NDSI值。当一个像素的NDSI≥ 40%时, 则该像素划分为积雪; 当NDSI< 40%, 并且第2波段(0.841 0-0.876 0 μ m)的反射率(R2)高于11%时, 则该像素划为水体; 当NDSI≥ 40%且R4> 10%时, 则该像素划为无积雪覆盖的地类[21]。此方法对祁连山冰沟流域的浅雪监测精度较差的原因除了地形复杂导致积雪分布破碎化(混合像元大量存在)外, 可能还有一部分原因要归于浅雪的光谱反射特征, 即由于部分可见光可以穿透浅雪层, 造成浅雪在可见光波段反射率较低, 而造成NDSI的阈值设置对浅雪提取没有好的适应性。所以, 本研究通过浅雪光谱特征分析来探讨其反射率对MODIS浅雪面积比例提取精度的影响。

2.1.2 MODIS数据 本研究所用MODIS数据行列号为h25v05, 为正弦曲线地图投影(sinusoidal projection), 包括MODIS逐日积雪产品MOD10A1、地表反射率产品MOD09GA, 以及植被指数产品MOD13Q1。其中, MOD10A1来源于美国国家冰雪数据中心NSIDC(National Snow and Ice Data Center)网站, 空间分辨率为500 m, 版本为V006, 时相为2016年11月2日。该产品包括积雪面积比例FRA(fractional snow cover)、雪被指数(NDSI)、积雪反照率SA(snow albedo)、质量评估QA(quality assessment)等多种数据。其中FRA数据编码值如表1所示。研究表明, MOD10A1所包含的积雪面积比例产品在研究区的适应性较差, 与同时相TM影像的积雪面积比例相差较大, 二者相关系数仅为0.740 0[22]。MOD09GA与MOD13Q1来源于美国国家航空航天局NASA(National Aeronautics and Space Administration)网站。地表反射率产品MOD09GA是在忽略了大气散射和吸收的情况下对地表光谱反射率的估算, 包括MODIS 1至7波段的反射率数据(ρ 17), 与MOD10A1同时相。MOD13Q1所包含的植被指数NDVI(mormalized difference vegetation index)数据的空间分辨率为250 m, 时间分辨率为16 d, 时相为2016年11月5日。

表1 MOD10A1积雪面积比例编码 Table 1 The fractional snow cover codes of MOD10A1

2.1.3 Landsat 8 OLI数据 本研究所用的Landsat 8陆地成像仪OLI数据来源于美国国家地质研究所USGS(United States Geological Survey)网站, 时相为2016年11月2日, 空间分辨率为30 m, 用于提取研究区积雪面积比例, 为了与MODIS标准产品进行区分, 这里用FSC表示, 并以其作为真值验证3种积雪面积比例反演模型的精度。

2.1.4 DEM高程数据 数字高程模型DEM(digital elevation model)是用一组有序的数值阵列形式表示地面高程信息的栅格数据。研究区90 m空间分辨率的DEM数据来源于USGS网站, 在ArcGIS 10.2.2软件中利用最近邻法将其空间分辨率重采样为500 m, 最邻近法是把原始图像中距离最近的像元值填充到新图像中, 然后将结果作为积雪面积比例BP神经网络模型的一个备选参数。

2.2 研究方法

2.2.1 数据预处理 数据预处理包括实测数据预处理与遥感数据预处理。1)实测数据预处理:对光谱数据与雪深数据进行一系列处理, 得到浅雪光谱反射率数据[23]; 然后分别利用ENVI 5.3软件与ArcGIS 10.2.2软件处理得到对应数码相片的积雪面积比例; 最后根据对应编号, 统计出不同面积比例浅雪的光谱反射率数据, 进行Savitzky-Golay卷积平滑处理, 制作出浅雪的光谱反射率曲线。2)遥感数据预处理:利用MODIS数据处理工具MRT(MODIS reprojection tool), 对MODIS产品(MOD10A1、MOD09GA和MOD13Q1)进行格式转换和坐标变换处理(正弦曲线投影转换为地理坐标, 椭球体选择WGS84, 重采样方法选用最邻近法, 图像文件转换为Geotiff格式)[13], 输出数据分别为NDSI、FRA、ρ 17和NDVI。同时, 对Landsat 8 OLI数据进行辐射定标和大气校正处理[24], 结合Landsat 8 OLI波段参数, 并参考SNOMAP算法[16], 利用OLI 影像的3、6波段计算NDSI, 根据研究区积雪特征, 设定NDSI的阈值为35%[25]。另外, 设定b5(第5波段)> 11%这个附加条件消除水体的干扰, 生成二值积雪分类图, 分辨率为30 m[26], 其中, 1代表积雪像元, 0代表非积雪像元。由于Landsat 8 OLI数据相对于MODIS具有更高的空间分辨率, 本研究假设Landsat 8 OLI积雪分类数据为地面真值, 并对MODIS积雪覆盖率数据进行验证。3)用研究区以及区域1和区域2的矢量图分别裁剪出对应区域的NDSI、FRA、ρ 17、NDVI、FSC和DEM。

2.2.2 线性回归模型法 Salomonson和Appel[27]提出了线性回归方法, 他们认为NDSI与FSC存在着某种线性关系, 通过建立二者之间的线性回归方程, 从而直接利用NDSI来推算FSC。该方法操作性强, 易于理解, 本研究基于上述原理分3步来实现研究区MODIS积雪面积比例模型的构建:1)计算每个MODIS格网(500 m× 500 m)的FSC。利用Landsat 8 OLI影像的二值积雪面积图与MODIS/NDSI图进行叠合分析, 利用ArcGIS 10.2.2软件的Fishnet工具统计出每个MODIS/NDSI网格内的FSC。其原理为先统计出每个MODIS/NDSI网格内所包含的Landsat 8 OLI分类图中值为1的积雪像元个数ns, 以及总像元数 nt26, 再利用式(2)计算出每个MODIS像元的FSC; 2)利用ArcGIS 10.2.2软件的Fishnet工具统计出每一个MODIS网格的NDSI的值; 3)最后利用统计法构建基于NDSI与FSC之间的线性回归模型。

FSC= nsnt(2)

式中:ns表示每个MODIS像元中的所含有的Landsat 8 OLI积雪分类图中值为1的积雪像元个数, nt表示总像元数。

2.2.3 全约束线性混合像元分解法 当具有不同波谱属性的物质出现在同一个像元内时, 就会出现混合像元。混合像元不完全属于某一种地物, 为了能让分类更加精确, 同时使遥感定量化更加深入, 需要将混合像元分解成一种地物占像元的百分含量(丰度), 即混合像元分解, 也叫亚像元分解。混合像元的存在, 是传统的像元级遥感分类和面积量测的精度难以达到使用要求的主要原因[28]。混合像元解混算法主要有线性混合像元分解法[29]、非负矩阵分解法[30]、稀疏回归解混算法[31]、贝叶斯分解法[32]、双线性分解算法[33]和梯度分解混合算法[34]等。

全约束最小二乘线性混合像元分解模型是解混算法中较常用的方法, 该模型认为像元在某一光谱波段的反射率是由各组成端元反射率与像元组分为权重系数的线性组合, 解混主要由两步构成, 第一步是提取“ 纯” 地物的光谱, 即端元提取; 第二步是用端元的线性组合来表示混合像元, 即混合像元分解。选取合适的端元是混合像元成功分解的关键。现有的端元提取算法大多以高光谱影像为基础, 主要分为凸面体分析与统计分析[35]。约束条件包括, 组分和为1, 且非负, 满足如下公式:

Yi= j=1mrijfj+ni (3)

j=1mrij=1 (4)

0≤ rij≤ 1 (5)

式(3)、(4)及(5)中, m为端元数, Yi是在i波段下的混合像元的总体反射率值, rij是第i个端元在i波段下的反射率值, fj是第j个端元在混合像元内所占的比例系数, ni是光谱在i波段的误差。

本研究基于上述原理, 利用连续最大角凸锥SMACC(squential mximum angle convex cone)算法从地表反射率产品MOD09GA影像中提取端元波谱以及丰度图像, 再利用全约束最小二乘混合像元解混方法进行MODIS混合像元分解[29], 提取研究区MODIS各像元内的积雪面积比例。采用均方根误差RMSE(root mean square error)对混合像元分解结果进行精度评价。

2.2.4 人工神经网络模型 BP神经网络是一种按误差逆向传播算法训练的多层前馈网络, 是人工神经网络的一种。最初由Werbos于1974年提出, 在1986年被Rumelhart等普及化[36], 是目前在医疗、环境和技术等领域应用最广泛的神经网络模型之一[37]。其学习规则是使用梯度下降法, 通过反向传播来不断调整网络的权值和阈值, 使网络的误差平方和达到设定值。BP神经网络模型拓扑结构包括输入层(input layer)、隐含层(hidden layer)和输出层(output layer), 它是包含多个隐含层的网络, 具备处理线性不可分问题的能力[38]

本研究以MODIS产品(MOD09GA的ρ 17, NDSI, NDVI, FRA)及DEM共11个自变量为网络输入参数, Landsat 8 OLI积雪面积比例真值(FSC)为网络输出参数, 基于LM(Levenberg-Marquardt)训练算法, 训练出了研究区积雪面积比例模拟的BP神经网络模型, 模型所用基本参数如表2所示。其中, BP神经网络模型的表现能力通过RMSE[式(6)]与拟合决定系数R2进行评价, 网络输入参数进行了偏差法标准化处理[式(7)]。

RMSE= i=1n(yi-yi')2n-1(6)

X= x-μσ(7)

式中:yi表示真值, yi'表示模拟值, i表示某一样本, n表示样本数量; μ 为均值, σ 为标准差。

为了确定研究区积雪面积比例模拟的最佳输入参数, 分别以11个自变量的16种组合方式(表3)进行模拟实验[39]。每次实验均用区域1进行BP神经网络模型构建, 然后在区域2对各种组合方式进行精度评价。

表2 积雪面积比例的BP神经网络模型参数 Table 2 BP artificial neural network model parameters of fractional snow cover
表3 BP神经网络模型的不同实验方案 Table 3 Different scheme of the BP artificial neural network model
3 结果与分析
3.1 浅雪光谱特征分析

图2是野外实测所得的不同面积比例的浅雪光谱反射率曲线(350-1 800 nm)。由于波长大于1 800 nm的波段积雪反射率受噪声干扰强烈, 也易受水分吸收影响, 考虑到积雪的光谱特征大多分布在可见光及近红外区域, 故在本研究中选取350-1 800 nm的波段, 分析其光谱特征。可以看出, 不同面积比例的浅雪光谱反射率曲线均符合积雪的光谱变化规律(图2)。其中, 积雪面积比例为5/6和6/6的浅雪光谱反射率在350-900 nm波段内出现了反射率大于1的情况, 造成这种结果的主要原因可能是:光谱仪探头视场范围内积雪面积比例较高, 加之测量时太阳光线较强, 积雪反射率整体较高, 从而导致光谱仪在此范围内达到了光饱和状态。不同积雪面积比例的浅雪光谱均存在多处反射“ 峰” 和吸收“ 谷” , 且位置大致相同。雪的粒径大小、雪花絮状分裂的形态和积雪的松紧程度不同都对雪被的光谱特性有明显的影响。受积雪光谱在可见光区的高反射率、短波近红外区的强吸收、以及积雪下垫面不同的影响, 在350-1 400 nm的波段内, 浅雪光谱反射率均随积雪面积比例的增大呈增高的趋势, 且反射率均在40%以上。在1 400-1 800 nm的波段内, 浅雪光谱反射率均随积雪面积比例的增大呈降低的趋势, 且反射率均在40%以下。裸地(主要为岩甸)的光谱反射曲线变化较平缓, 反射率介于9%~23%。整体上, 浅雪的光谱在可见光波段具有较高的反射率, 并没有出现因为太阳光容易穿透浅雪层而导致其在可见光波段反射率较低的情况。从而可以说明, 浅雪的光谱反射率对基于NDSI阈值法的MODIS全球积雪产品在该地区的浅雪提取精度几乎没影响, 精度差的主要原因是该地区复杂的地形导致的积雪分布破碎化, 即混合像元的大量存在。

图2 冰沟流域不同积雪面积比例的光谱反射曲线Fig. 2 The spectral reflectance of different fractional snow cover in Binggou Basin

3.2 积雪面积比例模型

本研究基于区域1全部像元做MODIS/NDSI与Landsat 8/FSC之间的回归分析。图3显示研究区MODIS/NDSI与Landsat 8 OLI 资料获取的地表真实积雪面积比例数据之间的散点图, 大量的数据点集中在趋势线的两侧, NDSI与FSC之间有很好的相关性, 拟合决定系数(R2)为0.601 7, 模型拟合结果较好。所以, 以此建立线性模型式(8), 应该能较好地反演出研究区的实际积雪面积比例。

FSC=0.813 3NDSI+0.252 2 (8)

式中:FSC表示Landsat 8 OLI 资料获取的地表真实积雪面积比例, NDSI表示雪被指数。

图3 雪被指数NDSI与Landsat 8 OLI积雪面积比例真值FCS之间的线性拟合Fig. 3 Scatter diagram of MODIS/NDSI versus ture value FSC from Landsat 8 OLI

结合研究区的土地覆盖类型, 通过SMACC端元提取算法从MOD09GA地表反射率数据中提取出了5种端元波谱(图4), 分别是积雪、植被、裸地、土壤和水域。可以看出, 5种端元光谱之间具有明显的差异, 变化趋势也基本与各地物光谱特征相似。端元提取之

后, 基于全约束偏最小二乘法线性混合像元光谱解混算法, 对研究区的地表反射率数据进行了混合像元分解, 得到了不同端元的丰度图像, 其中研究区积雪比例的范围为0~0.738 6, 平均值为0.242 4。

图4 人工选取的5类地物端元光谱曲线Fig. 4 The spectral curve of five surface features

表4列出了BP神经网络模型输入参数的组合方式的模拟评价结果, 并按照拟合决定系数(R2)进行了排序(从高到低)。结果表明, 16种组合方式构建的BP神经网络模型, 实测值与模拟值之间的R2介于0.339 2~0.553 4, 平均值为0.453 4, RMSE介于0.218 2~0.283 7, 平均值为0.236 0。方案8训练的模型[输入参数为:(ρ 17)+NDSI+DEM], 相较于其它输入参数训练的模型具有较高的R2(0.553 4)和较低的RMSE(0.218 2), 故以此参数组合的BP神经网络模型为预测研究区积雪面积比例的最佳模型。图5显示了最佳模型的回归拟合结果。可以看出, 4类数据集(训练集、验证集、测试集、全部数据集)的目标输出结果与仿真输出结果之间均具有很高的相关性, 且RMSE均较小。从验证集(图6)的泛化性能(R2=0.700 2, RMSE=0.205 5)和测试集(图6)的预测能力(R2=0.630 4, RMSE=0.161 3)来看, 均具有理想的模拟精度, 可以用该BP神经网络模型来反演整个研究区的积雪面积比例。

表4 输入参数的16种组合模型的模拟评价 Table 4 Evaluation of input parameters’ 16 combinations of BP

图5 BP神经网络模型最佳输入参数的模拟结果
注:(a)、(b)、(c)和(d)分别表示训练集、验证集、测试集和全部数据集的积雪面积比例真值与其对应的模拟值之间的拟合结果。
Fig. 5 The simulation results of BP artificial neural network model with optimal parameters
Note:(a), (b), (c) and (d) indicate the fitting results between true values and simulated values of the training data set, the validation set, the test set , the full data set, respectively.

3.3 模型精度验证

以整个研究区作为验证区域, 分别用3种MODIS积雪面积比例反演模型(线性回归、线性混合像元分解及BP神经网络)模拟得到研究区的积雪面积比例(图6)。可以看出, 线性回归模型与BP神经网络模型反演得到的积雪面积比例图与Landsat 8 OLI获取的真值积雪图更接近, BP神经网络模型模拟结果尤其在积雪边缘严重破碎部分, 呈现出更好的面积比例变化趋势; 而线性混合像元分解模型对于积雪破碎区的信息提取结果较差, 漏分现象很严重。

为了进一度定量地验证模型精度, 将3种模型的反演结果分别与Landsat 8 OLI的FSC进行统计分析, 并使用平均雪盖率、平均绝对误差MAE(mean absolute error)、正向平均误差(Ea)、负向平均误差(Eb)、RMSE以及R2等评价指标对上述模型进行精度评价(表5)。可以看出, 3种积雪面积比例反演模型中, 线性回归模型与BP神经网络模型在5个精度评价指标(平均积雪面积比例、MAE、Ea、RMSER2)上远优于线性混合像元分解模型, 说明线性混合像元分解模型(基于SMACC端元提取算法和全约束偏最小二乘法线性混合像元分解算法)在小尺度、地形复杂区域(如冰沟流域)的积雪面积比例提取方面具有一定的不足之处, 后续研究应该有针对性的对此方法做进一步研究; 线性回归模型的负向平均误差较大(Eb=-0.145 1), 表明线性回归模型倾向于在局部地区提取更多的雪盖信息, 而正向平均误差(Ea=0.171 7)小于线性混合像元分解模型(Ea=0.230 1)和BP神经网络模型(Ea=0.177 7), 在一定程度上, 可以说明线性回归模型在局部地区的漏分现象较严重; 线性回归模型与BP神经网络模型的积雪面积比例反演精度都较高, 但BP神经网络模型的均方根误差(0.061 6)远小于线性回归模型(0.164 3), 其稳定性有较大改进, R2也相对最大, 达0.650 9, 更适合于地形复杂地区积雪面积比例的监测。

图6 Landsat 8 OLI积雪图及MODIS 积雪面积比例图
注:(a)、(b)、(c)分别为由线性回归模型、线性混合像元分解模型、BP神经网络模型反演得到的积雪面积比例图, 分辨率为500 m; (d)是Landsat 8 OLI真实积雪面积比例图, 分辨率为30 m。
Fig. 6 Fractional snow cover from Landsat 8 OLI image and MODIS image
Note:(a), (b) and (c) indicate fractional snow cover retrieval map using linear regression model, linear unmixing model and BP artificial neural network model, respectively, with resolution of 500 m; (d) indicate the true map from Landsat 8 OLI with resolution of 30 m.

表5 3种MODIS积雪面积比例模型的误差统计 Table 5 Error statistics between the true value and three kinds of simulation value of MODIS fractional snow cover
4 讨论

BP神经网络模型反演积雪面积比例的最佳输入参数组合为(ρ 17)+NDSI+DEM, 而不是(ρ 17)+NDSI+NDVI+FRA+DEM。造成这种结果的主要原因可能是:模型引入更多输入信息的同时, 可能导致模型输入与结果之间的数据冗余增强, 反而使得输入信息数量不足, 在接下来的研究中应该有针对性地对此问题进行探讨。

线性混合像元分解模型的反演精度低于线性回归模型, 是3种模型里中精度最低的, 前人针对该方面的研究精度高于本研究[21]。造成这种结果的主要原因可能是:研究区太小, 所对应的MODIS像元没有代表性; MODIS数据的波段数较少, 而线性混合像元分解模型主要是针对高光谱的物质识别, 从而导致积雪提取精度差。在线性回归模型的研究中, 接下来可以考虑将分辨率更高的影像作为真值的载体, 如用分辨率可达厘米级的小型无人机拍摄的照片, 经过处理作为地面积雪面积比例的真值, 这样构建的模型可能精度更高。

BP神经网络模型的反演精度与稳定性相较于其它两种模型均较高, 是因为其可以在综合考虑多种因素的基础上, 通过大样本数据的深层次学习, 挖掘输入参数与输出参数之间的内在联系, 为研究区积雪面积比例的提取提供一种理想的方法, 为融雪径流估算、积雪的气候响应等研究提供更加丰富和准确的信息。

5 结论

本研究在分析研究区浅雪光谱特征的基础上, 探索了浅雪光谱反射率对MODIS浅雪面积比例提取精度的影响。然后在混合像元存在的事实下, 运用线性回归、线性混合像元分解及BP神经网络3种反演模型对研究区MODIS积雪面积比例进行了反演, 并对模型精度进行了验证, 主要获得了以下结论:

1)不同面积比例的浅雪光谱反射率曲线均符合积雪的光谱变化规律, 基于NDSI阈值法的MODIS浅雪提取精度几乎不受浅雪光谱反射率的影响, 由于该地区地形复杂, 从而导致积雪分布破碎化, 即混合像元的大量存在。

2)BP神经网络模型反演积雪面积比例的输入参数的最佳组合为(ρ 17)+NDSI+DEM。以其作为神经网络的输入参数时, 模型表现能力最好, 具有较为理想的预测能力和泛化能力。

3)线性混合像元分解模型在小尺度、地形复杂区域(如冰沟流域)的积雪面积比例提取方面精度较低, 线性回归模型和BP神经网络模型均具有较好的反演精度, 但是在模型的稳定性(基于Ea、Eb、RMSE等评价指标)方面, BP神经网络模型更胜一筹。

参考文献
[1] 王建. 卫星遥感雪盖制图方法对比与分析. 遥感技术与应用, 1999, 14(4): 29-36.
Wang J. Comparison and analysis on methods of snow cover mapping by using satellite remote sensing data. Remote Sensing Technology & Application, 1999, 14(4): 29-36. (in Chinese) [本文引用:1]
[2] 师银芳. 基于TM影像的祁连山冰沟实验区积雪信息提取研究. 兰州: 西北师范大学硕士学位论文, 2012.
Shi Y F. The study on extraction of snow information of the area using TM image in Binggou Watershed, Qilian Mountains. Master Thesis. Lanzhou: Lanzhou University, 2012. (in Chinese) [本文引用:2]
[3] 陈乾, 陈添宇. 祁连山区季节性积雪资源的气候分析. 地理研究, 1991, 10(1): 24-38.
Chen Q, Chen T Y. Climatical analysis of seasonal snow resources in QILIAN MT. Geographical Research, 1991, 10(1): 24-38. (in Chinese) [本文引用:1]
[4] 陈文年, 吴彦, 吴宁, 罗鹏. 3种高山植物的物候和种群分布格局在融雪梯度上的变化. 植物研究, 2011, 31(2): 206-212.
Chen W N, Wu Y, Wu N, Luo P. Variation in phenology and population distribution pattern of three alpine species along the snowmelt gradient. Bulletin of Botanical Research, 2011, 31(2): 206-212. (in Chinese) [本文引用:1]
[5] 刘琳, 孙庚, 吴彦, 何奕忻, 吴宁, 张林, 徐俊俊. 季节性雪被对青藏高原东缘高寒草甸土壤氮矿化的影响. 应用与环境生物学报, 2011, 17(4): 453-460.
Liu L, Sun G, Wu Y, He Y X, Wu N, Zhang L, Xu J J. Effect of seasonal snow cover on soil nitrogen mineralization in an alpine meadow on the eastern Tibetan Plateau. Chinese Journal of Applied & Environmental Biology, 2011, 17(4): 453-460. (in Chinese) [本文引用:1]
[6] 赵哈林, 周瑞莲, 赵悦. 雪生态学研究进展. 地球科学进展, 2004, 19(2): 296-304.
Zhao H L, Zhou R L, Zhao Y. Advance in snow ecology study in the world. Advance in Earth Sciences, 2004, 19(2): 296-304. (in Chinese) [本文引用:1]
[7] 唐志光, 王建, 彦立利, 李弘毅, 梁继. 基于MODIS的青藏高原亚像元积雪覆盖反演. 干旱区资源与环境, 2013, 27(11): 33-38.
Tang Z G, Wang J, Yan L L, Li H Y, Liang J. Estimating sub-pixel snow cover from MODIS in Qinghai-Tibet plateau. Journal of Arid Land Resources & Environment, 2013, 27(11): 33-38. (in Chinese) [本文引用:1]
[8] 孟宝平, 陈思宇, 崔霞, 冯琦胜, 梁天刚. 基于多源遥感数据的高寒草地生物量反演模型精度——以夏河县桑科草原试验区为例. 中国水稻科学, 2015, 32(11): 1730-1739.
Meng B P, Chen S Y, Cui X, Feng Q S, Liang T G. The accuracy of grassland vegetation biomass estimated model based on multi-source remote sensing data ——As a case of experimental area in Sangke grassland in Xiahe County. Pratacultural Science, 2015, 32(11): 1730-1739. (in Chinese) [本文引用:1]
[9] 柏延臣, 冯学智. 积雪遥感动态研究的现状及展望. 遥感技术与应用, 1997(2): 60-66.
Bai Y C, Feng X Z. Introduction to some research work on snow remote sensing. Remote Sensing Technology & Application, 1997(2): 60-66. (in Chinese) [本文引用:1]
[10] Lee S, Klein A G, Over T M. A comparison of MODIS and NOHRSC snow-cover products for simulating streamflow using the snowmelt runoff model. Hydrological Processes, 2010, 19(15): 2951-2972. [本文引用:1]
[11] Tekeli A E, Akyürek Z, Şorman A A, Şensoy A, Şorman A Ü. using MODIS snow cover maps in modeling snowmelt runoff process in the eastern part of Turkey. Remote Sensing of Environment, 2005, 97(2): 216-230. [本文引用:1]
[12] Sirguey P, Mathieu R, Arnaud Y. Subpixel Monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the southern Alps of New Zealand : Methodology and accuracy assessment. Remote Sensing of Environment, 2009, 113(1): 160-181. [本文引用:1]
[13] 郑有飞, 范旻昊, 张雪芬, 吴荣军. 基于MODIS遥感数据的混合像元分解技术研究和应用. 大气科学学报, 2008, 31(2): 145-150.
Zheng Y F, Fan M H, Zhang X F, Wu R J. Pixel unmixing technology of MODIS remote sensing data. Journal of Nanjing Institute of Meteorology, 2008, 31(2): 145-150. (in Chinese) [本文引用:2]
[14] 陈梦蝶, 黄晓东, 王玮, 梁天刚. CIVCO 地形校正模型对青藏高原地区积雪判别的有效性检验. 中国水稻科学, 2014, 31(2): 209-218.
Chen M D, Huang X D, Wang W, Liang T G. Validation test of the CIVCO topographic correction model on snow mapping algorithm in Tibetan Plateau. Pratacultrual Science, 2014, 31(2): 209-218. [本文引用:1]
[15] 郝晓华, 王建, 李弘毅. MODIS雪盖制图中NDSI阈值的检验——以祁连山中部山区为例. 冰川冻土, 2008, 30(1): 132-138.
Hao X H, Wang J, Li H Y. Evaluation of the NDSI threshold value in mapping snow cover of MODIS——A case study of snow in the middle Qilian Mountains. Journal of Glaciology & Geocryology, 2008, 30(1): 132-138. (in Chinese) [本文引用:1]
[16] Hall D K, Riggs G A, Salomonson V V. Development of methods for mapping global snow cover using moderate resolution imaging spectroradiometer Data. Remote Sensing of Environment, 1995, 54(2): 127-140. [本文引用:2]
[17] 杨志怀, 杨针娘, 王强. 祁连山冰沟流域径流分析与估算. 冰川冻土, 1992, 14(3): 251-257.
Yang Z H, Yang Z N, Wang Q. Runoff analysis and estimation of the Binggou Basin in the Qilian mountain. Journal of Glaciology & Geocryology, 1992, 14(3): 251-257. (in Chinese) [本文引用:1]
[18] 杨针娘, 杨志怀, 梁凤仙, 王强. 祁连山冰沟流域冻土水文过程. 冰川冻土, 1993, 15(2): 235-241.
Yang Z N, Yang Z H, Liang F X, Wang Q. Permafrost hydrological processes in Binggou basin of Qilian mountains. Journal of Glaciology and Geocryology, 1993, 15(2): 235-241. (in Chinese) [本文引用:1]
[19] 黄晓东, 张学通, 李霞, 梁天刚. 北疆牧区MODIS积雪产品MOD10A1和MOD10A2的精度分析与评价. 冰川冻土, 2007, 29(5): 722-729.
Huang X D, Zhang X T, Li X, Liang T G. Accuracy analysis for MODIS snow products of MOD10A1 and MOD10A2 in northern Xinjiang area. Journal of Glaciology & Geocryology, 2007, 29(5): 722-729. (in Chinese) [本文引用:1]
[20] Hall D K, Riggs G A, Salomonson V V, Digirolamo N E, Bayr K J. MODIS snow-cover products. Remote Sensing of Environment, 2002, 83(1-2): 181-194. [本文引用:1]
[21] 梁天刚, 高新华, 黄晓东, 张学通. 新疆北部MODIS积雪制图算法的分类精度. 干旱区研究, 2007, 24(4): 446-452.
Liang T G, Gao X H, Huang X D, Zhang X T. Study on the accuracy of MODIS snow cover mapping algorithm in northern Xinjiang. Arid Zone Research, 2007, 24(4): 446-452. (in Chinese) [本文引用:2]
[22] 张颖, 黄晓东, 王玮, 梁天刚. MODIS逐日积雪覆盖率产品验证及算法重建. 干旱区研究, 2013, 30(5): 808-814.
Zhang Y, Huang X D, Wang W, Liang T G. Validation and algorithm redevelopment of MODIS daily fractional snow cover products. Arid Zone Research, 2013, 30(5): 808-814. (in Chinese) [本文引用:1]
[23] 胡远宁, 崔霞, 孟宝平, 杨淑霞, 梁天刚. 甘南高寒草甸主要毒杂草光谱特征分析. 中国水稻科学, 2015, 32(2): 160-167.
Hu Y N, Cui X, Meng B P, Yang S X, Liang T G. Spectral characteristics analysis of typical poisonous weeds in Gannan alpine meadow. Pratacultural Science, 2015, 32(2): 160-167. (in Chinese) [本文引用:1]
[24] 王敏, 高新华, 陈思宇, 冯琦胜, 梁天刚. 基于Land sat8遥感影像的土地利用分类研究——以四川省红原县安曲示范区为例. 中国水稻科学, 2015, 32(5): 694-701.
Wang M, Gao X H, Chen S Y, Feng Q S, Liang T G. The land use classification based on Land sat 8 remote sensing image——A case study of Anqu demonstration community in Hongyuan County of Sichuan Province. Pratacultural Science, 2015, 32(5): 694-701. (in Chinese) [本文引用:1]
[25] 王玮. 青藏高原牧区积雪监测研究. 兰州: 兰州大学硕士学位论文, 2011.
Wang W, Monitoring snow cover in pastoral areas on Qinghai-Tibetan plateau. Master Thesis. Lanzhou: Lanzhou University, 2011. (in Chinese) [本文引用:1]
[26] 周强, 王世新, 周艺, 王丽涛. MODIS亚像元积雪覆盖率提取方法. 中国科学院大学学报, 2009, 26(3): 383-388.
Zhou Q, Wang S X, Zhou Y, Wang L T. Algorithm for MODIS subpixel snow fraction. Journal of the Graduate School of the Chinese Academy of Sciences, 2009, 26(3): 383-388. (in Chinese) [本文引用:1]
[27] Salomonson V V, Appel I. Estimating fractional snow cover from MODIS using the normalized difference snow index. Remote Sensing of Environment, 2004, 89(3): 351-360. [本文引用:1]
[28] 童庆禧, 张兵, 郑兰芬. 高光谱遥感. 北京: 高等教育出版社, 2006.
Tong Q X, Zhang B, Zheng L F. Hyperspectral Remote Sensing-Theory, Technology and Application. Beijing: Higher Education Press, 2006. (in Chinese) [本文引用:1]
[29] Heinz D C, Chein-I-Chang. Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery. IEEE Transactions on Geoscience & Remote Sensing, 2001, 39(3): 529-545. [本文引用:2]
[30] Lee D D, Seung H S. Learning the parts of objects with nonnegative matrix factorization. Nature, 1999, 401: 788-791. [本文引用:1]
[31] 王天成. 高光谱图像解混算法研究. 哈尔滨: 哈尔滨工业大学硕士学位论文, 2016.
Wang T C. Research on unmixing algorithms for hyperspectral images. Master Thesis. Harbin: Harbin Institute of Technology, 2016. (in Chinese) [本文引用:1]
[32] 李熙, 关泽群, 沈体雁. 基于贝叶斯网络的混合像元分解模型. 光电子·激光, 2008, 19(8): 1121-1126.
Li X, Guan Z Q, Shen T Y. Sub-pixel analysis based on Bayesian network. Journal of Optoelectronics Laser, 2008, 19(8): 1121-1126. (in Chinese) [本文引用:1]
[33] 李春芝. 基于高光谱影像的抗噪模型及解混算法研究. 上海: 华东师范大学博士学位论文, 2014.
Li C Z. Hyper-spectral unmixing algorithms based on an anti-noise model. PhD Thesis. Shanghai: East China Normal University, 2014. (in Chinese) [本文引用:1]
[34] 汤仪平, 金福江. 最速下降法和共轭梯度的混合算法及全局收敛. 华侨大学学报: 自然版, 2007, 28(2): 124-126.
Tang Y P, Jin F J. A Hybrid algorithm of the steepest descent method and the conjugate gradient method and its global convergence. Journal of Huaqiao University: Natural Science Edition, 2007, 28(2): 124-126. (in Chinese) [本文引用:1]
[35] 郝晓华, 王杰, 王建, 黄晓东, 李弘毅, 刘艳. 积雪混合像元光谱特征观测及解混方法比较. 光谱学与光谱分析, 2012, 32(10): 2753-2758.
Hao X H, Wang J, Wang J, Huang X D, Li H Y, Liu Y. Observations of snow mixed pixel spectral characteristics using a ground-based spectral radiometer and comparing with unmixing algorithms. Spectroscopy and Spectral Analysis, 2012, 32(10): 2753-2758. (in Chinese) [本文引用:1]
[36] Rumelhart D E, Hinton G E, Williams R J. Learning Internal Representations by Error Propagation. Cambridge: MIT Press, 1986: 318-362. [本文引用:1]
[37] 王德明, 王莉, 张广明. 基于遗传BP神经网络的短期风速预测模型. 浙江大学学报: 工学版, 2012(5): 837-841.
Wang D M, Wang L, Zhang G M. Short-term wind speed forecast model for wind farms based on genetic BP neural network. Journal of Zhejiang University: Engeering Science, 2012(5): 837-841. (in Chinese) [本文引用:1]
[38] 陈明. MATLAB神经网络原理与实例精解. 北京: 清华大学出版社, 2013.
Chen M. MATLAB Neural Network Principle and Instances of Precision Solution. Beijing: Tsinghua University Press, 2013. (in Chinese) [本文引用:1]
[39] Hou J, Huang C. Improving mountainous snow cover fraction mapping via artificial neural networks combined with MODIS and ancillary topographic data. IEEE Transactions on Geoscience & Remote Sensing, 2014, 52(9): 5601-5611. [本文引用:1]