楊曉峰，魏浩翰，張 強(qiáng)，劉朝海

(1.南京林業(yè)大學(xué)土木工程學(xué)院，南京 210037； 2.江蘇萊特北斗信息科技有限公司，江蘇 常州 213100)

0 引言

長(zhǎng)期、高精度的水位監(jiān)測(cè)可為防洪抗旱和水利建設(shè)提供重要的參考資料，且對(duì)于農(nóng)業(yè)生產(chǎn)及水資源的合理運(yùn)用具有重要意義。通常情況下，水庫(kù)會(huì)建立水位監(jiān)測(cè)站或設(shè)置標(biāo)尺進(jìn)行水位數(shù)據(jù)的讀取，但是這些方法浪費(fèi)人力物力、效率低、出錯(cuò)率高。因此，使用一種高效且精確的水位監(jiān)測(cè)方法代替?zhèn)鹘y(tǒng)的水位數(shù)據(jù)收集方法十分重要。

近年來，基于單天線的GNSS干涉反射測(cè)量(GNSS Interferometric Reflectometry，GNSS-IR)技術(shù)持續(xù)受到關(guān)注，可用于近地表環(huán)境參數(shù)(包括土壤濕度、積雪深度、植被、潮汐、水位等)反演[1-5]。對(duì)于已建立全球衛(wèi)星導(dǎo)航系統(tǒng)(Global Navigation Satellite System，GNSS)變形監(jiān)測(cè)網(wǎng)的大壩庫(kù)區(qū)，GNSS-IR技術(shù)無需增加額外設(shè)備，僅利用現(xiàn)有GNSS接收設(shè)備和大量免費(fèi)的GNSS信號(hào)得到的信噪比(Signal-to-Noise Ratio，SNR)數(shù)據(jù)即可反演得到水庫(kù)水位，利用這種新的水庫(kù)水位反演方法可代替?zhèn)鹘y(tǒng)方法實(shí)現(xiàn)水位讀取的準(zhǔn)確化、自動(dòng)化，并能提高GNSS信號(hào)的利用率，節(jié)約水庫(kù)水位監(jiān)測(cè)成本，對(duì)于水庫(kù)范圍內(nèi)的災(zāi)害監(jiān)測(cè)與預(yù)警具有獨(dú)特的優(yōu)勢(shì)和廣泛的應(yīng)用前景[6]。

國(guó)內(nèi)外學(xué)者針對(duì)GNSS-IR反演水位開展了大量研究。蘇曉容、WANG等在海潮反演過程中進(jìn)行小波去噪，消除復(fù)雜環(huán)境下SNR信號(hào)中的噪聲，使得GNSS-IR技術(shù)反演潮位的精度提高到分米級(jí)[7-8]；呂錚等利用GPS L1、L2和BDS B1、B2信噪比數(shù)據(jù)，在水庫(kù)水位反演過程中，分析信號(hào)頻率、高度角范圍及弧段長(zhǎng)度等因素對(duì)大壩水位反演結(jié)果的影響，并認(rèn)為GPS L1和BDS B1信號(hào)更適用于大壩水位反演，在此基礎(chǔ)上利用中值法進(jìn)行GNSS信號(hào)多頻融合，得到了0.1～0.13m的水庫(kù)水位反演精度[6]；王杰等引入潮波函數(shù)對(duì)復(fù)雜環(huán)境下GNSS-IR反演的海潮潮位進(jìn)行改正，使反演精度提高至0.13m[9]；王瑞芳利用經(jīng)驗(yàn)?zāi)B(tài)分解方法提取海港水面反射信號(hào)，使GNSS-IR反演結(jié)果的均方根誤差(Root Mean Square Error，RMSE)提高至0.15m[10]；張一等利用NARX回歸神經(jīng)網(wǎng)絡(luò)反演海面有效波高，RMSE為0.505m[11]；陳昊晟等基于GNSS-IR技術(shù)利用均值法融合GPS L1、GLONASS L1 L2和BDS B1頻段的SNR信號(hào)，獲得了0.162m的水位反演精度[12]。以上研究大多數(shù)是針對(duì)海洋潮汐的反演，針對(duì)內(nèi)陸水體水位的反演研究目前主要集中在信號(hào)頻率、高度角范圍及信號(hào)弧段長(zhǎng)度對(duì)反演結(jié)果的影響[6,13]，反演精度僅達(dá)到分米級(jí)。為了提高GNSS-IR反演水庫(kù)水位的精度，本文基于多頻多模GNSS信號(hào)融合反演的策略，選取位于南水北調(diào)山東省境內(nèi)雙王城水庫(kù)GNSS變形觀測(cè)站2017年10月1日—12月26日共87天的原始SNR數(shù)據(jù)為研究對(duì)象，首先通過設(shè)置高度角和方位角范圍選取有效的SNR觀測(cè)量；其次分別利用最小二乘擬合法和小波分解法對(duì)原始SNR序列進(jìn)行去趨勢(shì)項(xiàng)，提取去趨勢(shì)項(xiàng)后的SNR殘差序列；然后對(duì)全球定位系統(tǒng)(Global Positioning System，GPS)和北斗衛(wèi)星導(dǎo)航系統(tǒng)(BeiDou Navigation Satellite System，BDS)各頻段SNR殘差序列進(jìn)行Lomb-Scargle頻譜分析(LS譜分析)，得到對(duì)應(yīng)的水位反演結(jié)果；最后分別利用均值算法、中值算法、隨機(jī)森林算法和BP神經(jīng)網(wǎng)絡(luò)算法對(duì)各頻段SNR水位反演結(jié)果進(jìn)行數(shù)據(jù)融合，并與實(shí)測(cè)水位數(shù)據(jù)進(jìn)行對(duì)比分析。

1 基本原理

1.1 GNSS-IR水庫(kù)水位高度反演基本原理

衛(wèi)星發(fā)射的微波信號(hào)在傳播過程中不免會(huì)產(chǎn)生多路徑效應(yīng)，其直射信號(hào)與經(jīng)由平靜水面以及地表反射的反射信號(hào)在干涉效應(yīng)下形成的復(fù)合信號(hào)被接收機(jī)所接收，經(jīng)處理生成復(fù)合信號(hào)數(shù)據(jù)，即SNR數(shù)據(jù)?；贕NSS-IR技術(shù)可以充分利用SNR數(shù)據(jù)來反演水庫(kù)水位高度，反演的基本原理如圖1所示[14]。

圖1 GNSS-IR技術(shù)反演水庫(kù)水位高度的原理圖

SNR主要由以下兩部分組成[15]

SNR=SNRd+dSNR

(1)

式中，SNRd為趨勢(shì)項(xiàng)信號(hào)；dSNR為直射信號(hào)和反射信號(hào)經(jīng)由干涉效應(yīng)所形成的殘差序列。

從信號(hào)振蕩的角度考慮，SNR序列也可以通過振蕩幅度和相位差來表示[5]

SNR2≈Ac2=Ad2+Ar2+2AdArcosψ

(2)

式中，Ad為直射信號(hào)分量振蕩幅度；Ar為多路徑反射信號(hào)分量振蕩幅度；ψ為兩者的相位差。在衛(wèi)星運(yùn)動(dòng)過程中，直射信號(hào)和反射信號(hào)的相位差ψ隨著衛(wèi)星高度角e的變化而變化，體現(xiàn)在受信號(hào)干涉的SNR值的增強(qiáng)和減弱。由于多路徑和天線增益，低高度角下SNR值的振蕩更加明顯。

由于測(cè)量型GNSS接收機(jī)天線有效抑制了地表反射的多路徑信號(hào)，直射信號(hào)的振幅遠(yuǎn)大于反射信號(hào)的振幅，即Ad>>Ar。為了得到包含反射信息的SNR殘差序列，通常采用低階多項(xiàng)式擬合[14]的方法提取出SNR序列的趨勢(shì)項(xiàng)。低高度角下的SNR殘差序列可表示為[5]

(3)

dSNR=Acos(2πft+φ)

(4)

對(duì)SNR殘差序列進(jìn)行LS譜分析[16]，即可得到振蕩頻率f，則等效天線高可表示為[5]

(5)

1.2 小波分解原理

在式(1)中，SNR包含直射信號(hào)分量和反射信號(hào)分量，需要從SNR序列中消除趨勢(shì)項(xiàng)以得到dSNR殘差序列，進(jìn)而通過LS譜分析得到等效天線高。常用的消除趨勢(shì)項(xiàng)的方法為二階多項(xiàng)式法[14]，本文采用小波分解[17-18]的方法，其原理如下所示。

設(shè)SNR觀測(cè)值G(t)為

(6)

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

(7)

(8)

式中，Aj為低頻信號(hào)的小波系數(shù)，j為分解層數(shù)；Dj為高頻信號(hào)的小波系數(shù)；G(t)為原始信號(hào)，t為歷元。以db5作為小波基對(duì)原始SNR序列進(jìn)行5層分解處理，將原始SNR信號(hào)減去第五層低頻信號(hào)，即可消除趨勢(shì)項(xiàng)獲得包含水面信息的反射信號(hào)，并與傳統(tǒng)的二階多項(xiàng)式法得到的結(jié)果進(jìn)行對(duì)比。

1.3 多星多頻數(shù)據(jù)融合方法

理論上，根據(jù)式(1)～式(5)，利用一顆衛(wèi)星一個(gè)頻率的SNR觀測(cè)量即可實(shí)現(xiàn)水庫(kù)水位反演。但是，受限于測(cè)站周邊環(huán)境、衛(wèi)星信號(hào)質(zhì)量、接收機(jī)軟硬件差異等多種因素影響，不同衛(wèi)星、不同頻率反演結(jié)果互有差異，因此要探索多頻多模數(shù)據(jù)融合的最優(yōu)方法。本文分別采用均值算法、中值算法、隨機(jī)森林算法[19]和BP神經(jīng)網(wǎng)絡(luò)算法[20]對(duì)GPS和BDS各頻段的數(shù)據(jù)反演結(jié)果進(jìn)行融合，進(jìn)而探尋最優(yōu)的融合算法。四種算法的原理如下：

均值算法和中值算法是常用的數(shù)據(jù)融合算法，其主要區(qū)別在于：均值算法是對(duì)觀測(cè)數(shù)據(jù)取平均值，而中值算法是直接選取數(shù)據(jù)的中間值作為輸出結(jié)果。

隨機(jī)森林[19]是目前機(jī)器學(xué)習(xí)算法中一種常用的算法，將多個(gè)決策樹通過集成的方法融合在一起，最終結(jié)果由每個(gè)決策樹的結(jié)果綜合得到。每一個(gè)決策樹的訓(xùn)練樣本是通過自助采樣法抽取的，即隨機(jī)從所有特征中抽取一個(gè)子集，其隨機(jī)性抽取可避免過擬合的問題，且訓(xùn)練抽取的隨機(jī)性增加了各個(gè)決策樹之間的差異，使得最終融合的模型具有較高的精度[19]。

BP神經(jīng)網(wǎng)絡(luò)[20]是一種由誤差反饋訓(xùn)練的多層感知器網(wǎng)絡(luò)，主要由輸入層、隱藏層和輸出層組成。在一個(gè)神經(jīng)網(wǎng)絡(luò)中可以有多個(gè)隱藏層，同一層的神經(jīng)元相互獨(dú)立，相鄰層之間的神經(jīng)元完全相連。這種網(wǎng)絡(luò)由一系列節(jié)點(diǎn)組成，包括正向傳播和誤差反向傳播兩個(gè)過程。其最核心的特點(diǎn)就是：信號(hào)是前向傳播，而誤差是反向傳播。前向傳播過程中，輸入信號(hào)經(jīng)由輸入層、隱藏層逐層處理，到輸出層時(shí)，如果結(jié)果未達(dá)到期望要求，則進(jìn)入反向傳播過程，將誤差信號(hào)原路返回，修改各層權(quán)重。其基本思想是通過迭代調(diào)整權(quán)重，使網(wǎng)絡(luò)的實(shí)際輸出值和期望輸出值之間的均方根誤差最小化。因此具有良好的自學(xué)習(xí)、自適應(yīng)、魯棒性和泛化能力[20]。

2 試驗(yàn)配置及流程

2.1 數(shù)據(jù)選取

GNSS觀測(cè)數(shù)據(jù)由位于南水北調(diào)東線山東境內(nèi)的雙王城水庫(kù)GNSS變形監(jiān)測(cè)網(wǎng)SW50測(cè)站提供，數(shù)據(jù)類型是GPS和BDS雙系統(tǒng)SNR原始數(shù)據(jù)；實(shí)測(cè)水位數(shù)據(jù)由現(xiàn)場(chǎng)水位觀測(cè)站提供。如圖2所示，SW50測(cè)站位于庫(kù)區(qū)西北角，以觀測(cè)墩的方式安置在大壩內(nèi)側(cè)，附近有水位觀測(cè)站提供每日實(shí)測(cè)水位數(shù)據(jù)。從圖2可以看出，當(dāng)信號(hào)方位角在200°～300°時(shí)，反射表面為測(cè)站西側(cè)的地表；當(dāng)方位角在30°～120°時(shí)，反射表面從遠(yuǎn)到近分別為水庫(kù)水體和大壩護(hù)坡。由此可見，方位角應(yīng)選擇30°～120°，且高度角應(yīng)選擇水面反射區(qū)域，即5°～15°，盡量避開大壩護(hù)坡的表面反射信號(hào)。

(a)SW50測(cè)站實(shí)景照片

2.2 水位反演策略

收集了SW50測(cè)站2017 年10月1日—12月26日(年積日DOY 274～360)共87天觀測(cè)數(shù)據(jù)和同期水位監(jiān)測(cè)站的實(shí)測(cè)水位數(shù)據(jù)，具體的水位反演策略為：

1)SNR數(shù)據(jù)篩選。根據(jù)測(cè)站位置(圖2(a))和菲涅爾反射區(qū)(圖2(b))確定水面反射區(qū)域，其中方位角范圍為30°～120°，高度角范圍為5°～15°，選取運(yùn)動(dòng)軌跡在此范圍內(nèi)的可用衛(wèi)星，提取其有效的SNR觀測(cè)量。

2)SNR去趨勢(shì)項(xiàng)。分別使用二階多項(xiàng)式擬合法及小波分解法對(duì)SNR序列進(jìn)行去趨勢(shì)項(xiàng)，得到包含水面高度信息的SNR殘差序列。

3)LS譜分析。對(duì)SNR殘差序列進(jìn)行LS譜分析，獲取不同頻段SNR信號(hào)反演的等效天線高，并與對(duì)應(yīng)頻段DOY274反演的等效天線高求差，從而得到DOY274以來反演的水位變化值；在此基礎(chǔ)上，與實(shí)測(cè)水位變化值進(jìn)行對(duì)比分析，選出較好的去趨勢(shì)項(xiàng)方法。

4)多頻多模SNR信號(hào)融合。分別利用均值算法、中值算法、隨機(jī)森林算法和BP神經(jīng)網(wǎng)絡(luò)算法對(duì)GPS和BDS多頻多模信號(hào)反演結(jié)果進(jìn)行融合，并分析其精度。

3 結(jié)果與分析

3.1 SNR信號(hào)去趨勢(shì)項(xiàng)方法對(duì)比分析

為了比較二階多項(xiàng)式法和小波分解法去趨勢(shì)項(xiàng)的效果，分別用這兩種方法對(duì)原始SNR觀測(cè)序列進(jìn)行去趨勢(shì)項(xiàng)處理。以2017年10月1日(DOY 274)G26衛(wèi)星觀測(cè)到的S1C頻段SNR數(shù)據(jù)為例，圖3所示為利用二階多項(xiàng)式擬合去趨勢(shì)項(xiàng)和頻譜分析的序列圖。其中，圖3(a)為去趨勢(shì)項(xiàng)前后的SNR序列圖，黑色線為原始SNR序列，紅色線為SNR的去趨勢(shì)項(xiàng)序列，藍(lán)色線為消除趨勢(shì)項(xiàng)信號(hào)后的SNR殘差序列；圖3(b)為SNR殘差序列LS譜分析曲線圖。從圖3可以直觀地看出，對(duì)于S1C頻段的SNR數(shù)據(jù)，傳統(tǒng)的二階多項(xiàng)式法去趨勢(shì)效果比較理想。

圖4所示為利用db5作為小波基對(duì)上述SNR序列進(jìn)行5層小波分解處理后的結(jié)果。其中，黑色線為原始SNR序列；深藍(lán)色線為小波分解的5層高頻信號(hào)，即細(xì)節(jié)信號(hào)，反映原始SNR序列的細(xì)節(jié)信息；淺藍(lán)色線和紅色線為小波分解下的5層低頻信號(hào)，也稱作逼近信號(hào)，是原始SNR序列緩慢變化的部分，紅色線為趨勢(shì)項(xiàng)。將趨勢(shì)項(xiàng)與原始SNR序列相減可提取SNR殘差序列，結(jié)果如圖5(a)所示，對(duì)所得到的dSNR殘差序列進(jìn)行LS譜分析即可求得等效天線高，圖5(b)所示為G26衛(wèi)星S1C頻段的SNR殘差序列LS譜分析曲線圖。結(jié)合圖3和圖5可知，對(duì)于S1C頻段的SNR數(shù)據(jù)，二階多項(xiàng)式法和小波分解法去趨勢(shì)項(xiàng)效果均比較理想。

(a)

圖4 小波分解的原始SNR序列高頻和低頻信號(hào)(2017, DOY274，G26衛(wèi)星S1C頻段信號(hào))

(a) 消除趨勢(shì)項(xiàng)后的SNR殘差序列

為了進(jìn)一步分析不同去趨勢(shì)項(xiàng)方法對(duì)各SNR頻段的適用性，圖6所示為2017年DOY 274～360期間，分別利用二階多項(xiàng)式和小波分解對(duì)GPS的S1C、S2L、S2W和S5Q這4個(gè)頻段SNR觀測(cè)值去趨勢(shì)項(xiàng)后采用LS譜分析反演得到的水位曲線，同時(shí)與實(shí)測(cè)水位曲線進(jìn)行對(duì)比。總體上看，GPS的S1C、S2L和S5Q這3個(gè)頻段SNR的水位反演結(jié)果與實(shí)測(cè)水位有較好的一致性。其中，S1C和S5Q頻段，如圖6(a)、(d)所示，利用二階多項(xiàng)式擬合和小波分解去趨勢(shì)項(xiàng)后反演的水位曲線高度一致甚至幾乎重合。此外，S2L頻段小波分解的效果明顯優(yōu)于二階多項(xiàng)式擬合的效果，這是由于試驗(yàn)所用的GPS衛(wèi)星還未完全實(shí)現(xiàn)現(xiàn)代化，所用S2L頻段是混合了L2C碼與L2P碼的混合信號(hào)，其中的L2P碼信號(hào)大大降低了反演精度[13]，甚至有些研究認(rèn)為S2L信號(hào)不適用于水位反演[6-7]。另外，對(duì)于S2W頻段，如圖6(c)所示，不管是用二階多項(xiàng)式擬合還是小波分解的方法，其水位反演結(jié)果均不夠理想，這是由于S2W信號(hào)采用Z跟蹤(Z-Tracking)技術(shù)導(dǎo)致噪聲增加[21-22]。綜上可知，對(duì)于SNR觀測(cè)值去趨勢(shì)項(xiàng)方法的選擇，小波分解法整體上優(yōu)于二階多項(xiàng)式擬合法，但是GPS-S2W頻段用于水位反演的效果不甚理想。

(a) S1C

需要特別注意的是，圖6(c)的S2W頻段反演水位結(jié)果顯示，在DOY310之前小波分解結(jié)果與實(shí)測(cè)值符合較好，DOY310～315之間反演結(jié)果出現(xiàn)了突變，在DOY315之后又與實(shí)測(cè)結(jié)果有良好的一致性。究其原因，是由于GPS S2W頻段信號(hào)采用的Z跟蹤技術(shù)會(huì)導(dǎo)致L2信號(hào)被L1信號(hào)污染，出現(xiàn)“雙峰”現(xiàn)象[21]——S2W頻段第一個(gè)峰值對(duì)應(yīng)L2信噪比的多徑振蕩，第二個(gè)峰值為L(zhǎng)1污染L2產(chǎn)生的噪聲信號(hào)。為了分析GPS S2W“雙峰”現(xiàn)象對(duì)水位反演造成的干擾，圖7繪制了S2W信號(hào)的LS譜分析圖和水位反演結(jié)果。其中，圖7(a)為2017年DOY 314的S2W頻段的 LS譜分析圖，可以看出，“雙峰”的振幅非常接近，難以進(jìn)行有效區(qū)分?；诖耍瑘D6(c)在DOY310～315之間才會(huì)出現(xiàn)反演結(jié)果突變現(xiàn)象。圖7(b)為分別利用“雙峰”反演的水位時(shí)間序列圖，其中第一個(gè)主頻峰值的水位反演結(jié)果與實(shí)測(cè)水位有較好的一致性，其均方根誤差、平均偏差和相關(guān)系數(shù)分別為0.1439m、0.1094m和0.9987，但是反演的精度仍不如GPS其他頻段；第二個(gè)主頻峰值的水位反演結(jié)果與實(shí)測(cè)水位有較大的誤差，其均方根誤差、平均偏差和相關(guān)系數(shù)分別為0.6792m、0.5578m和0.9420。

(a)

為了進(jìn)一步驗(yàn)證上述結(jié)果，表1列出了GPS各頻段信號(hào)反演的水位相對(duì)于實(shí)測(cè)水位的精度指標(biāo)。

從表1可以明顯看出，與實(shí)測(cè)水位相比，無論采用二階多項(xiàng)式法還是小波分解方法，S2W頻段的反演結(jié)果的均方根誤差達(dá)到甚至超過了1m，再次驗(yàn)證了該頻段SNR觀測(cè)量用于水位反演效果不佳的結(jié)論。此外，小波分解和二階多項(xiàng)式方法在S1C頻段的反演精度相當(dāng)，而在S2L、S5Q頻段，小波分解方法明顯優(yōu)于二階多項(xiàng)式法。限于篇幅，BDS各頻段SNR去趨勢(shì)項(xiàng)后的水位反演結(jié)果僅用表格的形式列于表2中。從表2可以看出，對(duì)于S1I、S6I和S7I這3個(gè)頻段，小波分解方法的反演精度略優(yōu)于二階多項(xiàng)式擬合法。綜合表1和表2的數(shù)據(jù)，剔除S2W頻段第二主頻峰值反演結(jié)果后,對(duì)所有頻段反演結(jié)果取算數(shù)平均值，發(fā)現(xiàn)小波分解反演的水位與實(shí)測(cè)水位對(duì)比得到的平均均方根誤差、平均偏差和平均相關(guān)系數(shù)分別為0.1062m、0.0304m和0.9986，而二階多項(xiàng)式擬合法的精度結(jié)果分別為0.2245m、0.1147m和0.9978。

表1 各頻段GPS-SNR小波分解法和二階多項(xiàng)式法去趨勢(shì)項(xiàng)后反演水位精度對(duì)比

表2 各頻段BDS-SNR小波分解法和二階多項(xiàng)式法去趨勢(shì)項(xiàng)后反演水位精度對(duì)比

從理論上來說，傳統(tǒng)的二階多項(xiàng)式擬合法雖然可用于處理大多數(shù)衛(wèi)星頻段的SNR數(shù)據(jù)，但其基本原理是基于最小二乘思想，即基于殘差的平方和尋找數(shù)據(jù)的最優(yōu)函數(shù)匹配，等權(quán)地處理所有數(shù)據(jù)。當(dāng)原始SNR序列中出現(xiàn)一部分質(zhì)量不佳的信號(hào)時(shí)，殘差平方和會(huì)對(duì)擬合結(jié)果產(chǎn)生比較明顯的影響，使擬合曲線偏移，導(dǎo)致擬合出的趨勢(shì)項(xiàng)與實(shí)際直射信號(hào)相比出現(xiàn)較大的偏差[23]。而小波分解是基于信號(hào)學(xué)理論，通過頻率變化的快慢將非平穩(wěn)的SNR信號(hào)分離成低頻信號(hào)和高頻信號(hào)，從而實(shí)現(xiàn)去趨勢(shì)項(xiàng)。即使原始SNR序列中出現(xiàn)一部分質(zhì)量不佳的信號(hào)，也不會(huì)影響高頻低頻信號(hào)的分離效果。因此，理論分析與試驗(yàn)驗(yàn)證均說明，小波分解法用于原始SNR信號(hào)去趨勢(shì)項(xiàng)的效果要優(yōu)于二階多項(xiàng)式擬合法。

綜上可知，在水庫(kù)水位高度反演消除SNR趨勢(shì)項(xiàng)的過程中，小波分解法總體上優(yōu)于傳統(tǒng)的二階多項(xiàng)式擬合法，且S2W頻段信號(hào)用于水位反演效果不佳。同時(shí)也可以發(fā)現(xiàn)，不同系統(tǒng)不同頻段數(shù)據(jù)之間存在不同的系統(tǒng)偏差，使得在同一水位基準(zhǔn)下不同系統(tǒng)不同頻段數(shù)據(jù)的反演結(jié)果隨著DOY的增加出現(xiàn)偏移現(xiàn)象，需要根據(jù)誤差調(diào)整權(quán)重，通過多星多模融合以進(jìn)一步減少系統(tǒng)誤差的影響。

3.2 多頻多模信號(hào)融合反演水位變化分析

雖然單一頻段SNR信號(hào)能夠用于水庫(kù)水位反演，但是可用衛(wèi)星信號(hào)數(shù)量有限，難以發(fā)現(xiàn)可能存在的粗差，且不同頻段之間可能存在系統(tǒng)性偏差，因此要考慮進(jìn)行多星多頻信號(hào)融合反演[24]。在利用小波分解對(duì)SNR去趨勢(shì)項(xiàng)和LS譜分析獲得其反演結(jié)果的基礎(chǔ)上，分別利用均值法、中值法、隨機(jī)森林算法和BP神經(jīng)網(wǎng)絡(luò)算法這四種信號(hào)融合算法對(duì)測(cè)站接收到的全部GPS、BDS頻段信號(hào)反演結(jié)果進(jìn)行融合，并與實(shí)測(cè)水位進(jìn)行比較，結(jié)果如表3所示。

由表3可知，將所有頻段的反演水位進(jìn)行多星多頻融合后，均方根誤差只達(dá)到了分米級(jí)，甚至不如某些單頻信號(hào)(S1C、S2L、S5Q、S7I)的反演結(jié)果，總體來說反演效果不太理想。這是由于融合過程中反演精度較低的頻段作為輸入項(xiàng)會(huì)影響整體的反演精度。因此，在信號(hào)融合前需要設(shè)置閾值將反演精度較差的頻段進(jìn)行剔除。在確保數(shù)據(jù)利用率的前提下，分別設(shè)置0.1m、0.15m和0.2m為閾值，將均方根誤差小于閾值的頻段設(shè)置為可取頻段，大于閾值的頻段設(shè)置為不可取頻段予以去除，經(jīng)過比較發(fā)現(xiàn)，選擇0.1m的閾值最為合適。由表1可知，SW50測(cè)站GPS的S1C、S2L和S5Q均為可取頻段；由表2可知，BDS的S1I和S6I頻段的均方根誤差均超過閾值0.1m，均為不可取頻段，只有S7I頻段的反演結(jié)果達(dá)到厘米級(jí)，為可取頻段。水位反演所用信號(hào)頻段和對(duì)應(yīng)衛(wèi)星編號(hào)如表4所示。

表4 水位反演所用信號(hào)頻段和對(duì)應(yīng)衛(wèi)星編號(hào)

將前50天各頻段單獨(dú)反演結(jié)果作為BP神經(jīng)網(wǎng)絡(luò)和隨機(jī)森林算法的訓(xùn)練樣本進(jìn)行訓(xùn)練，輸出單元為后37天的水位變化結(jié)果。為了避免網(wǎng)絡(luò)收斂過慢，對(duì)輸入數(shù)據(jù)進(jìn)行歸一化處理。其中BP神經(jīng)網(wǎng)絡(luò)采用37-1-1型網(wǎng)絡(luò)結(jié)構(gòu)，采用L-M算法，各參數(shù)設(shè)置為[20]：性能函數(shù)mse，輸入傳遞函數(shù)tansig，輸出傳遞函數(shù)purelin，目標(biāo)誤差0.0003，學(xué)習(xí)速率0.01，最大訓(xùn)練次數(shù)1000。

將均值法、中值法、隨機(jī)森林和BP神經(jīng)網(wǎng)絡(luò)四種融合算法得到的水位變化結(jié)果與實(shí)測(cè)水位變化進(jìn)行對(duì)比分析，結(jié)果如表5所示。

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

由表5可知，對(duì)于水面較平靜的環(huán)境，四種算法的均方根誤差和平均偏差均達(dá)到了厘米級(jí)精度，且相關(guān)系數(shù)高于0.999，說明均具有良好的反演效果。其中，從均方根誤差和平均偏差2個(gè)指標(biāo)來衡量，BP神經(jīng)網(wǎng)絡(luò)算法的水位反演效果明顯優(yōu)于其他三種算法。從理論上來看，均值法和中值法對(duì)各頻段信號(hào)的反演結(jié)果分配的權(quán)重相同，因此對(duì)反演效果不佳的頻段不能有效抑制其反演誤差；隨機(jī)森林算法中雖然每個(gè)決策樹的訓(xùn)練樣本是隨機(jī)抽取的，但是每個(gè)樹的權(quán)重相等，其本質(zhì)也是一種特殊的等權(quán)算法；而BP神經(jīng)網(wǎng)絡(luò)可以通過每次訓(xùn)練的誤差來修正權(quán)重，使網(wǎng)絡(luò)的實(shí)際輸出值(反演水位結(jié)果)和期望輸出值(實(shí)測(cè)水位結(jié)果)之間的均方根誤差最小化。因此，理論分析和實(shí)測(cè)驗(yàn)證均說明BP神經(jīng)網(wǎng)絡(luò)算法更優(yōu)。

4 結(jié)論

以雙王城水庫(kù)SW50測(cè)站2017 年10月1日—12月26日的GNSS SNR原始觀測(cè)量為研究對(duì)象，全面分析了GPS和BDS雙系統(tǒng)多頻SNR信號(hào)的變化特征，在此基礎(chǔ)上提出了聯(lián)合小波分解與BP神經(jīng)網(wǎng)絡(luò)法反演水庫(kù)水位變化的策略，并得到以下結(jié)果：

1)在SNR去趨勢(shì)項(xiàng)方法的選擇上，小波分解算法優(yōu)于傳統(tǒng)的二階多項(xiàng)式擬合算法。

2)利用多頻多模SNR信號(hào)反演水庫(kù)水位變化，需要提前設(shè)置閾值以消除粗差、提高精度，并且在水庫(kù)、大壩等水面較為平靜的反演環(huán)境下，均值算法、中值算法、隨機(jī)森林算法和BP神經(jīng)網(wǎng)絡(luò)算法均能達(dá)到厘米級(jí)精度，其中BP神經(jīng)網(wǎng)絡(luò)算法用于水位反演的效果更優(yōu)。

同時(shí)，本研究首次對(duì)GPS S2W信號(hào)雙峰現(xiàn)象用于水位反演進(jìn)行探討，發(fā)現(xiàn)由于采用Z跟蹤技術(shù)，GPS S2W信號(hào)所呈現(xiàn)的“雙峰”現(xiàn)象對(duì)水位變化反演產(chǎn)生了明顯干擾，在后續(xù)工作中需要對(duì)此現(xiàn)象進(jìn)一步研究以準(zhǔn)確提取有用的信號(hào)。此外，由于所選測(cè)站硬件設(shè)備的限制，BDS數(shù)據(jù)不包含北斗3號(hào)的衛(wèi)星信號(hào)，未來需對(duì)更多的北斗衛(wèi)星反射信號(hào)反演厘米級(jí)水位的方法進(jìn)行探討。