利用多系统GNSS干涉反射测量估计长江巴东水位变化

陈昊晟 金双根，2 业明达 郭孝祖

0 引言

全球导航卫星系统(GNSS)，是目前所有全球导航卫星系统和区域增强系统的统称[1]，包括美国的全球定位系统(GPS)、俄罗斯的格洛纳斯系统(GLONASS)、中国的北斗导航卫星系统(BDS)、欧盟的伽利略导航卫星系统(GALILEO)，以及日本准天顶卫星系统(QZSS)和印度区域导航卫星 (IRNSS)等区域增强系统.GNSS持续地发射L波段信号，为用户提供高精度的定位、导航、授时信息，因而被广泛应用于国土测绘、交通运输、公共安全等领域[2].随着GNSS应用的不断拓展，GNSS反射信号可以提供地物反射表面的特征和地球物理参数，即GNSS反射测量技术(Global Navigation Satellite System Reflectometry,GNSS-R).传统的GNSS-R技术是一种通过大地测量型GNSS接收机获取反射信号来反演周围环境参数的技术，现在已经被作为一种全新的遥感手段，应用于测量海面高、土壤湿度、积雪厚度等[3].

水位信息是水文、气象等重要基础数据，其精确实时监测也事关人民生活财产安全以及国家社会安定.传统的水位监测主要通过安装水位计来获得特定位置的水位数据，但该方法需要消耗大量的人力、物力和财力，且在人烟稀少的偏远位置更受到限制.GNSS-R技术为水位监测提供了一种新的思路.1988年，欧州航天局(ESA) 提出GPS的L波段信号可以作为散射计进行海洋遥感研究[4].Martin-Neira等提出利用GNSS反射信号进行海洋测高的理念，使用原本被视为噪声的反射信号能够反演海面粗糙程度，利用从海面反射的反射信号和直射信号相结合来估计海面高度[5-7].2013年，Larson等[8]提出一种使用SNR(信噪比)数据的测高方法，称之为GNSS干涉反射测量(GNSS-IR)，该技术通过处理包含GNSS直射信号与反射信号干涉信息的信噪比数据，反演水位或其他环境参数.2014年， Löfgren等[9]对GPS与GLONASS信号分别运用GNSS-IR和相位分析两种方法来反演水位，两种方法的精度都达到了厘米级，且GNSS-IR反演水位的结果在粗糙海面条件下的表现要好于相位延迟分析的结果，在风速大的情况下水位反演结果仍能保证较高的精度.2019年，Chen等[10]利用小波去噪技术，基于GNSS-IR方法计算的海平面高度估计值重建潮汐波形，与初始估计结果相比，均方根误差由32 cm降低到19 cm，相关系数从0.83提升到0.95.2020年，南阳等[11]在西安灞河悬挑平台上使用3组不同的仪器进行GNSS-IR测高实验，分别用来探究频段、高度角区间、扼流圈天线对 GNSS-IR水面测高精度的影响，得到了最佳RMS结果为1.04 cm的精度，实验表明L1频段的信噪比数据更适合测高反演，同时给出了高度角区间选取的准则.

随着全球导航卫星系统应用的发展，GNSS-IR能够高精度、连续地监测水位，具有广泛的应用前景,且利用GNSS-IR监测水位的观测成本低、精度高，对于全球海平面变化的相关研究以及重点湖泊河流水位信息的监测具有重要意义.本文利用多系统GNSS的信噪比反演长江上游巴东站水位变化，并与实测结果进行对比和分析，探究GNSS-IR技术对长江水位监测的可行性，以期为未来GNSS-IR测站对于河流湖泊沿岸以及海岸海平面的水位观测提供参考.

1 GNSS-IR反演水位原理与方法

GNSS-IR测高技术利用GNSS系统多路径信号SNR数据来反演水位高度.GNSS接收机主要接收直射信号，当GNSS测站位于水面附近时，接收机能同时接收来自天空的直射信号以及经水面反射的反射信号.GNSS信号经过水面发生反射后，沿镜面反射的相干分量与其他方向不规则的散射分量共同组成反射信号，最终和直接信号一起被接收机所获取.

图1 GNSS-IR原理示意Fig.1 GNSS-IR schematic diagram

如图1所示，GNSS信号经过平静水面后发生反射，绿色部分为接收机获取的直接信号，橘红色部分为反射信号，两种信号会发生干涉，它们到达天线所走过路程差值可用天线与反射水面的相对高度h和卫星高度角e表示：

δ=2hsine,

(1)

将其转换为相位差可表示为

(2)

由式(2)可知，反射信号与直射信号的相位差ψ是载波波长λ、卫星高度角e、天线到反射面的距离h的函数，并和高度角的正弦值成正比[11].

信噪比观测值SNR是常规测量型GNSS接收机的附带观测值[12]，通常被用作表征载波信号质量好坏的指标：

(3)

PR表示接收信号的功率，N表示噪声功率.

在高度角较低的情况下，多路径效应加剧，从而引起SNR的变化[13-14]，其具体表现为低高度角区域的振荡现象.因此提取低高度角的SNR值，能够更好地利用多路径效应来分析反射信号所包含的水面信息.由于以dB-Hz为单位的SNR观测值对卫星高度角的变化响应较小，振幅波动起伏变化较小，一般将其线性化转换:

SNRvolts/volts=10SNRdB-Hz/20.

(4)

当GNSS接收机放置于水边时，SNR观测值可以表示为

(5)

式(5)中，SNR观测值是直射信号分量幅度Ad、多路径反射信号幅度Ar和两者之间相位差ψ的函数.由式(2)和式(5)可知，当卫星上升(或下降)时，直射信号和反射信号的相位差ψ随着高度角e而变化，从而导致受信号干涉影响的SNR值时而加强、时而减弱.由于多路径和天线增益的原因，低高度角下SNR值震荡明显.

由于测量型接收机的天线被设计为右旋圆极化，会接收到更多的右旋圆极化直接信号，而抑制多路径效应下的左旋圆极化的反射信号[15]，因此，直射信号的振幅远大于反射信号的振幅，即Ad≫Ar.通过低阶多项式来拟合SNR中的直射信号，再通过去除趋势项来获取包含反射信号信息的SNR残差序列[16].结合式(2)和式(5)，低高度角的残差序列可表示为

(6)

δSNR=Acos(2πft+φ).

(7)

通过信噪比的残差序列分析获得其频率f，从而进一步获得天线至水面的高度h.接收机记录的SNR数据通常会在时间上均匀采样，但仰角的正弦值分布不均匀，无法采用传统的快速傅里叶变换(Fast Fourier Transform，FFT)进行频谱分析，而Lomb-Scargle频谱分析法(简称为LSP谱分析)可以处理不均匀间隔的样本，因而采用LSP谱分析，结果如图2所示.

图2 Lomb-Scargle频谱分析 (DOY:143；PRN:G28；频段：GPS L1)Fig.2 Lomb-Scargle spectrum analysis (DOY:143，PRN:G28,Frequency：GPS L1)

对SNR残差序列进行LSP频谱分析，可以得到式(7)中频率f，根据

(8)

得到天线相位中心到水面的相对高度h，再结合实际测量和计算得到的天线相位中心的绝对高度，从而计算出水面的绝对高度，实现GNSS-IR技术对水面高度的测量.

具体操作流程如图3所示.

图3 GNSS-IR水位监测处理流程Fig.3 GNSS-IR water level monitoring & processing

2 巴东站GNSS-IR水位反演与分析

本次实验数据采集位于长江水利委员会水文局长江三峡水文水资源勘测黄陵庙分局巴东水位站(110.383 574 8°E ，31.052 923 7°N)，该站坐落于湖北恩施州巴东县信陵镇沿江大道38号.监测时间段为2021年年积日143日至172日(2021-05-22—2021-06-21)，连续观测30 d，接收机采样间隔设置为1 s.

接收机天线所接收的信号频段如表1所示.

本次实验获得了GPS、BDS以及GLONASS的SNR观测数据，水位计实测数据作为验证数据，该水位数据由长江三峡水文水资源勘测局提供.图4为巴东站GNSS-IR测站现场图.

表1 华测N72接收机天线所接收信号

图4 巴东站GNSS-IR测站Fig.4 Appearance of GNSS-IR station at Badong tation

首先提取观测文件中各个系统的5°～25°仰角的SNR数据，由于在低高度角的反射信号会引起SNR数据的振荡，因此可以通过去除SNR趋势项的方法获取SNR残差序列.再利用Lomb-Scargle频谱分析法提取残差序列中的主频率值，以此计算天线相位中心到水面的反射高度值，再根据天线相位中心的高度，最终计算出巴东站水面的实际高度，实现整个GNSS-IR的水位反演.

将不同系统的反演结果与巴东站实测水位数据进行对比和分析.后文出现的所有水位值结果均做了距平处理，使用的平均值为这30 d内最高水位和最低水位的平均，画图时所用的水位值都减去了该平均值.图5a—d中，红色线是由水位站水位计测量的真实水位结果，而蓝色点表示由GNSS-IR反演的水位结果.横轴表示UTC时下的年积日，纵轴表示水位的相对变化.

图5 GNSS-IR水位反演结果与水位站实测结果对比Fig.5 Comparison of water level between GNSS-IR results and in-sit measurements

通过图5可以发现，本次实验由GNSS-IR反演的水位结果和实测值具有比较好的一致性.和此前的GNSS-IR研究相比，此次实验监测的水面相对变化较大，但仍获得较为准确的水位监测结果.由图5a—d可以看出，GPS的L1频段、GLONASS的L1频段与BDS的B1频段都获得较好的水位反演结果，GLONASS的L2频段的反演结果相比之下稍差，其反演结果普遍比实测水位值要高.

对各系统反演水位的绝对误差情况进行分析(图6)，可以发现GPS L1、GLONASS L1以及BDS B1基本没有系统误差，而GLONASS L2频段反演结果与实测值相比偏大.

图7a—d为对实验所反演的水位结果与实际水位站测量结果的相关性分析，图中横轴表示的是由实验反演的水位相对变化值，而纵轴为由实测数据利用一维线性插值所得的水位相对值，蓝色圆点表示的是水位反演结果与实测值对应位置，红色的线表示散点的一维线性回归线.为了更直观地比较以上不同频段在水位测量的精度情况，对GNSS-IR不同频段反演结果计算了均方根误差(Root Mean Square Error，RMSE)以及决定系数(R2)，参考图7.

图6 GNSS-IR各系统反演水位绝对误差情况Fig.6 Absolute errors of water level retrieved by GNSS-IR systems

通过计算得到两者的回归系数R2与RMSE的值，可以发现4个频段的R2都在0.98以上；RMSE方面，除了GLONASS的L2频段达到了0.3 m，其他的GPS L1频段、GLONASS L1频段以及BDS的B1频段的RMSE值都在厘米级.综合比较，实验结果不仅可以反映水位变化的趋势，同时也具有较好的精度.

图7 不同系统水位反演结果与水位站实测结果相关性Fig.7 Correlation between the water level results retrieved by GNSS-IR systems and the in-situ measurements

综合比较散点图和相关性分析图(图8a、图8b)，GPS的L1频段、GLONASS的L1频段以及BDS的B1频段在水位反演实验获得的结果相关系数均大于0.99，且RMSE的值都在厘米级；GLONASS的L2频段相比之下稍差一点，但是相关系数也达到了0.982，RMES为0.303 m.联合所有系统的反演结果，相关系数约为0.995，RMSE为0.162 m(表2).

总体而言，巴东站GNSS-IR反演与实测水位数据吻合较好.由于弧段要求7°～23°，导致每日各系统的观测量受到限制，但是综合所有频段信号每日还是能获得最少6次，最多15次的水位观测.

GPS的L1频段、GLONASS的L1频段以及BDS的B1频段都达到了厘米级精度且精确度相当，而GLONASS的L2频段的反演结果明显不及其他3个频段的信号.这可能是由它们信号的频率范围造成的，除GLONASS的L2频段以外，其他3个频段频率均在1 550 MHz以上，而GLONASS L2频段的频率低于1 260 MHz.高频信号波长较小，更容易受到多路径效应的影响，所以SNR中反射信号的信息更多，对水位的反演也就越准确.

图8 GNSS-IR各系统联合反演水位反演结果Fig.8 Comparison between GNSS-IR joint inversion results and in-situ measurements a.water level;b.relative changes of water level

表2 GNSS-IR不同频段观测结果

3 结论

本文基于GNSS-IR观测数据对长江上游巴东水文站进行了水位反演，通过不同GNSS系统不同频段信号的水位反演结果与实测数据进行对比，获得以下结论：

1)GNSS-IR能够获得厘米级的水位变化精度，RMSE最低为6.43 cm；

2)GLONASS、GPS和BDS L1频段的估计结果较好，而GLONASS L2频段受频段信号频率范围的影响，其估计精度低于其他频段；

3)联合不同GNSS系统反演水位变化，提高了GNSS-IR反演水位的时间分辨率.

后续将进一步考虑GNSS-IR各卫星观测弧段高度角范围与GNSS-IR反演精度和频次的关系，通过对弧段的限制，平衡反演水位精度与反演结果的数据量，最终筛选出合适的弧段范围，保证在满足水位反演精度的条件下，尽可能多获得更好的水位观测，为实现未来GNSS-IR对于河流湖泊沿岸以及海岸海平面的水位观测提供支撑.

致谢：感谢巴东水位站陈文重主任、刘天成总工和刘峥鹏同志为本次研究提供的帮助.