基于多源遥感数据的疏勒河上游山区流域VIC-CAS模型积雪模拟效果评估

2021-05-24郭佳锴张世强

冰川冻土 2021年2期

郭佳锴，李哲，李飞，张世强

（1.西北大学陕西省地表系统与环境承载力重点实验室，陕西西安 710127； 2.西北大学城市与环境学院，陕西西安 710127）

0 引言

在全球水循环中，积雪的积累和消融过程调节着水的年内再分配，是干旱半干旱地区春季最重要的淡水资源［1］。积雪作为冰冻圈的重要组成部分，是气候系统中最敏感的变化因子之一［2］。积雪是全球气候变化过程中的重要因素，积雪变化对地表反照率、全球和地区能量平衡有重要影响［3-4］。同时，积雪分布和深度是天气数值预报模式的重要参数［5］。

积雪覆盖和雪深可以由地面观测、遥感反演和模型模拟获得［4-5］。地面观测雪深资料是目前研究积雪长期变化特征最广泛、最可靠的资料，具有时间尺度长，连续性好的优点，但也存在着空间分辨率低，空间分布不均匀，观测投入较大的不足［6］。近

年来，遥感积雪数据的不断发展，为寒区水文模型提供了丰富的数据源，改善了以往地面观测资料分辨率低，成本高的缺点［7］。目前，常见的积雪覆盖遥感产品有Landsat 和SPOT［8］，AVHRR［9］，MODIS 积雪产品［10］及SMMR、SSM/I［11-12］和AMSR-E［13-14］等微波积雪产品。MODIS 数据因其较高的时空分辨率和光谱分辨率在积雪研究与动态变化监测领域中广泛应用［15-16］。张颖等［17］研究发现MODIS 逐日积雪覆盖率产品在青藏高原地区精度高于MODIS 标准积雪产品。在积雪深度监测方面，被动微波遥感数据被认为是监测大尺度积雪深度时空变化的有效手段。现有雪深产品在青藏高原地区表现出较大的不确定性［18］。Dai 等［19］基于青藏高原的台站雪深观测数据和近年来对地面雪深的大规模调查资料，对利用SSM I/S 和AMSR-E 亮温数据反演雪深的算法进行了校准，进而生产了中国雪深数据集，这是目前青藏高原地区精度最高的雪深数据产品。

积雪的积累和消融过程是寒区水文模型的重要组成部分。寒区水文模型模拟积雪一般可分为度日因子和能量平衡两种方法［20］。其中，度日因子法是基于冰雪消融和气温之间的关系计算［21］，目前SRM［22］和SWAT［23］融雪模块均采用过这种方法。能量平衡模型是模拟积雪-大气，积雪-土壤及积雪内的物质和能量平衡［21］，常用的模型有ISNOBAL［24］、UEB［25］等，一些分布式水文模型也出现了融雪模块，如VIC［26］、DHSVM［27］等。近年来，已经出现了多个利用卫星遥感获得的积雪分布来驱动水文模型提高精度的成功案例［7，28-30］，表明利用卫星遥感积雪数据有助于改善以往水文模型仅仅利用出山口径流资料进行率定和检验的不确定性。如赵军等［28］在疏勒河流域上游将SRM 融雪径流模型与MODIS 积雪产品结合，取得了较好的模拟效果。Andreadis等［29］利用EnKF方法将遥感积雪面积数据和雪水当量数据同化到VIC 模型中，对中低海拔，融雪期和浅层积雪的积雪变量模拟精度取得显著提升。Che 等［30］利用EnKF 方法将被动微波遥感雪深数据同化到陆面模型中，提高了积累期的雪深估算精度。VIC-CAS［31］分布式水文模型中很好地考虑了单条冰川的融水过程和冰川变化，其模拟的疏勒河流域的冰川变化与遥感监测的单条冰川变化［32］具有很好的可比性，其径流模拟也取得了很好的效果［33］，为评估模型中对积雪的积累和消融过程的合理性提供了很好的基础。

因此，本研究选取了在降水等观测资料较为丰富的疏勒河上游山区流域，基于中分辨率成像光谱仪（Moderate Resolution Imaging Spectroradiometer，简称MODIS）逐日积雪覆盖资料和中国长序列遥感反演雪深资料，选取对冰川融水模拟较好的VICCAS 模型模拟的积雪覆盖和雪深进行评估，并对比模型在不同高程带的模拟效果差异，为进一步改进模型中的积雪过程提供支撑。

1 研究区与数据

1.1 研究区概况

图1 疏勒河上游山区流域气象站、雨量计和高程带分布及VIC-CAS模型中的182个子流域单元Fig.1 Spatial distribution of meteorological stations and rain gauges，elevation zones，and 182 sub-basins of VIC-CAS in mountainous upper reach of the Shule River basin

疏勒河是河西走廊三大内陆河之一，发源于祁连山腹地的疏勒南山北坡，其出山口昌马水文站以上为疏勒河上游（96.6°～99.0°E，38.2°～40.0°N）［34-35］（图1）。疏勒河上游地区地形主要由疏勒南山、托勒南山和疏勒河谷地组成，山区地势高峻、地形陡峭，谷地地形相对低缓。疏勒河上游面积约1.14×104km2，属于高原大陆性气候，研究区海拔介于2 100～5 750 m 之间，平均海拔为3 900 m［36-37］，其中2 000～3 000 m 高程带所占区域占整个研究区面积的7%，主要是河流出山口的河谷区；3 000～4 000 m高程带所占区域占整个研究区面积的45%，主要为疏勒河上游谷地中段，河谷两侧的山前缘区；4 000～5 000 m高程带所占区域占整个研究区面积的47%；5 000～6 000 m 高程带所占区域占整个研究区面积的2%［37］。研究区分布大量冻土，积雪和冰川［38］，多年平均气温约-4 ℃，多年平均降雨量约378.4 mm，降水主要集中在5—9月，降水量约占全年总降雨量的90%以上，冬季和春季以降雪为主［39］。

1.2 数据及预处理

1.2.1 MODIS积雪覆盖数据

本研究中的积雪覆盖观测数据利用美国国家雪冰中心（National Snow and Ice Data Center，NSIDC）网站下载的TERRA MODIS MOD09GA地表反射率产品计算获得。MOD09GA 的时间分辨率为日，空间分辨率为500 m。本研究收集了2002—2013年全年共计4 383 天的MOD09 产品1～7 通道的数字影像。在选择不同版本的MODIS积雪面积产品时，考虑到在青藏高原区域MOD09GA数据经过处理后的积雪覆盖与同时相Landsat TM/ETM+影像获取的积雪覆盖的相关系数r达到0.85，高于MOD10A1的0.74，最终选择MOD09GA作为检验数据［17］。

对MOD09GA 影像的预处理包括对卫星影像进行坐标变换。将正弦曲线投影转换为地理坐标，椭球体选为WGS84。选用双线性法重采样，将图像文件转换为GeoTIFF 格式，并利用研究区边界进行剪裁。根据雪的反射率特点，利用NDSI 算法对MODIS数据对每个像元的积雪覆盖度进行计算［40］，其公式为

式中：Band4、Band6 分别为MODIS 波段4 和波段6的反射率。

根据高扬等［41］在青藏高原地区不同土地覆盖类型NDSI 阈值积雪判别的研究结果，选取0.33 作为最佳阈值，并根据波段2 和波段4 排除水体和暗物质干扰。本研究利用MOD10A1产品中的云覆盖数据确定含云像元，在进行对比时对含云像元不作对比。不含云的各像元的积雪覆盖度，采用Salomoson 等［42］提出的NDSI 与真实亚像元积雪覆盖度之间的经验公式获得。

式中：FSC 为每个像元的积雪覆盖度；NDSI 为每个像元的NDSI 值。对逐日数据进行计算，最终得到2002—2013 年疏勒河上游山区流域MODIS 各子流域的日积雪覆盖度。

1.2.2 雪深数据

本研究中2002—2013 年疏勒河上游山区流域的雪深数据从国家冰川冻土沙漠科学数据中心网站（www.crensed.ac.cn）下载的中国雪深长时间序列数据集（the Long-term Snow Depth Dataset of China）中提取。该数据集中2002—2007 年的雪深数据用AMSR-E 反演获得，2008—2016 年的雪深采用SSMI/S传感器的亮度温度反演，其时间分辨率为日，空间分辨率均为0.25°。用于反演该雪深数据集的原始数据来自美国国家雪冰数据中心（NSIDC）处理的AMSR-E（2002—2007 年）和SSMI/S（2008—2019 年）逐日被动微波亮温数据，其首先通过对不同传感器的亮温进行交叉定标提高亮温在时间上的一致性，然后利用在针对中国地区修正的Chang算法进行雪深反演［11］。

本研究将下载的ASCII文件用Python语言批量转成栅格，利用研究区边界剪裁得到2002—2013年疏勒河上游山区流域逐月和逐年平均的积雪深度栅格文件。根据遥感观测雪深分别将最接近平均值的2004 年作为平雪年的代表，雪深最大的2008年作为多雪年代表，雪深最小的2013年作为少雪年代表。

1.2.3 VIC-CAS模型模拟的积雪覆盖度和雪深

VIC-CAS 模型是在VIC-3L 基础上改进的分布式水文模型［43］，其计算每个网格单元的能量和水分平衡，并考虑到积雪、融雪和土壤冻融过程［26］。VIC模型可以由日或日内尺度气象数据驱动，包括降水量、最高和最低气温以及风速数据，模型通过内置的插值程序生成各格网的气象驱动数据，模型参数包括土壤和植被参数，以及数字高程模型等［44］。VIC-CAS 模型增加了冰川模块，将单条冰川作为子流域，进而将单条冰川细分为100 m 间距高程带，利用度日因子逐高程带计算冰川消融量和物质平衡，并考虑了冰川面积的长期变化［45］。VIC-CAS 通过分别采用单层和双层模型计算冠层和地面积雪过程和融雪过程，其中在双层模型中考虑了积雪的积累和升华过程，并考虑了不同高程带积雪积累和消融的空间异质性［46］。VIC-CAS 积雪算法中的雪深与雪龄有关，考虑了积雪的密实化过程以及新雪的影响［47］。VIC-CAS 模型假定子流域内的积雪均匀分布，通过阈值分段考虑了积雪覆盖面积与雪深之间的关系，从而考虑了薄雪覆盖率较低的特征［29］。

Zhang等［33］利用VIC-CAS模型模拟了疏勒河上游山区流域的冰川径流和河川径流，与遥感获得的不同时期的单条冰川变化［32］和出山口观测径流对比表明，VIC-CAS 取得了很好的模拟效果。本研究采用Zhang 等［33］在疏勒河上游山区率定和验证的VIC-CAS 模型参数，与遥感获取的积雪覆盖度和积雪深度对比，分析了VIC-CAS 模型在疏勒河上游山区流域模拟的积雪覆盖度和积雪深度模拟效果。模型中的年降水梯度为根据1 139～4 156 m 的降水观测计算所得的14.654 mm·（100m）-1，模型的关键参数、率定过程、冰川模拟结果和河流径流结果详见文献［33］。本研究中直接选取了2002—2013 年疏勒河上游山区流域182个子流域的积雪覆盖度和雪深日模拟结果，其中，每个子流域的高程通过子流域中各高程带所占面积比率与高程的乘积累加得到。随后用属性表的子流域编号作为关键字，对VIC-CAS 模型结果数据文件中的积雪覆盖度和积雪深度两个变量的数据与疏勒河上游山区流域子流域矢量边界进行连接，从而将MODIS 积雪覆盖数据和中国雪深长时间序列数据与子流域编号链接统计。

考虑到遥感反演的雪深具有较大的不确定性，本研究中对模型的评价重点放在积雪覆盖度指标，对雪深指标主要对比其一致性，不作为评价重点。

1.3 统计分析

1.3.1 相关分析

采用相关系数（r）评价VIC-CAS 模型模拟与观测的积雪覆盖度和积雪深度值之间的相关性。

式中：Oi为VIC-CAS 模型模拟的积雪覆盖度和积雪深度；为VIC-CAS 模型模拟的日平均积雪覆盖度和日积雪深度的平均值；Mi为MODIS 观测的积雪覆盖度和中国雪深长时间序列数据集的积雪深度；为MODIS 影像积雪覆盖度和中国雪深长时间序列数据集积雪深度的平均值；n为总月数。

根据文献［48］，将r划分为6 个区间：r≤0 为负相关，0＜r≤0.2 为不相关或极弱相关，0.2＜r≤0.4 为弱相关，0.4＜r≤0.6为中相关，0.6＜r≤0.8为强相关，0.8＜r≤1为极强相关。

1.3.2 均方根误差

采用均方根误差（RMSE）分析方法评价VICCAS 模型模拟与观测的积雪覆盖度和积雪深度之间的离散程度。

式中：Mi为MODIS 影像积雪覆盖度和中国雪深长时间序列数据集的积雪深度；Oi为VIC-CAS模型模拟的积雪覆盖度和积雪深度的平均值；n为总月数。

1.4 评估流程

疏勒河上游山区VIC-CAS 模拟和遥感观测的积雪覆盖和雪深评估流程如图2所示。数据预处理包括对VIC-CAS 模型输出文件分别计算各子流域的月平均和年平均积雪覆盖度和积雪深度，对中国雪深长时间序列数据集进行栅格转换，提取研究区月平均和年平均积雪深度。对MOD09GA 数据的处理过程包括NDSI 计算，NDSI 与积雪覆盖度之间的经验线性回归方程计算，重采样，提取研究区范围。数据分析包含利用ArcPy语言对处理后的数据使用以表格显示分区统计工具统计182个子流域中积雪覆盖度和积雪深度的平均值，并计算r和

RMSE。

图2 疏勒河上游山区VIC-CAS模型模拟和遥感观测的积雪覆盖度和雪深评估流程Fig.2 Flow chart of the evaluation on snow coverage and snow depth simulated by VIC-CAS model and observed by remote sensing in mountainous upper reach of the Shule River basin

2 结果与分析

2.1 积雪覆盖度评估

2.1.1 平雪年、多雪年和少雪年积雪覆盖度

在疏勒河上游山区流域VIC-CAS 模型模拟和MODIS观测的月平均积雪覆盖度在平雪年，多雪年和少雪年的对比如图3 所示。r在多雪年最大，为0.67，属于强相关；少雪年为0.52，属于中相关；平雪年为0.37，属于弱相关。RMSE 在三个年份相差不大，多雪年为0.12，少雪年为0.09，平雪年为0.13。总体来看，VIC-CAS 模型在多雪年的模拟效果最好，少雪年和平雪年精度相对较低，且存在大量模型模拟为0 但MODIS 积雪观测有值的点，表明薄雪的模拟结果可能较差。

图3 VIC-CAS模型模拟和MODIS观测的月平均积雪覆盖度在疏勒河上游山区流域不同降雪年份的对比Fig.3 Comparisons in monthly averaged snow coverage between simulated by VIC-CAS model and observed by MODIS in mountainous upper reach of the Shule River basin in normal snow year（a），more snow year（b）and less snow year（c）

疏勒河上游山区流域模拟和观测的年平均积雪覆盖度在平雪年，多雪年和少雪年的空间分布如图4 所示。从空间分布看，VIC-CAS 模型模拟的积雪覆盖度偏小，多雪年模拟数据与MODIS 产品表现出最为相似的空间分布特征，即研究区中部高山区积雪覆盖度高，西北部河谷区积雪覆盖度低，表明其对积雪的空间分布模拟较好，少雪年和平雪年积雪空间分布模拟的特征相似性较差。

图4 VIC-CAS模型模拟和MODIS观测的年平均积雪覆盖度在不同降雪年份的空间分布Fig.4 Spatial distribution of yearly averaged snow coverage between simulated by VIC-CAS model and observed by MODIS in different years（a，b，c are simulated by VIC-CAS model；d，e，f are observed by MODIS）

2.1.2 不同高程带积雪覆盖度

VIC-CAS 模型模拟和MODIS 观测的月平均积雪覆盖度在不同高程带2 000～3 000 m，3 000～4 000 m，4 000～5 000 m 的对比如图5 所示。VICCAS模型在海拔4 000～5 000 m的模拟效果最好，其在平雪年，多雪年和少雪年中r分别为0.41，0.66和0.60，RMSE 分别为0.15，0.12，0.11，模拟精度较高。海拔2 000～3 000 m 模拟效果不佳，在平雪年出现了负相关-0.10。总的来看，VIC-CAS 模型在占比最高的4 000～5 000 m 高程带模拟精度最高，3 000～4 000 m 其次，在较低海拔的2 000～3 000 m积雪覆盖度模拟精度较低。

图5 VIC-CAS模型模拟和MODIS观测的月平均积雪覆盖度在疏勒河上游山区流域不同高程带的对比Fig.5 Comparisons in monthly averaged snow coverage between simulated by VIC-CAS model and observed by MODIS at different elevation zones in mountainous upper reach of the Shule River basin

2.2 雪深对比

2.2.1 平雪年、多雪年和少雪年雪深

VIC-CAS模型模拟和中国雪深长时间序列数据集观测的月平均积雪深度在平雪年，多雪年和少雪年的对比如图6 所示。VIC-CAS 模型在多雪年的r为0.44，属于中相关；在平雪年和少雪年相关程度较低，其中平雪年出现负相关-0.22。总体来看，VIC-CAS模型在不同降雪年份雪深一致性较低。平雪年和少雪年存在大量模拟积雪深度为0但积雪产品中深度不为0的情况，说明薄雪仍然是其主要差异所在。

图6 VIC-CAS模型模拟和中国雪深长时间序列数据集观测的月平均积雪深度在疏勒河上游山区流域不同降雪年份的对比Fig.6 Comparisons in monthly averaged snow depth between simulated by VIC-CAS model and observed by the Long-term Snow Depth Dataset of China in mountainous upper reach of the Shule River basin in normal snow year（a），more snow year（b）and less snow year（c）

模拟和观测的年平均积雪深度在平雪年，多雪年和少雪年的空间分布如图7 所示。可以看出VIC-CAS 模型模拟的雪深总体偏小，但很好地反映了雪深随高程变化的特点；观测的雪深数据与高程带的关系刻画不明显，主要是较粗的空间分辨率所致。

图7 VIC-CAS模型模拟和中国雪深长时间序列数据集观测的年平均积雪深度在不同降雪年份的空间分布Fig.7 Spatial distribution of yearly averaged snow depth between simulated by VIC-CAS model and observed by the Long-term Snow Depth Dataset of China in different years（a，b，c are simulated by VIC-CAS model；d，e，f are observed by the Long-term Snow Depth Dataset of China）

2.2.2 不同高程带雪深

VIC-CAS 模型模拟和中国雪深长时间序列数据集观测的月平均积雪深度在2 000～3 000 m、3 000～4 000 m、4 000～5 000 m 高程带的对比如图8所示。可以看出，VIC-CAS 模拟值与观测的积雪深度的r与高程带无关，其均在多雪年呈现中相关，少雪年无相关性，平雪年呈现负相关的特点，其可能仍与观测的积雪深度分辨率较粗有关。

3 讨论

图8 疏勒河上游山区流域VIC-CAS模型模拟和中国雪深长时间序列数据集观测的月平均积雪深度在不同高程带下的对比Fig.8 Comparisons in monthly averaged snow depth between simulated by VIC-CAS and observed by the Long-term Snow Depth Dataset of China at different elevation zones in mountainous upper reach of the Shule River basin

传统积雪产品大多采用地面气象站观测值空间插值和人工观测获得，如Bulygina 等［49］利用856个观测站的积雪深度数据研究俄罗斯的积雪覆盖和雪深的年际变化。这种方法需要耗费大量人力物力，且难以获得连续的，大面积的积雪信息［50］。在疏勒河上游山区流域不存在有关积雪的气象台站数据，无法提供有效的地面实测资料进行交叉验证的数据支持。利用遥感获取积雪覆盖和雪深能够有效地获取积雪的空间分布［1］。因此，本研究证明了采用遥感反演的积雪对整个研究区的积雪模拟效果进行评估是有效可行的。

从不同积雪覆盖的遥感产品的精度看，张颖等［17］研究表明，MOD10A1 产品对于积雪破碎区的信息提取较差，山体阴影的漏分现象更为严重，较小区域的积雪识别误差大，精度低。王雪璐等［51］研究表明，MOD10A1 在青海省的积雪分类精度明显低于在其他研究区的精度验证结果，经过NDSI 阈值调整后的MOD09GA 产品的精度高于同时期MOD10A1 产品。本研究中采用了这一算法，因此积雪覆盖产品应具有较高的精度。

目前来看，雪深数据产品空间分辨率仍较粗，难以反映高程分布的影响。Dai 等［18］在2016 年，基于SSMI/S（F17）对所有传感器的亮温数据进行校正，最终得到长期的中国积雪深度数据，并在2019年［52］开发了一种适用于被动微波的积雪识别方法，提升了产品的精度。尽管本研究中所采用的积雪产品是目前在中国精度最高的雪深数据产品，但依旧存在空间分辨率低和不确定性较大的问题，所以本研究中的评估以积雪覆盖结果为主，雪深结果重在对比两者之间的一致性。

需要指出的是，遥感反演的积雪覆盖和雪深也具有较大的不确定性。MODIS 积雪覆盖数据的不确定性主要在于MOD09GA 产品进行归一化处理的时候受到云，水体等影响以及处理过程中NDSI与积雪覆盖度的线性回归算法的误差。Salomonson等［42］建立了NDSI 与真实亚像元积雪覆盖度之间的关系，该回归方程是在西伯利亚地区积雪观测的基础上建立的，对其他地区的适用性可能存疑。本研究中存在大量观测值为0 的数据，可能是由于线性回归算法得到的积雪覆盖度剔除了小于0.1的数据产生的误差［42］。遥感雪深反演数据的不确定性主要在于产生中国地区雪深长时间序列数据集的被动微波卫星遥感资料空间分辨率很低，反演过程时积雪的空间异质性问题难以解决。

从空间分布对比看，VIC-CAS 模型模拟的积雪覆盖度和雪深偏小，自2009年起在研究区内开始布设较多的降水观测设备，获取了较为准确的降水梯度，模型的冰川融水和径流模拟效果均较好［33］，为对比积雪的相关结果提供了很好的基础。本研究发现VIC-CAS 模型和遥感反演的薄雪覆盖度的差异较大，在不同高程带对积雪的模拟精度差别不大，但在2 000～3 000 m 模拟精度较低，说明在低海拔区存在较大不确定性，这可能与模型中对积雪再分布和风吹雪的算法和参数化方案有关，特别是在薄雪情况下，其可能存在较大的不确定性。因此，需要在流域尺度上进一步加强对风吹雪分布、过程和参数化方案的观测和模拟优化。