基于NoSREx实验的积雪演化及微波辐射分析研究

周静恬，邱玉宝，黄 琳，Juha LEMMETYINEN，石利娟， 李青寰， 施建成

（1.中国科学院空天信息创新研究院数字地球重点实验室，北京 100094；2.可持续发展大数据国际研究中心，北京 100094；3.中国科学院大学，北京 100049；4.中国科学院空天信息创新研究院-芬兰气象研究所北极观测联合研究中心（JRC-AO），芬兰索丹屈莱FI-99660；5.芬兰气象研究所北极空间中心，芬兰赫尔辛基FI-00560；6.中国科学院国家空间科学中心，北京 100190）

0 引言

积雪是冰冻圈最为活跃的季节性要素，是全球气候变化的灵敏指示器［1］。针对积雪的监测具有重要的现实意义，当前被动微波遥感在积雪监测方面具有不可替代的作用。积雪微波辐射传输模型模拟可开展过程和机理研究，仿真环境可重现积雪演变、积雪与微波相互作用等重要过程，是被动微波积雪参数——雪深和雪水当量反演算法发展的基础。

针对积雪的物理特性是影响空间监测的重要因素，其演变特性存在较大的空间异质性［2-3］，积雪物理特性的演化影响微波辐射［4-5］，从而影响雪深/雪水当量反演算法［6-7］。当前计算雪深/雪水当量的算法主要是半经验线性算法［8-9］，其不确定性来源于雪物理特性变化及下垫面变化和大气变化影响［4］，其中雪的物理特性演变过程对算法影响较大［10］。由于积雪反演算法在时间上和空间上还存在很大的不确定性，导致精度受到影响［7，11-12］。通过改进算法参数来更好捕捉动态积雪演变过程，或耦合发射率模型和陆表模型，或耦合发射率模型和由天气驱动的积雪演变模型将很可能提高雪深/雪水当量反演算法准确性［13］。结合实验测量的雪深等辅助信息［14］或了解积雪的先验知识［15］有助于减少雪深/雪水当量反演算法的不确定性。

为应对上述挑战，论文采用在北欧地区的2009—2013年积雪地面和地基微波实验观测数据集［16-17］，研究不同积雪期的分期判断方法，探讨了芬兰积雪实验场积雪特性的时间序列变化以及地基微波辐射计观测亮温差与雪深变化关系，最后采用MEMLS模型［18-20］开展积雪微波辐射的模拟、对比和分析研究。

1 实验数据、模型和方法

1.1 实验数据

北欧积雪实验（Nordic Snow Radar Experiment，NoSREx）是ESA Earth Explorer 7候选任务CoReH2O［21］（寒冷地区水文学高分辨率观测）中Phase A研究组成部分，该实验目的是在北极寒带森林地区的代表性地点提供整个冬季连续时间序列的地面、主被动微波积雪观测。实验区域位于芬兰气象研究所北极研究中心Sodankylä的观测站，地理位置为67.3618°N，26.6338°E（图1）。实验采用了包括10.65、18.7、37、90 GHz四个频率的微波辐射计系统，其中在2011—2012年实验期间频段21 GHz替换了90 GHz，四个入射角分别为30°、40°、50°、60°，获得了水平和垂直极化下的观测亮温以及标准差。地面观测数据包括积雪参数的人工和自动测量，人工测量是指每隔3～5天进行雪坑（snow pit）观测，其参数包括雪深、密度和雪水当量、雪密度剖面、雪温度剖面、雪层评估，观测天数在2009—2013四年依次为91、31、23、32天。自动测量参数包括温度、露点温度、风速、气压、雪深、地表湿度和温度、气温、雪水当量等，每天共测量8次，每次间隔3小时，观测天数在2009—2013年依次为365、365、366、365天［16］。

图1 NoSREx IOA的网络摄像头图像以及主要微波仪器，SnowScat散射计和SodRad辐射计系统的照片［16］Fig.1 Webcam image of NoSREx IOA and photographs of main microwave instruments，the SnowScat scatterometer and the SodRad radiometer system［16］

该数据集已用于积雪相关模型模拟［22-23］、反演积雪参数［24-25］、改进算法［26］、耦合雪物理和辐射传输模型［27-29］、分析观测数据［29-31］等研究。

1.2 多层积雪微波辐射模型

多层积雪微波辐射传输模型（MEMLS）是由Wiesmann等［18-19］开发的一个针对多层积雪的被动微波辐射传输模型，适用频率范围为5～100 GHz。它以辐射传输为基础，将积雪根据物理特性分为多层，利用六流近似理论来描述每个雪层内部的多次散射与吸收，同时考虑了雪层之间的界面散射。由于模型在特定环境下开发，在不同环境下的适应性不同，为了将其有效应用到芬兰积雪试验场，在该研究中需要调整模型的输入参数地表均方根高度参数的值。

1.3 研究方法

1.3.1 积雪积累期、稳定期和消融期确定

根据气温和积雪厚度的变化，把积雪期划为3个雪期，其中积雪“积累期”为出现降雪后积雪逐渐累积的过程，表现为雪深随时间呈现显著增加趋势，气温基本已降低至0℃以下；积雪“稳定期”为气温连续低于0℃时，雪深随时间增加趋势减小，即积雪压实的过程，表现为雪深出现增加明显变缓或甚至下降趋势；积雪“消融期”为随着温度升高且连续5日大于0℃，积雪进入融化或消融，可能伴随再冻结过程，表现为雪深迅速减小，气温基本维持在0℃以上。

针对芬兰Sodankylä地区，根据其地面观测数据的温度和雪深组合特征（图2），研判3个时期开始和结束时间，方法为：积累期开始时间为雪深一周内连续大于5 cm，气温处于从0℃以上过渡到0℃以下阶段，此前气温不低于-20℃；积累期结束时间也是稳定期开始时间，为连续3个周内雪深变化范围为±5 cm以内，气温至少连续5天低于0℃；稳定期结束时间也是消融期开始时间，为连续5天气温大于0℃，雪深处于减小状态，一般为50～70 cm；消融期结束时间为雪深一周内连续小于5 cm，气温处于0℃以下过渡到0℃以上阶段。

图2 2009—2013年气温和雪深时序变化图Fig.2 Changes of temperature and snow depth from 2009 to 2013

1.3.2 模型输入确定

MEMLS模型输入包括频率、入射角、天空背景亮度温度、散射系数、积雪参数、地表温度以及积雪-土壤界面反射率［19］。天空背景的亮度温度所需的大气透射率参照Pulliainen等［32-33］的提出的统计反演方法，天空背景亮度温度所需的上下行大气亮度温度由Aschbache［34］提出的公式估算。散射系数参照Matzler等［18］的工作，采用玻恩近似计算。输入的积雪参数包括积雪层数、积雪温度、积雪湿度、积雪密度、积雪盐度、积雪厚度，指数相关长度，指数相关长度采用Durand等［35］的经验公式通过积雪粒径估算，其他参数由NoSREx实验的地面测量获得［16］。陆地积雪的盐度通常认为是0［36-37］，本次模拟也将积雪盐度设为0。地表温度由NoSREx实验的地面测量数据获得。水平极化和垂直极化的积雪-土壤界面反射率，本研究是利用Wang&Choudhury（1981）半经验模型（简称QHN）计算［38-39］。QHN模型中所需的土壤介电常数采用Dobson模型计算［40］。土壤-雪界面的粗糙程度通过土壤粗糙度来反映［41-42］。土壤粗糙度通常用均方根高度和表面相关长度两个统计变量表示，其具体定义可见文献［43-44］。土壤粗糙度较难直接测量，通常采用最小化成本函数（cost function）［45-48］的方法估算。本研究选用的代价函数如下：

式中：CF为代价函数的值；m为模拟次数；TM为测量亮温；TS为模拟亮温。

2 积雪演化特征分析

2.1 不同积雪期的判定

气温和雪深随时间变化如图2所示。气温日波动性和月波动性较大，年波动性相似，一般在11—12月降低到0℃以下，在4—5月上升到0℃以上；每年雪深的变化趋势相似，一般从10—11月逐渐累积，在3—4月达到最大值，在4月下旬至5月迅速融化。根据图2以及本文1.3.1章节的分类标准，判定芬兰Sodankylä的积雪期分为积累期、稳定期和消融期，积累期一般为10月—次年2月，时长约为4～5个月；稳定期为2—4月，时长约为2个月；消融期为4—5月，时长约为1个月，其具体结果时间如表1。

表1 不同雪期分类的时间（年-月-日）Table1 Time classification of different snow peroid（yyyy-mm-dd）

2.2 雪颗粒形状变化

图3展示了自然积雪的分层特性。如图可知，融化状态（Melt Forms，MF）主要出现在积累期早期（10—12月）和消融期（4—5月）。降水粒子（Precipitation Particles，PP）和分解碎片降水粒子（Decomposing and Fragmented precipitations particles，DF）主要出现在近雪表层（图中均为绿色），且主要在4月前出现。在积累期后期和稳定期，积雪分层颗粒形状出现的种类主要是圆形颗粒（Rounded Grains，RG）、片状颗粒（Faceted Crystals，FC）、深霜（Depth Hoar，DH）。在消融期主要类别为融化状态。图3（c）可见深蓝色较少，表明2011—2012年缺乏深霜层，这与该年较其他年份更温和的气温有关，色带参考季节性陆面积雪的国际分类［49］。

2.3 雪粒径变化

结合图3～4可知：（1）在整个雪季，研究区域的垂直剖面上每个积雪层粒径不同，一般底层积雪粒径最大且为深霜层，即深霜层粒径最大，底层粒径会从小变大再变小，表层粒径一直较小且变化不大。（2）在积累期早期，10—12月，雪深较浅时，雪的粒径也较小，一般在1.5 mm以内。（3）随着雪深的积累，在积雪底层的粒径逐渐增大。在积累期早期10—12月从较小的0.25～1 mm，在积累期后期（1—2月）增长为1～3 mm。（4）随着雪深增加速度减缓，在积累期后期和稳定期（1—4月），近地表的粒径维持在2 mm左右，均出现了增加到至少2.5 mm的现象。（5）在融雪期4—5月，观测雪层的颗粒形状大多数为融雪状态，积雪融化，粒径相较稳定期的表层大底层小，为1 mm左右。（6）靠近积雪表面的粒径值始终较小，基本维持在1.5 mm以内，但少数情况下会出现较大值。（7）观测的平均粒径的最大值出现在每年的稳定期，2—3月，值为2.5～4 mm，均出现在近地表层。

图3 2009—2013年自然积雪的分层特性Fig.3 Stratification characteristics of snow cover from 2009 to 2013

3 微波辐射亮温模拟分析

3.1 微波亮温差对雪深变化的依赖分析

为研究线性亮温梯度算法（18 GHz和37 GHz，V和H）在北欧实验区域的适用性，分析了微波亮温差对雪深变化的依赖性（图5）。积累期（12月—次年2月）和稳定期前期（2—3月），雪深和亮温差（18～37 GHz）整体来看具有一致的变化趋势，局部来看具有相反的波动性，即雪深小幅度减小亮温差反向增加，如图5中矩形框中的部分：黑线和锯齿形状的红色和蓝色曲线反向波动，这可能是由于融雪使得亮温差减小，而后积雪重新冻结亮温差增大。在稳定期后期（3—4月），由于3月雪深均大于60 cm，亮温差对雪深增加不再敏感趋于饱和，体现了算法的局限性［21］；局部反向现象较之前不明显，这可能是由于早期融雪形成的雪壳结构逐渐松弛为更典型的冬末降雪［31］，雪粒径和雪颗粒形状均改变。在消融期（4—5月），随着积雪迅速融化雪深减小，亮温差值波动性较大基本维持在±15 K以内，亮温差与雪深线性相关性不明显。

图5 2009—2010年雪深和亮温差（18 GHz和37 GHz）在50°入射角条件下随时间的变化Fig.5 Variation of snow depth and brightness temperature difference（18 GHz and 37 GHz）with time at 50°incident angle from 2009 to 2010

雪深的变化和亮温差变化在积累期和稳定期总体上具有相似性，但局部来看没有明显的一致性，这可能是由于积雪亮温还会受到如含水量、积雪粒径以及雪颗粒形状等的影响，而这些参数均随着时间变化。由上一章可知整个雪季积雪演变较大，因此需在算法中考虑随时间演化的积雪特征，已有研究考虑积雪动态变化的算法［8］。微波亮温差对雪深和线性函数在不同时期依赖性不同，进一步表明提高反演算法精度需考虑参数的动态变化。

3.2 MEMLS模型参数值调整

模型模拟所需要的土壤粗糙度参数较难直接测量，且其对模型模拟结果有一定影响，在模拟前需要确定土壤粗糙度参数的值。地表均方根高度的值由土壤本身的性质决定，不同情况下该数值可以在数毫米到几十毫米之间变化［39］。采用1.3.2章节所述方法，得出的不同频率和入射角的最优值（表2），大多数情况下，在值为0.001 m时均方根误差最小，因此在模拟时将垂直和水平极化下的地表粗糙度值均设为0.001 m。潘金梅等［18-20］的研究表明Sodankylä地区的地表均方根高度大约为0.001 m，与本研究的估算结果吻合。

表2 地表均方根误差参数最优值Table 2 Optimal values of surface root mean square error parameters

图4 2009—2013年雪季积雪粒径演变Fig.4 Evolution of snow grain size in snow season from 2009 to 2013

3.3 微波辐射亮温模拟结果验证

调整参数值后，基于MEMLS模型模拟2009—2013年的微波辐射亮温，用观测亮温和模拟亮温的均方根误差评估模拟的准确性。如表3所示，在10.65 GHz频率下，入射角度为50°和60°，极化方式为垂直时模拟结果最好，RMSE分别为4.79 K和4.72 K。其次模拟结果较好的是18.7 GHz频率下，入射角为50°和60°，极化方式为垂直极化，RMSE分别为7.13 K和7.25 K。频率为37 GHz时，模拟结果较差，当水平极化入射角为60°时RMSE达到21.76 K。模拟结果最差的是在90 GHz下，极化方式为垂直极化，入射角为50°和60°，RMSE分别为22.37 K和22.14 K。因此在低频波段的垂直极化的模拟结果较好，高频模拟结果较差。

表3 模拟和观测亮度温度的均方根误差Table 3 Root mean square errors of simulated and observed brightness temperatures

比较不同积雪期的模拟准确度（表4），在10.65 GHz且垂直极化下，除入射角为40°，稳定期的模拟亮温的RMSE均小于其他两时期；在18.7 GHz且垂直极化下，在入射角为40°和60°，稳定期的模拟亮温的RMSE小于其他两时期；在37 GHz且垂直极化下，MEMLS模型在稳定期的模拟效果较积累期和消融期更好；在90 GHz的模拟结果较差，模拟亮温和观测亮温的RMSE均大于20 K。将2009—2013年模拟出的亮温和观测亮温进行了对比分析（图6）。总体上，水平极化较垂直极化的观测和模拟亮温波动性更大；10.65 GHz和18.7 GHz相较37 GHz和90 GHz的模拟亮温更为接近观测亮温；稳定期较积累期和消融期的模拟亮温更为接近观测亮温。

图6 2009—2013年入射角50°、垂直（V）和水平极化（H）下10.65、18.7、37 GHz和90 GHz的测量和模拟亮度温度Fig.6 Measured and simulated brightness temperature at 10.65，18.7，37 GHz and 90 GHz at incident angle 50°，vertical（V）and horizontal polarization（H）from 2009 to 2013

表4 不同积雪期模拟和观测亮度温度的均方根误差Table 4 Root mean square errors of simulated and observed brightness temperatures in different snow period

4 结论与讨论

依据2009—2013年北欧积雪观测数据，分析了芬兰Sodankylä研究区域的雪季的不同时期的积雪物理演化过程，并利用MEMLS模型对积雪的微波辐射亮温进行了模拟，分析认为：基于温度和雪深的变化情况，可将整个雪季分为积累期（10月—次年2月）、稳定期（2—4月）、消融期（4—5月）。

地面观测数据分析表明，积雪物理演化过程有如下3点特征：1）气温和雪层温度在稳定期达到最低，雪深在稳定期末积累到最厚；2）积累期和稳定期颗粒形状主要为圆形粒径、片状颗粒和深霜，深霜粒径较大；3）积累期早期会出现融化状态和降水粒子，消融期积雪颗粒形状主要为融化状态。

从3个积雪期与雪深、亮温差的关系来看，积累期（10月—次年2月）和稳定期前期（2—3月），雪深和亮温差（18～37 GHz）整体来看变化趋势一致，在稳定期后期（3—4月），雪深均超过60 cm，亮温差趋于饱和。在消融期（4—5月），积雪迅速融化雪深减小，亮温差（18和37 GHz）基本维持在±15 K以内且波动性较大，积雪微波辐射亮温差和雪深两者的关系随着雪季的不同时期变化，难以用静态半经验反演算法进行描述。这可能是由于积雪物理特性本身（从而微波辐射）随时间具有较大的变异性，后续可考虑针对不同积雪期的积雪演化动态过程来改进积雪反演算法。

MEMLS模型输入的地表均方根高度参数设为0.001 m［20］较为适合芬兰实验环境，校正后的MEMLS模型在较低频率（10.65 GHz、18.7 GHz），较高入射角（50°和60°）且垂直极化下能更好的模拟地基辐射计观测亮温。3个雪期均在低频（10.65 GHz和18.7 GHz）的垂直极化下，模拟结果较好。在稳定期的10.65 GHz、50°入射角且垂直极化下模拟结果最佳，RMSE最低为2.49 K。对于37 GHz且垂直极化下，稳定期的模拟效果较积累期和消融期更好，这表明根据气温和雪深变化来分类积雪期，可更好的表征积雪演化对微波辐射信号的影响。模拟结果在低频段（10.65 GHz和18.7 GHz）较好，可考虑结合较低频段来改进反演算法。

MEMLS模型模拟和观测亮温的误差是由模型本身的局限性以及观测参数误差等因素造成的。模拟结果受各个输入参数的影响，不同参数对模拟结果影响程度不同［50］。驱动MEMLS模型模拟的参数通过积雪剖面调查获取，该调查为了尽可能详细地捕捉积雪变化，由专家每隔3～5天且连续四个冬季开展雪坑（snow pit）测量。根据该高质量的积雪剖面调查结果进行模拟，模拟和观测亮温仍存在误差，这表明一些影响积雪辐射传输过程的重要变量可能在模型中并没有得到有效的体现，因此有必要开展新的积雪模型研究。此外，地面观测参数在模型模拟和积雪反演中能够起到重要的指导作用，但为模型和反演算法确定合理的积雪结构等效参数从而让模型更好的反映真实情况同样重要。虽然辐射计观测点和积雪刨面调查点的空间距离非常近，但积雪结构的空间变异性依然会导致模拟的结果和辐射计观测结果出现差异，这在一定程度上反映出将传统地面调查技术获取的积雪结构信息直接运用于模型模拟的局限性。而且，积雪剖面调查获取积雪结构信息的过程中会存在调查者主观因素的影响（例如如何确定积雪中每一层的边界以及如何选取粒径的测量位置），这也会导致辐射计观测结果和模拟结果出现偏差。一般而言，高频波段比低频波段对积雪粒径更敏感，而水平极化比垂直极化对积雪分层结构更敏感。本次模拟结果低频优于高频，且低频（10.65 GHz和18.7 GHz）的垂直极化优于水平极化，高频（90 GHz）的水平极化优于垂直极化，这说明积雪结构信息在模型模拟中的重要性。一方面，积雪剖面调查取样过程中的细微偏差足以对模拟结果产生显著影响；另一方面，模型模拟中的理想化条件和自然条件下积雪的状况存在区别，物理理论和实际观测的有效连接还迫切的需要更深入的研究。