本文内容来源于《测绘学报》2022年第11期（审图号GS京（2022）1229号）
多模多频GNSS-IR水位反演中的频间偏差分析及改正王笑蕾, 何秀凤, 宋敏峰, 陈殊, 牛紫瑾
河海大学地球科学与工程学院, 江苏 南京 211100
基金项目：江苏省自然科学基金(BK20190496);国家自然科学基金(42004018;41830110);中央高校基本科研业务费专项资金(B200202015)
摘要：准确的水位监测对于水资源调控、水灾监控及气候气象研究十分重要。近年来, 随着全球导航卫星系统的不断发展, 一种GNSS干涉遥感水位反演技术被提出。目前, GNSS-IR水位反演中主要有3类误差: 高度变化误差、高度角弯曲误差和对流层延迟误差, 相关的改正方法也被陆续提出。本文研究了4个GNSS跟踪站(HKQT、SW50、SW51和SW52)的多模多频反演结果, 发现除了已经发现的3类误差外, 还存在一种明显的频间偏差。现有研究对于GNSS-IR频间偏差的研究很少, 尚未达成将其归为GNSS-IR误差源的共识。为了进一步挖掘频间偏差的特性、加深对该误差的认识, 本文计算了不同信号的有效高度(reflector height, RH), 发现不同频率信号的RH反演值间存在差距; 将该差距与波长进行对比分析, 发现该差距与信号波长存在显著的线性关系(相关系数>95%)。该偏差量级在分米级, 表现出了明显的频间偏差系统特性。针对发现的频间偏差及其相关特性, 本文提出了相关的误差改正方法, 结果表明: 顾及频间偏差的改正结果比未顾及频间偏差的反演结果的均方根误差高1.5~12 cm。同时, 改正后的融合反演值的精度较未改正的独立信号的反演值提高了30%~80%;精度的提高得益于多模多频信号提供的大量冗余数据以及对各类误差(包括系统误差、粗差和随机误差)的正确处理。
关键词：GNSS-IR 水位反演 频间偏差 误差改正
引文格式：王笑蕾, 何秀凤, 宋敏峰, 等. 多模多频GNSS-IR水位反演中的频间偏差分析及改正[J]. 测绘学报，2022，51(11)：2328-2338. DOI: 10.11947/j.AGCS.2022.20210461 WANG Xiaolei, HE Xiufeng, SONG Minfeng, et al. Analysis of inter-frequency bias in multi-mode multi-frequency GNSS-IR water level retrieval and correction method[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(11): 2328-2338. DOI: 10.11947/j.AGCS.2022.20210461
阅读全文：http://xb.chinasmp.com/article/2022/1001-1595/20221109.htm准确的水位监测对于水资源调控、水灾监控及气候气象研究十分重要。传统水位计监测技术存在造价较高或需人工维护的缺点，而且需要附加技术才能实现基准统一[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误差源的共识。为了进一步挖掘频间偏差的特性、加深对该误差的认识，本文分析此类频间偏差的特性，并提出相关的误差改正方法。
当仅存在一次多路径反射条件下，直射信号与反射信号之间的路径差为D=2h(sin e)[3-5]；其中，h为反射面与天线中心之间的垂直距离，称为有效高度(reflector height, RH)；e为高度角。因此，二者之间的相位差Δφ为[3-5]式中，ωφ为角频率。如果忽略输入数据对应时间中h的变化，即在反射面在一定时间内静止的假设下，式(2)通过求导可写为式中，为静态假设下获得的RH。式(3)简单变形，可得=(λf)/2。对输入SNR数据进行趋势项剔除、高度角和方位角截取后，利用Lomb-Scargle Periodogram(LSP)方法对残余SNR弧段进行频谱分析可以获得频率f，再根据=(λf)/2获得。由于易获得，因此反演处理中，一般先通过频谱分析获得，再从值中扣除相关误差，最后通过站点坐标、WGS-84和水位基准间的转换，即可获得水位反演值[3-4]。根据已有研究，本文汇总了已有共识的3类GNSS-IR系统误差，并给出相关的误差改正或削弱方法。由于水面是处在变化之中的，假设反射面静止会引入误差。考虑反射面处于变化状态中，引入高度变化速率。式(2)重新求导，可得式中，，为高度角变化速率。式(3)和式(4)分别是在静态假设和非静态假设下求导，二式左边项相同，故可得简单变形后，可得和h间的关系。因此，通过在反演序列中扣除改正项，即可改正高度变化误差。改正的关键是ḣ的准确求解：一般是先对序列进行曲线拟合，再对拟合曲线求导获得，详见文献[4, 16—17]；也可将ḣ作为未知参数进行求解获得，详见文献[23]。此类误差主要与水位变化速率有关，在水位变化速率较大的区域的误差量级可达分米级，甚至米级；除了水面较为平静的状态(ḣ < 10-6 m/s)，一般都必须改正该误差。大气折射效应会引起信号弯曲，从而导致计算高度角与实际高度角之间的偏差，而这一高度角偏差又会进一步导致反演结果的偏差。文献[18]发现该高度角弯曲误差，并提出利用大气折射改正公式(6)改正高度角偏差δe的方法式中，T为温度，单位为℃；P为气压，单位为mb。也可利用其他的大气折射公式进行高度角偏差改正[27]。利用改正后的高度角进行反演计算，即可削弱此误差。该误差与高度角有关，高度角越小，误差越大；该误差会引起反演结果的尺度误差及均值误差。对于5°~20°高度角区间，反演误差量级在亚厘米至毫米级，可视精度要求情况进行改正。大气折射效应会引起信号延迟，由于反射信号比直射信号传播路径更长，因而受到的延迟影响更大，从而造成反演结果的误差。文献[19, 24—26]利用GPT2w模型和VMF1映射模型来修正该误差。先计算反射信号和直射信号间的相对延迟值τT式中，Δτhz=τhz(-h)-τhz(0)，Δτwz=τwz(-h)-τwz(0)，Δτhz和Δτwz分别为高度h处与高度0处的相对干延迟模型值和相对湿延迟模型值；mh(e)和mw(e)分别为干映射值和湿映射值。大气延迟对反演结果的影响值为式中，ΔhT为延迟误差改正量。该误差与RH大小及高度角有关，该误差会引起反演结果的尺度误差及均值误差。一般情况下(有效高度10 m以内，最低高度角选择5°)，影响量级在厘米级至毫米级；可视精度要求情况进行改正。对于某些RH较大的站点(例如安置于灯塔上部、高坝顶部、高楼顶部等的站点)或者极低高度角的情况，该误差必须改正。本文选取了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附近有一个水位监测站，提供每天一次的水库水位测量值。 |
| Fig. 1 GNSS net of Shuangwangcheng Reservoir |
|
Tab. 1 Frequency, wavelength, code and SNR type of the used signals表选项
|
| 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可提供实测的潮位数据。 |
| Fig. 3 HKQT site surroundings |
|
利用1.1节所述原理，先剔除SNR序列的趋势线；再根据选定的高度角区间和方位角区间截取SNR弧段；然后根据1.2.2节的公式改正高度角弯曲偏差，获得改正后的高度角值；最后对SNR弧段进行LSP分析并获得对应的RH值。根据1.1节所述原理，可以通过站点坐标及基准转换参数获得特定常数，然后通过常数减去RH就可以得到水位反演值。然而，由于双王城水库的水位观测结果并不固定在特定框架下，即其起算的基准面无法确定，因此无法通过基准转换参数获得常数值。在文中，通过RH与水位观测值之间的对应关系，计算获得对应的常数。在计算对应常数的过程中，笔者发现不同的信号对应的常数不同，且差距较大，即不同信号的反演值间存在均值偏差。本文统计了不同信号对应的常数、反演个数、反演值与实测值间的均方根误差(root mean squared error, RMSE)和相关系数(root mean squared error, CORR)。同时为了更好地表现信号间的频间偏差，本文统计了各信号常数与GPS L1信号常数之间的差距，相关结果记录见表 2。各个信号的反演序列如图 5所示。
Tab. 2 Statistics of RHs of Shuangwangcheng sites
表选项
|
| Fig. 5 Retrievals of sites of Shuangwangcheng Reservoir |
|
图 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]，这些硬件造成的误差，理论上与频率线性相关。经过与2.2节所述相同的处理后，获得了HKQT站各个信号的RH值。根据1.1节所述原理，可以通过站点坐标及基准转换参数获得特定常数，然后通过常数减去RH就可以得到潮位反演值。IOC(Intergovernmental Oceanographic Commission, http://www.ioc-sealevelmonitoring.org/)提供了HKQT站的三维坐标结果以及相应框架转换参数，根据相关数据计算获得用于估计潮位的常数值为7.797 m。然而，根据2.2节的分析，频间偏差是存在的。因此，通过RH与潮位观测值之间的对应关系，计算获得对应的常数。在计算对应常数的过程中，发现不同的信号对应的常数不同，且差距较大。本文统计了不同信号对应的常数、反演个数、反演值与实测值间的RMSE和CORR。同时为了更好地表现信号间的频间偏差，统计了各信号常数与L1信号常数之间的差距，相关结果记录见表 3。各个信号的反演序列及反演序列与实测序列的偏差如图 6所示。
Tab. 3 Statistics of RHs of site in Quarry Bay, Hong Kong表选项
|
| Fig. 6 Retrievals and bias of HKQT |
|
图 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。考虑高度变化误差、对流层延迟误差及频间偏差，则和h之间的关系为本文使用了窗口长度为2 h、窗口步长为10 min的滑动窗口将序列进行切分。对于属于时间窗口i，信号j，轨迹l的反演值，可写为i, j, l，则窗口内的方程组可写为式中，ei, j, l、ΔhTi, j, l、ti, j, l分别为i, j, l对应的高度角、对流层延迟误差及时刻；λj为信号j的波长，λ1为L1对应的信号波长；Ns为测站接收的信号类型的个数；Ni, s, t为时间窗口i内最后一个信号类型的轨迹个数；下标为i、Ns、Ni, s, t的参数代表的是时间窗口i内最后一个信号最后一个轨迹对应的参数；hi和ḣi为对应时间窗口i所输出的RH估值和RH变换速率的估值。为了保证解算的正确，本文更推荐在获得各个独立信号RH的情况下，根据频间偏差推算获得对应的a，建立下述方程组式中，为式(10)或式(11)等号左边的向量；Ai为系数矩阵；，为未知参数向量。本文使用抗差估计方法来解算式(12)，获得对应的未知数估值及其不确定度，即中误差。该方法能够较好地处理RH反演序列中的粗差，很适合处理GNSS-IR反演数据中粗差较多的情况[29]。利用3.1节所述方法，对第2节获得的各个信号的RH结果进行处理，获得各个窗口的hi输出结果，并将其转换至水位基准，获得对应的水位序列及其不确定度。值得注意的是，由于双王城站点的GPS L2P信号有双波峰现象，因此相关RH反演值不参与改正方程组的建立。本文相关站点误差改正后的多模多频融合水位反演结果如图 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 |
|
本文分析了频间偏差的相关特性，发现了频间偏差与波长之间存在线性关系这一规律；同时，本文顾及高度变化误差、对流层延迟误差及频间偏差，提出了相关的误差改正方法。本文结果表明：顾及频间偏差的改正结果比未顾及频间偏差的反演结果的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-sealevelmonitoring.org/)提供的实测潮位数据。
作者简介第一作者简介：王笑蕾(1991—), 女, 博士, 副教授, 研究方向为GNSS遥感。E-mail: chd_wxl@qq.com
通信作者：何秀凤, E-mail：xfhe@hhu.edu.cn