库区GNSS-R 精细化反演水面高度及其验证研究

桑文刚，刘迎春，何秀凤，王昭然

(1.山东建筑大学 测绘地理信息学院,济南 250101；2.河海大学 地球科学与工程学院,南京 211100)

0 引言

水位信息对湖泊、大型水库的水资源优化调度与智能化管理具有极其重要的作用.近年来，水位监测方法正逐步由传统的人工水尺向浮子式水位计、压力式水位计、超声波水位计等升级换代，在自动化程度、测量精度方面有了较大地提高，可以满足一般水位监测要求[1].但其监测时空分辨率较低，尤其是在库区高边坡等复杂地形以及恶劣天气条件下，其精度、可靠性及可用性难以保障.同时，利用全球卫星导航反射信号(GNSS-R)开展水面高度反演方面取得了一定的研究成果[2-7].一方面，GNSS 在覆盖范围、获取数据量、所需成本上具有其他监测方法无法比拟的优势[8].另一方面，国内如三峡等大型工程，在库岸高危边坡及坝体上建立了数量众多的GNSS为主的自动化变形监测系统.若充分利用现有的GNSS 测站，则能够在不增加投入成本的基础上，有效弥补现有手段时空分辨率低精度低的问题，极大地提高了GNSS 数据的利用率，对于大坝、湖泊等水位监测具有重要意义[9].

MARTIN-NEIRA 等[10]在1993 年首次提出利用GNSS-R 反演海平面的高度变化，开辟了利用GNSS-R 监测测站周围环境这个新兴的研究领域.2012 年，LARSON 等[11]利用瑞典Onsala 空间天文台的GPS 信噪比(SNR)观测数据，成功反演了海平面变化的时间序列，且与周围的验潮站资料进行对比分析，发现较差均方根(RMS)为4.8 cm，精度能达到10 cm.在国内，采用GNSS-R技术进行水位反演还处于起步阶段，2012 年，吴继忠等[12]给出一种利用单台测量型 GPS 接收机反射信号测量水面高度的方法，利用SNR 数据对平静湖面进行了测高试验，结果表明利用SNR 反演水面高度在较好的条件下精度能达到±3 cm.2018 年，匡翠林等[13]基于数据的分析软件，反演了国际GNSS 服务(IGS)站SC02 长达14 年的海平面变化时间序列，通过与该站附近的验潮站数据进行对比，发现两者具有很好的一致性，反演精度为8 cm.2019 年，SONG 等[14]使用GNSS-R技术来监测双王城水库水位的变化，通过GPS SNR 测量估算的水位变化的每日时间序列可以达到厘米级的精度.2020 年，张弛等[15]利用平静湖面上的 GNSS 静态观测数据进行水面高度的反演，确定了GNSS-R技术反演水面高度的能力，验证了普通大地测量型接收机应用于监测水面的有效性.程文[16]利用SNR 数据成功反演了Calcasieu Pass GPS 站的水位变化，与附近验潮站监测结果表现出良好的一致性，均方根误差(RMSE)在8 cm 左右，同时证明了岸基GPS 站可对水位连续监测.同年，南阳等[17]在西安灞河悬挑平台上使用3 组不同的仪器进行测高实验，指出目前L1 频段的SNR 更适合测高反演.

为了充分利用我国大中型水库、大坝为数众多的GNSS 变形监测系统提供的丰富的数据资源，发挥GNSS-R 反演水位在该领域的技术优势，本文详细推导了GNSS SNR 数据反演水面高度的原理，给出相应的数据处理方法，并选择济南卧虎山水库开展实验，通过在坝体及泄洪闸口不同位置设站观测，证明该方法的精度及可靠性.

1 GNSS-R 反演水面高度原理及数据处理方法

1.1 基于SNR 水面测高原理

由图1 可知，位于水系岸边的GNSS 设备所接收的是卫星的直射信号与水面反射信号叠加的混合信号.其中反射信号成分所占比重会在SNR 得到直观反映.因此，SNR 通常被用来评判信号质量，其值越高表示接收到的信号质量越好.直射信号和反射信号之间的干扰会在SNR 数据中产生特征变化，而该特征变化取决于接收机天线相位中心与反射介质表面之间的距离[9].

图1 GNSS-R 反演水面的原理

由图1 可知，接收机天线同时接收卫星直射及水面的反射信号，h是天线相位中心到瞬时水平面的垂直反射距离，θ 是瞬时卫星的直射信号和水平面的夹角，即卫星高度角.红色部分为直射信号与反射信号之间的路径差D，可表示为

其中，D1表示为

依据图1 所示几何关系，由D1可以推导出D2

由式(2)和式(3)可得，直射信号与反射信号之间的路径差D为

根据反射信号与直射信号之间的路径差可以得到相位差φ为

式中，φ是信号的载波波长 λ、卫星高度角 θ、垂直反射距离h的函数，并且与高度角的正弦值 sin θ 成正比.

根据接收机载波相位跟踪环SNR 和振幅之间的相关关系[18]

式中：AH是合成信号的振幅；Ad是直射信号振幅；Ar是反射信号振幅；cos θ 是直射信号和反射信号夹角的余弦值.

在合成信号中，直射信号占据主要成分[19]，且直射信号的振幅Ad远大于反射信号的振幅Ar，为了有效分析反射信号成分的影响，需采用二次多项式去除趋势项Ad得到 SNR 残差序列，用 δ(SNR) 来表示.

结合式(5)和式(6)，SNR 残差序列可用式(7)表示：

式中，A、φ是 δ(SNR) 中反射信号的振幅和相位.令x=sin θ，f=，式(7)可表示为

利用频谱分析获取f，进而可计算垂直距离h.考虑到卫星高度角正弦值变化是不均匀的，根据式(7)可知 δ(SNR) 相对于高度角正弦值是非等间隔采样的，因此需要采用Lomb-Scargle 周期算法代替传统的快速傅里叶变换(FFT)方法，以解决非等间隔采样问题，来获得残差序列频谱能量的主频[20].

对于残差序列 δSNR(tj)，j=1，2，3···，N，其功率谱可表达为频率f的函数

式中：Pδ(f)是频率f为的周期信号的功率；δSNR(tj)是残差序列信号；tj是残差序列信号所对应的时间；N为数据统计量；τ 为时间平移不变量，用来调节非等间隔采样的数据引起的相移[21-22].

故可以通过Lomb-Scargle 周期算法获取f，进而得到有效高度值[23]，即反射面到天线相位中心的垂直反射距离h

1.2 数据处理方法

基于以上SNR 反演水位的原理，在GNSS-R技术反演水面高度的数据处理过程中，精细化获取包含有水面反射信号的SNR 数据，需要对采集的原始数据进行预处理.在测站站心坐标系下，计算观测时段各卫星的高度角和方位角.根据接收机与水域相对位置，初步判定反射点落入该水域卫星的伪随机噪声码(PRN)及其对应时段.从观测文件中筛选出所需的SNR 数据后，进一步对其精细化提取.在卫星信号频率选择上，SNR 数值较小会导致反演精度较低，经对比，S1 均值大于S2 均值，因此选用L1 频率的SNR数据[9].在卫星高度角选择上，在相对平静的水库湖泊，选取5°～20°高度角区间卫星，SNR 数据因更易受到多路径效应的影响，周期性震荡更明显，有利于水面高度的反演.通过绘制高度角与SNR 关系曲线图，通常现形式有上升弧段、下降弧段以及包含完整的上升下降弧段.由于频谱分析时，高度角的正弦值应单调递增.倘若两个弧段都有记录，则需将SNR 数据划分为两段，如若只含有上升或者下降单个弧段，则直接处理.完成可用于反演的SNR 信号精细化提取后，进而用二次多项式去除SNR 趋势项，获得SNR 残差序列，对该残差序列进行Lomb-Scargle 算法分析获得主频率，最后根据式(10)实现水面高度的反演.数据处理技术流程如图2 所示.

图2 GNSS-R 数据处理技术流程图

2 GNSS-R 水位监测实验及结果分析

2.1 数据来源

为了验证上述GNSS-R技术监测水位变化理论方法的有效性以及探究GNSS-R 测高的可靠性和适用性，本文利用布设在水库岸边的两个测站的观测数据进行实验分析.实验地点为济南市卧虎山水库，如图3 所示，实验所用的接收机为Novatel OEM615 接收机，可以同时接收GPS 信号和GLONASS 信号.本文只使用GPS L1 波段的SNR 数据，接收机的采样间隔设置为1 s，实验观测时间为北京时间2019 年12 年1 日.实验期间同时实测一个水位数据，用于检验GNSS-R 的水位反演结果.

图3 测站位置、观测环境及仪器布设

因本实验对反射信号的要求最大程度上来自水面，在测站的地理位置选择上应避免其他反射信号的干扰.其中测站1 位于卧虎山水库大坝上，该位置附近的斜坡可以在一定程度上阻挡来自岸边的反射信号干扰；测站2 选取的位置是在水库闸门附近一个竖直的平台上，最大程度保证了接收机记录的反射信号全部来自水面的反射信号.图3 中(a)、(b)、(c)是测站1 和测站2 的具体位置和周围环境以及仪器布设.

2.2 实验分析

为了验证GNSS-R技术应用于水库水面的可靠性，将数据预处理软件对数据进行分析，基于该软件可得到观测时段落入反射区的卫星的情况，对测站1 落入反射区卫星G5 的L1 频率所有历元时间的有效SNR 数据可视化分析.经对比L1 频率的噪声小且精度高，因此采用GPS L1 波段的SNR 数据.由双纵轴图4 可知，卫星高度角呈单调下降趋势，因此不需再进行弧段划分.经过去除趋势项，得到SNR 残差序列，即PRN 5 号卫星在下降阶段去除趋势项后的SNR 残差序列.最后对得到的残差序列进行频谱分析，通过LSP 频谱分析得到主频率，然后根据式(10)得到测站1 对应的垂直反射距离h，具体如图5 所示.

图4 G5 卫星SNR 和高度角时序图

图5 测站1 垂直反射距离h

为了验证在不同测站GNSS-R 测高方法的适用性，在测站2 进行了同样的测高实验.实验数据选用测站2 落入反射区卫星G15 的L1 频率所有历元时间的有效SNR 数据.同时由双纵轴图6 可知，卫星高度角是下降的，所以只含有单个下降弧段，因此不需要再进行划分.图7 给出了PRN 15 号卫星在下降阶段去除趋势项、重采样后的SNR 残差序列，截至高度角5°～20°，横坐标为随高度角变化的重采样时间点，纵坐标为 SNR 值的线性变化值伏特.最后对残差序列进行频谱分析，根据式(10)得到对应的垂直反射距离h，如图8 所示.

图6 G15 卫星L1 SNR 和高度角时序图

图7 去除趋势项后的SNR 残差序列图

图8 测站2 垂直反射距离h

2.3 结果综合分析

由表1 可知，将通过GNSS-R技术反演得到天线相位中心到水平面的距离与实测距离对比分析可得，垂直反射距离h的反演精度可以达到厘米级，测站1 较差为3 cm，测站2 较差为5 cm.同时由实验结果可知，按照最大程度使得反射信号来自于水面的需求选择的测站1 和测站2 均可获得有效水面高度.即无论是在卧虎山水库大坝上的测站1，还是在水库闸门附近竖直平台上的测站2，均能取得较好的水位反演结果.

表1 反演水面高与实测高度对比表 m

3 结束语

为了充分发挥GNSS-R 在水位监测上的优势，有效利用大坝已建成的GNSS 变形监测系统测站的数据，本文在详细推导GNSS-R技术监测水面高度的原理基础上，以SNR 为主要观测量，分析精细化提取用于反演卫星及其有用时段方法，进而进行去除趋势项及LSP 频谱分析，并详细介绍了数据处理方法与技术流程.通过在济南卧虎山水库坝体及泄洪闸不同位置进行实测数据反演、分析，并与现场实测水位高度记录数据进行对比，结果表明通过精细化选择卫星及其有效观测时段SNR，按照本文理论方法均得到有效反演高度，通过与实测水面高度比较，两测站较差分别为3 cm、5 cm.反演精度达到厘米级，满足一般水位监测的需要.充分验证了利用现有大坝和库区已建立的GNSS 变形监测系统，进行GNSS-R 水位监测可行性.随着我国北斗三号(BDS-3)投入使用，可用导航定位卫星增加，大大提高水位监测时空分辨率，使得水位监测的智慧化运营管理更加灵活、准确、多样.