小波与BP神经网络联合反演GNSS-IR高精度水库水位变化

杨晓峰，魏浩翰，张 强，刘朝海

(1.南京林业大学土木工程学院，南京 210037； 2.江苏莱特北斗信息科技有限公司，江苏 常州 213100)

0 引言

长期、高精度的水位监测可为防洪抗旱和水利建设提供重要的参考资料，且对于农业生产及水资源的合理运用具有重要意义。通常情况下，水库会建立水位监测站或设置标尺进行水位数据的读取，但是这些方法浪费人力物力、效率低、出错率高。因此，使用一种高效且精确的水位监测方法代替传统的水位数据收集方法十分重要。

近年来，基于单天线的GNSS干涉反射测量(GNSS Interferometric Reflectometry，GNSS-IR)技术持续受到关注，可用于近地表环境参数(包括土壤湿度、积雪深度、植被、潮汐、水位等)反演[1-5]。对于已建立全球卫星导航系统(Global Navigation Satellite System，GNSS)变形监测网的大坝库区，GNSS-IR技术无需增加额外设备，仅利用现有GNSS接收设备和大量免费的GNSS信号得到的信噪比(Signal-to-Noise Ratio，SNR)数据即可反演得到水库水位，利用这种新的水库水位反演方法可代替传统方法实现水位读取的准确化、自动化，并能提高GNSS信号的利用率，节约水库水位监测成本，对于水库范围内的灾害监测与预警具有独特的优势和广泛的应用前景[6]。

国内外学者针对GNSS-IR反演水位开展了大量研究。苏晓容、WANG等在海潮反演过程中进行小波去噪，消除复杂环境下SNR信号中的噪声，使得GNSS-IR技术反演潮位的精度提高到分米级[7-8]；吕铮等利用GPS L1、L2和BDS B1、B2信噪比数据，在水库水位反演过程中，分析信号频率、高度角范围及弧段长度等因素对大坝水位反演结果的影响，并认为GPS L1和BDS B1信号更适用于大坝水位反演，在此基础上利用中值法进行GNSS信号多频融合，得到了0.1～0.13m的水库水位反演精度[6]；王杰等引入潮波函数对复杂环境下GNSS-IR反演的海潮潮位进行改正，使反演精度提高至0.13m[9]；王瑞芳利用经验模态分解方法提取海港水面反射信号，使GNSS-IR反演结果的均方根误差(Root Mean Square Error，RMSE)提高至0.15m[10]；张一等利用NARX回归神经网络反演海面有效波高，RMSE为0.505m[11]；陈昊晟等基于GNSS-IR技术利用均值法融合GPS L1、GLONASS L1 L2和BDS B1频段的SNR信号，获得了0.162m的水位反演精度[12]。以上研究大多数是针对海洋潮汐的反演，针对内陆水体水位的反演研究目前主要集中在信号频率、高度角范围及信号弧段长度对反演结果的影响[6,13]，反演精度仅达到分米级。为了提高GNSS-IR反演水库水位的精度，本文基于多频多模GNSS信号融合反演的策略，选取位于南水北调山东省境内双王城水库GNSS变形观测站2017年10月1日—12月26日共87天的原始SNR数据为研究对象，首先通过设置高度角和方位角范围选取有效的SNR观测量；其次分别利用最小二乘拟合法和小波分解法对原始SNR序列进行去趋势项，提取去趋势项后的SNR残差序列；然后对全球定位系统(Global Positioning System，GPS)和北斗卫星导航系统(BeiDou Navigation Satellite System，BDS)各频段SNR残差序列进行Lomb-Scargle频谱分析(LS谱分析)，得到对应的水位反演结果；最后分别利用均值算法、中值算法、随机森林算法和BP神经网络算法对各频段SNR水位反演结果进行数据融合，并与实测水位数据进行对比分析。

1 基本原理

1.1 GNSS-IR水库水位高度反演基本原理

卫星发射的微波信号在传播过程中不免会产生多路径效应，其直射信号与经由平静水面以及地表反射的反射信号在干涉效应下形成的复合信号被接收机所接收，经处理生成复合信号数据，即SNR数据。基于GNSS-IR技术可以充分利用SNR数据来反演水库水位高度，反演的基本原理如图1所示[14]。

图1 GNSS-IR技术反演水库水位高度的原理图

SNR主要由以下两部分组成[15]

SNR=SNRd+dSNR

(1)

式中，SNRd为趋势项信号；dSNR为直射信号和反射信号经由干涉效应所形成的残差序列。

从信号振荡的角度考虑，SNR序列也可以通过振荡幅度和相位差来表示[5]

SNR2≈Ac2=Ad2+Ar2+2AdArcosψ

(2)

式中，Ad为直射信号分量振荡幅度；Ar为多路径反射信号分量振荡幅度；ψ为两者的相位差。在卫星运动过程中，直射信号和反射信号的相位差ψ随着卫星高度角e的变化而变化，体现在受信号干涉的SNR值的增强和减弱。由于多路径和天线增益，低高度角下SNR值的振荡更加明显。

由于测量型GNSS接收机天线有效抑制了地表反射的多路径信号，直射信号的振幅远大于反射信号的振幅，即Ad>>Ar。为了得到包含反射信息的SNR残差序列，通常采用低阶多项式拟合[14]的方法提取出SNR序列的趋势项。低高度角下的SNR残差序列可表示为[5]

(3)

dSNR=Acos(2πft+φ)

(4)

对SNR残差序列进行LS谱分析[16]，即可得到振荡频率f，则等效天线高可表示为[5]

(5)

1.2 小波分解原理

在式(1)中，SNR包含直射信号分量和反射信号分量，需要从SNR序列中消除趋势项以得到dSNR残差序列，进而通过LS谱分析得到等效天线高。常用的消除趋势项的方法为二阶多项式法[14]，本文采用小波分解[17-18]的方法，其原理如下所示。

设SNR观测值G(t)为

(6)

式中，t为历元。小波分解的算法为[17]

(7)

(8)

式中，Aj为低频信号的小波系数，j为分解层数；Dj为高频信号的小波系数；G(t)为原始信号，t为历元。以db5作为小波基对原始SNR序列进行5层分解处理，将原始SNR信号减去第五层低频信号，即可消除趋势项获得包含水面信息的反射信号，并与传统的二阶多项式法得到的结果进行对比。

1.3 多星多频数据融合方法

理论上，根据式(1)～式(5)，利用一颗卫星一个频率的SNR观测量即可实现水库水位反演。但是，受限于测站周边环境、卫星信号质量、接收机软硬件差异等多种因素影响，不同卫星、不同频率反演结果互有差异，因此要探索多频多模数据融合的最优方法。本文分别采用均值算法、中值算法、随机森林算法[19]和BP神经网络算法[20]对GPS和BDS各频段的数据反演结果进行融合，进而探寻最优的融合算法。四种算法的原理如下：

均值算法和中值算法是常用的数据融合算法，其主要区别在于：均值算法是对观测数据取平均值，而中值算法是直接选取数据的中间值作为输出结果。

随机森林[19]是目前机器学习算法中一种常用的算法，将多个决策树通过集成的方法融合在一起，最终结果由每个决策树的结果综合得到。每一个决策树的训练样本是通过自助采样法抽取的，即随机从所有特征中抽取一个子集，其随机性抽取可避免过拟合的问题，且训练抽取的随机性增加了各个决策树之间的差异，使得最终融合的模型具有较高的精度[19]。

BP神经网络[20]是一种由误差反馈训练的多层感知器网络，主要由输入层、隐藏层和输出层组成。在一个神经网络中可以有多个隐藏层，同一层的神经元相互独立，相邻层之间的神经元完全相连。这种网络由一系列节点组成，包括正向传播和误差反向传播两个过程。其最核心的特点就是：信号是前向传播，而误差是反向传播。前向传播过程中，输入信号经由输入层、隐藏层逐层处理，到输出层时，如果结果未达到期望要求，则进入反向传播过程，将误差信号原路返回，修改各层权重。其基本思想是通过迭代调整权重，使网络的实际输出值和期望输出值之间的均方根误差最小化。因此具有良好的自学习、自适应、鲁棒性和泛化能力[20]。

2 试验配置及流程

2.1 数据选取

GNSS观测数据由位于南水北调东线山东境内的双王城水库GNSS变形监测网SW50测站提供，数据类型是GPS和BDS双系统SNR原始数据；实测水位数据由现场水位观测站提供。如图2所示，SW50测站位于库区西北角，以观测墩的方式安置在大坝内侧，附近有水位观测站提供每日实测水位数据。从图2可以看出，当信号方位角在200°～300°时，反射表面为测站西侧的地表；当方位角在30°～120°时，反射表面从远到近分别为水库水体和大坝护坡。由此可见，方位角应选择30°～120°，且高度角应选择水面反射区域，即5°～15°，尽量避开大坝护坡的表面反射信号。

(a)SW50测站实景照片

2.2 水位反演策略

收集了SW50测站2017 年10月1日—12月26日(年积日DOY 274～360)共87天观测数据和同期水位监测站的实测水位数据，具体的水位反演策略为：

1)SNR数据筛选。根据测站位置(图2(a))和菲涅尔反射区(图2(b))确定水面反射区域，其中方位角范围为30°～120°，高度角范围为5°～15°，选取运动轨迹在此范围内的可用卫星，提取其有效的SNR观测量。

2)SNR去趋势项。分别使用二阶多项式拟合法及小波分解法对SNR序列进行去趋势项，得到包含水面高度信息的SNR残差序列。

3)LS谱分析。对SNR残差序列进行LS谱分析，获取不同频段SNR信号反演的等效天线高，并与对应频段DOY274反演的等效天线高求差，从而得到DOY274以来反演的水位变化值；在此基础上，与实测水位变化值进行对比分析，选出较好的去趋势项方法。

4)多频多模SNR信号融合。分别利用均值算法、中值算法、随机森林算法和BP神经网络算法对GPS和BDS多频多模信号反演结果进行融合，并分析其精度。

3 结果与分析

3.1 SNR信号去趋势项方法对比分析

为了比较二阶多项式法和小波分解法去趋势项的效果，分别用这两种方法对原始SNR观测序列进行去趋势项处理。以2017年10月1日(DOY 274)G26卫星观测到的S1C频段SNR数据为例，图3所示为利用二阶多项式拟合去趋势项和频谱分析的序列图。其中，图3(a)为去趋势项前后的SNR序列图，黑色线为原始SNR序列，红色线为SNR的去趋势项序列，蓝色线为消除趋势项信号后的SNR残差序列；图3(b)为SNR残差序列LS谱分析曲线图。从图3可以直观地看出，对于S1C频段的SNR数据，传统的二阶多项式法去趋势效果比较理想。

图4所示为利用db5作为小波基对上述SNR序列进行5层小波分解处理后的结果。其中，黑色线为原始SNR序列；深蓝色线为小波分解的5层高频信号，即细节信号，反映原始SNR序列的细节信息；浅蓝色线和红色线为小波分解下的5层低频信号，也称作逼近信号，是原始SNR序列缓慢变化的部分，红色线为趋势项。将趋势项与原始SNR序列相减可提取SNR残差序列，结果如图5(a)所示，对所得到的dSNR残差序列进行LS谱分析即可求得等效天线高，图5(b)所示为G26卫星S1C频段的SNR残差序列LS谱分析曲线图。结合图3和图5可知，对于S1C频段的SNR数据，二阶多项式法和小波分解法去趋势项效果均比较理想。

(a)

图4 小波分解的原始SNR序列高频和低频信号(2017, DOY274，G26卫星S1C频段信号)

(a) 消除趋势项后的SNR残差序列

为了进一步分析不同去趋势项方法对各SNR频段的适用性，图6所示为2017年DOY 274～360期间，分别利用二阶多项式和小波分解对GPS的S1C、S2L、S2W和S5Q这4个频段SNR观测值去趋势项后采用LS谱分析反演得到的水位曲线，同时与实测水位曲线进行对比。总体上看，GPS的S1C、S2L和S5Q这3个频段SNR的水位反演结果与实测水位有较好的一致性。其中，S1C和S5Q频段，如图6(a)、(d)所示，利用二阶多项式拟合和小波分解去趋势项后反演的水位曲线高度一致甚至几乎重合。此外，S2L频段小波分解的效果明显优于二阶多项式拟合的效果，这是由于试验所用的GPS卫星还未完全实现现代化，所用S2L频段是混合了L2C码与L2P码的混合信号，其中的L2P码信号大大降低了反演精度[13]，甚至有些研究认为S2L信号不适用于水位反演[6-7]。另外，对于S2W频段，如图6(c)所示，不管是用二阶多项式拟合还是小波分解的方法，其水位反演结果均不够理想，这是由于S2W信号采用Z跟踪(Z-Tracking)技术导致噪声增加[21-22]。综上可知，对于SNR观测值去趋势项方法的选择，小波分解法整体上优于二阶多项式拟合法，但是GPS-S2W频段用于水位反演的效果不甚理想。

(a) S1C

需要特别注意的是，图6(c)的S2W频段反演水位结果显示，在DOY310之前小波分解结果与实测值符合较好，DOY310～315之间反演结果出现了突变，在DOY315之后又与实测结果有良好的一致性。究其原因，是由于GPS S2W频段信号采用的Z跟踪技术会导致L2信号被L1信号污染，出现“双峰”现象[21]——S2W频段第一个峰值对应L2信噪比的多径振荡，第二个峰值为L1污染L2产生的噪声信号。为了分析GPS S2W“双峰”现象对水位反演造成的干扰，图7绘制了S2W信号的LS谱分析图和水位反演结果。其中，图7(a)为2017年DOY 314的S2W频段的 LS谱分析图，可以看出，“双峰”的振幅非常接近，难以进行有效区分。基于此，图6(c)在DOY310～315之间才会出现反演结果突变现象。图7(b)为分别利用“双峰”反演的水位时间序列图，其中第一个主频峰值的水位反演结果与实测水位有较好的一致性，其均方根误差、平均偏差和相关系数分别为0.1439m、0.1094m和0.9987，但是反演的精度仍不如GPS其他频段；第二个主频峰值的水位反演结果与实测水位有较大的误差，其均方根误差、平均偏差和相关系数分别为0.6792m、0.5578m和0.9420。

(a)

为了进一步验证上述结果，表1列出了GPS各频段信号反演的水位相对于实测水位的精度指标。

从表1可以明显看出，与实测水位相比，无论采用二阶多项式法还是小波分解方法，S2W频段的反演结果的均方根误差达到甚至超过了1m，再次验证了该频段SNR观测量用于水位反演效果不佳的结论。此外，小波分解和二阶多项式方法在S1C频段的反演精度相当，而在S2L、S5Q频段，小波分解方法明显优于二阶多项式法。限于篇幅，BDS各频段SNR去趋势项后的水位反演结果仅用表格的形式列于表2中。从表2可以看出，对于S1I、S6I和S7I这3个频段，小波分解方法的反演精度略优于二阶多项式拟合法。综合表1和表2的数据，剔除S2W频段第二主频峰值反演结果后,对所有频段反演结果取算数平均值，发现小波分解反演的水位与实测水位对比得到的平均均方根误差、平均偏差和平均相关系数分别为0.1062m、0.0304m和0.9986，而二阶多项式拟合法的精度结果分别为0.2245m、0.1147m和0.9978。

表1 各频段GPS-SNR小波分解法和二阶多项式法去趋势项后反演水位精度对比

表2 各频段BDS-SNR小波分解法和二阶多项式法去趋势项后反演水位精度对比

从理论上来说，传统的二阶多项式拟合法虽然可用于处理大多数卫星频段的SNR数据，但其基本原理是基于最小二乘思想，即基于残差的平方和寻找数据的最优函数匹配，等权地处理所有数据。当原始SNR序列中出现一部分质量不佳的信号时，残差平方和会对拟合结果产生比较明显的影响，使拟合曲线偏移，导致拟合出的趋势项与实际直射信号相比出现较大的偏差[23]。而小波分解是基于信号学理论，通过频率变化的快慢将非平稳的SNR信号分离成低频信号和高频信号，从而实现去趋势项。即使原始SNR序列中出现一部分质量不佳的信号，也不会影响高频低频信号的分离效果。因此，理论分析与试验验证均说明，小波分解法用于原始SNR信号去趋势项的效果要优于二阶多项式拟合法。

综上可知，在水库水位高度反演消除SNR趋势项的过程中，小波分解法总体上优于传统的二阶多项式拟合法，且S2W频段信号用于水位反演效果不佳。同时也可以发现，不同系统不同频段数据之间存在不同的系统偏差，使得在同一水位基准下不同系统不同频段数据的反演结果随着DOY的增加出现偏移现象，需要根据误差调整权重，通过多星多模融合以进一步减少系统误差的影响。

3.2 多频多模信号融合反演水位变化分析

虽然单一频段SNR信号能够用于水库水位反演，但是可用卫星信号数量有限，难以发现可能存在的粗差，且不同频段之间可能存在系统性偏差，因此要考虑进行多星多频信号融合反演[24]。在利用小波分解对SNR去趋势项和LS谱分析获得其反演结果的基础上，分别利用均值法、中值法、随机森林算法和BP神经网络算法这四种信号融合算法对测站接收到的全部GPS、BDS频段信号反演结果进行融合，并与实测水位进行比较，结果如表3所示。

由表3可知，将所有频段的反演水位进行多星多频融合后，均方根误差只达到了分米级，甚至不如某些单频信号(S1C、S2L、S5Q、S7I)的反演结果，总体来说反演效果不太理想。这是由于融合过程中反演精度较低的频段作为输入项会影响整体的反演精度。因此，在信号融合前需要设置阈值将反演精度较差的频段进行剔除。在确保数据利用率的前提下，分别设置0.1m、0.15m和0.2m为阈值，将均方根误差小于阈值的频段设置为可取频段，大于阈值的频段设置为不可取频段予以去除，经过比较发现，选择0.1m的阈值最为合适。由表1可知，SW50测站GPS的S1C、S2L和S5Q均为可取频段；由表2可知，BDS的S1I和S6I频段的均方根误差均超过阈值0.1m，均为不可取频段，只有S7I频段的反演结果达到厘米级，为可取频段。水位反演所用信号频段和对应卫星编号如表4所示。

表4 水位反演所用信号频段和对应卫星编号

将前50天各频段单独反演结果作为BP神经网络和随机森林算法的训练样本进行训练，输出单元为后37天的水位变化结果。为了避免网络收敛过慢，对输入数据进行归一化处理。其中BP神经网络采用37-1-1型网络结构，采用L-M算法，各参数设置为[20]：性能函数mse，输入传递函数tansig，输出传递函数purelin，目标误差0.0003，学习速率0.01，最大训练次数1000。

将均值法、中值法、随机森林和BP神经网络四种融合算法得到的水位变化结果与实测水位变化进行对比分析，结果如表5所示。

表5 不同方法反演水位的精度

由表5可知，对于水面较平静的环境，四种算法的均方根误差和平均偏差均达到了厘米级精度，且相关系数高于0.999，说明均具有良好的反演效果。其中，从均方根误差和平均偏差2个指标来衡量，BP神经网络算法的水位反演效果明显优于其他三种算法。从理论上来看，均值法和中值法对各频段信号的反演结果分配的权重相同，因此对反演效果不佳的频段不能有效抑制其反演误差；随机森林算法中虽然每个决策树的训练样本是随机抽取的，但是每个树的权重相等，其本质也是一种特殊的等权算法；而BP神经网络可以通过每次训练的误差来修正权重，使网络的实际输出值(反演水位结果)和期望输出值(实测水位结果)之间的均方根误差最小化。因此，理论分析和实测验证均说明BP神经网络算法更优。

4 结论

以双王城水库SW50测站2017 年10月1日—12月26日的GNSS SNR原始观测量为研究对象，全面分析了GPS和BDS双系统多频SNR信号的变化特征，在此基础上提出了联合小波分解与BP神经网络法反演水库水位变化的策略，并得到以下结果：

1)在SNR去趋势项方法的选择上，小波分解算法优于传统的二阶多项式拟合算法。

2)利用多频多模SNR信号反演水库水位变化，需要提前设置阈值以消除粗差、提高精度，并且在水库、大坝等水面较为平静的反演环境下，均值算法、中值算法、随机森林算法和BP神经网络算法均能达到厘米级精度，其中BP神经网络算法用于水位反演的效果更优。

同时，本研究首次对GPS S2W信号双峰现象用于水位反演进行探讨，发现由于采用Z跟踪技术，GPS S2W信号所呈现的“双峰”现象对水位变化反演产生了明显干扰，在后续工作中需要对此现象进一步研究以准确提取有用的信号。此外，由于所选测站硬件设备的限制，BDS数据不包含北斗3号的卫星信号，未来需对更多的北斗卫星反射信号反演厘米级水位的方法进行探讨。