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

陳 殊 何秀鳳 王笑蕾 宋敏峰

1 河海大學(xué)地球科學(xué)與工程學(xué)院，南京市佛城西路8號(hào)，211100

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

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

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

1 方法與原理

1.1 GNSS-MR反演潮位原理

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

(1)

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

直射信號(hào)與反射信號(hào)相位差為：

(2)

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

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

(3)

(4)

1.2 小波分析提取瞬時(shí)潮位

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

(5)

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

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

(6)

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

圖1 小波分析頻譜

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

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

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

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

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

2.1 站點(diǎn)介紹

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

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

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

以BRST站為例，本文選取BRST站2020年年積日197～211的SNR數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，其中信噪比數(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ù)實(shí)驗(yàn)中未采用這兩種信號(hào)。為獲取質(zhì)量更高的SNR數(shù)據(jù)，本文分別從GNSS不同SNR類(lèi)型數(shù)據(jù)LSP反演結(jié)果的精度及其小波譜表現(xiàn)來(lái)評(píng)定SNR數(shù)據(jù)的質(zhì)量等級(jí)。

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

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

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

圖4 BRST站GNSS不同類(lèi)型SNR序列小波譜

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

3 實(shí)驗(yàn)與結(jié)果

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

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

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

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

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

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

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

為進(jìn)一步驗(yàn)證基于小波分析的多模多頻瞬時(shí)潮位反演的有效性，本文又選取HKQT站2020年年積日211～225的SNR數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。其中信噪比數(shù)據(jù)GPS有S1C、S2W、S2X、S5X，GLONASS有S1C、S1P、S2C、S2P，Galileo有S1X、S5X、S7X、S8X，Beidou有S2I、S7I。同樣進(jìn)行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í)測(cè)數(shù)據(jù)具有良好的對(duì)應(yīng)關(guān)系。但LSP方法得到的結(jié)果較少，而小波分析的反演結(jié)果可以清晰反映潮位起伏變化。表3為小波分析結(jié)果和LSP結(jié)果RMSE值、相關(guān)系數(shù)值、日均反演點(diǎn)數(shù)量和時(shí)間間隔統(tǒng)計(jì)。

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

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

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

4 結(jié) 語(yǔ)

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

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