基于MODIS反演雪深的融雪径流模拟

陈智梁，王 娟，李春红，王冉旋，王 奕，马志贵

（1.国家能源集团新疆吉林台水电开发有限公司，新疆伊犁835000；2.南瑞集团有限公司/国网电力科学研究院，南京211106）

新疆地区河流大多以冰雪融水为主要补给水源，其春、秋季来水以冰雪融水为主，夏季洪水多为融雪与暴雨的叠加。对于此类地区而言，融雪径流是河流的基础性来水，准确的融雪径流预报可以提高水资源利用效率，有效预防洪涝灾害［1-3］。

SRM 融雪径流模型因结构简单、应用效果较好已广泛应用于融雪径流预报［4-18］，它以积雪覆盖面积作为模型驱动之一，计算逐日融雪径流过程。随着遥感影像时空分辨率和光谱分辨率的逐渐提高，遥感已成为有效的积雪监测手段［19］，其中MODIS 遥感数据因具有较高的时空分辨率，其积雪覆盖产品通常被作为SRM 模型的输入［4-15］。NSIDC（National Snow and Ice Data Center）发布的基于MODIS 传感器的积雪覆盖产品有两类：全球逐日积雪覆盖产品MOD10A1/MYD10A1 和全球8日合成的积雪覆盖产品MOD10A2/MYD10A2，空间分辨率均为500 m［16-18］。逐日积雪覆盖产品的积雪分类精度受云层影响严重，在多云或阴天时，平均积雪识别率仅为17.8%［18，19］，因此融雪径流预报通常使用合成产品MOD10A2。在春秋积雪不稳定季节，积雪覆盖面积在短时间内就会发生明显变化，譬如一场薄薄的新雪会带来较大的积雪覆盖面积，又会在短短几天内融化殆尽，因此8日合成的遥感积雪产品无法反映积雪覆盖信息的逐日变化过程，逐日融雪径流计算也将产生较大的偏差。

SRM 模型假定雪盖区在单位计算时段内有充足的积雪可融化，雪深不设限，仅根据流域的积雪覆盖率、温度和度日因子计算逐日融雪径流过程。若每日的积雪覆盖率准确，则该假定对融雪径流计算结果的影响不大；但如前所述，可用的遥感积雪覆盖产品为8日合成，无准确的逐日积雪覆盖面积，因此要保证融雪径流计算精度，需要获取相对准确的逐日积雪覆盖信息。积雪覆盖面积是积雪的二维空间体现，精细化空间分布的雪深既可计算雪盖面积也反映了积雪量，同时也是融雪径流计算的中间变量，因此将空间分布的雪深信息作为SRM 模型的输入，进行雪盖面积、积雪量、融雪量间的转换计算是一种可行的解决途径。首先由空间分布的雪深信息计算雪盖面积、雪盖率，并将其代入SRM 模型计算融雪径流；而后用融雪径流计算所得的融雪深度计算雪深变化，进而得到雪深变化后的雪盖率，此雪盖率即为SRM 模型下一时段的输入。依次循环，即可计算逐日的融雪量、融雪深度和积雪覆盖面积。实测融雪径流和雪深空间分布信息是计算过程中模型参数调整校核的依据。

综上，为了减少因积雪信息不足造成的融雪径流预报误差，本文采用MODIS 遥感信息反演积雪深度，用雪深信息作为SRM 模型的输入，并在模型中增加逐日雪深和雪盖率计算模块，最终以实测径流和雪深数据作为融雪径流模拟的校核依据，得到雪深、积雪覆盖和径流的逐日变化过程。该方法在新疆某河A水库以上流域（以下简称A流域）进行了试验应用。

1 研究区概况

A 流域为狭长形，集水面积2 576 km2，高程范围1 891～4 614 m，流域上游高海拔区分布有永久性积雪和冰川；洪水成因主要是河源区高山带冰川及永久性积雪的融冰雪水，年内最大洪峰主要集中在6-7月，洪峰过程一般为5～10 d。夏季降雨多集中在上游河段，洪水过程可分为消融型洪水和消融与暴雨叠加型的混合洪水；全年过程多为一个融雪型主峰，在此基础上叠加若干短时段的暴雨洪水。

图1 A流域DEMFig.1 DEM of A basin

流域内布设有水文自动测报系统，包括5 个雨量气温站（A1～A5）、2个雨量气温雪深站（S1、S2）和1个水库站，雨量气温站大多布设在河谷地区，测站分布见图2，图中测站旁标注数字表示测站位置的高程。采集的水文气象数据为2016-2020年，数据时段长为日；其中仅2019、2020年有实测雪深数据。

图2 A流域测站分布图Fig.2 Station distribution in A basin

2 原理方法

2.1 MODIS反演雪深

已发布的MODIS 逐日和合成积雪资料均为积雪覆盖信息［16-19］，无雪深数据。为获取空间分布的雪深信息，本文基于傅华等［20］研制的MODIS积雪遥感监测系统，集合地表不同覆盖物的光谱特征，结合季节、地形、下垫面等积雪深度的影响因素，利用MODIS 高光谱、多波段反射率和测站观测雪深，采用数学统计方法，建立MODIS积雪深度回归模型。从而在计算积雪面积、覆盖率的基础上，获取空间分布的雪深、雪水当量信息，制作500 m 分辨率的积雪深度产品。因受云层影响，雪深为晴空日或10日合成产品［20，21］。

积雪深度的反演主要依赖于积雪深度与可见光波段的反射率间存在较好的相关性［21-24］：当积雪深度小于20 cm 时，雪面反射率随积雪深度的增加而增加；当积雪深度大于20 cm 时，雪面反射率随积雪深度增加的趋势减缓，当积雪达到一定深度后，雪面反射率趋于饱和。

当积雪深度相同时，下垫面类型不同，反射率也存在差别，因此进行积雪深度反演时，必须要考虑下垫面的差异，根据积雪深度反演原理，用数理统计方法（逐步回归判别、贝叶斯判别等），建立积雪深度与影响因子之间的线性关系［21］：

式中：S为积雪深度；x1，...，xn为影响因子，根据雪深反演原理定义为MODIS 多个特定通道反射率的函数；a1，...，an+1为回归系数，通过历史资料的统计分析确定。因影响因子较多，本文仅针对A 流域山地地形，依据土地覆盖类型（森林、灌木、草原、荒漠等）将其进行分类处理，每一类单独进行回归计算；计算结果结合坡度、坡向等对积雪分布及反射率的影响进行订正，从而提高反演精度［25，26］。

2.2 基于雪深的SRM模型

因不同垂向高程带的融雪径流特征差异明显，SRM 模型计算中通常将流域划分为若干区（高程带），每区单独计算融雪径流过程，再汇流至流域出口断面进行叠加。SRM 模型的大致原理为度日法［26］，即基于积雪覆盖面积和度日因子逐日计算融雪径流和降水径流，并将它们叠加到当日的退水流量上，得到每日的日径流量［4-14］。本文将积雪覆盖率替换为雪深，则基于雪深的SRM融雪径流模型的计算公式如下：

式中：Q为日平均径流量，m3/s；C是径流系数，其中Cs、Cr分别为融雪和降雨的径流系数；α为度日因子，表示单位时间、温度的融雪径流深，cm/℃·d；T为度日数，℃·d；△T为计算区域平均高程与气温站高程不同而产生的气温调整值，℃·d；P为日降水量，cm；K是退水系数，表示在没有融雪或降水时期逐日径流的下降率；n为径流计算的日数序列；常数项是径流深到径流量的转换系数［25］。遥感反演雪深产品的空间分辨率为500 m 网格，As为计算区域内雪深>0 的所有网格面积，A为计算区域总面积，As与A的比值为无量纲数，代表了区域的积雪覆盖率。

对于10日合成雪深产品而言，假定某产品为第n日的雪深，则下一个合成产品为第n+9日雪深，两个产品的间隔期内无雪深数据，因此需要以合成雪深为初始值，依据降水、气温等信息，采用公式计算得到间隔期内逐日雪深。

雪深计算首先需判别当日（第n日）气温是否超过临界温度：低于临界温度时，降水通过雪水当量转化为雪深，无融雪，雪深增加；高于临界温度时，积雪消融，雪深减小。积雪深度计算公式为：

式中：Hn为第n日雪深，cm；Hm，n为第n日积雪消融的深度，cm；Hs，n为第n日降雪深度，cm。

降雪深度和消融的积雪深度均通过相应的水量换算得到，计算公式为：

式中：Rm，n为降水量或径流量，m3，为SRM 模型计算量；ρw、ρs分别为水、雪密度，g/cm3。

通常，新下的松软的雪密度为0.04～0.1 g/cm3，融雪期的雪密度为0.6～0.7 g/cm3，雪密度的影响因素较多，新疆地区雪密度时空分布可参见黄慰军等研究成果［28］。因融雪径流量、积雪深度互为相关计算值，因此需借助合成的反演积雪深和实测径流量校核，由融雪径流模型迭代计算得到逐日的雪深和径流过程。

3 结果分析

3.1 流域分区

因气温随高程变化，而气温是融雪径流的重要影响因素，因此需要对垂直高程大于500 m的流域进行高程带划分。基于ASTER GDEM-Global Elevation Data 的DEM 数据（空间分辨率30 m），采用ArcGIS 软件的Arc Hydro Tools 工具包生成数字流域，提取相应的地形信息；在此基础上采用ArcGIS 的空间分析工具条Contour 进行高程带划分，每500 m 左右分为一带。A 流域依据500 m等高线进行整合，分为4个高程带，见表1、图3。

图3 A流域高程分带图Fig.3 Elevation distribution in A basin

表1 A流域高程带信息表Tab.1 Elevation information in A basin

3.2 MODIS反演雪深分析

MODIS 反演的积雪深度产品为二进制的无符号短整型格式（unsigned short），头文件后缀为*.HDR；数据文件采用等经纬度投影；A 流域反演数据行列数为：columns=1 518，lines=548；数据范围为0～255，其中0～49 代表雪深值，cm；250～252 代表水体（雪深为0），253、254代表云层（无效信息）。将雪深产品与研究区进行空间叠置，可得流域内500 m分辨率的雪深分布信息，图4为2019-01-08 A流域雪深空间分布图。

图4 遥感反演雪深分布图（2019-01-08）Fig.4 Distribution of retrieved snow depth（2019-01-08）

在具备实测雪深数据的2019、2020年，选取积雪期内的MODIS 反演雪深产品，从中提取雪深站位置的数据，依据积雪等级和判识标准划分表［24］，将其与图2 中S1、S2 雪深站的实测雪深进行对比分析，结果见表2。

表2 可见，2019、2020年遥感反演雪深数据的准确率均高于85%，遥感反演的年积雪日平均雪深与实测平均值的误差均小于10%，2020年误差在2%以内，表明雪深反演精度较好。

表2 遥感反演与实测雪深对比结果统计Tab.2 Comparison of the snow depth between RS retrieved and measured

在遥感反演雪深的基础上，对雪深的空间分布信息进行处理，将空间网格积雪数据与流域分区叠置，计算各分区的积雪相关信息，包括雪深面积分布、平均雪深、积雪覆盖率等。此处的平均雪深指有积雪覆盖区域的平均雪深，无积雪区不统计在内，因此雪量应表示为平均雪深和雪盖率的函数，对于非稳定积雪期的3-6月，仅依据平均雪深无法确定分区的雪量和积雪融水量。处理后流域各分区2016-2020年平均积雪深度和雪盖率变化过程见图5、图6，由图可看出，①雪深和雪盖率呈现明显的年际变化，其中雪盖率的年际变化规律更明显；②雪深、雪盖率与高程均呈现一定的相关性，其中雪盖率随高程变化的相关性更强；③融雪期的雪盖率为与高程相关的缓慢变化的过程，而积雪期的雪盖率变化迅猛，且与高程关系不大；④较雪盖率而言，积雪区内不同高程带的平均雪深增加了随时间跳变得无序性。

图5 各区2016-2020年平均雪深变化过程图Fig.5 Variation process of average snow depth in 2016-2020

图6 各区2016-2020年雪盖率变化过程图Fig.6 Variation process of snow cover in 2016-2020

3.3 模型参数确定

基于雪深的SRM 模型参数与原基于雪盖的SRM 模型基本相同，参数确定方式也类似，部分参数采用RS、GIS 手段或直接通过历史资料推求得到［1，4-9］。

气温直减率可表示为气温随高程变化的函数。在对应的高程区域内，根据不同高程测站的历史气温数据，采用γ =公式推求。其中γ 为气温直减率，表示每100 m高程气温下降的度数，单位为℃/100 m；h1、h2和T1、T2分别为两个测站的高程和历史平均气温。依据2016-2019年实测气温数据计算各高程带的气温直减率，当1 个高程带包含2 个以上测站时，首先计算每相邻高程两测站间的气温直减率，再依据高程差进行加权计算，最终计算所得各高程带气温直减率见表3。

表3 不同高程带气温直减率Tab.3 The lapse rate of air temperature in elevation bands

无降雨和融雪径流情况下，退水系数可表示为径流量随时间变化的函数。退水系数为第n天的实测径流量，x、y为需要确定的常数，可根据退水段的历史实测流量资料，构建基于Qn和Qn+1的双对数散点图推求得到［4，5］。依据A站2016-2019年退水期日径流绘制退水过程散点图（图7），根据拟合公式选取两点，Q1=100 m3/s，Q2=250 m3/s 代入上式计算得x=906，y=-0.003 8。

图7 A站2016-2019年退水过程散点图Fig.7 Scatter plot of recession flow for A Station，2016-2019

基于雪深的SRM 模型中的融雪径流系数CS、降雨径流系数CR和度日因子α无法通过历史资料直接确定，需要将降水、气温、蒸发和雪深的时间序列数据作为模型输入，在设定初始参数的情况下，通过模型计算得到各高程带逐日雪深和流域出口的流量过程，经过与实测雪深和径流对比，进行参数的不断调试，使计算和实测流量过程拟合最优，从而确定最终参数［29，30］。

3.4 模拟精度分析

以MODIS 反演雪深为初始值，采用基于雪深的SRM 模型计算逐日雪深和径流过程。以径流误差、Nash系数R2作为计算径流与实测径流吻合程度的评价指标，以反演雪深作为计算雪深的校核标准，对A流域2016-2020年径流过程进行模拟，其中2016-2019年为模型率定期，2020年为检验期。2016-2020年A 流域径流模拟结果见表4（2020年数据截至12月10日），2019、2020年径流模拟过程见图8。结果表明，基于雪深的融雪径流模拟精度较高，基本反映了流域径流过程，其中2016-2020年模拟年径流误差均小于7%，各年的径流拟合Nash系数都大于0.85。

表4 径流模拟结果Tab.4 Simulation results of snowmelt runoff

图8 模拟径流与实测径流对比图Fig.8 Comparison of simulated runoff and measured runoff

融雪径流计算中雪深和雪盖率的准确性直接影响融雪径流的模拟精度，因雪深空间分布的不均匀性难以刻画，本文依据雪深的空间分布计算区域雪盖率，并在融雪径流模拟中依据遥感反演信息对计算雪盖率进行校核。图9为A流域不同高程带2019、2020年3-6月融雪期的遥感反演和模拟计算雪盖率变化过程，该图清晰反映了高程对积雪覆盖率的影响，同时可看出经校核的逐日计算雪深较遥感雪深更精细地反映雪盖率随时间的变化。

图9 各区计算雪盖率与反演雪盖率对比图Fig.9 Comparison of calculated snow cover and retrieved snow cover

4 结 论

为了减少MODIS 合成雪盖产品因时间分辨率不足造成的融雪径流预报误差，本文以新疆A流域为研究对象，在建立MO⁃DIS积雪深度回归模型、反演积雪深度的基础上，将空间分布的雪深信息作为SRM 模型的输入，并在模型中增加逐日雪深和雪盖率计算模块，最终以实测径流和雪深数据作为融雪径流模拟校核的依据，模拟计算A流域2016-2020年雪深、积雪覆盖和径流的逐日变化过程，确定SRM模型参数。

A流域MODIS反演雪深与2个实测站雪深相比，2019、2020年准确率分别为85.29%和88.00%，遥感反演的日平均雪深与实测雪深的误差小于10%。基于雪深的融雪径流模拟在2016-2020年模拟年径流误差均小于7%，最高的为2020年检验期，误差小于1%；同时每年的径流拟合Nash系数都大于0.85。A 流域试验应用结果表明，MODIS 反演雪深和基于雪深的融雪径流模拟精度较高，MODIS 反演雪深有助于准确刻画流域积雪的时间、空间分布情况，采用基于雪深的SRM 模型计算，可更精细地反映积雪随时间的变化过程，同时促进融雪径流预报精度的提升。

积雪信息的时空分辨率和准确性是影响融雪径流预报的关键，因受云层影响，基于MODIS 信息反演的积雪产品时间分辨率较低，而基于被动微波遥感的积雪产品空间分辨率低，将两者进行时空尺度上的融合可为融雪径流预报提供更精确的积雪信息数据源，将是遥感和融雪径流预报的发展方向。 □