一种GPS多路径信号的雪深估测方法

晋亚斐 ，周 威，李厚朴

(1.交通运输部天津水运工程科学研究所，天津 300456；2.海军工程大学 导航工程系，武汉 430033)

地球冰冻圈内现存的大量积雪是影响全球气候变化的关键因素之一，其不仅可用于评估水资源储量，同时对地表能量平衡进行反馈调节[1]。由于全球气温上升，表层积雪逐渐消融，融化的雪水一方面为河流补充水流量，另一方面造成海平面高度的缓慢上升。因此，如何快速、精准监测地表的雪深变化量并分析变化规律是一项重要的科学研究。近年来，GNSS-R (Global Navigation Satellite System Reflectometry) 技术探测地表环境特征参数已成为研究热点之一。经过数年发展，GNSS-R应用领域由海洋遥感(海面高度、海风和有效波高等[2-4])扩展到陆地遥感(土壤湿度[5]和积雪厚度[6]等)。由于GPS卫星系统的建设成熟，卫星数量多且轨迹覆盖全球，其发射的L波段电磁波不仅能减弱信号穿透大气层时的衰减，同时敏感反映反射层含水量的变化，因而GPS可作为一种理想信号源[7]。

针对GNSS-R的发展，国内外已取得了较多的研究成果。现有的GNSS-R技术分为双天线和单天线两种模式，前者需采用特制的GNSS-R接收机，而后者使用常规测地型接收机。由于连续运行参考网络 (Continuously Operating Reference Stations, CORS) 逐步完善，基于CORS站探测地表雪深的研究具有更重要的实际意义。Larson等[8-9]首次展示如何利用测地型GPS接收机接收的多路径反射信号，证明反射信号的幅度和相位会随天线周围土壤湿度的变化而变化。接着，Larson等[10]建立GPS反射信号的频率和振幅的数学模型，并正式确立单频反射分量的信噪比观测量与多路径误差之间的映射关系，此为GPS-R模型提供了理论基础。根据GPS-R土壤湿度反演模型，Larson团队[11]利用该模型遥测科罗拉多地区的暴风雪厚度，并与实际测量值对比，验证了测量型GPS接收机同样适用于雪深探测应用。考虑到天线极化以及天线与地面之间的响应特性，Nievinski等人[12]构建了一种新的基于物理的多路径正演雪深模型，估测值与实测值的相关系数达到0.97。黄良珂等人[13]针对数据源单一问题，探索GPS卫星新型的L5载波在雪深监测方面的可靠性，获取了良好的结果。但是，GNSS-R探测雪深仍处于理论阶段，模型变量参数的设置、卫星波段的选择以及测量范围等仍需更多的研究工作。

根据前期研究可知，采用GPS卫星L2C载波[14]和 L5载波[15]反演雪深的精度较高，但是部分CORS接收机仅能接收旧的L1C/A和L2波段信号，因此基于GPS L1C/A信号的GPS-R技术可探测更长时间的雪深序列。本文详细阐述了信噪比观测量与反射高度的关系，以瑞典KIRU测站附近区域为研究对象，验证GPS L1C/A波段的信噪比数据在反演地表雪深方面的适用性。

1 基本理论

表层反射的多路径效应是由相同频率的直射信号与反射信号的相位不同导致的，其相位会随着卫星位置的变化而变化，从而两者发生相干震荡，这一现象将以振幅的形式反映于接收的信噪比数据[16]。在低卫星高度角时，多路径效应的影响将作为信号的主导因素，故可用反射信号的信噪比来构建多路径反演模型。图1为GPS-R技术监测雪深的示意图。设H为接收机天线与无雪时地表的垂直距离，h为接收机天线与积雪表层的垂直距离，hS为积雪厚度，θ为卫星信号的高度角。假设积雪内部均匀且电磁波穿透较少，则大部分信号在积雪表层反射，此为GPS-R遥感理论成立的前提条件。

图1 GPS-R技术监测地表雪深示意图Fig.1 Schematic diagram of GPS-R technology for monitoring surface snow depth

信噪比是一种评估卫星的信号质量和噪声特性的指标，其主要受到天线增益、接收机噪声以及多路径效应的影响。当直射信号分量Sd和反射信号分量Sr同时被接收机的天线接收时，相位的延迟导致两者在接收机内部会发生相干而产生合成信号Sn，其与SNR的关系可用式(1)表达[17]

(1)

低高度角时，相干效应引起的信号波动较明显。信号趋势表现出类抛物线状，中间平稳两端震荡。为了获取包含了环境特征信息的反射信号，一般采用两阶或三阶多项式拟合信号的趋势项，并将两者进行作差。除去趋势部分(主要包含直射信号分量)后，此处仅仅保留含震荡信息的反射信号分量。图2表示03号和14号卫星信噪比的干涉效应。图2中，2-a和2-c的横轴表示卫星高度角，2-b和2-d的横轴表示高度角的正弦值，纵坐标均表示线性化后的SNR，其单位为Volts。此处仅绘出低高度角(5°～25°)范围内的SNR残差序列，其余的反射分量可忽略不计。去趋势后的反射信号SNR的振幅[18]可表示为

(2)

(3)

式中：A为振幅大小；λ为波长；ψ为与反射面有关的相位偏移量。由图3可知，SNR残差序列会随着高度角θ的变化显示出类三角函数的波形。尽管卫星高度角θ是随时间变化的已知量，但是由于时间间隔是不均匀的，长时间的观测序列中无法保证其被整周期地截断，一般不使用快速傅里叶变换(FFT)方法[19]。因此，Lomb-Scargle (L-S)谱分析法[20]用于分析该反射波形的干涉测量度量(相位、振幅、频率等)。最后通过式(4)可求得到天线与反射层之间的垂直高度h。

图3 GPS L1C/A载波在不同高度角的菲涅尔反射区Fig.3 First fresnel zone of GPS L1 C/A at different elevation angles

(4)

式中：f为信号的振荡频率。考虑到每颗卫星在一个年积日内拥有重复的运行轨迹，最多可获取4段dSNR序列，因此每日雪深值需求解众多反演结果的平均值来确定，从而实现反演地表雪深信息的目的。

利用GPS-R技术监测雪深，仅需将大地测量型接收机安置在平坦的陆地上，保证有足够的卫星信号经积雪表层反射后进入天线内部。目前国际上大量气象监测站附近布设了GNSS跟踪站，观测数据可用于近实时、连续的雪深监测。综上，基于GPS-R的雪深遥感技术的基本流程如下：(1)根据卫星观测站周围的地形、植被分布以及卫星接收数据的范围来筛选卫星数据；(2)根据最大化利用GPS观测数据的原则，筛选在较低高度角范围(5°～25°)之内的信噪比数据；(3)通过低阶多项式拟合分离直射分量和反射分量，获取反射信号的高精度dSNR数据序列；(4)利用L-S谱分析法对SNR序列进行频谱分析，获取有效垂直距离，并进一步获取雪深值。

2 实例论证

本试验采用IGS服务中心 (International GNSS Service, http:www.igs.org) 提供的瑞典KIRU跟踪站(纬度：67.86°N，经度：20.97°W，高程：335 m)连续6个月共184个年积日的GPS观测数据，测站高度为1 m。的ESRANGE气象监测站提供了时间分辨率为1 d的实际雪深测量数据 (https:gis.ncdc.noaa.gov)。KIRU跟踪站与ESRANGE气象站位于同一区域(相距约5 km)，可认为是并址观测，对比两者的试验结果可进行算法验证。

美国NCDC中心(National Climatic Data Center)为了确定测量区域的范围且各卫星的运行轨迹不同，需要根据反射区域的阈值大小确定卫星高度角和方位角的覆盖区域，从而筛选出可用来反演雪深的信号。根据菲涅尔反射区理论[21]可知，电磁波的总能量是不变的，在不考虑介质层吸收的情况下，第一菲涅尔区的反射信号最强，能量损失率最少，接收天线可获得介质层的最大信号能量。通常认为，地表反射信号主要来自第一菲涅尔区。由于卫星高度角逐渐变大，GPS信号的菲涅尔反射区的椭圆面积逐渐变小，不考虑地表粗糙度的影响。图3为不同高度角下GPS L1C/A载波的菲涅尔反射区。结合文献[22]的结果可得，反射区域为若干个以GNSS观测站为中心的椭圆。因此，本实验提取了高度角为0°～20°，方位角为30°～330°的信噪比，最大化地保证了实验所需的数据。

随机选取2015年地表雪深具有明显差异的2个年积日的GPS卫星数据进行雪深分析，对应的实测雪深为0 m和0.46 m，结果如图4。图4-a和图4-b分别表示处理后的SNR序列及对应的L-S谱分析。图4-a的横轴表示卫星高度角的正弦值，纵轴表示去除趋势项后的SNR序列。图4-b的横轴表示天线到积雪层的垂直高度，纵轴表示SNR的频谱振幅强度。经由L-S谱分析后，频谱振幅强度的峰值对应的距离(横轴数值)即为垂直反射高度h，由于已知YEL2跟踪站的高度，由hs=H-h可求得某一时段内的估测雪深值hs。

由图4-a可知，随着地表积雪厚度增加，SNR序列的间隔距离增加，即信号的振荡频率(红色横线划出的部分)逐渐减弱，但整体的波动范围(±4 V内)无明显变化。同时，对于一个SNR信号序列来说，SNR的数值会随卫星高度角增大而逐渐减小。该现象出现的原因是积雪层增加使表层的粗糙度增大，信号部分被吸收并发生散射，多路径反射的信噪比随之降低。图4-b表明：(1)两种不同雪深情况下，波峰的幅度值的特征是保持一致的，随着积雪厚度的变化而变化较少；(2)当雪深增加时，使得频谱振幅的最高峰值左移，即垂直反射高度变小。此外，图4-b 中最高波峰右侧出现了一个类似的较高波形，此时可能是因为地面不平整，信号的反射面高度低于测站地面的高度，从而导致垂直反射高度超过了天线高度。上述结果表明雪深变化量与低高度角反射分量信噪比的能量频谱变化关联密切，同时使得信噪比产生了一定的规律性变化。

由于地表积雪厚度是随时间逐渐变化的长序列，利用GPS多路径反射信号估测雪深的可靠性需进一步得到验证，图5-a给出了KIRU观测站2015年连续6个月共184个年积日的雪深估测结果及其测量精度分析。图中，GPS L1C/A反演雪深的整体点位趋势与实测雪深数据在长时间序列上保持了较高的一致性。其中，有两次明显的强降雪过程，分别为2015年11月中旬和2015年12月中下旬，积雪厚度快速增加。对于雪深的急剧变化，GPS反演结果均对此现象表现出强正相关。结合图5-b中的估测偏差值，对不同时段的雪深进行分析，可得几个初步的结果：(1)当地表雪深不足10 cm时，由于GPS-R模型估测的误差影响较大，估测结果不够精准，难以反映出地表积雪的实际情况。其原因是测站周围地表环境复杂(包含各类其他的地物)，地表无积雪(或少积雪)使得地表介质的物理属性差异较大，信噪比数据受到其干扰导致产生的偏差增大；(2)雪深处于15 cm以上，反演值与实测值的一致性较好，该结果较好地体现了雪深的变化趋势。积雪厚度达到40 cm左右时，估测偏差有所降低，GPS信号的探测敏感度更佳。同时，对反演值与实测值的相关性进行比较，结果如图5-c所示。其中，直线斜率为1，绝大部分的雪深点都分布于直线附近。此时，反演值和实测值的均方根误差为4.5 cm，相关系数r为0.96。直线斜率与相关系数十分接近，说明基于GPS卫星较旧的L1C/A信噪比同样可用于地表雪深探测的应用。

图5 GPS卫星估测雪深与实测雪深的精度对比Fig.5 Comparison of GPS L1 C/A-estimated snow depth and in situ snow depth

3 结束语

高精度的地表积雪空间变化等数据的长期监测对全球水资源和气候研究具有十分重要的意义。本文介绍了一种以普通的大地测量型GPS接收机为接收平台，基于信噪比数据的GPS-R雪深遥感技术，该技术能有效节省资金且扩展常规CORS观测网络的应用领域。本文首先阐述了基于GPS L1C/A信噪比反演雪深的观测模式、理论方法以及实例验证，并给出了处理信号的具体流程。结果表明：(1) 单天线模式下，陆基GPS-R技术能以较高的时空分辨率探测地表的雪深变化数据；(2) 随着积雪的不断积累，SNR残差序列的振荡频率变化较缓，且SNR数值会随着卫星高度角增大而减小；(3) 对于地表有积雪和无积雪两种情况下，前者的测量精度较优于后者，从而验证了多路径反射信号与积雪厚度变化之间的强相关。同时，本实验仍存在较多问题亟待解决与完善，如：实验的数据量有限、未考虑反射面粗糙度的影响、多路径信号的有效提取方法等，因此更多的研究工作才能使得GPS-R算法在测量精度和稳定性方面得到长足进步。

致谢：感谢国际IGS中心提供的GPS观测数据和美国国家气象数据中心提供的实地雪深监测数据。