基于小波分析的多模多頻GNSS-MR潮位反演

陳 殊 何秀鳳 王笑蕾 宋敏峰

1 河海大學地球科學與工程學院，南京市佛城西路8號，211100

全球氣候變暖導致海平面不斷升高，這不僅會影響氣候變化，也對人類生命安全造成嚴重威脅。因此，監(jiān)測海平面高度變化意義重大[1]。近年來，GNSS多路徑反射測量(GNSS multipath reflectometry, GNSS-MR)被用于海面高度探測，目前利用GNSS-MR技術進行潮位監(jiān)測在國內(nèi)外已有一定的研究基礎。Larson等[2]利用大地測量型接收機獲取的GPS實測數(shù)據(jù)對潮位進行反演，實驗表明GPS-MR反演結(jié)果與驗潮站實測數(shù)據(jù)具有較好的一致性。L?fgren等[3]分別使用GPS和GLONASS系統(tǒng)L1和L2頻段信號進行潮位反演，證實GLONASS系統(tǒng)也可用于潮位反演。Roussel等[4]基于最小二乘方法，對GPS和GLONASS系統(tǒng)反演的潮位值進行融合。張雙成等[5]利用不同時間段的GPS信噪比數(shù)據(jù)進行海面高度反演，研究表明GPS-MR反演結(jié)果與實測潮位較差約為10 cm。Jin等[6]首次使用Beidou系統(tǒng)SNR(signal-to-noise ratio)和相位組合方法估計海面變化。Wang等[7]利用基于穩(wěn)健估計的數(shù)據(jù)融合方法進行4個系統(tǒng)(GPS、GLONASS、Galileo和Beidou)聯(lián)合潮位反演，結(jié)果表明多模多頻GNSS-MR可以有效提高潮位反演精度和時間分辨率。從國內(nèi)外研究來看，GNSS-MR技術可以用于監(jiān)測潮位變化。

目前GNSS-MR潮位反演的原理和方法已經(jīng)相對成熟，即利用LSP[8]方法提取SNR中頻率信息，再將該信息轉(zhuǎn)換為垂直反射高度RH值，最后統(tǒng)一到潮位基準上[9]。但是，基于LSP的潮位反演方法在一段SNR序列中通常只能得到一個時刻的潮位值，這樣可能會出現(xiàn)反演點缺失且難以刻畫潮位起伏變化的情況。雖然多模多頻觀測可以提供更多有效衛(wèi)星弧段的SNR數(shù)據(jù)，但存在個別信號質(zhì)量較差，從而影響潮位反演結(jié)果的精度。針對GNSS-MR技術反演時間分辨率不足的問題，Wang等[10]提出利用小波分析方法提取SNR序列瞬時頻率，進而得到潮位值，可大幅提高反演點的時間分辨率。但是小波分析方法過于依賴SNR數(shù)據(jù)質(zhì)量，使得該方法反演潮位的精度有時比LSP方法低。

本文針對小波分析精度偏低的問題，利用多模多頻GNSS數(shù)據(jù)，分析各個頻段SNR數(shù)據(jù)質(zhì)量，選取質(zhì)量更好的SNR數(shù)據(jù)進行小波分析潮位反演，分別以法國BRST站和香港HKQT站15 d的數(shù)據(jù)進行實驗，并與LSP反演結(jié)果進行分析比較。

1 方法與原理

1.1 GNSS-MR反演潮位原理

當GNSS測站位于海邊時，衛(wèi)星直射信號與海水表面反射信號相干的合成SNR可表示為[11]：

(1)

式中，Am和Ad分別為直射信號和反射信號振幅，ψ為直射信號與反射信號之間的相位偏差。

直射信號與反射信號相位差為：

(2)

式中，D為反射信號與直射信號的路徑差，λ為信號波長，e為衛(wèi)星高度角，h為垂直反射距離。

假設反射面靜止，根據(jù)式(2)可知：

(3)

(4)

1.2 小波分析提取瞬時潮位

小波分析是一種多分辨率時頻分析方法，通過伸縮和平移形成一系列靈活窗口對信號進行多尺度細化分析[14]。在小波分析中，分析函數(shù)為一個小波ψ，小波通常被定義為[15]：

(5)

式中，m、n分別為縮放因子和平移因子，二者屬于整數(shù)集合Z={0,±1,±2,…}，t為時間。對于去除趨勢項后的SNR序列在尺度a(>0)和位置b的小波變換可表示為[10]：

W(a,b;δS(sine),ψ(sine))=

(6)

式中，δS(sine)為去除趨勢項后的SNR序列，*表示復共軛，W(a,b)為小波變換系數(shù)。MATLAB中‘centfrq’和‘scal2frq’兩個函數(shù)可以在給定小波尺度情況下計算尺度與采樣頻率之間的關系。

圖1 小波分析頻譜

對于PBAY站，3～11 m是能夠反映潮位變化的有效高度區(qū)間[12]，區(qū)間以外的能量聚集部分則認為是噪聲。在3～11 m范圍內(nèi)，有兩個明顯的能量聚集區(qū)，分別位于橫軸sine=0.15、縱軸6.1 m和橫軸sine=0.41、縱軸6.3 m處。找出每個sine對應的最大振幅值(圖2(a))，再確定其對應的有效高度，即可將sine、振幅值、有效高度的三維圖轉(zhuǎn)變成sine和垂直反射距離的二維圖(圖2(b))。

圖2 小波分析頻譜二維表示

根據(jù)以上步驟即可提取SNR序列瞬時頻率，進而得到瞬時潮位反演值，但由于SNR數(shù)據(jù)質(zhì)量原因，即便選取每個歷元能量最大值對應的有效高度，依然會存在反演值出現(xiàn)跳變的情況。為剔除這種粗差，本文根據(jù)LSP方法獲取的RH值計算標準偏差σ和均值μ，剔除在[μ-3σ,μ+3σ]區(qū)間外的瞬時垂直反射距離，同時利用小波分析得到的瞬時潮位反演值也需要采用潮汐調(diào)和函數(shù)進行海面動態(tài)改正[10]，技術流程如圖3所示。

圖3 基于小波分析反演潮位流程圖

2 多模多頻GNSS站點和數(shù)據(jù)分析

2.1 站點介紹

BRST站(48.4°N，4.5°W)位于法國西海岸Brest海港岸邊，距離BRST站292 m的Brest驗潮站可提供實測數(shù)據(jù)[16]。BRST站配備有Trimble NetR9大地測量型接收機，可提供GPS、GLONASS、Galileo、Beidou和SBAS衛(wèi)星觀測數(shù)據(jù)，采樣間隔為30 s。為獲取來自海面的反射信號，方位角區(qū)間設置為130°～270°，有效高度角區(qū)間為5°～30°。

HKQT站(22.3°N，114.1°E)位于香港北邊海岸，安裝有Trimble NetR5大地測量型接收機，可提供GPS、GLONASS、Galileo、Beidou、QZSS和SBAS衛(wèi)星觀測數(shù)據(jù)，采樣間隔為1 s、5 s和30 s[17]。該站點海面反射信號所對應區(qū)域的方位角為-60°～105°，有效高度角為4°～9°。距離HKQT站2 m的Quarry Bay實測潮位數(shù)據(jù)可用于對比分析。

2.2 數(shù)據(jù)分析

以BRST站為例，本文選取BRST站2020年年積日197～211的SNR數(shù)據(jù)進行實驗，其中信噪比數(shù)據(jù)GPS有4種(S1C、S2W、S2X、S5X)，GLONASS有4種(S1C、S1P、S2C、S2P)，Galileo有5種(S1X、S5X、S7X、S8X、S6X)，Beidou有5種(S1X、S5X、S2I、S7I、S6I)。其中，Galileo中S6X和Beidou中S7I數(shù)據(jù)記錄有缺失，因此在后續(xù)實驗中未采用這兩種信號。為獲取質(zhì)量更高的SNR數(shù)據(jù)，本文分別從GNSS不同SNR類型數(shù)據(jù)LSP反演結(jié)果的精度及其小波譜表現(xiàn)來評定SNR數(shù)據(jù)的質(zhì)量等級。

表1為BRST站4個衛(wèi)星系統(tǒng)每種SNR類型數(shù)據(jù)LSP反演結(jié)果的日均反演點數(shù)量、均方根誤差(RMSE)和相關系數(shù)(CORR)統(tǒng)計結(jié)果，以評估采樣率及反演精度。BRST站GPS、GLONASS、Galileo和Beidou系統(tǒng)中反演精度最高的 SNR類型依次為S5X、S2P、S7X和S6I。不同SNR類型數(shù)據(jù)質(zhì)量等級為：對于GPS，S5X>S2X>S2W>S1C；對于GLONASS，S2P>S2C>S1P>S1C；對于Galileo，S7X>S1X>S5X>S8X；對于Beidou，S6I>S5X >S2I>S1X。

表1 BRST站每種SNR類型的日均反演點數(shù)量、RMSE和CORR

圖4為BRST站2020年年積日201的GPS PRN27、GLONASS PRN16、Galileo PRN33和Beidou PRN21衛(wèi)星不同類型SNR序列小波譜圖。從圖中可以看出，小波譜圖不僅可以確定每一時刻不同垂直反射距離RH的功率，而且還可以反映SNR序列的質(zhì)量。質(zhì)量較好的SNR序列具有明顯的功率集中區(qū)且由噪聲引起的散亂區(qū)較少，而質(zhì)量較差的SNR序列功率分布較離散且噪聲多[9]。從小波譜圖可以看出，每種SNR類型數(shù)據(jù)質(zhì)量等級為：對于GPS，S5X>S2X>S2W>S1C；對于GLONASS，S2P>S2C>S1C>S1P；對于Galileo，S7X>S5X>S1X>S8X；對于Beidou，S6I>S2I >S5X>S1X。

圖4 BRST站GNSS不同類型SNR序列小波譜

比較根據(jù)LSP反演結(jié)果對SNR數(shù)據(jù)質(zhì)量進行評定可以發(fā)現(xiàn)，質(zhì)量最好的SNR類型兩種方法評定結(jié)果一致，其他類型存在個別差異。小波譜的表現(xiàn)可直接決定瞬時反演結(jié)果的準確性，因此需要篩選質(zhì)量較好的SNR類型進行反演以保證精度。

3 實驗與結(jié)果

3.1 BRST站GNSS-MR瞬時潮位結(jié)果

根據(jù)上節(jié)分析，本文最終選擇GPS系統(tǒng)S5X、S2X，GLONASS系統(tǒng)S2P，Galileo系統(tǒng)S7X，Beidou系統(tǒng)S6I數(shù)據(jù)進行GNSS-MR瞬時潮位反演，同時與LSP方法得到的潮位反演結(jié)果進行對比。BRST站LSP和小波分析反演結(jié)果如圖5所示。

圖5 BRST站LSP和小波分析方法潮位反演值

由圖5可知，BRST站利用LSP方法和小波分析方法獲得的潮位結(jié)果均與實際潮位對應良好。其中，LSP方法得到的結(jié)果較少，而小波分析得到的瞬時反演結(jié)果基本能描繪潮位的變化趨勢，滿足潮位監(jiān)測需求。為對兩種方法進行評定并驗證小波分析方法獲取潮位的精度是否可靠，本文將小波分析結(jié)果以及LSP結(jié)果與驗潮站實測數(shù)據(jù)進行比較，表2為RMSE值、相關系數(shù)值、日均反演點數(shù)量和時間間隔統(tǒng)計。

從表2可以看出，小波分析方法和LSP方法相關系數(shù)均優(yōu)于98%，小波分析方法的反演精度比LSP方法低3%～10%，但日均反演點數(shù)量卻比LSP方法提高約12倍，且平均時間間隔約為0.2 h，表明小波分析方法可大大提高反演結(jié)果的時間分辨率。

表2 BRST站小波分析和LSP方法反演結(jié)果統(tǒng)計

3.2 HKQT站GNSS-MR瞬時潮位結(jié)果

為進一步驗證基于小波分析的多模多頻瞬時潮位反演的有效性，本文又選取HKQT站2020年年積日211～225的SNR數(shù)據(jù)進行實驗。其中信噪比數(shù)據(jù)GPS有S1C、S2W、S2X、S5X，GLONASS有S1C、S1P、S2C、S2P，Galileo有S1X、S5X、S7X、S8X，Beidou有S2I、S7I。同樣進行LSP反演結(jié)果精度分析和小波譜圖質(zhì)量分析，最終選擇GPS系統(tǒng)S5X，GLONASS系統(tǒng)S2P，Galileo系統(tǒng)S5X，Beidou系統(tǒng)S7I數(shù)據(jù)。LSP和小波分析反演結(jié)果如圖6所示。

由圖6可知，HKQT站利用LSP方法和小波分析方法獲得的潮位結(jié)果均與實測數(shù)據(jù)具有良好的對應關系。但LSP方法得到的結(jié)果較少，而小波分析的反演結(jié)果可以清晰反映潮位起伏變化。表3為小波分析結(jié)果和LSP結(jié)果RMSE值、相關系數(shù)值、日均反演點數(shù)量和時間間隔統(tǒng)計。

圖6 HKQT站LSP和小波分析方法潮位反演值

表3 HKQT站小波分析和LSP方法反演結(jié)果統(tǒng)計

從表3可以看出，HKQT站LSP方法反演結(jié)果的相關性優(yōu)于小波分析方法，小波分析方法的反演精度比LSP方法低8%～17%，但日均反演點數(shù)量卻比LSP方法提高約13倍，除Beidou外小波分析方法的平均時間間隔也為0.2 h左右。

4 結(jié) 語

利用小波分析方法從SNR序列中獲取潮位值可以大大提高GNSS-MR技術反演潮位的時間分辨率，但小波分析方法的反演精度易受SNR數(shù)據(jù)質(zhì)量的影響。本文針對該問題，利用多模多頻數(shù)據(jù)從LSP反演精度和小波譜圖質(zhì)量兩個方面分析評定4種系統(tǒng)不同SNR類型數(shù)據(jù)的質(zhì)量，進而選擇最優(yōu)的SNR數(shù)據(jù)類型進行瞬時頻率潮位反演。從BRST站和HKQT站的反演結(jié)果來看，基于小波分析的多模多頻GNSS-MR潮位反演方法能夠在精度損失較小的情況下大幅提高反演點數(shù)量和數(shù)據(jù)采樣率。這種基于多模多頻 SNR數(shù)據(jù)的小波分析方法，可充分利用海面多路徑反射信號，實現(xiàn)對SNR數(shù)據(jù)的最大利用。

但是小波分析方法無法脫離LSP方法單獨使用，目前其粗差還需要根據(jù)LSP反演結(jié)果的標準偏差σ和均值μ計算的淘汰區(qū)間進行剔除。同時利用小波分析得到的瞬時潮位反演值也需要采用LSP反演的時間序列進行潮波函數(shù)擬合來確定海面動態(tài)變化情況。未來需要更加合理的粗差剔除方法，從小波譜圖的表現(xiàn)來對粗差進行探測和去除，進一步提高小波分析方法用于潮位反演的精度。