雪深和雪水当量被动微波反演及应用进展

王泽坤， 甘甫平， 闫柏琨， 李贤庆， 李和谋

(1.中国矿业大学(北京)煤炭资源与安全开采国家重点实验室，北京 100083； 2.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083； 3.中国自然资源航空物探遥感中心，北京 100083)

0 引言

雪深和雪水当量不仅是积雪观测中最重要的要素，也是冰冻圈研究的重要组成部分，同时还是地球系统中非常重要的变量[1-2]。其时空分布特点对全球变化、地球系统科学、水循环与水资源以及雪灾和洪涝监测起着非常重要的作用[3-4]。

目前对雪深和雪水当量等积雪要素的观测和反演主要采用野外现场实地测量、地面台站长期观测和卫星遥感区域观测等手段。全球积雪区域主要位于中高纬度地区、南北两极以及高山地区，空间异质性较强。传统观测手段和地面站点获取的数据具有“点尺度”的特点而存在较大的局限性[5]。卫星观测具有长时间、大区域、动态监测全球范围积雪的能力[6-7]，即“面尺度”的特点。对 “点尺度” 与“面尺度”信息的协同观测使区域性的雪深和雪水当量反演以及定量分析成为可能[4-5]。卫星遥感积雪监测方式主要有光学雪覆盖监测及微波反演雪深和雪水当量2种，其中被动微波遥感是现阶段最有效的雪深和雪水当量反演方式[8-11]。

基于地面观测、卫星观测等积雪数据，国内外学者对积雪的微波辐射特性与雪深、雪密度、雪粒径、积雪结构、液态水含量和下垫面特性等内在关系和规律进行深入研究，建立了诸多积雪微波辐射传输模型[12-16]，发展了雪深和雪水当量的多种反演方法[8,17-22]。由静态到动态、简单到复杂的递进反演模式发展出可以应用于气候变化、水文监测等领域的雪深和雪水当量产品，其中作为三大稳定积雪区之一的青藏高原地区一直是研究的重点。

本文从数据源获取、算法和模型反演、主要应用等方面对被动微波反演雪深和雪水当量的发展动态进行了总结，并对其未来趋势做出展望。在综合分析前人研究和应用的基础之上，归纳了积雪数据获取的主要方式及观测的主要参数，为雪深和雪水当量反演提供数据支撑； 梳理了前人对雪深和雪水当量反演算法的研究，评价了不同方法的适用条件，为不同地域的应用提供了算法选择； 以青藏高原为例，探讨了雪深和雪水当量反演主要应用进展，并为今后积雪被动微波反演及应用的深入开展提供了科学的参考建议。

1 积雪观测方式与主要参数

根据目前积雪数据的3种获取方式，野外现场实地测量、地面台站长期观测得到的参数不仅可以直接用于积雪观测应用与理论研究，而且可以用于遥感反演算法建立、真实性检验与结果评价，还可以作为可靠的积雪特性输入信息进行模型模拟。卫星遥感区域观测获得的数据包括光学遥感数据和微波遥感数据，可以提供区域或全球的积雪信息。不同观测方式得到的主要参数见表1。

表1 积雪观测主要参数和相关仪器手段[6,23-24]Tab.1 Main parameters and related instrumentsof snow observation[6,23-24]

1.1 野外现场实地测量

对野外积雪实地观测，主要是结合实际应用需求，采用区域与路线相结合的方式，利用相关仪器在实地现场对积雪物理参数进行实时测量，为积雪特征研究、卫星-地面结合参量反演、真实性检验以及参量反演结果评价等提供同步或准同步数据。野外实地观测，可分别在点、线、面尺度上进行以积雪物理参数观测为主的野外综合试验，建立更加规范统一的积雪特性数据集和时空动态变化数据库，从而系统地开展相关的理论、方法等研究[23]。但由于受野外条件等限制，野外现场实地测量难以获取长时间动态的监测数据，观测区域尺度相对受限。

1.2 地面台站长期观测

地面台站是长期、定点且直接获取积雪信息的有效观测手段。利用站点仪器对积雪物理属性进行不间断的自动化观测，为积雪观测与应用提供连续动态的实时数据。中国气象局在全国765个气象基准站实施了大范围雪深观测[25]，此后陆续在2 400个气象站上开展积雪观测业务[26]，推动了积雪观测与应用的研究和要素反演技术的发展。但受地形地貌、交通、人类生活条件以及设备运行等限制，地面台站观测积雪数据时空分布不连续、不均匀，部分数据区域代表性有限，难以全面反映积雪的整体分布特征[27-28]。

1.3 卫星遥感区域观测

卫星遥感技术可以提供大范围的积雪信息监测。相较于光学遥感，微波遥感显示出了其在积雪监测中的绝对地位[6,24]，其中被动微波遥感尤甚。表2是国内外雪水当量反演相关被动微波传感器特

表2 相关被动微波传感器特征参数[29]Tab.2 Related passive microwave sensorcharacteristic parameters[29]

征参数的汇总[29]。被大量使用的有1978年发射的Nimbus-7卫星携带的多通道扫描微波辐射计(scanning multichannel microwave radiometer，SMMR)[17]提供的微波亮温数据，以及微波辐射成像仪(special sensor microwave imager，SSM/I)[30]、微波成像仪/探测器(special sensor microwave imager/sounder，SSMIS)[31]、高级微波辐射计(advanced microwave scanning radiometer for the earth observing system，AMSR-E)[22]和第二代高级微波辐射计(advanced microwave scanning radiometer-2，AMSR-2)[32]等传感器获取的数据。风云三号(FY-3)A星搭载的微波成像仪(microwave radiation imager，MWRI)[33]是我国首个业务化运行的星载被动微波辐射计，并在积雪监测等方面取得了较好的应用。

2 雪深和雪水当量被动微波反演算法

基于野外现场实地测量、地面台站长期观测和卫星遥感区域观测等数据，在对积雪物理特征和辐射传输模型等研究的基础之上，研发了卫星-地面结合的半经验算法、基于积雪辐散射特性的物理模型和机器学习算法等雪深和雪水当量被动微波反演算法。不同类型的算法适用于不同的积雪环境，其反演的效果、精度也不一样，目前还没有发展出卫星像元尺度上普遍适用的雪深和雪水当量反演算法。

2.1 半经验算法

半经验算法是基于雪深和雪水当量与亮温梯度之间的线性关系拟合得来的，由于系数获取的方式不同，可以分为静态算法和动态算法。

2.1.1 静态算法

静态算法是国际上应用最早的积雪反演算法，主要介绍Chang算法、NASA96算法、用于中国区域修订的Chang算法等。后面2种算法是在Chang算法的基础之上，针对指定的应用目标、覆被类型和特定区域进行修订而成。

2.1.1.1 Chang算法

Chang等[17]利用均质积雪辐射传输模型计算了积雪覆盖的向上微波辐射亮温，得到针对特定载荷SMMR数据的反演算法。通过假设积雪粒径为0.3 mm，积雪密度为0.3 g/cm3，建立亮温模拟结果与站点实测结果的线性回归方程为：

SD=1.59(TB,18H-TB,37H)，

(1)

式中：SD为雪深；TB,18H和TB,37H分别为18 GHz和37 GHz水平极化亮温。

总体而言，该算法对于大尺度雪深和雪水当量估算效果较好，但有时会低估积雪。这是由于当SD估计值小于2.5 cm时考虑到SMMR 25 km×25 km像元尺度内的异质性，则认为地表没有积雪；SD大于100 cm时会出现亮温差饱和，即亮温差不再随雪深变厚而增加，所以该算法普遍适用于SD在(2.5，100)cm范围的情况[11]。

2.1.1.2 NASA96算法

植被的存在会削弱积雪的微波辐射，同时其自身的辐射又会叠加到总的辐射亮温中。如果忽略植被的影响，反演算法将会低估森林地区的雪水当量[17,34]。Foster等[18]在Chang算法的基础上引入森林覆盖度来提高雪深的反演精度，修正后的NASA96算法为：

SD=a(TB,18H-TB,37H)/(1-f)，

(2)

式中：a为常数，通常在欧亚大陆取0.78，在北美和其他地区取1.59；f为森林覆盖度。由于森林在18 GHz和37 GHz发射率相似，亮温差仅由无森林覆盖的区域(1-f)贡献。

由于该算法忽略了雪粒径和积雪类型的差异，Foster等[35]又将森林覆盖率和雪粒径这2个动态变化参数引入到反演算法中。新算法能很好地反映雪的积累和消融时相，但是针对高山型、海洋型积雪进行反演时的误差依然较大。

2.1.1.3 中国区域反演算法

由于中国区域的积雪密度相对较小，青藏高原等地区的积雪下覆深霜层较为发育，所以在欧美地区发展的反演算法并不完全适用于我国。前期车涛[19]基于SMMR和SSM/I传感器对Chang算法进行了初步修订，发展了针对中国区域雪深反演的新算法，即

SDSMMR=0.78(TB,18H-TB,37H)+b，

(3)

SDSSM/I=0.66(TB,19H-TB,37H)+b，

(4)

式中：TB,19H为19 GHz水平极化亮温；b为补偿值，随月份变化。

之后Che等[20]通过剔除湿雪、冻土、降水和寒漠等干扰像元，对中国区域的SMMR和SSM/I雪深算法系数做了进一步修正，得到精度更高的算法，计算公式为：

SDSMMR=[0.78(TB,18H-TB,37H)+b]/(1-f) ，

(5)

SDSSM/I=[0.66(TB,19H-TB,37H)+b]/(1-f) 。

(6)

考虑到积雪覆盖度对雪深反演的影响，Jiang等[21]利用2002—2009年的AMSR-E微波亮温和地面台站数据建立中国区域不同下垫面的雪深反演算法，即

SD=fgrassSDgrass+fbarrenSDbarren+fforestSDforest+ffarmlandSDfarmland，

(7)

式中：fgrass，fbarren，fforest和ffarmland分别为每个像元内草地、裸地、森林和农田的覆盖度；SDgrass，SDbarren，SDforest和SDfarmland分别为相对“纯”像元(主要地物覆盖度大于或等于85%)建立的雪深统计反演算法。

该算法对浅雪反演效果较好，但往往会低估深雪。目前已成为我国FY-3卫星雪深和雪水当量反演的业务化算法。

2.1.2 动态算法

积雪物理属性的时空特征每时每刻都在发生变化，而静态算法却忽略了这种变化。将此影响引入到微波亮温的计算中，便可得到更符合实际积雪情况的动态反演算法。主要包括温度梯度指数(temperature gradient index，TGI)动态算法和Kelly(2009)动态算法等。

2.1.2.1 TGI动态算法

Josberger等[8]利用TGI来描述积雪上温度梯度的累积效应，表示了积雪粒径增长的一个累积变化指标，即

(8)

式中：TGI为温度梯度指数； C为临界温度梯度值，是一个比例常数，20 ℃/m；Tground为地表-雪界面的温度，即雪土界面温度；Tair为空气-雪界面温度，即雪表面温度；SD(t)为t时刻积雪深度。

研究发现去除季节末的湿雪后，亮温梯度(spectral gradient，SG)与TGI之间存在很好的线性关系，即

SG=(TB,19H-TB,37H)=αTGI+β，

(9)

式中：SG为亮温梯度；α和β为经验系数，由SG和TGI的关系回归得到。

通过对上式求导，则求得t时刻雪深SD(t)为：

(10)

虽然TGI算法考虑了积雪物理属性的变化过程，但是它只有在Tground-Tair>0时才有意义，并且需要连续的温度梯度和光谱梯度数据来确保积分的正常运算。

2.1.2.2 Kelly(2009)动态算法

该算法主要用于AMSR-E和AMSR-2雪深和雪水当量产品生产中[22]，除了考虑不同频率对不同厚度积雪的敏感性之外，还考虑了森林的影响，采用线性混合方程求解森林区和非森林区混合像元的雪深。公式为：

SD=fSDf+(1-f)SD0，

(11)

式中SDf和SD0分别代表森林区和非森林区雪深。

森林区和非森林区的雪深计算公式分别为：

(12)

(13)

式中：fd为森林密度；pol19和pol37分别为19 GHz和37 GHz的亮温极化差；TB,10V，TB,19V和TB,37V分别为10 GHz，19 GHz和37 GHz垂直极化亮温。

该算法忽略了雪粒径对反演的影响，Tedesco等[36]在此基础上进一步将雪粒径函数作为动态系数，引入冻土分布数据来表征其对亮温和雪粒径的影响，反演效果大大改进。

总体而言，半经验算法是目前全球业务化被动微波雪深反演的常见算法，它实质上是一种线性亮温梯度法。但是该类算法缺乏关于积雪特性参数(积雪粒径、积雪密度等)对微波信号干扰的研究，存在一定误差。一般只在特定区域上表现稳定且精度良好，而在不同研究区内会得到不同的反演结果[17-18,20-22,37]，在实际应用中，需要有针对性地选择不同的算法。

2.2 基于物理模型的反演算法

雪密度、雪粒径等积雪特性和下垫面的变化都会影响积雪微波辐射[17,38-39]，微波信号对积雪层具有较高的穿透力，并且随着雪深增加微波亮温表现出减小的趋势[40]。在微波低频波段，下垫面特性是影响干雪辐射的主要因素，而在高频波段，积雪颗粒的体散射则起了关键作用。这种随着频率增加散射作用增强、积雪辐射减弱的特点，可以用来探测地表雪深和雪水当量。积雪特性参数的时空演化对雪水当量反演精度影响很大[41]，许多学者为了降低这种演化带来的误差，发展了基于物理模型的反演算法。Jiang等[42]通过建立积雪辐射模拟数据库，发现雪水当量与积雪辐射和衰减信号之间有较好的统计回归关系，进而利用该关系反演雪水当量。该算法是基于纯积雪物理模型发展而来的，积雪微波辐射传输模型可以准确地模拟积雪覆盖地表的微波辐射散射特性，是解释观测数据和提取积雪参数的重要途径。

HUT(Helsinki University of technology)模型[12]是一种零阶半经验模型，其散射相函数被简化为只有前向散射，可采用一个经验系数q=0.96来度量被散射到前向的辐射能量比例，考虑了森林覆盖度和大气的影响。MEMLS(microwave emission model for layered snowpacks)模型[13]是一种可以模拟多层积雪辐射特性的模型，采用六流技术描述积雪层内多次散射和吸收过程，适用于5～100 GHz频率范围。Pan等[43]发现MEMLS模型在深雪时的模拟效果总体上要优于HUT模型，还有研究发现该模型可以模拟出雪深50 cm时的亮温差饱和现象。Tsang等[14]提出致密介质传输理论模型(dense media radiative transfer，DMRT)，应用其计算积雪中电磁波的有效传播常数，并据此修正积雪的消光系数、散射系数和反照率。之后Du等[15]提出了DMRT-AIEM-MD(dense media radiative transfer-advanced integral equation model-matrix doubling)微波辐射模型，改进了粗糙界面的建模过程。

积雪模型可以反映积雪参数变化引起的亮温改变，同时较好地模拟积雪层内的微波辐射传输过程[44]。但由于它在一定程度上对积雪几何结构做了简化处理，且依赖的观测样本相对有限，从而限制了积雪模型在全球尺度上雪深和雪水当量反演中的应用。伴随着计算机技术的飞速发展，Bicontinuous模型[16]等需要大量计算资源的数值计算模型将有望成为未来积雪模拟的主流发展趋势[45]。

2.3 机器学习算法

由于简单的半经验线性模型简化了雪水当量的实际反演，而复杂机理模型又缺乏直接可操作性。国内外学者开始利用机器学习的方法来描述微波亮温与积雪参数之间的非线性函数关系，旨在提高反演精度。

2.3.1 人工神经网络(artificial neural network，ANN)

ANN模型是人类神经系统一阶近似的数学模型，目前已被广泛应用于地表参数估算研究中。Forman等[46]使用ANN对北美积雪区不同频率和极化方式下的亮温进行估算，结果显示该算法在大陆尺度数据演化中较为高效。在实际反演中，ANN对先验知识的依赖性很小，能较好地模拟微波亮温与积雪参数之间的非线性关系，但其训练算法和结构对结果影响很大，要求输入的参数必须相对独立，限制了该算法进一步的应用[47]。

2.3.2 支持向量机(support vector machine，SVM)

SVM在图像分类中应用较为广泛，很多学者研究发现，基于SVM方法的雪深反演算法精度较高[30]，对积雪属性的变化更敏感，可以在很大程度上代替神经网络算法[47]。Xiao等[48]运用该算法建立了基于不同植被类型和不同积雪期情况下的雪深反演模型，发现该算法不仅提高了反演精度，还在一定程度上减小了“积雪饱和效应”。

2.3.3 随机森林(random forest，RF)

RF模型是由Breiman和Cutler 在2001年提出的一种基于决策树的算法，它通过对大量决策树的汇总提高了模型的预测精度。RF不需要顾虑一般回归分析面临的多元共线性的问题，不用做变量选择，有清楚的规则并且可以训练得更快[49]。Bair等[50]用ParBal能量平衡模型进行模拟，把雪水当量和亮温值输入到RF模型训练，利用RF模型实时反演山区雪水当量，相对神经网络算法精度有所提高。

2.3.4 马尔科夫链蒙特卡洛(Markov chain Monte Carlo，MCMC)

MCMC是实现多频率联合反演雪水当量的一种方法。Durand等[51]和Pan等[43]发展了基于贝叶斯理论的雪水当量反演算法BASE-PM(Bayesian-based algorithm for SWE estimation-passive microwave)，引入了雪粒径、密度和温度等相关先验知识，也考虑了多层积雪的情况。虽然该方法计算量十分庞大，但是对算法开发者要求较低，有望成为未来遥感参数反演的有力工具。

机器学习方法使用范围广、反演精度高，在一定程度上避免了线性算法在不同地区应用系数不同的局限[52]。然而，机器学习的准确性完全依赖于训练样本是否具有代表性[30]，缺乏物理模型的解释[38,49]，难以在全球算法中推广应用。

3 青藏高原雪深和雪水当量遥感反演

青藏高原作为世界上平均海拔最高的区域，降雪主要集中在当年10月—次年4月[53]，有着丰富的积雪和冰川资源，被誉为“世界第三极”和“亚洲水塔”。由于受多种因素影响，青藏高原的被动微波积雪观测存在很大局限性，已有的雪深和雪水当量反演算法在该地区误差较大[17]。国内外学者进行了一系列算法改进，得到了相应的雪深和雪水当量积雪产品，使其更好地应用于青藏高原，为该地区可持续发展提供重要的科学依据，同时有助于提升对全球气候变化研究和水文现象监测的能力。

1)改进了雪深和雪水当量等积雪要素的反演算法。由于地形、冻土和雪粒径等影响，观测数据的空间分辨率、时间分辨率的限制，以及卫星频谱的观测有效性的制约，需要有针对性地发展青藏高原雪深和雪水当量反演算法。针对站点观测值难以代表被动微波亮温数据像元(大约25 km×25 km)区域的积雪情况，构建一套有效的基于单位像元上地面点尺度观测转换为卫星面尺度观测的方法，是发展青藏高原积雪参数反演算法和真实性检验的重中之重[54]。李长春等[55]为了提高山区积雪深度的反演精度，基于FY-3B-MWRI亮温数据、土地利用和数字高程模型等辅助数据，分析了坡度、坡向等地形因子对微波遥感反演雪深的影响，构建基于地形效应的微波辐射模型，提出基于地形校正的山区积雪深度反演算法； 李晓静等[56]利用SSM/I数据对我国积雪判识方法做了进一步改进，不仅达到有效判识积雪的效果，还能给出相应积雪状态的定性信息，在很大程度上减弱了高原冻土对积雪识别的影响； Jiang等[21]引入了对浅雪敏感的高频通道89 GHz，用4个频率组合发展了FY-3B-MWRI雪深产品的业务化算法。此外，还建立了不同粒径的积雪端元数据库对提高青藏高原积雪产品反演精度也具有指导意义[57]。

2)发展了一系列雪深和雪水当量积雪产品。比如，Liang等[58]、Gao等[59]和Huang等[60]将MODIS 500 m×500 m每日积雪产品与AMSR-E每日雪水当量数据进行融合，大大提高了产品的精度。通过对被动微波数据进行降尺度，可以有效地解决混合像元带来的不确定性问题。Wang等[61]针对青藏高原雪深高估现象，将地理辅助数据(位置和高程)、光学遥感积雪日数(MODIS)和被动微波遥感亮温(AMSR-2)进行融合发展了一套综合多种因素方法的空间降尺度雪深反演算法，该算法使空间分辨率由原来的10 km降尺度为500 m，产品均方根误差(root mean square error，RMSE)由4 cm提高到2 cm。

3)得到了在气候变化、水文监测等领域的广泛应用。青藏高原雪水当量的时空变化对全球气候演变具有一定的指示作用[62]。Huang等[63]基于气象台站资料对中国三大稳定积雪区近70 a间的雪深变化进行了分析，发现整体上雪深有增加的趋势，而最大雪深往往出现在青藏高原地区。青藏高原积雪对亚洲气候的影响是通过季风实现的，总体来说呈负相关关系。Senan等[64]发现青藏高原地区在厚雪年会推迟印度洋夏季季风开始的日期，并造成印度地区降水减少。此外融雪径流蕴藏着体量丰富的淡水，直接影响第二年春季河川径流量变化[65]。随着全球气候持续变暖，青藏高原的积雪融水贡献率逐年递增[66]，雪深和雪水当量产品实时监测的直接反馈可以为洪水预警等提供宝贵的数据资料。

4 结论与展望

目前雪深和雪水当量等积雪观测已经实现了地面实测手段精准化、站点观测自动化和卫星观测多样化，为积雪参数的获取奠定了坚实的技术支持。基于被动微波遥感观测，研发了半经验算法、物理模型算法以及机器学习算法等方法，实现对雪深和雪水当量等积雪观测要素的区域性反演，支撑对全球气候变化、两极和青藏高原地区积雪特征认识、水资源管理等方面的研究。

在遥感研究呈现综合化和跨学科趋势的背景下，为了给积雪反演提供精度可靠的定量化遥感产品和数据支撑，综合分析已有的遥感理论成果和反演方法是当务之急。除了获取像元尺度的地面观测真值外，还要利用光学、微波传感器等多源遥感资料获取积雪特性的空间差异信息，建立积雪参数综合数据库，生成时空连续、更高精度的遥感产品，为雪深和雪水当量的反演提供更精确的卫星遥感观测数据。未来可以考虑从以下几个方面着手：

1)构建一套有效的点面升尺度转换方法，加强像元尺度上的地面连续观测与真实性检验。由于现有地面台站在空间上分布极不均匀导致其所测积雪数据的代表性也不强，而且进行野外人工积雪观测时不同国家和地区之间没有统一规范的积雪观测标准，难以形成比较系统的全球地面积雪观测数据，同时对卫星观测雪深和雪水当量的精度评估又存在偏颇性。因此，在充分利用先验知识的基础上，通过克里金等空间插值手段发展点面转换方法，结合采样点的观测数据实现最优估计，为几十千米级的像元尺度反演与验证提供方法支撑。

2)充分利用多源遥感数据，依靠积雪数据同化系统，发展雪深和雪水当量反演的最优算法。模型模拟和遥感观测都具有各自的优势和局限，如果将其结合并融合，依靠模型提供先验知识与相关理论，依靠多源(光学、微波等)遥感观测提供不同知识的要素，修改或发展模型与反演算法，使二者优势互补。在同化系统框架内融合积雪过程模型、积雪站点观测和多源卫星遥感观测等数据，提升反演的精度，预测未来雪深和雪水当量等积雪要素演变规律。

3)开发更高精度的雪深和雪水当量反演产品，提高积雪时空变化动态监测能力。在长时间序列产品生产中，ANN和SVM等机器学习算法的应用越来越广泛，为高级产品的生产提供坚实的技术支持。