符平貴，匡翠林，楚彬

(1.中南大學(xué)地球科學(xué)與信息物理學(xué)院,長(zhǎng)沙 410083；2.湖南省測(cè)繪科技研究所,長(zhǎng)沙 410007)

0 引言

北斗衛(wèi)星導(dǎo)航系統(tǒng)(BDS)作為中國(guó)自主建設(shè)、獨(dú)立運(yùn)行的衛(wèi)星導(dǎo)航系統(tǒng)，與其他衛(wèi)星系統(tǒng)相比，提供了更多的衛(wèi)星信號(hào)，使得該系統(tǒng)具有全球衛(wèi)星導(dǎo)航反射測(cè)量(GNSS-R)應(yīng)用的巨大潛力.Jin 等[1]首次利用BDS 三頻信號(hào)(B1I、B2I 和B3I)的信噪比(SNR)數(shù)據(jù)與三頻相位組合反演海平面高度，與驗(yàn)潮站結(jié)果相比，相關(guān)系數(shù)為0.87～0.91，均方根誤差(RMSE)小于0.6 m.Wu 等[2]研究了BDS 中地球靜止規(guī)道(GEO)衛(wèi)星的水位反演性能，實(shí)驗(yàn)結(jié)果表明GEO 衛(wèi)星的時(shí)延多普勒和相位測(cè)高具有較高的反演精度，適用于監(jiān)測(cè)海平面變化小的區(qū)域，且時(shí)間分辨率較高.陳發(fā)德等[3]基于MAYG 站的B1C、B2a 和B2 SNR 數(shù)據(jù)對(duì)水位進(jìn)行了反演，實(shí)驗(yàn)結(jié)果表明BDS的三頻數(shù)據(jù)反演精度略低于其他系統(tǒng).Zheng 等[4]分析了多頻多模GNSS-R 反演水位精度，結(jié)果表明北斗二號(hào)衛(wèi)星導(dǎo)航系統(tǒng)(BDS-2) 地球同步軌道衛(wèi)星(GEO)不適用于海岸水位反演，BDS-2 中軌道衛(wèi)星(MEO)的監(jiān)測(cè)精度優(yōu)于BDS-2 傾斜地球同步軌道(IGSO)，北斗三號(hào)衛(wèi)星導(dǎo)航系統(tǒng)(BDS-3)MEO 的監(jiān)測(cè)精度與BDS-2 MEO 相當(dāng).Liu 等[5]對(duì)北斗系統(tǒng)SNR 反演水位高度的潛力進(jìn)行了研究，但并未對(duì)B2b 和B2a+B2b等新信號(hào)進(jìn)行水位反演.上述研究主要集中于如何提高反演精度和不同軌道高度衛(wèi)星的測(cè)高性能，并未對(duì)北斗系統(tǒng)全頻段的水位反演性能進(jìn)行全面分析.

本文主要評(píng)估北斗系統(tǒng)全頻段GNSS-R 水位反演的性能.首先描述GNSS-R 水位反演方法，然后對(duì)比頻譜分析與非線性擬合算法的反演算法精度，在選擇合適的反演算法基礎(chǔ)上，引入Savitzky-Golay 濾波對(duì)北斗全頻段反演結(jié)果進(jìn)行濾波算法去噪，提升水位觀測(cè)精度.

1 GNSS-R 水位反演方法

1.1 GNSS-R 水位反演原理

圖1 為 GNSS-R 反演海平面高度原理圖.圖中h為天線相位中心(APC)到瞬時(shí)海平面的距離，θ 為直射信號(hào)與瞬時(shí)海平面的夾角.

圖1 GNSS-R 反演海平面高度原理

為了有效地抑制測(cè)站周?chē)瓷湮镆鸬亩嗦窂叫?yīng)，直射信號(hào)和反射信號(hào)的振幅值需要滿(mǎn)足如下的關(guān)系：

干涉合成信號(hào)振幅的總體趨勢(shì)取決于衛(wèi)星發(fā)射的直射信號(hào)的振幅Ad，而多路徑引起的反射信號(hào)的振幅Am則表現(xiàn)為低高度角下的周期性波動(dòng).SNR 與干涉后合成信號(hào)振幅之間滿(mǎn)足如下關(guān)系：

式中：Ac為GNSS 測(cè)站天線接收的合成信號(hào)振幅；cos ψ為直射信號(hào)與反射信號(hào)夾角的余弦值.

結(jié)合式(1)可知Ad與Am在數(shù)值上相差較大，通常利用低階多項(xiàng)式去除SNR 的趨勢(shì)項(xiàng)Ad.把圖2(a)中的趨勢(shì)項(xiàng)移除后，剩余項(xiàng)即可認(rèn)為是受多路徑影響所產(chǎn)生的SNR 振蕩序列，如圖2(b)所示，可利用其反演反射體的環(huán)境參數(shù).SNR 振蕩序列可表達(dá)為

圖2 SNR 觀測(cè)序列

式中：A為 振幅值；h為反演的距離；λ 為波長(zhǎng)；θ 為衛(wèi)星高度角；φ為延遲相位.

由式(3)可知，水位反演就是依據(jù) δ SNR 觀測(cè)值求解h.Lomb-Scargle 譜分析和非線性擬合是兩種常用求解h的方法，下面分別給予描述.

1.1.1 Lomb-Scargle 譜分析方法

設(shè)t=sin θ，f=2h/λ，則可將式(3)簡(jiǎn)化成標(biāo)準(zhǔn)的余弦函數(shù)形式：

由式(5)中可知，只需要求解出SNR 振蕩序列的頻率f即可反算出高度h.

SNR 振蕩序列是以高度角正弦值為變量的不等間距信號(hào)，而且接收到的信號(hào)中存在噪聲干擾.在應(yīng)用傅里葉變換時(shí)，這些噪聲和不均勻分布會(huì)對(duì)頻譜分析產(chǎn)生影響，導(dǎo)致頻譜圖與SNR 信號(hào)的實(shí)際結(jié)果有很大偏差，甚至?xí)霈F(xiàn)虛假峰.Lomb-Scargle 頻譜分析方法不僅可以有效地獲得弱周期信號(hào)的頻譜估計(jì)值、解決非均勻采樣信號(hào)的頻譜估計(jì)等問(wèn)題，還可以克服數(shù)據(jù)缺失和序列長(zhǎng)度不足等困難.下面簡(jiǎn)要闡述Lomb-Scargle 頻譜分析方法.

對(duì)于時(shí)序序列X(tj)，j=1,2,3,··,N，其功率譜可定義為關(guān)于頻率的函數(shù)：

式中：Px(f) 為 頻率f的周期信號(hào)功率；X(tj) 為離散實(shí)驗(yàn)數(shù)據(jù)；tj為離散實(shí)驗(yàn)數(shù)據(jù)的時(shí)間；N為實(shí)驗(yàn)數(shù)據(jù)統(tǒng)計(jì)量；τ 為時(shí)間平移不變量.

1.1.2 非線性擬合SNR 信號(hào)方法

Strandberg 等[6]提出利用非線性最小二乘法求解海面高度，將式(3)添加一項(xiàng)衰減因子得到下式:

式中，γ 與反射體的性質(zhì)有關(guān).

為了保證有穩(wěn)定的數(shù)值解，將式(7)轉(zhuǎn)換為式(8)：

式中：C1與C2共同決定式(9)中的振幅A與式(10)的相位φ：

同理，基于SNR 振蕩序列按照式(8)進(jìn)行非線性擬合，即可求解反射高.

1.2 Savitzky-Golay 濾波去噪

受測(cè)站環(huán)境因素影響，GNSS-R 水位反演結(jié)果中不可避免會(huì)存在一些誤差，導(dǎo)致反演精度下降.本文擬通過(guò)Savitzky-Golay 濾波算法[7]降低各種誤差帶來(lái)的影響.Savitzky-Golay 濾波是一種時(shí)間域范圍內(nèi)基于局部多項(xiàng)式的最小二乘擬合濾波方法，它能實(shí)現(xiàn)對(duì)數(shù)據(jù)特征的無(wú)衰減平滑.

下面對(duì)Savitzky-Golary 濾波算法進(jìn)行描述.取某點(diǎn)附近的 2m+1 個(gè) 連續(xù)數(shù)據(jù)進(jìn)行p階多項(xiàng)式擬合，其中p≤2m+1，記該多項(xiàng)式為

式中：F(t) 為擬合后的多項(xiàng)式；bk為多項(xiàng)式系數(shù)；tp為多項(xiàng)式自變量的冪次方.

最小二乘法通過(guò)最小化RMSE 之和進(jìn)行曲線擬合：

式中：S為RMSE；f(i) 為擬合前離散時(shí)間序列在第i點(diǎn)處的值.

為了使該誤差平方項(xiàng)達(dá)到最小值，令誤差表達(dá)式S對(duì)各系數(shù)求導(dǎo)數(shù)并使之為零，即：

因此，給定需要擬合的點(diǎn)數(shù)m和多項(xiàng)式階數(shù)p，求解方程組(14)，就可以求出系數(shù)序列(bp0,bp1,bp2,··,bpp)，從而確定多項(xiàng)式F(t) .

使用Savitzky-Golay 濾波算法對(duì)反演結(jié)果序列進(jìn)行去噪時(shí)，根據(jù)反演結(jié)果序列的特征選擇合適的多項(xiàng)式階數(shù)及窗口長(zhǎng)度.本文選取窗口大小為21，曲線擬合階數(shù)為3.

2 實(shí)驗(yàn)及結(jié)果分析

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

實(shí)驗(yàn)采用國(guó)際GNSS 服務(wù)(IGS)SAS2 站的北斗系統(tǒng)SNR 數(shù)據(jù)進(jìn)行水位反演，圖3(a)為SAS2 站點(diǎn)環(huán)境，圖3(b)為該站方位角80°～360°、高度角5°～30°的菲涅爾區(qū)，可看出GNSS 接收機(jī)距離水面較近，周邊存在一些非水面反射.為客觀評(píng)價(jià)北斗全頻段的反演效果，同時(shí)選取距站點(diǎn)約8 m 的驗(yàn)潮站Sassnitz(http://www.ioc-sealevelmonitoring.org/station.php?code=sass)的實(shí)測(cè)水位作為驗(yàn)證數(shù)據(jù).

圖3 SAS2 站點(diǎn)環(huán)境及菲涅爾區(qū)

2.2 頻譜分析與非線性擬合方法比較分析

本實(shí)驗(yàn)分別使用頻譜分析與非線性擬合對(duì)SAS2 站2021 年年積日1～365 天的北斗全頻段SNR數(shù)據(jù)求取反演結(jié)果并與該站GPS L5 的反演結(jié)果作對(duì)比，研究表明GPS L5 可取得較好的反演精度[8].通過(guò)將反演結(jié)果與驗(yàn)潮站實(shí)測(cè)水位進(jìn)行對(duì)比分析以驗(yàn)證頻譜分析與非線性擬合算法的準(zhǔn)確性.實(shí)驗(yàn)選擇了RMSE、相關(guān)系數(shù)、反演點(diǎn)數(shù)量(所有SNR 振蕩信號(hào)擬合出的反演高度h的數(shù)量)作為GNSS-R 水位反演性能的指標(biāo).

上述實(shí)驗(yàn)結(jié)果相關(guān)統(tǒng)計(jì)指標(biāo)如表1 所示.可以看出：頻譜分析方法中B2a 頻段反演效果相對(duì)于其他頻段較差，B3 頻段反演效果最優(yōu)；非線性擬合方法中B1 頻段反演效果相對(duì)于其他頻段較差，B3 反演效果最優(yōu)；非線性擬合方法的反演結(jié)果都相差不大，其中RMSE 最大與最小的差值為7 cm，而頻譜分析方法RMSE 最大與最小的差值達(dá)到了9.6 cm.北斗全頻段中通過(guò)非線性擬合的反演結(jié)果中，除B1C、B1、B2a 的RMSE 相對(duì)GPS L5 的RMSE 平均提高了0.87 cm，其余頻段的RMSE 都優(yōu)于GPS L5 頻段.而頻譜分析方法的B2a 和B2b 的RMSE 相對(duì)GPS L5 的RMSE 平均提高了2 cm.雖然非線性擬合方法得到的B1C、B1、B3 的反演精度略低于頻譜分析，但其反演時(shí)間分辨率均高于頻譜分析，且其余頻段反演精度均高于頻譜分析，尤其B2a 頻段的RMSE 相對(duì)于頻譜分析提高了17.05%.

表1 北斗全頻段水位反演結(jié)果統(tǒng)計(jì)信息

2.3 水位反演結(jié)果Savitzky-Golay 去噪

受測(cè)站環(huán)境因素的影響，水位反演值中不可避免地會(huì)帶來(lái)誤差，本文使用Savitzky-Golay 去噪算法降低誤差影響.對(duì)2.2 節(jié)中每個(gè)頻段的非線性擬合反演結(jié)果進(jìn)行去噪處理，去噪前后的反演結(jié)果如圖4 所示.圖中僅展示了SAS2 第1 天到第30 天濾波前后水位數(shù)據(jù)與驗(yàn)潮站實(shí)測(cè)數(shù)據(jù)的對(duì)比，可以看出去噪后的反演結(jié)果明顯與驗(yàn)潮站實(shí)測(cè)水位數(shù)據(jù)更加符合.圖5 展示了濾波前后水位數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的殘差圖，可以看出,濾波前的誤差圖中反演誤差變化范圍較大；而使用Savitzky-Golay 去噪后，反演誤差變化范圍小于濾波前的誤差變化范圍.

圖4 水位反演結(jié)果與驗(yàn)潮站數(shù)據(jù)對(duì)比

圖5 水位反演結(jié)果相對(duì)驗(yàn)潮站數(shù)據(jù)的偏差

去噪前后的結(jié)果與驗(yàn)潮站實(shí)測(cè)水位的統(tǒng)計(jì)結(jié)果如表2 所示.可以看出，經(jīng)過(guò)Savitzky-Golay 去噪處理后的水位反演精度均有很大提升，去噪后北斗全頻段的RMSE 相對(duì)去噪前的RMSE 平均減少了7.4 cm，其中提升效果最顯著的是B2a 頻段，其RMSE 減少了9 cm，相關(guān)系數(shù)平均提升了16.40%.

表2 GNSS-R 濾波前后精度統(tǒng)計(jì)

為提高水位監(jiān)測(cè)的時(shí)間分辨率和精度，考慮到北斗不同頻率之間的互補(bǔ)性，將北斗系統(tǒng)全頻段的反演水位組合在一起.圖6(a)是北斗全頻段反演水位值組合，可以看出反演點(diǎn)數(shù)量明顯增加，且時(shí)間分辨率相對(duì)于單一頻率的時(shí)間分辨率更高，其RMSE 為5.4 cm，優(yōu)于B1C、B1 以及B2a 頻段.圖6(b)為反演值與驗(yàn)潮站相關(guān)性分析，可以發(fā)現(xiàn)相關(guān)系數(shù)優(yōu)于B1C、B1以及B2a 等頻段，與驗(yàn)潮站實(shí)測(cè)數(shù)據(jù)有良好的一致性.總的來(lái)說(shuō)，直接將北斗系統(tǒng)不同頻率的反演結(jié)果組合，可達(dá)到提高精度和時(shí)間分辨率的目的.

圖6 北斗全頻段水位反演組合結(jié)果與驗(yàn)潮站比較

3 結(jié)束語(yǔ)

本文系統(tǒng)評(píng)估了我國(guó)北斗全頻段SNR 數(shù)據(jù)的GNSS-R 水位反演性能.首先比較了兩種常用的反演水位算法的精度，實(shí)驗(yàn)表明非線性擬合方法綜合效果優(yōu)于頻譜分析方法，雖個(gè)別頻段反演精度略低于頻譜分析方法，但非線性擬合方法的反演點(diǎn)數(shù)量多于頻譜分析方法.然后，針對(duì)反演結(jié)果中存在隨機(jī)誤差較大，本文通過(guò)Savitzky-Golay 算法對(duì)反演的水位結(jié)果進(jìn)行了濾波去噪，數(shù)值表明濾波后反演水位的RMSE 極大降低，達(dá)到厘米級(jí)，且提升了與潮位數(shù)據(jù)的相關(guān)性.整體上，北斗全頻段的反演性能與GPS L5 相當(dāng).