多模多频GNSS-IR水位反演中的频间偏差分析及改正

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

河海大学地球科学与工程学院，江苏 南京 211100

准确的水位监测对于水资源调控、水灾监控及气候气象研究十分重要。传统水位计监测技术存在造价较高或需人工维护的缺点，而且需要附加技术才能实现基准统一[1]。同时，卫星测高技术无法获得近岸水面高度且时间分辨率较低[2]。随着全球导航卫星系统(GNSS)的不断发展与完善，一种GNSS干涉遥感(GNSS-interferometry reflectomety,GNSS-IR)技术被发现可以用来进行水位(包括潮位)[3-7]、雪深[8-11]、土壤湿度[12-14]及地表冻融[15]等环境参数的监测。其中，GNSS-IR水位监测技术只基于沿岸的测量型接收机，根据其信噪比(signal-to-noise ratio,SNR)或其他观测数据中的干涉特性，便可完成水位监测[3-4]。利用新兴的GNSS-IR技术，可以实现全自动、小成本、高时间分辨率、长期连续的水位监测，并且监测结果自动固定在稳定框架下[5]。

文献[3—4]提出GNSS-IR技术，并由此发展出一套经典的GNSS-IR水位反演理论。目前，地基GNSS-IR技术主要有两类误差源：反射面高度变化和大气折射效应。文献[4]发现海面起伏会引起反演结果中的高度变化误差，并提出了相应的改正算法。文献[16]利用5个GNSS站进行了潮位反演，并且利用经验潮波系数来改进误差改正算法。文献[17]基于潮波分析模型，提出了一种不需要经验潮波系数的海面高度拟合方法，进一步完善了高度变化误差的改正算法。大气折射效应则会引起两类误差：一类是由于大气折射引起的信号弯曲带来的高度角弯曲误差；另一类是由于大气折射引起的对流层延迟误差。文献[18]针对大气折射引起的高度角弯曲误差，利用折射改正公式来改正高度角偏差。文献[19]利用VMF1映射函数模型和GPT2w对流层延迟模型来改正对流层延迟误差。

2020年6月23日，北斗三号全球导航卫星系统星座建设完成。同时，GPS和GLONASS正在推进现代化，Galileo在不断发展，各个区域系统也在持续改进和建设。多系统发展使得GNSS-IR技术拥有了丰富的信号源，而四系统的多模多频信号已经均被证明可以监测水位。文献[3]证明了GPS信号的水位可监测性。文献[10]证实了GLONASS系统的水位可监测性。文献[21]证明了北斗数据的水位可监测性。文献[22]利用GPS、GLONASS和北斗信号进行了水位监测试验。文献[23]利用GPS、GLONASS、Galileo和北斗系统进行了水位反演，分析了四系统的监测性能，并提出了一套多模多频融合反演算法。多模多频信号的引入，提升了时间采样，也为提高精度提供了可能。但是多模多频数据也带来了更多的误差，需要对这些误差进行正确认识及合理处理，以获得最优反演性能。

为了更好地挖掘GNSS-IR水位反演误差，本文研究4个GNSS连续跟踪站(HKQT、SW50、SW51和SW52)的多模多频反演结果，发现除了上述已经发现的3类误差，还存在一种明显的频间偏差。而目前对GNSS-IR频间偏差的研究极少，尚未达成将其归为GNSS-IR误差源的共识。为了进一步挖掘频间偏差的特性、加深对该误差的认识，本文分析此类频间偏差的特性，并提出相关的误差改正方法。

1 GNSS-IR原理及系统误差

1.1 GNSS-IR原理

当仅存在一次多路径反射条件下，直射信号与反射信号之间的路径差为D=2h(sine)[3-5]；其中，h为反射面与天线中心之间的垂直距离，称为有效高度(reflector height,RH)；e为高度角。因此，二者之间的相位差Δφ为[3-5]

Δφ=2πD/λ=4πh(sine)/λ

(1)

式中，λ为信号波长。相位差中隐藏了一个频率信息f

(2)

式中，ωφ为角频率。如果忽略输入数据对应时间中h的变化，即在反射面在一定时间内静止的假设下，式(2)通过求导可写为

(3)

1.2 GNSS-IR系统误差

根据已有研究，本文汇总了已有共识的3类GNSS-IR系统误差，并给出相关的误差改正或削弱方法。

1.2.1 高度变化误差

(4)

(5)

1.2.2 高度角弯曲误差

大气折射效应会引起信号弯曲，从而导致计算高度角与实际高度角之间的偏差，而这一高度角偏差又会进一步导致反演结果的偏差。文献[18]发现该高度角弯曲误差，并提出利用大气折射改正公式(6)改正高度角偏差δe的方法

(6)

式中，T为温度，单位为℃；P为气压，单位为mb。也可利用其他的大气折射公式进行高度角偏差改正[27]。利用改正后的高度角进行反演计算，即可削弱此误差。该误差与高度角有关，高度角越小，误差越大；该误差会引起反演结果的尺度误差及均值误差。对于5°～20°高度角区间，反演误差量级在亚厘米至毫米级，可视精度要求情况进行改正。

1.2.3 对流层延迟误差

大气折射效应会引起信号延迟，由于反射信号比直射信号传播路径更长，因而受到的延迟影响更大，从而造成反演结果的误差。文献[19,24—26]利用GPT2w模型和VMF1映射模型来修正该误差。先计算反射信号和直射信号间的相对延迟值τT

(7)

(8)

式中，ΔhT为延迟误差改正量。该误差与RH大小及高度角有关，该误差会引起反演结果的尺度误差及均值误差。一般情况下(有效高度10 m以内，最低高度角选择5°)，影响量级在厘米级至毫米级；可视精度要求情况进行改正。对于某些RH较大的站点(例如安置于灯塔上部、高坝顶部、高楼顶部等的站点)或者极低高度角的情况，该误差必须改正。

2 频间偏差分析

2.1 站点介绍

本文选取了4个GNSS连续跟踪站。其中3个站点(SW50、SW51和SW52)来自山东双王城水库的GNSS大坝监测系统，另一个站点(HKQT)来自香港卫星定位参考站网(satellite positioning reference station network，SatRef)。

双王城水库位于山东潍坊寿光市，水库大坝轴线总长为9.636 km，坝高约12.5 m，最大库容量达到6150万m3，是南水北调东线胶东干线工程的重要调蓄水库。为了监测水库大坝的稳定性，在大坝周围建立了3个监测站(SW50、SW51和SW52)和1个基准站(SW43)，如图1所示。该GNSS网配备CHC N72接收机及CHCC220GR天线，记录每15 s一次的GPS(L1和L2)和北斗(B1、B2和B3)观测值。相关信号对应的频率、波长及SNR类型见表1。监测站SW50、SW51和SW52均可接收到来自水库水面的反射信号。本文选取5°～15°高度角区间SNR弧段进行反演。同时， 根据图2中的反射区分布，选定对应的水域方位角为：SW50的方位角区域为60°～165°，SW51的方位角区域为130°～230°，SW52的方位角区域为270°～350°。算例选取2017-11-22—2018-01-10期间的数据进行分析。在SW50附近有一个水位监测站，提供每天一次的水库水位测量值。

图1 双王城水库GNSS监测网Fig.1 GNSS net of Shuangwangcheng Reservoir

图2 双王城站点信号反射区Fig.2 FFZs of sites of Shuangwangcheng Reservoir

HKQT站位于香港鲗鱼涌，属于香港卫星定位参考站网，站点环境如图3所示。HKQT站配有Trimble NetR9接收机及Trimble 59800.0型天线。本文使用其接收到的四系统——GPS(L1C/A、L2P、L2C和L5)、GLONASS(G1C、G1P、G2C和G2P)、Galileo(E1、E5、E7和E8)和北斗(B1和B2)5 s采样的观测数据进行反演。相关信号对应的频率、波长及SNR类型见表1。算例选取5°～15°高度角区间SNR弧段进行反演。同时，根据图4中的反射区分布，选定对应的海域方位角为-60°～105°。算例选取2017年DOY 224—DOY 244期间的数据进行分析。距离HKQT站点2m处有一隶属Intergovern-mental Oceanographic Commission(IOC,http:∥www.ioc-sealevelmonitoring.org/)的验潮站Quarry Bay可提供实测的潮位数据。

图3 HKQT站点环境Fig.3 HKQT site surroundings

2.2 双王城站点反演序列

图4 HKQT站点信号反射区Fig.4 FFZs of HKQT

表1 使用信号对应的频率、波长、码及SNR类型

图5中，L1、B1、B2及B3的反演与水位测量对应极好；而L2的反演序列表现出双序列现象，这是由于L2信号在LSP谱图中表现出的双波峰现象导致[18]。由表2可知，L1、B1、B2及B3的反演精度很好，RMSE在1.74～11.63 cm，且大部分集中在3～6 cm。这种高精度结果主要是由于水库表面平静、反射条件好、SNR多径特性明显导致。表2中，对于3个监测站，L2信号的常数比L1信号大0.19 m，B1的常数与L1相似，B2的常数比L1信号大0.20～0.21 m，B3的常数比L1信号大0.13～0.15 m；这意味着，L2信号的RH比L1信号整体小0.19 m，B1的RH与L1相似，B2的RH比L1信号小0.20～0.21 m，B3的RH比L1信号小0.13～0.15 m。比较各个信号间的偏差值及表1中各个信号的波长值，可以发现偏差大小与波长呈显著线性关系(CORR=98.82%)。因此，偏差可表示为a×δλ，其中a为常数。根据拟合，双王城站点3个站反演结果频间偏差的a=3.6。这种偏差量级在分米级，明显不是相位中心偏移或是随机误差导致；而且其呈现了明确的系统特性(与波长相关)。因此，可以确定有一类与波长相关的系统误差确实存在。而波长λ等于光速与信号频率的比值，即该偏差也与频率呈现系统特性，故其为一种频间偏差。这种频间偏差可能是由电磁偏差、菲涅耳反射系数的组合相互作用及天线方向图和两个圆极化之间的功率偏移引起的[28]，这些硬件造成的误差，理论上与频率线性相关。

图5 双王城水库站点反演序列 Fig.5 Retrievals of sites of Shuangwangcheng Reservoir

表2 双王城站点RH反演结果统计表

2.3 香港鲗鱼涌站点反演序列

经过与2.2节所述相同的处理后，获得了HKQT站各个信号的RH值。根据1.1节所述原理，可以通过站点坐标及基准转换参数获得特定常数，然后通过常数减去RH就可以得到潮位反演值。IOC(Intergovernmental Oceanographic Commission,http:∥www.ioc-sealevel monitoring.org/)提供了HKQT站的三维坐标结果以及相应框架转换参数，根据相关数据计算获得用于估计潮位的常数值为7.797 m。然而，根据2.2节的分析，频间偏差是存在的。因此，通过RH与潮位观测值之间的对应关系，计算获得对应的常数。在计算对应常数的过程中，发现不同的信号对应的常数不同，且差距较大。本文统计了不同信号对应的常数、反演个数、反演值与实测值间的RMSE和CORR。同时为了更好地表现信号间的频间偏差，统计了各信号常数与L1信号常数之间的差距，相关结果记录见表3。各个信号的反演序列及反演序列与实测序列的偏差如图6所示。

表3 香港鲗鱼涌站点RH反演结果统计表

图6中，GPS、GLONASS、Galileo和北斗各个信号的反演序列与潮位实测序列对应良好。根据表3，HKQT站各个信号的反演精度在低分米级，RMSE在15～23 cm。表3中，波长与L1近似的信号，其对应的常数与L1信号的常数相同或近似；而与L2波长近似的信号，其对应的常数较L1的常数低0.12 m。这意味着，与L1波长近似的信号的RH与L1信号的RH近似，L2信号及与其波长近似的信号的RH比L1信号整体小0.12 m。比较表3中各个信号间的偏差值及表1中各个信号的波长值，可以发现：偏差大小与波长呈显著线性关系(CORR=96.97%)。偏差可表示为a×δλ，其中a为常数。根据拟合，HKQT站点对应的a=2.3。

图6 HKQT站点反演序列及偏差Fig.6 Retrievals and bias of HKQT

3 频间偏差改正

3.1 误差改正方法

(9)

(10)

(11)

式(10)或式(11)可简写为

(12)

3.2 改正结果

利用3.1节所述方法，对第2节获得的各个信号的RH结果进行处理，获得各个窗口的hi输出结果，并将其转换至水位基准，获得对应的水位序列及其不确定度。值得注意的是，由于双王城站点的GPS L2P信号有双波峰现象，因此相关RH反演值不参与改正方程组的建立。本文相关站点误差改正后的多模多频融合水位反演结果如图7所示。

图7 各站误差改正后的多模多频融合水位反演序列Fig.7 Multi-GNSS combined retrievals with error correction of sites

经过误差改正处理后，各个测站的水位反演结果与实测水位结果的对应关系较好(CORR>99%)。各个测站的误差改正后得出多模多频融合水位反演序列的RMSE，SW50为1.56 cm，SW51为0.67 cm，SW52为0.99 cm，HKQT为6.98 cm。与第2节所述的各个信号的反演结果的RMSE相比，改正后的融合反演值的精度提高了30%～80%。同时，本文计算了未顾及频间偏差， 只顾及其他误差的改正结果，得出相关的多模多频融合水位反演序列的RMSE，SW50为3.15 cm，SW51为5.90 cm，SW52为2.24 cm，HKQT为18.97 cm。与顾及频间偏差的改正结果相比，未顾及频间偏差的反演结果的RMSE高1.5～12 cm。具体而言，SW50和SW52站点改正频间偏差后提高了约1.5 cm，SW51站提高了约5 cm，而HKQT站提高了约12 cm。这是由于SW51站的北斗反演值个数显著大于SW50和SW52站，因此，在反演结果改正中，频间偏差对SW51的误差影响更大；而HKQT站，由于其多系统多频信号均可用，因此在反演结果改正中，频间偏差对融合反演结果的误差影响远大于双王城各站点。图8显示了HKQT站未顾及频间偏差的改正结果，可以看出，未改正频间偏差时，融合反演结果的精度变差、粗差增多。

图8 HKQT站未顾及频间偏差的多模多频融合水位反演结果及其不确定度Fig.8 Multi-GNSS combined retrievals without consideration of inter-frequency bias and corresponding uncer-tainties at HKQT

4 结论与讨论

本文分析了频间偏差的相关特性，发现了频间偏差与波长之间存在线性关系这一规律；同时，本文顾及高度变化误差、对流层延迟误差及频间偏差，提出了相关的误差改正方法。本文结果表明：顾及频间偏差的改正结果比未顾及频间偏差的反演结果的RMSE高1.5～12 cm；同时，改正后的融合反演值的精度较未改正的独立信号的反演值提高了30%～80%；精度的提高得益于多模多频信号提供的大量冗余数据及对各类误差(包括系统误差、粗差和随机误差)的正确处理。值得注意的是，改正方法中的相关方程组是经典平差函数模型的表达方式，可借鉴测量数据处理方法进行模型扩充和解算方法扩充；在有新的误差被发现后，也可以根据误差特性进行扩展。第3.1节所述误差改正模型具有很好的扩展性。

除了系统误差，本文还表现出了一类L2P信号独有的反演结果特性。双王城站点中，L2的反演序列表现出双序列现象，这是由于L2信号在LSP谱图中表现出的双波峰现象导致[18]。而HKQT站则并未表现出L2P双波峰现象导致的反演结果很差的情况。对于此类表现在L2P上的独特情况——有时出现、有时不出现“双波峰”的情况，其根本原因还没有定论，猜测与P码调制、L1干扰或对L2的半无码跟踪有关。这类独特现象还需要更深一步的研究。

此外，本文只研究了GNSS-IR水位反演的情况，后续可以对雪深、冰厚等参数进行GNSS-IR反演，以探讨在其他参数反演中是否存在频间偏差。同时，对于形成频间偏差的具体原因还未被发现；而该误差成因对于更好地认识及改正频间偏差非常重要。第3.1节中，利用a×δλ来表示频间偏差大小，而表2和表3中，某些信号并不完全符合a×δλ的数值，有约1～3 cm的差距。文献[28]指出，频间偏差的产生可能与电磁偏差、菲涅耳反射系数的组合相互作用及功率偏移等接收机硬件有关，而这些硬件造成的误差理论上与频率线性相关。因此笔者认为，对于特定站点，a取决于接收机硬件，它对于不同信号是唯一的，相关差距应该是随机误差导致的。但是，该结论仅基于现有研究推断而来，a与信号特征是否有关依然待研究。

目前，GNSS-IR技术仍处在快速发展阶段，还有许多的误差特性还未理清，例如地形引起的误差[30]，不同高度角SNR质量引起的误差[31]等，也有许多的误差现象还没有被发现成因。需要后续研究者对此进行深入研究，以实现更好的GNSS-IR监测性能，推进GNSS-IR技术的实际应用进程。

致谢:感谢香港卫星定位参考站网提供的GNSS数据，感谢IOC(Intergovernmental Oceano-graphic Commission, http:∥www.ioc-sealevel monitoring.org/)提供的实测潮位数据。