宋依娜，肖鹏峰,2,3，张学良，卓 越，马 威

1.南京大学地理与海洋科学学院，江苏 南京 210023；2.自然资源部国土卫星遥感应用重点实验室，江苏 南京 210023；3.江苏省地理信息技术重点实验室，江苏 南京 210023

积雪是冰冻圈和水圈重要的组成部分，在全球变化研究中具有重要地位[1]。一方面，积雪作为一种特殊的地表覆盖物，通过反射大部分太阳辐射以及在消融过程中吸收大量能量，对全球气候系统产生了显著的降温作用[2-3]。同时，积雪融水作为全球水资源的重要组成部分，为全球1/6人口提供了水源[4]，并进一步产生了水力发电、农业灌溉等其他用途[5]。另一方面，积雪具有时空不稳定性，过量或非适时的积雪也会引发自然灾害，给人类生命和财产安全带来巨大的威胁和损失[6]。而积雪深度作为积雪的重要因子，表征局部气候环境特征与水资源条件，是水文预测、气候模拟、雪灾监测与评价等模型的重要输入参数[7]。获取高精度雪深空间分布信息可为水资源管理、气候变化研究和防灾减灾工程等提供科学支撑。

合成孔径雷达(synthetic aperture radar,SAR)因其有较高的空间分辨率[8]且对积雪特性有着较高的敏感度等优势[9-10]，在流域尺度的雪深反演中具有较强的适用性[11-12]。极化SAR通过回波获取关于积雪的散射特征[13]，被广泛应用于积雪相关研究，如干湿雪的分类、雪湿度和雪密度等参数的反演[14-16]，但是目前的研究对极化SAR的极化相位信息的分析较少[17]。极化SAR计算同一天线、不同极化方式的极化相干，可获取极化相干系数和极化相位信息，包括同极化相位差(co-polarized phase difference,CPD)。积雪为各向异性介质，导致HH和VV极化信号在雪层中具有不同的传播速度，因此产生CPD。文献[18]根据此原理建立积雪CPD模型，定量分析了CPD与积雪各向异性结构、雪深的关系。但是，根据CPD模型反演干雪深度需要冰粒的各向异性结构(以冰粒轴间比表示)和积雪密度两个参数，其中冰粒轴间比参数比较抽象，难以通过野外观测获得。而且，目前国际上基于CPD模型进行雪深反演的研究主要使用X波段，暂未见使用C波段极化SAR数据的研究报道。文献[19]假设研究区积雪各向异性介电常数不变，使用一个实测站点处的积雪密度0.07 g/cm3和假定轴间比1.5作为雪深反演算法的输入，利用全极化TerraSAR-X数据反演整个影像的雪深，并以该实测站点的雪深进行验证，获得较高的反演精度。但是，由于积雪的各向异性、密度空间分布不均，仅使用一个实测点的积雪密度和假定轴间比反演整个影像的雪深，将给反演结果带来较大误差。文献[20]根据实测点处CPD的正负，结合该点处实测雪深、雪密度等参数，代入CPD正演模型，求得冰粒轴间比，再借助积雪密度反演算法[21]获取全局雪深分布，但此反演方法较为复杂，且积雪密度反演算法可能会带来新的误差。

本文基于CPD模型提出了新的雪深反演方法，将反演所需要的参数减少为遥感获取的CPD数据以及实测雪深数据，根据全部实测点的雪深与CPD的关系，拟合得到最适合研究区积雪特性的反演模型系数，用它们反演研究区的雪深。为了验证反演方法在C波段、不同积雪条件下的应用能力，使用GF-3全极化数据计算CPD，以新疆阿尔泰山南坡克兰河上游为研究区，根据实测雪情差异划分深雪区与浅雪区进行雪深反演，分析反演结果的差异，以期为山区雪深反演提供新思路。

1 研究区与数据

1.1 研究区概况

克兰河流域位于新疆北部，阿尔泰山南麓，在北温带大陆性气候条件影响下，所形成的积雪为典型的“干寒型”积雪[22]，年平均雪深在40 cm以上。克兰河属融雪径流补给为主的河流，年平均径流量6.0×108m3[23]。积雪融水对当地农业灌溉及畜牧业发展具有重要影响，研究该流域的积雪深度具有重要意义。

研究区范围如图1所示，位于克兰河流域上游东支小东沟源头区域。研究区西南部为戈壁，地形起伏较为平缓；东北部是阿尔泰山南麓侵蚀中山山地，海拔较高，地形起伏较大，是典型的山区环境[24]。

1.2 研究数据

(1)SAR数据。CPD为HH和VV极化通道进行极化相干得到复相干系数的相位角，SAR卫星必须具有HH和VV极化通道极化相干成像的能力才能计算CPD，目前常见的可以计算CPD的SAR卫星有Radarsat-2、TerraSAR-X/TanDEM-X、ALOS-PALSAR等。研究表明C波段是CPD对雪深变化仍有较高敏感性的频率下限[25]，本文研究选取我国的GF-3卫星数据，探讨基于CPD模型使用C波段极化SAR数据反演雪深的能力。GF-3卫星于2016年8月10日发射，最高分辨率为1 m，设计有12种成像模式[26]。综合考虑GF-3卫星数据的成像时间和图像覆盖区域、成像质量、同步观测试验可行性等，选取2018年1月17日获取的阿勒泰地区GF-3卫星全极化条带数据作为实验数据，影像范围为47°49′—48°07′N,88°03′—88°24′E。研究所用的SAR影像相关参数见表1。

(2)野外实测数据。CPD反演雪深的可行性分析、模型参数的确定与反演结果的精度验证均需地面实测数据的支持，因此于2018年1月17日对研究区进行星地同步积雪观测，设计了积雪样带观测和积雪属性剖面观测相结合的观测方案。此时为研究区的积雪稳定期，观测期间最高气温在0℃以下，未发生降水。积雪样带观测共计44个样点，主要测量雪深、积雪密度、雪层温度等参数。积雪剖面观测共计8个剖面点，主要测量雪深、雪层密度、雪层温度、雪层介电常数、液态水含量等参数。雪深直接由钢制直尺测量，取观测点周围3 m的范围内3次测量的平均值，测量精度为0.1 cm；积雪密度由钢制直尺、托板和自制的量雪筒测量，测量精度为1 g/cm3；自积雪表层向下每10 cm划分雪层，测得积雪参数代表该雪层积雪参数，其中雪层温度利用经过标定的针式温度计进行测量，测量精度为0.1℃，雪层介电常数、密度和液态水含量通过雪特性分析仪SnowFork测量3次取平均值得到[27]。

表1 所用的GF-3影像的成像参数Tab.1 Parameters of the used GF-3 image

雪深测量结果如图2所示，可见研究区雪深空间分布不均，西南部戈壁区域积雪较浅，东北部山地区域积雪较深。这是因为来自西南方向的暖湿气团在东移过程与阿尔泰山成近正交角，地形对气流的抬升作用显著，有利于山区降水，使得山麓的降水多于河谷平原地区[28]。为了探究模型在不同积雪条件下的适用性，将研究区分为深雪区与浅雪区两个区域，分别进行反演与分析。浅雪区的雪深观测点16个，平均雪深27.5 cm；深雪区的雪深观测点36个，平均雪深72.6 cm。

深雪区观测剖面5个，平均积雪密度0.23 g/cm3，平均体积含水量0.64%；浅雪区观测剖面3个，平均积雪密度0.17 g/cm3，平均体积含水量0.27%。雪层含水量都接近于0，可认为积雪为干雪，满足CPD模型条件。分层积雪密度与液态水含量的观测结果如图3所示。可以看出浅雪区密度与含水量分层垂直差异小，深雪区分层积雪密度存在明显的垂直差异，雪层密度呈中间大、底部和顶部小的特征。深雪区底层积雪的含水量高于表层积雪，垂直分布表现出较大的差异。图4为研究区雪粒图像，可以看出深雪区较浅雪区雪粒垂直分异更明显，积雪重结晶变质作用更强。

图1 研究区位置Fig.1 Location of the study area

图2 观测雪深点的空间分布Fig.2 Position of the measured snow depth points

图3 积雪剖面观测结果Fig.3 Distribution of snow depth,snow density,and snow wetness of the measured snow samples

图4 各雪层的雪粒观测结果Fig.4 Observed results of snow particle of each layer

(3)其他数据。根据阿勒泰国家基准气候站提供的气温和雪深数据，在星地同步观测前，1月13日前后发生了一次较大的降雪，自1月13日至2月底无较大降雪天气过程，平均雪深约为14.5 cm。此外，采用国家基础地理信息中心发布的全球30 m地表覆盖数据GlobeLand 30掩膜森林等复杂散射机理的区域，使其不参与雪深反演；采用30 m分辨率的SRTM DEM数据作为数字高程数据。

2 本文研究方法

2.1 积雪各向异性结构

积雪的物理特性与其微观几何结构有关。在理想状态下，干雪可以看作是冰粒和空气的二相混合物[29]。为了模拟雪的各向异性结构，假定干雪冰粒为规则的椭球体，以椭球体中心为原点，建立三维笛卡儿坐标系(x,y,z)，z轴为冰粒晶体对称轴，平行于重力方向，x和y平面平行于平坦地表。ax、ay和az是椭球体在x、y、z方向上的半轴长，假设ax=ay，冰粒轴间比r=ax/az，表示冰粒的形状。r=1时，为各向同性随机结构的新降雪；r>1时，积雪在自身重力压实作用下迅速密实化[30]，逐渐转变为各向异性、圆盘状水平排列结构的雪粒[31]，r越大，各向异性越强；r<1时，积雪受温度梯度作用发生重结晶和变质作用，形成各向异性、针状垂直排列的变质雪粒[32-33]，r越趋近于0，各向异性越强。由新降雪到发生重结晶与变质作用变成陈雪的过程，伴随冰粒晶体形状及轴间比r的演变如图5所示。

图5 积雪冰粒演化过程Fig.5 Evolution of snow particle shape

2.2 同极化相位差模型

当电磁波穿透干雪层时，积雪粒子的大小远小于电磁波波长，干雪可认为是均质的、体散射可忽略不计的介质[34-35]，因此后向散射信号主要来源于下垫面，因雪层折射率而延迟。文献[36]基于Maxwell-Garnett理论将冰川冰的各向异性结构与观测到的负相位差联系起来，建立了CPD模型。文献[37]进一步基于干雪微结构及其微波极化特性，将CPD模型与雪的微观结构联系起来。

如图6所示，当电磁波以θ入射角穿透一个深度为SD的干雪层，冰粒在电磁波极化电场作用下极化，依据Laplace方程及其边界条件，计算出3个方向的退极化因子Ni，i∈{x,y,z}[38]

(1)

Nx=Ny=0.5(1-Nz)

(2)

图6 极化信号穿透雪层产生相位差示意图(改编自文献[37])Fig.6 Electromagnetic wave penetrates the snow layer and results in the phase difference (adapted from reference[37])

根据Maxwell-Garnett等效介质理论计算冰粒各向异性介电常数εeff,i，i∈{x,y,z}

(3)

式中，Ni为退极化因子；εice为冰的等效介电常数；εair为空气的等效介电常数；fvol为积雪与冰的密度比，故各向异性介电常数可由积雪密度ρsnow、冰粒轴间比r计算得到。

因为折射率n与介电常数εe存在εe=n2的关系，故HH极化和VV极化的折射率nH和nV与3个轴向上的等效介电常数εeff,i、入射角θ满足以下关系

(4)

(5)

HH和VV极化信号在雪层中的折射率不同，导致后向散射信号所传播的路径和距离不同，从而产生了相位差异。利用波长与相位的关系，以及极化信号传播路径和雪层的几何关系，得到不同雪深引起的同极化相位差

(6)

2.3 雪深反演方法

根据式(6)可计算厚度为SD的积雪导致的相位差，由此可以推导雪深反演模型为

(7)

根据nV、nH与εeff,i的关系、εeff,i与积雪密度ρsnow、Ni的关系以及Ni与冰粒轴间比r的关系可知，在确定SAR数据计算CPD之后，输入积雪密度和积雪冰粒轴间比参数即可得到雪深。

本文研究在文献[19]的基础上进一步深入，通过模拟CPD与雪深的关系可知，在假定积雪各向异性介电常数不变的情况下，CPD与雪深呈线性关系，如图7所示。当研究区有多个实测点时，可根据实测点雪深数据与从SAR数据计算的CPD拟合得到最适合研究区积雪特性的线性模型经验系数a、b，从而得到该区域全部像元的雪深。即

(8)

(9)

式中，γc为绝对相干系数；SVV和SHH分别为VV和HH极化单视复数影像；〈〉为滤波器。本文借助干涉雷达数据处理的软件Gamma计算CPD，包括数据预处理、极化相干、空间滤波、地理编码等步骤。考虑实测样点的空间分布，设置地理编码CPD的空间分辨率为10 m。

图7 积雪各向异性介电常恒定下CPD与雪深的关系Fig.7 Relationship between CPD and snow depth based on constant snow anisotropic relative permittivity

3 反演结果分析

3.1 反演精度分析

深雪区与浅雪区积雪属性差异大，区域内部积雪属性差异小。为获得更精确的反演结果，分区域进行雪深反演。对于同一实测点，实测雪深可以看作一个常数值，而不同滤波窗口大小计算得到的CPD不同，则反演得到的拟合关系、反演精度不同。本文使用留P法进行交叉验证，以有效避免过拟合和欠拟合状态，评价不同滤波窗口大小下深雪区、浅雪区的反演模型优劣。为了充分利用所有实测点数据，P取1，即每次假定一个点的雪深未知，使用其余点处实测雪深值与CPD拟合得到模型系数，反演得到该点雪深，然后使用相关系数R和均方根误差RMSE来分析所有点上雪深实测值与反演结果的误差，如图8所示。

结果表明，浅雪区的最优滤波窗口为59×59像元，相关系数R为0.83，均方根误差RMSE为2.72 cm；深雪区的最优滤波窗口为33×33像元，相关系数R为0.54，均方根误差RMSE为11.69 cm。随着滤波窗口大小的增加，浅雪区和深雪区的R均呈现出先上升再下降的趋势，RMSE呈现先大幅上涨波动、然后减小、最后稳定的趋势。在滤波窗口25×25像元以下的结果，R与RMSE出现较大波动，这是因为滤波窗口较小则不能有效抑制相干斑噪声，相位的波动使相干性降低；当滤波窗口过大，则会平滑CPD的变化，削弱图像的细节，也导致相干性降低。

图8 反演结果R和RMSE与滤波窗口大小的关系Fig.8 Relationship between R,RMSE and the window size of filter

不同区域最优滤波窗口大小下雪深的空间分布如图9所示。图9(a)显示深雪区反演得到的雪深自西向东逐渐减少，南部和北部的雪深小于中部区域。图9(c)为浅雪区雪深反演的结果，可见除西北部反演雪深明显大于其他区域以外，其他区域反演雪深呈相间分布。掩膜森林、水体等具有复杂的散射机制及受人类活动影响严重的人造地表区域的雪深结果如图9(b)、(d)所示。

3.2 反演不确定性分析

地形通过影响水热条件从而影响积雪空间分布，地形起伏带来雷达构像几何关系的变化从而影响雷达散射机制，对雪深的反演结果产生影响[39]。为分析地形对反演精度的影响，选取坡度作为变量与反演误差进行相关性分析。如图10所示，雪深反演误差与坡度显著相关，随着坡度的增加，雪深的反演误差呈现增加趋势，且浅雪区的相关性强于深雪区。由浅雪区实测点坡度分布范围可知，浅雪区的地形起伏较深雪区更大，由于构建CPD模型时假定冰粒的z轴平行于重力方向，复杂的地形会导致极化信号穿透雪层的几何关系发生改变，影响模型的适用性，同时电磁波回波由于受到地形起伏干扰作用会引入过多的噪声，从而增大反演误差。

图9 最优滤波窗口下的雪深反演结果Fig.9 Spatial distribution ofretrieved snow depth with optimal window size of filter

CPD模型以雪层均质、忽略体散射为基础，然而严格来说，自然界中的媒介质都是不均匀的，不可避免会引入误差。由实测结果可知，深雪区积雪密度较浅雪区大，由图4实测雪粒图像也可以看出深雪区经过重力压实与变质作用导致粒子粗糙度更大，雪层均质性更差，极化信号穿透雪层时体散射增强，此时，后向散射信号不仅仅取决于下垫面，也由雪层中不同的散射体决定，有效散射中心向雪层上表面移动，因此同极化相干性受到体去相关的影响而降低[18]；根据电磁波在介质中的穿透深度公式[29]，针对本文所用的GF-3数据，对不同湿度下的积雪穿透深度进行了模拟，结果如图11所示。可以看出C波段穿透深度随着积雪湿度增加而显著下降。由图3的剖面雪深与湿度的关系可以看出，由于深雪区较强的土壤热传导作用导致积雪底层含水量较大，部分底层含水量超过3%，水的介电常数在很大程度上改变了雪的介电特性，大幅降低了极化信号对雪层的穿透能力，导致极化信号不能到达下垫面，因此深雪区雪深反演精度低于浅雪区。

图10 反演误差与坡度的关系Fig.10 Correlation between slope and inversion errors

此外，为了探讨不同入射角观测几何条件下CPD模型反演能力的差异，使用入射角发生改变时单位雪深引起的CPD变化来表征，单位雪深引起的CPD变化越大则代表CPD对该入射角下的电磁波越敏感，反演能力越强。模拟CPD与不同入射角之间的关系如图12所示。CPD随SAR电磁波入射角的增加而增加，单位雪深引起的CPD变化也随着SAR电磁波入射的增加而增加，但二者是非线性关系。当入射角增大时，HH极化与VV极化的折射率差增大，所引起的CPD变化也增大，即CPD对入射角较大的SAR电磁波更加敏感。故不同入射角观测几何条件下，本文方法的反演能力存在差异，它更适合具有大入射角的SAR卫星。

图11 C波段极化信号积雪穿透深度模拟Fig.11 Simulation the penetration depth of C band polarimetric signal

图12 模拟CPD与不同入射角之间的关系Fig.12 Simulation the relationship between CPD and incidence angle

3.3 反演方法比较

为验证本文所提出的反演方法的优势，与已有方法进行对比。根据文献[19]提出的雪深反演方法，使用研究区平均积雪密度和假定轴间比1.5代入反演模型，求得雪深；根据文献[20]提出的雪深反演方法，结合不同滤波窗口下的实测点处CPD、实测雪深、雪密度等参数求得轴间比，再结合该点处的实测积雪密度，输入模型反演该点处的雪深。使用留P法(P=1)进行交叉验证来评价3种反演方法优劣，比较结果见表2。根据文献[19]提出的雪深反演方法得到掩膜后的深雪区、浅雪区的雪深分布结果分别如图13(a)、13(c)所示，根据文献[20]提出的雪深反演方法得到掩膜后的深雪区、浅雪区的雪深分布结果分别如图13(b)、13(d)所示。

表2 本文方法与已有方法的雪深反演精度比较Tab.2 Comparison of the proposed method with the existing methods

图13 已有雪深反演方法的反演结果Fig.13 Spatial distribution of retrieved snow depth based on the Patil’s and Majumdar’s methods

可以看出，对于本文研究所用GF-3数据，本文方法比另外两种方法精度更高，且将反演所需要的参数减少为遥感获取的CPD数据和进行模型训练的实测雪深数据。比较结果证明了本文方法具有较好的反演积雪深度的能力。

4 结 论

本文以新疆阿勒泰克兰河上游地区为研究区，在获取C波段全极化GF-3数据及地面同步观测数据的基础上，基于CPD模型构建线性的半经验雪深反演模型，根据实测积雪属性差异，分深雪区、浅雪区反演雪深，探讨了反演方法在不同积雪条件下的适用性，对反演不确定性进行了分析并与已有方法进行比较，主要结论包括：

(1)在积雪各向异性介电常数视为恒定的理想情况下，CPD仅是雪深的函数，可根据全部实测点的雪深数据与从SAR数据获取的CPD数据，拟合得到最适合研究区积雪特性的模型经验系数，反演研究区雪深。反演精度的高低与计算CPD过程中使用的滤波器的窗口大小有关，浅雪区的最优滤波窗口为59×59像元，此时R为0.83,RMSE为2.72 cm，深雪区的最优滤波窗口为33×33像元，此时R为0.54,RMSE为11.69 cm。

(2)雪深反演误差与坡度显著相关，随着坡度的增加，雪深的反演误差呈现增加的趋势，雪深反演不确定性受雪层变质程度、含水量及卫星入射角观测几何条件影响，反演方法对于干燥、雪层变质结晶程度低、均质的积雪，以及具有大入射角的SAR卫星有更好的适用性。

(3)对比已有基于CPD模型的雪深反演方法，本文所提出的反演方法精度更高，并且将反演所需要的参数减少为遥感获取的CPD数据和进行模型训练的实测雪深数据，能够可靠地反演积雪深度。