姜桂婷， 李水艷

(河海大學(xué) 理學(xué)院， 江蘇 南京 211100)

0 引 言

潮位的精確檢測(cè)對(duì)于我國(guó)海洋事業(yè)的發(fā)展起著不可或缺的作用. 目前， 有成千上萬(wàn)的全球定位導(dǎo)航系統(tǒng)(Global Navigation Satellite System， GNSS)[1]正在全世界運(yùn)行， 多路徑原理[2]即信號(hào)中同時(shí)包含直射信號(hào)以及周?chē)乇淼姆瓷湫盘?hào)信息， 從而利用現(xiàn)成的大量商用GNSS提供其周?chē)h(huán)境的幾何信息[3-4].

GNSS干涉遙感(GNSS-Interferometr Reflectomety, GNSS-IR)技術(shù)基于常規(guī)的大地測(cè)量型接收機(jī)， 利用SNR序列中的干涉振蕩特性， 便可完成對(duì)站點(diǎn)周?chē)h(huán)境的監(jiān)測(cè)[5-7]. Larson等[8]首次將地基(岸基)GPS 的信噪比(Signal-to-Noise Ratio， SNR) 序列用于潮位監(jiān)測(cè)中. 此后， Larson等[9]在 Kachemak 灣利用 PBAY 站的 SNR 信號(hào)進(jìn)行潮位反演， 并且提出改正動(dòng)態(tài)海面偏差的方法. 隨后， L?fgren等[10]又進(jìn)一步研究改進(jìn)此方法.

目前， GNSS-IR 潮位反演技術(shù)已經(jīng)有一套較為完善的經(jīng)典潮位反演原理及技術(shù)流程： 即利用 LSP(Lomb-Scargle Periodogram)提取 SNR 中的頻率信息， 并將該信息轉(zhuǎn)換為天線相位中心與反射面之間的垂直距離(Reflector Height， RH)， 簡(jiǎn)稱(chēng)為有效高度， 再將其統(tǒng)一到潮位基準(zhǔn)上. LSP 方法對(duì)于一段輸入的 SNR 序列只能輸出一個(gè)頻率值. 海面的時(shí)刻變化特性表明對(duì)于具有一定時(shí)間長(zhǎng)度的SNR來(lái)說(shuō)， 與海面高度相關(guān)的頻率也隨之變化， 由此可知LSP方法的不合理性， 且其損失了反演點(diǎn)的分辨率. 針對(duì)這個(gè)問(wèn)題， 前人提出Window Lomb-Scargle譜分析(winLSP)方法， 通過(guò)滑動(dòng)窗口， 用LSP方法提取由窗口調(diào)整的每個(gè)段的頻率. 該方法雖然能夠輸出更多的反演點(diǎn)值， 但其結(jié)果精度遠(yuǎn)不及 LSP方法.

小波變換常作為數(shù)據(jù)信噪比分析研究的一種時(shí)頻分析方法. Bilich等[11]利用小波分析提取的SNR序列的光譜含量來(lái)確定造成顯著多路徑錯(cuò)誤的因素. Wang 等[12-13]利用連續(xù)小波變換來(lái)提取多模多頻 SNR 序列中的瞬時(shí)頻率， 以期望提高 GNSS-IR 潮位反演中的數(shù)據(jù)分辨率， 與此同時(shí)證明了該頻率在潮位檢索中的潛力.

對(duì)于多模多頻站點(diǎn)潮位分析， 利用連續(xù)小波變換可以得到不同SNR序列的小波譜圖. 本文基于系統(tǒng)不同質(zhì)量的SNR序列分析， 給出一定的正則化降噪模型， 從而獲得一個(gè)高質(zhì)量的RH序列來(lái)進(jìn)行潮位分析.

1 背景知識(shí)

1.1 GNSS-IR分析原理

信號(hào)接收器接收來(lái)自衛(wèi)星的直射信號(hào)與反射信號(hào). 在只有一次多路徑反射時(shí)， 與直射信號(hào)相比， 反射信號(hào)傳播了額外的距離D=D1-D2， 如圖 1 所示.

圖1 一次多路徑反射原理圖Fig.1 Single multipath reflection schematic diagram

由圖 1 中的幾何關(guān)系容易得到

(1)

根據(jù)式(2)計(jì)算出直射信號(hào)與反射信號(hào)之間的相位差φ

(2)

式中：λ為信號(hào)的波長(zhǎng).

在靜態(tài)或準(zhǔn)靜態(tài)的情況下， 即不考慮海面和衛(wèi)星的運(yùn)動(dòng)時(shí)， 有

(3)

式中：f為正弦線SNR序列的頻率； d表示微分符號(hào).

由此得

(4)

考慮到海面時(shí)刻都在變化的特性， 特別是在潮波變化較大的地方， 靜態(tài)假設(shè)變得不合理. 高度速率誤差可以通過(guò)增加改變海面高度的一階修正項(xiàng)來(lái)修正， 即

修正后的RH時(shí)間序列h為

(5)

最終海平面高度為

s=C-h(t),

(6)

式中：s為海平面高度；C為恒定常數(shù)， 與站高、 GNSS系統(tǒng)的框架基準(zhǔn)和驗(yàn)潮站的高度基準(zhǔn)有關(guān).

1.2 連續(xù)小波變換

與傅里葉變換一樣， 連續(xù)小波變換使用內(nèi)積來(lái)度量信號(hào)f和分析函數(shù)之間的相似性， 然而連續(xù)小波變換不再使用三角函數(shù)基， 而是使用隨時(shí)間急速衰減的小波基ψ， 通過(guò)平移和伸縮， 得到小波系數(shù)矩陣

(7)

式中：a為尺度參數(shù)， 控制小波的伸縮；b為平移參數(shù)， 控制小波的平移.

連續(xù)小波變換同時(shí)從時(shí)域和頻域來(lái)分析信號(hào)， 其頻率信息包含在尺度中， 尺度減小， 中心頻率提高， 反之亦然. 通過(guò)連續(xù)小波變換可以了解信號(hào)的頻率成分， 與此同時(shí)還能確定它在時(shí)域上存在的具體位置. Abry在1997年給出了在給定尺度上具有指定采樣周期的小波與頻率之間的映射.

2 數(shù)據(jù)處理過(guò)程

2.1 連續(xù)小波變換提取瞬時(shí)頻率

對(duì)于δS(sine)序列， 其連續(xù)小波變換為

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

(8)

式中： *代表復(fù)共軛；δS(sine)為通過(guò)二次擬合去除直射信號(hào)后的SNR序列.

利用a和b的值， 可以得到小波變換系數(shù)矩陣W(a,b). 同時(shí)， MATLAB小波工具箱軟件提供了2個(gè)函數(shù)： centfrq 和scal2frq， 給出了指定尺度小波和實(shí)際頻率的近似關(guān)系. 小波分析頻譜中， 每個(gè)歷元對(duì)應(yīng)多個(gè)頻率(即有效高度)， 本文選擇每一歷元能量最大值對(duì)應(yīng)的RH值， 認(rèn)為該 RH 值是該歷元的最顯著頻率對(duì)應(yīng)的RH值.

2.2 頻率域加權(quán)

基于期望得到的新的RH序列具有一定的連續(xù)性， 可通過(guò)計(jì)算每個(gè)RH序列異常點(diǎn)的個(gè)數(shù)來(lái)確定其加權(quán)的權(quán)重， 異常點(diǎn)越多的序列所占權(quán)重越小.

若得到的RH序列分別為ri(i=1,2…)， 稱(chēng)滿足條件|r(i+1)-r(i)|>0.5的點(diǎn)為異常點(diǎn)， 其中r(i)代表r序列的第i個(gè)元素.

設(shè)其對(duì)應(yīng)異常點(diǎn)個(gè)數(shù)為si(i=1,2…)， 選取權(quán)重為

(9)

根據(jù)經(jīng)驗(yàn)可得， 頻域率的加權(quán)比時(shí)間域的加權(quán)更加適合， 因此， 進(jìn)一步加權(quán)從SNR序列中得到的小波系數(shù)矩陣Wi(i=1,2，…)， 可以得到

(10)

取W每列小波系數(shù)最大處對(duì)應(yīng)的RH值， 得到所需要的新r序列.

2.3 正則化去噪

正則化方法最常用的一個(gè)應(yīng)用領(lǐng)域就是去噪. 正則化過(guò)程是在已知r序列的基礎(chǔ)上尋找一個(gè)更好的估計(jì). 由于希望能在保證連續(xù)平滑性的基礎(chǔ)上減少噪音， 因此， 得到的最小二乘問(wèn)題為

(11)

則此優(yōu)化問(wèn)題的解為

(12)

3 數(shù)值實(shí)驗(yàn)與比較

3.1 數(shù)據(jù)來(lái)源

站點(diǎn)選取了能夠接收到多模多頻GNSS信號(hào)的Multi-GNSS Experiment (MGEX)站： MAYG站.

MAYG站位于印度洋附近的馬約特島(12.8°S, 45.3°E)， 安裝Trimble TRM59800.00天線(不帶天線罩)， 該天線與Trimble NETR9 GNSS 接收機(jī)連接， 記錄采樣間隔30 s的觀測(cè)數(shù)據(jù). 該站可以接收到多模多頻信號(hào). 在位于MAYG 站幾米處政府間海洋學(xué)委員會(huì)(Intergovernmental Oceanographic Commission) 建立了Dzaoudzi驗(yàn)潮站， 該潮位站以1 min的時(shí)間分辨率收集潮位數(shù)據(jù)， 以此得到了實(shí)測(cè)潮位. MAYG站的海域方位角為20°～170°.

本文選取MAYG站2017年第250年積日GPS PRN 07, GLONASS PRN 23, Galileo PRN 09衛(wèi)星的SNR 數(shù)據(jù)進(jìn)行分析. 具體使用的SNR類(lèi)型如表 1 所示.

表 1 SNR數(shù)據(jù)類(lèi)型Tab.1 SNR data type

表 1 中， 對(duì)于GPS系統(tǒng)， 使用3種SNR 數(shù)據(jù)類(lèi)型： S1C， S2W， S2X； 對(duì)于GLONASS系統(tǒng)， 使用3種SNR數(shù)據(jù)類(lèi)型： S1C， S1P， S2C； 對(duì)于Galileo系統(tǒng)， 使用3種SNR數(shù)據(jù)類(lèi)型： S1X， S5X， S7X.

3.2 數(shù)據(jù)實(shí)驗(yàn)

首先對(duì)GPS系統(tǒng)3種SNR序列進(jìn)行實(shí)驗(yàn)研究， 3種SNR序列如圖 2 所示.

利用連續(xù)小波變換提取SNR序列的瞬時(shí)頻率， 得到如圖 3 所示的小波頻譜圖.

(a) S1C

(b) S2W

(c)S2X圖3 3種信號(hào)對(duì)應(yīng)的小波頻譜圖Fig.3 Wavelet spectrum corresponding to three signals

圖 3(a) 是由GPS系統(tǒng)S1C類(lèi)型SNR序列得到的小波時(shí)頻圖， 圖 3(b) 是由GPS系統(tǒng)S2W類(lèi)型SNR序列得到的小波時(shí)頻圖， 圖 3(c)是由GPS系統(tǒng)S2X類(lèi)型SNR序列得到的小波時(shí)頻圖. 小波譜圖不僅給出了每個(gè)歷元不同RH(或頻率)的功率(或幅度)， 而且還反映了SNR序列的質(zhì)量. 質(zhì)量較好的SNR系列具有邊界清晰的功率集中區(qū)域和較少的由噪聲引起的散亂區(qū)域. 根據(jù)圖中小波分析譜圖的表現(xiàn)， 可以看出每個(gè)SNR類(lèi)型的質(zhì)量等級(jí)為： S2X> S2W> S1W.

對(duì)于3幅小波時(shí)頻圖選取最大能量處對(duì)應(yīng)的RH值， 得到如圖 4 所示的RH序列.

應(yīng)用式(9)和式(10)得到r序列， 需要先對(duì)序列的異常點(diǎn)進(jìn)行基本預(yù)處理

r(i+1)=1.001*r(i).

(a) S1C

(b) S2W

(c) S2X圖4 3種信號(hào)對(duì)應(yīng)的RH序列Fig.4 RH sequence corresponding to three signals

對(duì)于得到的新r序列， 應(yīng)用模型(11)進(jìn)行降噪， 得到如圖 5 所示的序列， 其中， 正則化參數(shù)矩陣對(duì)角元素為

圖5 正則化降噪后的RH序列Fig.5 RH sequence after regularization and noise reduction

3.3 潮位對(duì)比

GPS系統(tǒng)3種不同類(lèi)型RH序列通過(guò)式(5)和式(6)得到的潮位數(shù)據(jù)與實(shí)測(cè)潮位對(duì)比如圖 6 所示.

(a) L1C/A 類(lèi)型

(b) L2W類(lèi)型

(c) L2C類(lèi)型圖6 未采用正則化去噪處理得到的潮位數(shù)據(jù)與實(shí)測(cè)潮位的對(duì)比Fig.6 Comparison between tidal level data obtained withoutregularization denoising and measured tidal level

圖 6 中黑色點(diǎn)代表實(shí)測(cè)的潮位數(shù)據(jù)， 灰色點(diǎn)代表由GPS不同類(lèi)型SNR數(shù)據(jù)推導(dǎo)出的潮位.

降噪處理后推導(dǎo)得到的潮位與實(shí)測(cè)潮位的對(duì)比如圖 7 所示.

通過(guò)對(duì)比圖 6 與圖 7， 能直觀地看到降噪處理后推導(dǎo)出的潮位與實(shí)際潮位更加接近. 為了使對(duì)比更具統(tǒng)計(jì)學(xué)的意義， 比較其均方根誤差(Root Mean Square Error， RMSE)， 該方法衡量了推導(dǎo)值與記錄值之間的誤差， 數(shù)值越小說(shuō)明誤差越小.

同時(shí)本文也用此方法對(duì)GLONASS和Galileo系統(tǒng)的不同類(lèi)型的SNR數(shù)據(jù)進(jìn)行了潮位分析， SNR數(shù)據(jù)類(lèi)型與表1相對(duì)應(yīng)， 由于篇幅有限， 就不展開(kāi)敘述， 得到具體結(jié)果的均方根誤差如表 2 所示.

圖7 降噪處理后的潮位與實(shí)測(cè)潮位對(duì)比Fig.7 Comparison between the tide level obtained bynoise reduction and the measured tide level

表 2 所示為GPS系統(tǒng)第7號(hào)衛(wèi)星， GLONASS第23號(hào)衛(wèi)星和Galileo第9號(hào)衛(wèi)星對(duì)應(yīng)3種SNR數(shù)據(jù)潮位分析結(jié)果的RMSE值， 以及加權(quán)降噪后的RMSE值， 數(shù)值對(duì)比發(fā)現(xiàn)RMSE值有明顯的降低. 由此可知本文方法利用同系統(tǒng)不同類(lèi)型SNR序列的信息， 在提高數(shù)據(jù)利用率的同時(shí)降低了誤差.

表 2 均方根誤差比較Tab.2 Comparison of RMSE

4 結(jié) 論

本文首先利用連續(xù)小波變換來(lái)提取SNR序列的瞬時(shí)頻率， 由于實(shí)際衛(wèi)星接收到的信號(hào)噪音很大， 因此得到的潮位數(shù)據(jù)誤差較大. 然后利用同系統(tǒng)不同類(lèi)型SNR序列的綜合信息進(jìn)行正則化的降噪處理， 與正則化方法處理前的SNR序列進(jìn)行的潮位數(shù)據(jù)相比， 至少能提高20%左右的精度， 具有一定的應(yīng)用和推廣意義.