0 引言

全球导航卫星系统（Global Navigation Satellite System，GNSS）可以为用户提供全天候高精度的导航、定位、授时信息．随着GNSS的快速发展，有学者发现GNSS反射信号中蕴含着地物反射表面信息和地球物理参数，由此衍生出一种新的测量技术，即GNSS反射测量技术（GNSS-R）^[1]．GNSS-R广泛应用于海面风速和风场、海冰厚度、土壤湿度、积雪厚度等参数的探测^[2-7]．

1993年，欧洲航天局（ESA，简称欧空局）提出利用GNSS反射信号进行海面高度测量^[8-10]，GNSS-R测高技术得以发展．相对于传统测高手段和雷达测高技术，GNSS-R测高技术具有监测范围广、受环境影响小、全天候、低成本的优点^[11-12]．文献^[13]利用广泛使用的Lomb-Scargle周期图（Lomb-Scargle Periodogram，LSP）法，对GNSS信噪比（Signal-to-Noise Ratio，SNR）数据的主要多路径频率进行估计，并反演海面高度．由于潮位数据时间分辨率的高低直接影响潮位变化信息的有效性，因此在保证一定精度条件下，潮位数据时间分辨率越高，潮位分析和海平面的监测精度越高，效果越好．但是基于LSP方法的潮位反演在一段SNR序列中通常只能得到一个时刻的潮位值，进行粗差后会出现反演点缺失、难以刻画潮位起伏的缺点．结合GNSS多模多频的优势，文献^[14]利用GPS系统和GLONASS系统L1和L2频段信号进行潮位反演，证明GLONASS系统也可以实现潮位反演．文献^[15]利用BDS系统的SNR数据，结合相位组合方法估计海潮位，通过多源数据结合提升反演结果的质量．文献^[16]提出一种多模多频数据联测的方法，通过提升观测数据量来提高反演值数量．文献^[17-18]利用最小二乘方法进行多系统联合潮位反演，提高了反演值数量．此外，文献^[19-20]利用小波分析方法提取SNR序列瞬时频率，进而得到潮位值，在精度损失较小的情况下大幅提高反演点数量和数据采样率，但因依赖LSP方法，无法单独使用．文献^[21]提出一种利用最优频段潮位反演结果改正其他频段反演结果的方法，即进行LSP频谱分析时，选择距离最优频段反演潮位范围最接近的峰值反演频率，不仅能够增加反演结果数量，而且也可以提高反演结果的精度，但其结果受限于最优频段的精度．还有学者提出窗口LSP（Windows LSP，WINLSP）分析法，利用时间窗口或高度角窗口截取SNR序列，每个窗口得到1个反演值，从而提高反演结果的时间分辨率．文献^[22]选择窗口长度为800个SNR数据截取SNR序列，可以将每天的潮位数量提升至80个左右．文献^[23]用窗口长度为5°的高度角窗口截取SNR序列，能够在十几分钟内得到50多个潮位值．但是，如果时间窗口选取太短，将会导致截取的SNR序列不能满足频谱分析的要求; 此外，用高度角窗口截取后，虽然反演点数量增多，但会导致反演精度降低，且有些测站反演潮位的有效高度角区间较小，致使这种方法无法使用．

目前，GNSS-R潮位反演的研究大多利用LSP方法对SNR数据的主要多路径频率进行估计，其结果往往会得到多个峰值．当频谱分析可靠性较高时，最高峰值往往高于次峰值很多倍^[21]，即峰噪比较大; 当峰噪较小时，SNR序列可能有多个主频率，最高峰值对应的频率可能未携带有效潮位信息^[22]，即存在较多的噪声数据．

Savitzky-Golay（SG）平滑滤波方法是基于最小二乘拟合的卷积方法提出的，它通过高阶多项式对滑动窗口内数据进行最小二乘拟合，用于平滑连续数据，并被广泛应用于信号平滑和数据降噪，包括雷达信号、光谱数据、地震监测数据^[24-26]等数据的降噪．在提取多频潮位反演最优值的过程中，当峰噪比较小时，潮位结果存在较多噪声数据，故最优值的提取也是一个噪声去除的过程．因此，本文从多路径频率的估计入手，提出利用SG平滑滤波方法优化提取SNR序列的最优潮位反演值．首先，获取LSP频谱分析的前4个峰值; 然后，利用SG平滑滤波方法筛选最优峰值，并反演得到最优潮位反演值; 最后，采用BRST站和MAYG站30 d的数据进行验证，并与LSP、WINLSP 两种方法的反演结果进行对比，结果表明本文所提方法能够在保证精度的前提下，增加反演结果的数量，提高数据的利用率和反演结果的时间分辨率．

1 GNSS-R潮位反演原理及SG平滑滤波

1.1 GNSS-R潮位反演原理

GNSS接收机在低高度角时会接收到直射信号和反射信号合成的干涉信号^[27]．GNSS-R潮位反演就是根据直射信号与反射信号之间的延迟，利用信号路径与反射面之间的几何关系，由SNR数据来反演潮位值．其基本原理如图1所示．

图1中，e表示卫星高度角，h表示接收机天线相位中心到水面之间的相对高度．GNSS信号经过水面发生镜面反射，反射信号与直射信号之间的路程差可以表示为δ:

由式（1）可以进一步求出其相位差φ:

式中，λ表示载波波长，φ表示相位差．由式（2）看出，相位差与高度角正弦值成正比．随着卫星运动，高度角不断变化，相位差也随之发生改变．

当GNSS接收机放置于水边时，接收到的干涉信号可表示为

式中，A_c表示合成信号的幅值，A_d表示直射信号的幅值，A_r表示多路径反射信号的幅值．

由式（2）和式（3）可知，当卫星位置变化时，直射信号和反射信号的相位差φ随高度角e变化，表示载波信号质量好坏的SNR数据也会随之发生变化．受多路径和天线增益的影响，在低高度角时，多路径效应加剧，导致SNR数据出现振荡现象^[15]．所以利用低高度角的SNR数据，能够更好地分析反射信号中的潮位信息．直射信号决定了干涉信号的整体变化趋势，即A_d远大于A_r．为了提取可用于反演地表参数的反射信号，需要使用低阶多项式去除SNR数据的趋势项．去除趋势项后的SNR残差可表示为

式中，A表示信号振幅，φ表示相位．若以t表示sin e，以f表示$\frac{2h}{\lambda }$，则式（4）可以化为一个标准的余弦函数:

式（5）中，t为非等间隔采样．使用LSP频谱分析方法SNR残差进行频谱分析^[13]，可以得到式（5）中的频率f．接收机天线相位中心到水面的相对高度h可由$\frac{\lambda f}{2}$求得，利用天线相位中心的实测高度减去h，便可以计算出相应坐标系下的水面高度，实现GNSS-R技术潮位反演．

1.2 SG平滑滤波

SG平滑滤波是一种卷积滑动窗口加权平均算法，在滤除噪声时能够保证信号的形状、宽度不变，使得到的数据波形尽可能逼近原数据波形^[24]，广泛应用于数据平滑降噪．设一个以x（i）为中心包含2M+1个数据点的窗口，构造一个p阶多项式q（n）拟合该数组^[28]，如下:

式中，a₀，a₁，···，a_m为拟合系数．经过最小二乘拟合得到残差C:

当残差C最小时，滤波效果最佳^[28]．首先，求得C最小时的多项式系数，得到拟合曲线; 然后，取数据中心点处的拟合值作为滤波后的值; 最后，通过移动窗口得到原数据的拟合点．本文获取LSP频谱分析的前4个峰值，即在同一时刻有4个反演值，选择反演值和拟合值差值最小的作为输出结果．

2 数据处理

2.1 站点介绍

BRST测站（48.4°N，4.5°W）位于法国西海岸的Brest海港岸边，安装的是Trimble NetR9大地测量型接收机，数据采样间隔为30 s．附近的Brest验潮站可以提供采样间隔为1 min的验潮数据．为收集来自海面的反射信号，实验中GNSS的方位角设定为130°~270°，有效高度角区间设定为5°~30°．

MAYG测站（12.78°S，45.26°E）位于印度洋西北方的马约特岛，安装的是Trimble NetR9大地测量型接收机和Trimble TRM59800.00天线，数据采样间隔为30 s．附近的Dzaoudzi验潮站可以提供采样间隔为1 min的验潮数据．实验中GNSS的方位角设定为20°~170°，有效高度角区间设定为5°~20°．

2.2 数据处理

SG平滑滤波提取潮位的数据处理流程如图2所示．下面以2021年年积日第157~172天，BRST测站G1号卫星观测数据为例进行具体介绍．

1）对G1号卫星观测数据进行分析，得到SNR残差序列．图3为G1号卫星2021年第160天的SNR残差序列．

2）对SNR残差进行LSP频谱分析，结果如图4所示．

3）提取频率值．图4中f₁~f₄为前4个最大频率振幅的频率值．从图4中可以发现f₁和f₂振幅值相近，即峰噪比较小．按照目前主流方法进行分析，将会舍弃由该段SNR数据得到的反演值，导致该时段内潮位值空缺．

4）首先，根据4个频率分别计算潮位高度; 然后，进行粗差剔除; 最后，利用SG平滑滤波进行最优潮位值提取．SG平滑滤波时，均选取窗口大小为5，进行二次曲线拟合，并取窗口中心点处的值作为拟合值，选择反演值和拟合值差值最小的作为输出结果．

5）为了保证滤波精度，根据预先设定的阈值，对该差值进行限制．如果满足限制条件，便得到反演结果; 反之，超过限制条件，则对优化反演结果重新进行滤波．由于目前基于GNSS-R技术进行潮位反演的精度主要在分米级和厘米级，故该阈值取值范围设为0~1 m，优选为0.1~0.5 m．

图5为多频反演结果与SG平滑滤波提取到的最优潮位．

图5中，h₁~h₄为f₁~f₄对应的反演高度值，黑色实线为每一个窗口的二次拟合函数，黑色虚线为SG平滑滤波后提取到的潮位值连线．

图6为接收机附近的Brest验潮站2017—2021年验潮数据．可以看出，BRST测站附近有效高度区间为0~8 m，区间外的值认为是误差．图5中第160天的h₂被作为误差剔除，滤波后输出的潮位值为h₁．

3 实验与结果

3.1 BRST站结果分析

选用BRST测站2021年年积日第152~182天（30 d）的GPS L1波段数据，分别用LSP、WINLSP，以及SG平滑滤波结合LSP、WINLSP方法反演潮位．数据处理中，WINLSP和SG-WINLSP的窗口长度设定为100个SNR数据（约50 min），步长设定为20个SNR数据．将4种方法的反演结果与验潮站实测数据进行对比，结果如图7所示．

图7a中，黑色实线代表验潮站潮位值，黑色圆点代表LSP反演结果，灰色圆点为SG-LSP反演结果; 图7b中，黑色圆点代表WINLSP反演结果，灰色圆点为SG-WINLSP反演结果．可以看出，4种方法获得的潮位结果均与验潮站数据对应良好．

4种方法的每日有效时段和单日数据量对比如图8所示．图8a为有效时段分布曲线，可以看出，SG-WINLSP方法反演结果的有效时段最多，其值最少为18 h，最多为23 h，平均每天可以提供20 h的潮位数据．LSP方法反演结果的有效时段最少，其值最少为6 h，最多为17 h，平均每天可以提供10 h的潮位数据．WINLSP和SG-LSP方法平均每天可以分别提供17 h和13 h潮位数据．

图8b为第155天4种方法反演结果数量对比，图中横轴表示小时，纵轴表示该天内每个小时内的反演值数量．可以看出，1、2、9、22时，无有效的反射信号可用，4种方法均没有反演结果．此外，LSP方法还有11个小时没有可用值; SG-LSP方法仅在8、15、16、21时没有可用值; WINLSP方法仅在7、18时没有可用值; SG-WINLSP均有可用值．因此，利用SG平滑滤波方法能够提高反演值的数量和每日有效时段．

4种方法都是对同一SNR弧段进行分析，当SNR质量不佳时，LSP方法无法提取到有效潮位信息．但WINLSP方法是将SNR序列分割为多个窗口后进行反演，剔除粗差后，符合要求的反演值能够被保留下来．SG-LSP和SG-WINLSP方法通过提取多频数据，在干扰较大的情况下依然能够提取SNR序列中的有效数据．为了进一步验证SG平滑滤波方法提取潮位的精度，将4种方法的反演结果与验潮站实测数据进行比较，结果如表1所示．

从表1可以看出:LSP方法和SG-LSP方法反演结果相关系数达到了0.98，SG-LSP方法日均反演点数量达到19.30个，比LSP方法提高了34.3%，最大时间间隔达到6.63 h，比LSP方法减少了43.2%；同WINLSP方法相比，SG-WINLSP方法反演结果的精度提高了7 cm，达到0.522 2 m，相关系数达0.94，日均反演点数量提高了24.2%，达到60.03个，最大时间间隔减少了25.4%，达到3.08 h．

4种方法反演结果的频率分布如图9所示，图中横轴为频谱分析提取到的频率值位次，纵轴为反演值数量．可以看出:SG-LSP方法的反演结果中f₁占67%、f₂占20%、f₃占7%、f₄占6%; SG-WINLSP方法的反演结果中f₁占68%、f₂占18%、f₃占8%、f₄占6%，且4个频率的反演值均超过了100．上述结果表明SG平滑滤波方法可以有效提取f₁~f₄中的潮位信息，提高了数据的利用率．

为了验证不同窗口和步长对WINLSP和SG-WINLSP两种方法反演结果的影响，设置5组不同的窗口和步长进行对比，结果如表2所示．可以看出:当窗口长度一定时，不同步长的反演结果RMSE值相差不大，随着步长的减小，反演结果数量增多; 当步长一定时，随着窗口长度的减小，反演结果的RMSE值增大，反演值数量增多．这一结果表明窗口长度主要影响反演值的精度，同时会影响反演值数量，而步长主要影响反演结果的数量，对反演结果的精度影响较小．

3.2 MAYG站结果分析

选择MAYG测站2022年年积日第110~140天（30 d）GPS L1波段数据进行分析，4种方法的反演结果如图10所示．这里，WINLSP和SG-WINLSP的窗口长度为80个SNR数据（约40 min），步长为10个SNR数据．

由图10可知，4种方法的潮位反演结果与验潮站数据对应良好．

4 种方法对应的每日有效时段和和单日数据量对比如图11所示．

图11a为4种方法反演结果的有效时段分布曲线．可以看出:SG-WINLSP方法反演结果有效时段最多，其值最少15 h，最多20 h，平均每天可以提供18 h的潮位数据; LSP方法反演结果的有效时段最少，其值最少6 h，最多16 h，平均每天可以提供11 h的潮位数据; WINLSP和SG-LSP方法平均每天分别可以提供16 h和13 h潮位数据．图11b为第133天4种方法反演结果数量对比细节．可以看出该天的5、6、11、13、17、18、22时，由于GPS卫星无有效的反射信号可用，4种方法均没有反演结果．除此之外，LSP方法在该天中还有8个时段没有可用值; SG-LSP方法通过多频提取后，在1、2、3、14、16、21时没有可用值; WINLSP方法通过分割SNR数据，仅在2、10、21时没有可用值; SG-WINLSP均有可用值，且数据量均居首位．

4 种方法的反演结果与验潮站实测数据的对比结果如表3所示

从表3可以看出:同LSP方法相比，SG-LSP方法的反演精度与其相当，达到0.220 9 m，相关系数为0.96，日均反演点数量达到15.23个，提高了19.6%，最大时间间隔达到7.07 h，减少了29.4%；同WINLSP方法相比，SG-WINLSP方法的反演精度较之提高了7 cm，达到0.314 7 m，相关系数为0.92，日均反演点数量达到71.17个，提高了45.9%，最大时间间隔达到4.50 h，减少了28.6%．这一结果表明SG平滑滤波方法可在保证反演精度的情况下，提高反演结果的时间分辨率．

图12为4种方法反演结果的频率分布，横轴为频谱分析提取到的频率值f₁~f₄，纵轴为反演值数量．可以看出:SG-LSP方法的反演结果中f₁占74%、f₂占13%、f₃占8%、f₄占5%; SG-WINLSP方法的反演结果中f₁占65%、f₂占19%、f₃占9%，f₄占7%，每个频率的反演值均超过100．表明SG平滑滤波方法可以有效提取f₁~f₄中的潮位信息，提高了数据的利用率．

4 结束语

本文利用SG平滑滤波提取最优潮位值，进而反演出潮位的变化情况．同时，用LSP、WINLSP 2种方法反演潮位结果作为对比，获得以下结论:

1）SG平滑滤波方法能够在保证精度的前提下，提高反演结果的时间分辨率．

2）在干扰因素较大的情况下，频谱分析中主频不明显，存在误差信息，SG平滑滤波方法能够提取到f₁~f₄携带的潮位信息，提高了数据的利用率．

3）与验潮站实测值相比，对于BRST站:LSP反演结果RMSE值为0.317 9 m，日均反演值14.37个; SG-LSP反演结果RMSE值为0.321 1 m，日均反演值19.30个; WINLSP反演结果RMSE值为0.592 5 m，日均反演值48.33个; SG-WINLSP反演结果RMSE值为0.522 2 m，日均反演值60.03个; 4种方法与实测值相关系数均优于0.92．对于MAYG站:LSP反演结果RMSE值为0.218 5 m，日均反演值12.73个; SG-LSP反演结果RMSE值为0.220 9 m，日均反演值15.23个; WINLSP反演结果RMSE值为0.386 9 m，日均反演值48.77个; SG-WINLSP反演结果RMSE值为0.314 7 m，日均反演值71.17个; 4种方法与实测值相关系数均优于0.88．可以发现:相对于LSP方法，WINLSP方法虽然提高了反演值数量，但是精度也大大减小，而SG-LSP方法能够在保证精度的前提下，提高反演值数量; 相对于WINLSP，SG-WINLSP方法不仅能够提升反演值数量，还能在精度上有一定提高；BRST站和MAYG站均提高了7 cm．整体来看，MAYG站的精度优于BRST站，但是MAYG站反演值与实测值相关系数却低于BRST站，可以认为是BRST站相对于MAYG站潮位变化范围更大、振荡更加剧烈导致的．

后续将进一步寻找适合潮位反演融合算法，发挥GNSS多模多频的优势，获取数量充足的SNR序列，在保证反演精度条件下，尽可能提高潮位反演值的时间分辨率．

参考文献