陈殊，何秀凤，王笑蕾，宋敏峰

河海大学 地球科学与工程学院,南京 211100

1 引言

沿海地区是世界上经济活动和城市化最集中的地区，目前世界上有一半人口生活在离海60 km以内地区，而且沿海地区的人口增长速度正在加快。然而由于风暴潮和海平面上升等自然因素，沿海地区将面临越来越多威胁，所以对海面变化进行实时监测并研究其变化规律具有重要意义（Melet 等，2018）。随着GNSS 技术的不断发展与完善，GNSS 信号不仅可以用于定位、导航和授时，利用GNSS反射信号还可以监测潮位变化（Jin等，2014）。

自Martin-Neira（1993）提出用GNSS反射信号反演海潮变化以来，GNSS-MR 技术逐渐成熟，成为沿海海平面测高的另一种可行方案（Geremia-Nievinski等，2020）。但是，GNSS-MR潮位反演技术在当前的实际应用依然存在精度不够高和分辨率低的两大问题（Larson 等，2017）。反演结果的单点精度因信号类型而异，质量好的信号反演的精度更高，采样率由卫星轨迹数量、signal-tonoise ratio（SNR）数据质量以及环境的几何形状决定。多模多频GNSS 可以观测到更多的卫星、获取质量更好的信号和更为开阔的信号感知区域，这为GNSS-MR 反演潮位提供了更多的反演值从而提高精度和时间分辨率（Larson 等，2017；Tabibi等，2017）。Löfgren 和Haas（2014）分别使用了GPS 和GLONASS 的L1 和L2 频段的信号进行潮位反演，证实了GLONASS 系统也可以用于潮位反演。Roussel等（2015）基于最小二乘方法，对GPS 和GLONASS 系统反演的潮位值进行融合。Strandberg等（2016）分别用L1、L2C 和L2P 频段的SNR 数据进行潮位反演，并对信号的质量进行分析。Jin等（2017）第一次用BeiDou系统L2、L6、L7 信号估计海面变化。Wang等（2019）用稳健估计的方法对4 个系统（GPS、GLONASS、Galileo 和BeiDou）潮位反演结果进行融合，验证了多模多频GNSS-MR 可以有效地提高潮位反演的精度和时间分辨率。何秀凤等（2020）利用多模多频GNSS信号反演了沿海台风风暴潮，证实了在极端天气情况下，GNSS-MR 技术依然可以探测到潮位的异常变化。

目前，多模多频GNSS-MR 融合算法主要使用“bisquare”定权方法（Wang等，2019），并将定权中涉及的常系数确定为1；这种定权方法对于偏差大于1 的观测值无法很好的处理。然而，由于GNSS-MR 潮位反演中存在各类误差，反演观测值经常大于1 m。考虑到IGGIII 定权模型（Yang 等，2002）可以对不同大小的偏差进行更细致更合理的处理，文章提出基于IGGIII 模型的稳健估计方法来进行多模多频GNSS-MR 融合，并验证多模多频GNSS-MR 融合方法在不同测站环境用于潮位监测的可行性。

2 方法与原理

2.1 GNSS-MR经典海面反演原理

多路径效应主要与反射面的结构和电介质参数有关，当GNSS 测站位于海边时，接收机在海域方向会接收到卫星直射信号与海水表面反射信号相干的合成信号，这种相干性表现在卫星观测文件SNR的变化上（金双根 等，2017）。SNR可以表示为（Larson等，2008）

式中，Ad和Am分别表示为直射信号和反射信号的振幅；ψ是直射信号和反射信号之间的相位偏差。直射信号与反射信号的相位差为

式中，D为反射信号与直射信号的路径差，λ为信号波长，e为卫星高度角，h为天线相位中心到反射面的垂直距离。

假设反射面时静止的，根据式（2）有

式中，f是SNR 中受干涉影响部分的信号频率，是静态假设下得到的反射高度。式（3）中给出了与f的关系，即对于f的提取，经典方法是先用低阶多项式来去除SNR 序列中的直射信号，再用LSP（Lomb-Scargle Periodogram）分析方法（Lomb，1976；Scargle，1982）从去除直射信号的残差序列δSNR中求出f。

2.2 GNSS-MR潮位监测误差改正

为了求解动态变化的高度h，首先是得到的时间序列，再由此来估算，从时间序列中扣除即可得到改正后的时间序列h。

卫星信号在经过对流层区域时会发生折射效应，这种对流层延迟会使求得的垂直反射距离出现偏差（Williams 和Nievinski，2017）。在GNSS 反射信号反演海面高的研究中，Anderson（2000）首先提出了大气折射引起的无线信号几何弯曲可能会影响反演结果精度。Santamaría-Gómez 和Watson（2017）在GPS-MR 海平面估计中使用了基于大气压力和温度的经验天文折射模型来消除对流层延迟。Williams 和Nievinski（2017）使 用GPT2W（Global Temperature and Pressure）和VMF1（Vienna mapping functions）映射函数，来消除对流层延迟的影响。对流层延迟τT计算公式为

式中，ΔτZ=τZ(-H) -τZ(0)是天线和地面位置的天顶延迟差，mh和mw分别为干延迟映射函数和湿延迟映射函数。于是由对流层延迟τT引起的反演高度偏差ΔhT为

这样经过海面动态改正和对流层延迟改正后即可得到GNSS 接收机到海面实际的垂直距离h，再用接收机的大地高减去垂直距离便是真实的潮位值。

2.3 多模多频GNSS-MR潮位反演值融合方法

GNSS 系统中各个频段的SNR 数据都可以按照GNSS-MR 经典海面反演理论及误差改正方法对潮位进行监测。在整个技术流程中，首先需要确定高度角和方位角区间，获得来自海面的反射信号，再对截取的残差序列δSNR用LSP 分析方法求得有效高度在LSP 分析过程中，需要进行质量控制，一般选取峰值与背景噪声的比值（peak-to-noise ratio）大于3 的反演值（Larson 等，2013a）。最后从海面动态变化和对流层延迟两个方面对垂直高度进行修正。对于多模多频GNSS-MR 潮位反演值，本文采用基于状态方程和稳健估计的方法，考虑海面动态变化以及对流层延迟对数据进行融合（Wang等，2019），实现了每10分钟获得一个潮位值。在建立方程之前，首先需要设定一个滑动窗口，窗口的宽度需要容纳足够多的潮位反演点以满足建立方程求解。设窗口长度为T，窗口的滑动步长为第i个窗口的时间则在第i个窗口内可建立状态方程

式中，j表示信号类型，l表示窗口内反演值的序号，其中|tl-ti|＜T/2。则第i个窗口内所有的反演值可以用方程组表示

可以简化为矩阵方程

但是不同系统和频段信号的SNR 不同，这导致了在LSP 分析过程中每个SNR 弧段的主导频率存在差异，进而反演精度也不相同（Tabibi 等，2017），所以不同的信号应对应不同的权重，本文采用基于IGGIII 模型的稳健估计方法进行自动赋权。

稳健估计的基本原理为：首先通过最小二乘法得到反演值的估计残差；根据残差确定的各反演值新的权因子，进行下一轮最小二乘解算；反复迭代，直到前后两次估计值的变化值小于限差（方兴 等，2018）。

根据式（11）进行第一次迭代

第一次迭代的残差为

在接下来的迭代过程中，精度较低的反演值将被重新定权，逐渐降低其在求解垂直距离hi过程中的作用，以便达到提高精度的目的。潮位反演值主体服从正态分布，可疑和显著异常反演值只是占少部分。在稳健估计中，应充分利用主体可靠信息，所以本文采用IGGIII 模型来自动赋权，图1 为BRST 站反演数据迭代过程中残差vi所对应的权重，从图1中可以看出权重大致分为三段，即正常段，保持原有的权不变；可疑段，降权，使其小于原始权；淘汰段，权值为0，淘汰掉异常值，这不但能够削弱粗差对结果的影响还可以充分的保留精确数据，避免了数据严重失真（Yang等，2002）。IGGIII模型表示为

图1 潮位反演值权重分布Fig.1 Weight distribution of sea level retrievals

在（k+1）次迭代中

3 GNSS站点和数据分析

本文选取BRST 和HKQT 两个测站来验证多模多频GNSS-MR 技术用于潮位变化监测的有效性，这两个测站均可接收来自四系统的卫星观测数据。

3.1 站点介绍

BRST站（48.4°N，4.5°W）位于法国西海岸的Brest 海港岸边，周围海平面日涨落幅在7 m 左右（Santamaría-Gómez等，2015）。图2（a）即为BRST站点的测站环境，BRST站配备了Trimble NETR9大地测量型接收机，提供了GPS、GLONASS、Galileo、BeiDou和SBAS的卫星观测数据，采样间隔为30 s。为了验证多模多频GNSS-MR 技术监测潮位的精度，在实验中利用距离BRST 站292 m 的Brest 验潮站实测数据进行对比分析。

图2 BRST站的观测环境和第一菲涅尔反射区情况Fig.2 Platform at BRST station and the first Fresnel zones for site BRST

由于利用GNSS-MR 技术监测潮位变化仅需获取来自海面的卫星反射信号，本文以5°，10°，20°，30°的高度角来绘制第一菲涅尔反射区（Roesler 和Larson，2018），BRST站海面信号反射区域如图2（b）所示。为了获取来自海面的反射信号，本文设置两个方位角和高度角区间来剔除陆地的干扰信号。在方位角为130°—165°的区域，高度角区间为5°—20°的反射信号被使用；在方位角为165°—330°的区域，高度角区间为12°—25°的反射信号被使用（Wang等，2019）。

HKQT站（22.3°N，114.1°E）位于香港北边海岸，安装了Trimble NETR5 大地测量型接收机，提供GPS、GLONASS、Galileo、BeiDou、QZSS和SBAS的卫星观测数据，采样间隔有1 s、5 s 和30 s。HKQT 测站环境及海面反射信号如图3 所示，该站点来自海面反射信号区域的方位角为-60°—105°，有效高度角为4°—9°。距离HKQT 站2 m 的Quarry Bay实测潮位数据可以用于实验进行对比分析。

图3 HKQT站的观测环境和第一菲涅尔反射区情况Fig.3 Platform at HKQT station and the first Fresnel zones for site HKQT

3.2 数据分析

以BRST 站为例，本文选取BRST 站2020 年年积日为197—253 的SNR 数据进行实验，其中信噪比数据GPS 有S1C、S2W、S2X、S5X 等4 种，GLONASS 有S1C、S1P、S2C、S2P 等4 种，Galileo有S1X、S5X、S7X、S8X、S6X 等5 种，BeiDou 有S1X、S5X、S2I、S7I、S6I等5种。

实验期间内来自各个系统每天可用的SNR 弧段数统计结果如图4所示，从图中可以看出GPS和GLONASS 系统实验期间可用的SNR 弧段数波动较小，信号质量较为稳定，每天使用的弧段数在110左右。Galileo 和BeiDou 系统可用的SNR 弧段数波动较大，信号质量不稳定，但是由于记录的信噪比类型都有5 种，平均每天可用的弧段数在120 左右。为了进一步分析各系统不同波段的SNR 信号质量及其反演能力，以2020年年积日为227的4个系统数据进行分析。图5 为5 个系统各信号SNR 序列和对应的LSP 分析结果，其中BeiDou 系统由于数据记录不完整，只接收到S1X、S5X、S2I、S6I等4种信号。

图4 GPS，GLONASS，Galileo和BeiDou系统每天的SNR弧段数Fig.4 Number of SNR arcs used per day for GPS，GLONASS，Galileo，and BeiDou signals

图5 四系统SNR序列和LSP分析结果Fig.5 SNR sequence and LSP analysis of quad-constellation

由图5（a）四系统SNR 序列可知，在低高度角情况下，GPS 中S2W 的能量明显低于其他三种信号，SNR 的最大差值在10 db 左右，GLONAS、Galileo 和BeiDou 的不同载波信号SNR 差值较小，整体趋势较为一致。图5右边为四系统SNR序列对应的LSP分析结果，其中波峰对应的横坐标即为垂直反射距离，可以看出四个系统不同信号反演结果差值在0.08 m左右。GPS中S2W的SNR序列能量最低，在LSP分析中整体振幅也是最低；S2W和S2X的LSP 分析结果有明显的次波峰，信号质量较差，而S5X在LSP分析中只有一个较大波峰，信号质量最好。GLONASS 和Beidou 中信号的SNR 序列有较好的一致性，其LSP分析结果振幅值比较稳定，都只有一个明显的波峰。Galileo中除了S1X在LSP中振幅值最低且最大波峰不明显，S5X、S7X、S8X、S6X 四种信号LSP 结果都没有虚假高峰，能够较好的得到反演结果。四系统不同信号的LSP表现出不同的结果，每个信号对应于不同的峰值，这表明同一系统不同类型的SNR 会导致精度水平不一样的结果。

4 单模单频GNSS-MR潮位反演序列

以BRST 测站为例，对去除趋势项的SNR 序列进行LSP分析，可以获得每种信号的垂直反射距离的时间序列。在设定LSP分析中峰值噪声比大于3后，4个系统潮位有效反演值如图6所示。

图6 BRST站四系统GNSS-MR潮位反演结果Fig.6 Sea levels at BRST estimated from quad-constellation

由图6可知，4个系统每个信号通过GNSS-MR获取的潮位反演结果都能够描述潮位变化的趋势，但是这些结果与验潮站潮位测量值存在一定偏差。为了进一步评估各个信号的反演效果，同时对海面动态变化误差和对流层延迟误差ΔhT进行改正，BRST 测站每个信号反演序列、-和--ΔhT的均方根误差统计结果如表1 所示，经过两次误差改正后的精度提高百分比如图7所示。

结合表1 和图7 可知通过海面动态改正和对流层延迟改正，4 个系统的反演精度大约提高了50%，GPS 中S5X 反演结果优于S1C、S2W 和S2X，精度最高；GLONASS 中S2C 和S2P 反演精度要优于S1C 和S1P；Galileo 的S5X 反演精度最高，S8X、S7X 和S6X 较次之，S1X 精度最低；BeiDou 系统整体反演精度要低于其他3 个系统，其中S7I 反演精度最低，S6I反演精度最高。

对于BRST 测站而言，经过误差改正后精度最好的SNR 数据分别是GPS 的S5X、GLONASS 的S2P、Galileo 的S5X 和BeiDou 的S61。分析可知随着信号的频率下降，反演的精度会随之提高，这可能是因为频率会影响天线增益模式和随机表面粗糙度，从而影响信号的质量（Nievinski 和Larson，2014）。

5 GNSS-MR多模多频联合潮位反演

5.1 BRST站联合反演结果

在BRST测站中，结合表1和图6可以看出单系统单频的反演结果精度低且时间分辨率差，而单系统多频又会有反演结果在某一时间段堆叠的情况。为了解决上述问题，本文在窗口长度为2 h，滑动步长为10 min 的时间窗口中建立状态方程，在稳健回归模型下实现四系统反演结果融合。滑动窗口长度的选择非常关键：长度过短，窗口中数据量过少，不足以进行可靠的抗差估计；窗口过长，则潮位的非线性运动造成的融合反演误差越大。文章从精度和数据损失率（损失数据量/窗口数量）两个方面分别比较0.5 h、1 h、1.5 h、2 h、2.5 h、3 h 窗口潮位反演差异（表2），最终确定最优窗口长度为2 h。

表2 BRST站不同窗口长度反演精度和数据损失率对比Table 2 Comparison of inversion accuracy and data loss rate of different window lengths in BRST station

图8（a）为四系统联合反演的结果，图中黑色实线表示验潮站实测数据，蓝色的点为联合反演结果，可以看出联合反演的结果较好的对应了验潮站的实测数据，同时实现了10 min 一个点的均匀采样，满足日常潮位监测的要求。改正后精度为12.43 cm，相较于各系统单信号的反演结果，精度大约提高了40%—60%。图8（b）为反演潮位与实测潮位的残差分布，散点大都在-0.5 m—0.5 m 范围内，在涨潮和落潮时残差波动较大。图8（c）是反演潮位与实测潮位的线性回归模型，相关系数为99.56%。

图8 BRST站多模多频GNSS联合潮位反演以及误差分析Fig.8 Multi-GNSS sea level retrieval series and error analysis of BRST

5.2 HKQT站联合反演结果

为了进一步验证基于稳健回归模型融合算法的有效性，本文选取HKQT站2020年年积日为203—232 期间GPS、GLONASS、Galileo、BeiDou 系统的SNR 数据进行实验。其中信噪比数据GPS 有S1C、S2W、S2X、S5X 等4 种，GLONASS 有S1C、S1P、S2C、S2P等4种，Galileo有S1X、S5X、S7X、S8X等4种，BeiDou有S2I、S7I两种。经过海面动态改正和对流层延迟改正后，四系统整体的潮位反演精度约为25.42 cm。

在HKQT 站中，同样从精度和数据损失率两个方面分别比较不同窗口潮位反演差异（表3），在保证精度和数据利用率的情况下最后选择2 h的时间窗口，滑动步长为10 min，然后在每个窗口中进行四系统反演结果融合。图9（a）为四系统联合反演的结果，图中黑色实线表示验潮站实测数据，红色的点为联合反演结果，联合反演的结果可以清晰的反映潮位的起伏变化。改正后精度为7.09 cm，相较于四系统整体的反演结果，精度提高了约72%。图9（b）为反演潮位与实测潮位的残差分布，散点大都集中在0 m 左右，整体表现更为稳定。图9（c）是反演潮位与实测潮位的线性回归模型，相关系数为99.06%。与BRST站的联合反演结果相比，基于稳健回归模型的数据融合算法在HKQT 站即潮位起伏较小的区域表现更好，精度提高更多。虽然BRST 站获取的信号种类多于HKQT 站，但是由于Galileo 和BeiDou系统存在个别信号质量低而影响了整体反演精度。

表3 HKQT站不同窗口长度反演精度和数据损失率对比Table 3 Comparison of inversion accuracy and data loss rate of different window lengths in HKQT station

图9 HKQT站多模多频GNSS联合潮位反演以及误差分析Fig.9 Multi-GNSS sea level retrieval series and error analysis of HKQT

6 结论

针对GNSS-MR 技术在实际用于潮位监测中存在精度不高和时间分辨率低两大问题，本文提出一种基于IGGIII 模型的稳健回归方法，在综合考虑海面动态变化和对流层延迟的情况下对多模多频GNSS-MR 反演结果进行融合，有效改善了潮位反演结果的精度和时间分辨率。本文主要得出如下结论：

（1）在潮位反演中多模多频数据相较于单系统单频数据可以大大提高反演结果的时间分辨率，但即使经过海面动态改正和对流层延迟改正，反演结果的精度依旧没有得到有效提高。在稳健回归模型下考虑海面动态变化误差和对流延迟，利用GNSS-MR 技术进行潮位反演，BRST站和HKQT站反演结果RMSE分别为12.43 cm和7.09 cm，与传统方法相比精度提高约40%—70%。

（2）传统方法获取潮位反演值由于SNR 数据采样时间不规律，潮位值也无法实现均匀采样。在基于稳健回归模型的数据融合算法中可以在特定的时间窗口内建立方程求解潮位值，BRST 站和HKQT站均可实现10 min的均匀采样，这大大提高了GNSS-MR 技术用于潮位监测的时间分辨率，满足了潮位监测需求。

利用这种基于IGGIII 模型的稳健回归方法对潮位反演值进行融合目前也存在一定局限性：（1）IGGIII模型也无法完全去除异常值对结果的影响，系统内存在的低质量信号会影响整体反演精度。（2）某一个时间点可能存在卫星数量少或者信号质量差的情况，这会导致这一时间段状态方程无法求解从而导致潮位数据缺失。（3）反演潮位时获得的潮位值是根据一段时间的SNR 序列获得，一般将其与SNR 序列的平均时刻对应；而验潮站测量结果也为短时间的平均值，因此反演值和观测值之间会存在时间尺度上的表示性误差。对于此类误差，需要在未来进一步研究。