小波與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號的衛星信號，未來需對更多的北斗衛星反射信號反演厘米級水位的方法進行探討。