陳殊,何秀鳳,王笑蕾,宋敏峰
河海大學(xué) 地球科學(xué)與工程學(xué)院,南京 211100
沿海地區(qū)是世界上經(jīng)濟(jì)活動和城市化最集中的地區(qū),目前世界上有一半人口生活在離海60 km以內(nèi)地區(qū),而且沿海地區(qū)的人口增長速度正在加快。然而由于風(fēng)暴潮和海平面上升等自然因素,沿海地區(qū)將面臨越來越多威脅,所以對海面變化進(jìn)行實時監(jiān)測并研究其變化規(guī)律具有重要意義(Melet 等,2018)。隨著GNSS 技術(shù)的不斷發(fā)展與完善,GNSS 信號不僅可以用于定位、導(dǎo)航和授時,利用GNSS反射信號還可以監(jiān)測潮位變化(Jin等,2014)。
自Martin-Neira(1993)提出用GNSS反射信號反演海潮變化以來,GNSS-MR 技術(shù)逐漸成熟,成為沿海海平面測高的另一種可行方案(Geremia-Nievinski等,2020)。但是,GNSS-MR潮位反演技術(shù)在當(dāng)前的實際應(yīng)用依然存在精度不夠高和分辨率低的兩大問題(Larson 等,2017)。反演結(jié)果的單點精度因信號類型而異,質(zhì)量好的信號反演的精度更高,采樣率由衛(wèi)星軌跡數(shù)量、signal-tonoise ratio(SNR)數(shù)據(jù)質(zhì)量以及環(huán)境的幾何形狀決定。多模多頻GNSS 可以觀測到更多的衛(wèi)星、獲取質(zhì)量更好的信號和更為開闊的信號感知區(qū)域,這為GNSS-MR 反演潮位提供了更多的反演值從而提高精度和時間分辨率(Larson 等,2017;Tabibi等,2017)。L?fgren 和Haas(2014)分別使用了GPS 和GLONASS 的L1 和L2 頻段的信號進(jìn)行潮位反演,證實了GLONASS 系統(tǒng)也可以用于潮位反演。Roussel等(2015)基于最小二乘方法,對GPS 和GLONASS 系統(tǒng)反演的潮位值進(jìn)行融合。Strandberg等(2016)分別用L1、L2C 和L2P 頻段的SNR 數(shù)據(jù)進(jìn)行潮位反演,并對信號的質(zhì)量進(jìn)行分析。Jin等(2017)第一次用BeiDou系統(tǒng)L2、L6、L7 信號估計海面變化。Wang等(2019)用穩(wěn)健估計的方法對4 個系統(tǒng)(GPS、GLONASS、Galileo 和BeiDou)潮位反演結(jié)果進(jìn)行融合,驗證了多模多頻GNSS-MR 可以有效地提高潮位反演的精度和時間分辨率。何秀鳳等(2020)利用多模多頻GNSS信號反演了沿海臺風(fēng)風(fēng)暴潮,證實了在極端天氣情況下,GNSS-MR 技術(shù)依然可以探測到潮位的異常變化。
目前,多模多頻GNSS-MR 融合算法主要使用“bisquare”定權(quán)方法(Wang等,2019),并將定權(quán)中涉及的常系數(shù)確定為1;這種定權(quán)方法對于偏差大于1 的觀測值無法很好的處理。然而,由于GNSS-MR 潮位反演中存在各類誤差,反演觀測值經(jīng)常大于1 m??紤]到IGGIII 定權(quán)模型(Yang 等,2002)可以對不同大小的偏差進(jìn)行更細(xì)致更合理的處理,文章提出基于IGGIII 模型的穩(wěn)健估計方法來進(jìn)行多模多頻GNSS-MR 融合,并驗證多模多頻GNSS-MR 融合方法在不同測站環(huán)境用于潮位監(jiān)測的可行性。
多路徑效應(yīng)主要與反射面的結(jié)構(gòu)和電介質(zhì)參數(shù)有關(guān),當(dāng)GNSS 測站位于海邊時,接收機(jī)在海域方向會接收到衛(wèi)星直射信號與海水表面反射信號相干的合成信號,這種相干性表現(xiàn)在衛(wèi)星觀測文件SNR的變化上(金雙根 等,2017)。SNR可以表示為(Larson等,2008)
式中,Ad和Am分別表示為直射信號和反射信號的振幅;ψ是直射信號和反射信號之間的相位偏差。直射信號與反射信號的相位差為
式中,D為反射信號與直射信號的路徑差,λ為信號波長,e為衛(wèi)星高度角,h為天線相位中心到反射面的垂直距離。
假設(shè)反射面時靜止的,根據(jù)式(2)有
式中,f是SNR 中受干涉影響部分的信號頻率,是靜態(tài)假設(shè)下得到的反射高度。式(3)中給出了與f的關(guān)系,即對于f的提取,經(jīng)典方法是先用低階多項式來去除SNR 序列中的直射信號,再用LSP(Lomb-Scargle Periodogram)分析方法(Lomb,1976;Scargle,1982)從去除直射信號的殘差序列δSNR中求出f。
為了求解動態(tài)變化的高度h,首先是得到的時間序列,再由此來估算,從時間序列中扣除即可得到改正后的時間序列h。
衛(wèi)星信號在經(jīng)過對流層區(qū)域時會發(fā)生折射效應(yīng),這種對流層延遲會使求得的垂直反射距離出現(xiàn)偏差(Williams 和Nievinski,2017)。在GNSS 反射信號反演海面高的研究中,Anderson(2000)首先提出了大氣折射引起的無線信號幾何彎曲可能會影響反演結(jié)果精度。Santamaría-Gómez 和Watson(2017)在GPS-MR 海平面估計中使用了基于大氣壓力和溫度的經(jīng)驗天文折射模型來消除對流層延遲。Williams 和Nievinski(2017)使 用GPT2W(Global Temperature and Pressure)和VMF1(Vienna mapping functions)映射函數(shù),來消除對流層延遲的影響。對流層延遲τT計算公式為
式中,ΔτZ=τZ(-H) -τZ(0)是天線和地面位置的天頂延遲差,mh和mw分別為干延遲映射函數(shù)和濕延遲映射函數(shù)。于是由對流層延遲τT引起的反演高度偏差ΔhT為
這樣經(jīng)過海面動態(tài)改正和對流層延遲改正后即可得到GNSS 接收機(jī)到海面實際的垂直距離h,再用接收機(jī)的大地高減去垂直距離便是真實的潮位值。
GNSS 系統(tǒng)中各個頻段的SNR 數(shù)據(jù)都可以按照GNSS-MR 經(jīng)典海面反演理論及誤差改正方法對潮位進(jìn)行監(jiān)測。在整個技術(shù)流程中,首先需要確定高度角和方位角區(qū)間,獲得來自海面的反射信號,再對截取的殘差序列δSNR用LSP 分析方法求得有效高度在LSP 分析過程中,需要進(jìn)行質(zhì)量控制,一般選取峰值與背景噪聲的比值(peak-to-noise ratio)大于3 的反演值(Larson 等,2013a)。最后從海面動態(tài)變化和對流層延遲兩個方面對垂直高度進(jìn)行修正。對于多模多頻GNSS-MR 潮位反演值,本文采用基于狀態(tài)方程和穩(wěn)健估計的方法,考慮海面動態(tài)變化以及對流層延遲對數(shù)據(jù)進(jìn)行融合(Wang等,2019),實現(xiàn)了每10分鐘獲得一個潮位值。在建立方程之前,首先需要設(shè)定一個滑動窗口,窗口的寬度需要容納足夠多的潮位反演點以滿足建立方程求解。設(shè)窗口長度為T,窗口的滑動步長為第i個窗口的時間則在第i個窗口內(nèi)可建立狀態(tài)方程
式中,j表示信號類型,l表示窗口內(nèi)反演值的序號,其中|tl-ti|<T/2。則第i個窗口內(nèi)所有的反演值可以用方程組表示
可以簡化為矩陣方程
但是不同系統(tǒng)和頻段信號的SNR 不同,這導(dǎo)致了在LSP 分析過程中每個SNR 弧段的主導(dǎo)頻率存在差異,進(jìn)而反演精度也不相同(Tabibi 等,2017),所以不同的信號應(yīng)對應(yīng)不同的權(quán)重,本文采用基于IGGIII 模型的穩(wěn)健估計方法進(jìn)行自動賦權(quán)。
穩(wěn)健估計的基本原理為:首先通過最小二乘法得到反演值的估計殘差;根據(jù)殘差確定的各反演值新的權(quán)因子,進(jìn)行下一輪最小二乘解算;反復(fù)迭代,直到前后兩次估計值的變化值小于限差(方興 等,2018)。
根據(jù)式(11)進(jìn)行第一次迭代
第一次迭代的殘差為
在接下來的迭代過程中,精度較低的反演值將被重新定權(quán),逐漸降低其在求解垂直距離hi過程中的作用,以便達(dá)到提高精度的目的。潮位反演值主體服從正態(tài)分布,可疑和顯著異常反演值只是占少部分。在穩(wěn)健估計中,應(yīng)充分利用主體可靠信息,所以本文采用IGGIII 模型來自動賦權(quán),圖1 為BRST 站反演數(shù)據(jù)迭代過程中殘差vi所對應(yīng)的權(quán)重,從圖1中可以看出權(quán)重大致分為三段,即正常段,保持原有的權(quán)不變;可疑段,降權(quán),使其小于原始權(quán);淘汰段,權(quán)值為0,淘汰掉異常值,這不但能夠削弱粗差對結(jié)果的影響還可以充分的保留精確數(shù)據(jù),避免了數(shù)據(jù)嚴(yán)重失真(Yang等,2002)。IGGIII模型表示為
圖1 潮位反演值權(quán)重分布Fig.1 Weight distribution of sea level retrievals
在(k+1)次迭代中
本文選取BRST 和HKQT 兩個測站來驗證多模多頻GNSS-MR 技術(shù)用于潮位變化監(jiān)測的有效性,這兩個測站均可接收來自四系統(tǒng)的衛(wèi)星觀測數(shù)據(jù)。
BRST站(48.4°N,4.5°W)位于法國西海岸的Brest 海港岸邊,周圍海平面日漲落幅在7 m 左右(Santamaría-Gómez等,2015)。圖2(a)即為BRST站點的測站環(huán)境,BRST站配備了Trimble NETR9大地測量型接收機(jī),提供了GPS、GLONASS、Galileo、BeiDou和SBAS的衛(wèi)星觀測數(shù)據(jù),采樣間隔為30 s。為了驗證多模多頻GNSS-MR 技術(shù)監(jiān)測潮位的精度,在實驗中利用距離BRST 站292 m 的Brest 驗潮站實測數(shù)據(jù)進(jìn)行對比分析。
圖2 BRST站的觀測環(huán)境和第一菲涅爾反射區(qū)情況Fig.2 Platform at BRST station and the first Fresnel zones for site BRST
由于利用GNSS-MR 技術(shù)監(jiān)測潮位變化僅需獲取來自海面的衛(wèi)星反射信號,本文以5°,10°,20°,30°的高度角來繪制第一菲涅爾反射區(qū)(Roesler 和Larson,2018),BRST站海面信號反射區(qū)域如圖2(b)所示。為了獲取來自海面的反射信號,本文設(shè)置兩個方位角和高度角區(qū)間來剔除陸地的干擾信號。在方位角為130°—165°的區(qū)域,高度角區(qū)間為5°—20°的反射信號被使用;在方位角為165°—330°的區(qū)域,高度角區(qū)間為12°—25°的反射信號被使用(Wang等,2019)。
HKQT站(22.3°N,114.1°E)位于香港北邊海岸,安裝了Trimble NETR5 大地測量型接收機(jī),提供GPS、GLONASS、Galileo、BeiDou、QZSS和SBAS的衛(wèi)星觀測數(shù)據(jù),采樣間隔有1 s、5 s 和30 s。HKQT 測站環(huán)境及海面反射信號如圖3 所示,該站點來自海面反射信號區(qū)域的方位角為-60°—105°,有效高度角為4°—9°。距離HKQT 站2 m 的Quarry Bay實測潮位數(shù)據(jù)可以用于實驗進(jìn)行對比分析。
圖3 HKQT站的觀測環(huán)境和第一菲涅爾反射區(qū)情況Fig.3 Platform at HKQT station and the first Fresnel zones for site HKQT
以BRST 站為例,本文選取BRST 站2020 年年積日為197—253 的SNR 數(shù)據(jù)進(jìn)行實驗,其中信噪比數(shù)據(jù)GPS 有S1C、S2W、S2X、S5X 等4 種,GLONASS 有S1C、S1P、S2C、S2P 等4 種,Galileo有S1X、S5X、S7X、S8X、S6X 等5 種,BeiDou 有S1X、S5X、S2I、S7I、S6I等5種。
實驗期間內(nèi)來自各個系統(tǒng)每天可用的SNR 弧段數(shù)統(tǒng)計結(jié)果如圖4所示,從圖中可以看出GPS和GLONASS 系統(tǒng)實驗期間可用的SNR 弧段數(shù)波動較小,信號質(zhì)量較為穩(wěn)定,每天使用的弧段數(shù)在110左右。Galileo 和BeiDou 系統(tǒng)可用的SNR 弧段數(shù)波動較大,信號質(zhì)量不穩(wěn)定,但是由于記錄的信噪比類型都有5 種,平均每天可用的弧段數(shù)在120 左右。為了進(jìn)一步分析各系統(tǒng)不同波段的SNR 信號質(zhì)量及其反演能力,以2020年年積日為227的4個系統(tǒng)數(shù)據(jù)進(jìn)行分析。圖5 為5 個系統(tǒng)各信號SNR 序列和對應(yīng)的LSP 分析結(jié)果,其中BeiDou 系統(tǒng)由于數(shù)據(jù)記錄不完整,只接收到S1X、S5X、S2I、S6I等4種信號。
圖4 GPS,GLONASS,Galileo和BeiDou系統(tǒng)每天的SNR弧段數(shù)Fig.4 Number of SNR arcs used per day for GPS,GLONASS,Galileo,and BeiDou signals
圖5 四系統(tǒng)SNR序列和LSP分析結(jié)果Fig.5 SNR sequence and LSP analysis of quad-constellation
由圖5(a)四系統(tǒng)SNR 序列可知,在低高度角情況下,GPS 中S2W 的能量明顯低于其他三種信號,SNR 的最大差值在10 db 左右,GLONAS、Galileo 和BeiDou 的不同載波信號SNR 差值較小,整體趨勢較為一致。圖5右邊為四系統(tǒng)SNR序列對應(yīng)的LSP分析結(jié)果,其中波峰對應(yīng)的橫坐標(biāo)即為垂直反射距離,可以看出四個系統(tǒng)不同信號反演結(jié)果差值在0.08 m左右。GPS中S2W的SNR序列能量最低,在LSP分析中整體振幅也是最低;S2W和S2X的LSP 分析結(jié)果有明顯的次波峰,信號質(zhì)量較差,而S5X在LSP分析中只有一個較大波峰,信號質(zhì)量最好。GLONASS 和Beidou 中信號的SNR 序列有較好的一致性,其LSP分析結(jié)果振幅值比較穩(wěn)定,都只有一個明顯的波峰。Galileo中除了S1X在LSP中振幅值最低且最大波峰不明顯,S5X、S7X、S8X、S6X 四種信號LSP 結(jié)果都沒有虛假高峰,能夠較好的得到反演結(jié)果。四系統(tǒng)不同信號的LSP表現(xiàn)出不同的結(jié)果,每個信號對應(yīng)于不同的峰值,這表明同一系統(tǒng)不同類型的SNR 會導(dǎo)致精度水平不一樣的結(jié)果。
以BRST 測站為例,對去除趨勢項的SNR 序列進(jìn)行LSP分析,可以獲得每種信號的垂直反射距離的時間序列。在設(shè)定LSP分析中峰值噪聲比大于3后,4個系統(tǒng)潮位有效反演值如圖6所示。
圖6 BRST站四系統(tǒng)GNSS-MR潮位反演結(jié)果Fig.6 Sea levels at BRST estimated from quad-constellation
由圖6可知,4個系統(tǒng)每個信號通過GNSS-MR獲取的潮位反演結(jié)果都能夠描述潮位變化的趨勢,但是這些結(jié)果與驗潮站潮位測量值存在一定偏差。為了進(jìn)一步評估各個信號的反演效果,同時對海面動態(tài)變化誤差和對流層延遲誤差ΔhT進(jìn)行改正,BRST 測站每個信號反演序列、-和--ΔhT的均方根誤差統(tǒng)計結(jié)果如表1 所示,經(jīng)過兩次誤差改正后的精度提高百分比如圖7所示。
結(jié)合表1 和圖7 可知通過海面動態(tài)改正和對流層延遲改正,4 個系統(tǒng)的反演精度大約提高了50%,GPS 中S5X 反演結(jié)果優(yōu)于S1C、S2W 和S2X,精度最高;GLONASS 中S2C 和S2P 反演精度要優(yōu)于S1C 和S1P;Galileo 的S5X 反演精度最高,S8X、S7X 和S6X 較次之,S1X 精度最低;BeiDou 系統(tǒng)整體反演精度要低于其他3 個系統(tǒng),其中S7I 反演精度最低,S6I反演精度最高。
對于BRST 測站而言,經(jīng)過誤差改正后精度最好的SNR 數(shù)據(jù)分別是GPS 的S5X、GLONASS 的S2P、Galileo 的S5X 和BeiDou 的S61。分析可知隨著信號的頻率下降,反演的精度會隨之提高,這可能是因為頻率會影響天線增益模式和隨機(jī)表面粗糙度,從而影響信號的質(zhì)量(Nievinski 和Larson,2014)。
在BRST測站中,結(jié)合表1和圖6可以看出單系統(tǒng)單頻的反演結(jié)果精度低且時間分辨率差,而單系統(tǒng)多頻又會有反演結(jié)果在某一時間段堆疊的情況。為了解決上述問題,本文在窗口長度為2 h,滑動步長為10 min 的時間窗口中建立狀態(tài)方程,在穩(wěn)健回歸模型下實現(xiàn)四系統(tǒng)反演結(jié)果融合?;瑒哟翱陂L度的選擇非常關(guān)鍵:長度過短,窗口中數(shù)據(jù)量過少,不足以進(jìn)行可靠的抗差估計;窗口過長,則潮位的非線性運動造成的融合反演誤差越大。文章從精度和數(shù)據(jù)損失率(損失數(shù)據(jù)量/窗口數(shù)量)兩個方面分別比較0.5 h、1 h、1.5 h、2 h、2.5 h、3 h 窗口潮位反演差異(表2),最終確定最優(yōu)窗口長度為2 h。
表2 BRST站不同窗口長度反演精度和數(shù)據(jù)損失率對比Table 2 Comparison of inversion accuracy and data loss rate of different window lengths in BRST station
圖8(a)為四系統(tǒng)聯(lián)合反演的結(jié)果,圖中黑色實線表示驗潮站實測數(shù)據(jù),藍(lán)色的點為聯(lián)合反演結(jié)果,可以看出聯(lián)合反演的結(jié)果較好的對應(yīng)了驗潮站的實測數(shù)據(jù),同時實現(xiàn)了10 min 一個點的均勻采樣,滿足日常潮位監(jiān)測的要求。改正后精度為12.43 cm,相較于各系統(tǒng)單信號的反演結(jié)果,精度大約提高了40%—60%。圖8(b)為反演潮位與實測潮位的殘差分布,散點大都在-0.5 m—0.5 m 范圍內(nèi),在漲潮和落潮時殘差波動較大。圖8(c)是反演潮位與實測潮位的線性回歸模型,相關(guān)系數(shù)為99.56%。
圖8 BRST站多模多頻GNSS聯(lián)合潮位反演以及誤差分析Fig.8 Multi-GNSS sea level retrieval series and error analysis of BRST
為了進(jìn)一步驗證基于穩(wěn)健回歸模型融合算法的有效性,本文選取HKQT站2020年年積日為203—232 期間GPS、GLONASS、Galileo、BeiDou 系統(tǒng)的SNR 數(shù)據(jù)進(jìn)行實驗。其中信噪比數(shù)據(jù)GPS 有S1C、S2W、S2X、S5X 等4 種,GLONASS 有S1C、S1P、S2C、S2P等4種,Galileo有S1X、S5X、S7X、S8X等4種,BeiDou有S2I、S7I兩種。經(jīng)過海面動態(tài)改正和對流層延遲改正后,四系統(tǒng)整體的潮位反演精度約為25.42 cm。
在HKQT 站中,同樣從精度和數(shù)據(jù)損失率兩個方面分別比較不同窗口潮位反演差異(表3),在保證精度和數(shù)據(jù)利用率的情況下最后選擇2 h的時間窗口,滑動步長為10 min,然后在每個窗口中進(jìn)行四系統(tǒng)反演結(jié)果融合。圖9(a)為四系統(tǒng)聯(lián)合反演的結(jié)果,圖中黑色實線表示驗潮站實測數(shù)據(jù),紅色的點為聯(lián)合反演結(jié)果,聯(lián)合反演的結(jié)果可以清晰的反映潮位的起伏變化。改正后精度為7.09 cm,相較于四系統(tǒng)整體的反演結(jié)果,精度提高了約72%。圖9(b)為反演潮位與實測潮位的殘差分布,散點大都集中在0 m 左右,整體表現(xiàn)更為穩(wěn)定。圖9(c)是反演潮位與實測潮位的線性回歸模型,相關(guān)系數(shù)為99.06%。與BRST站的聯(lián)合反演結(jié)果相比,基于穩(wěn)健回歸模型的數(shù)據(jù)融合算法在HKQT 站即潮位起伏較小的區(qū)域表現(xiàn)更好,精度提高更多。雖然BRST 站獲取的信號種類多于HKQT 站,但是由于Galileo 和BeiDou系統(tǒng)存在個別信號質(zhì)量低而影響了整體反演精度。
表3 HKQT站不同窗口長度反演精度和數(shù)據(jù)損失率對比Table 3 Comparison of inversion accuracy and data loss rate of different window lengths in HKQT station
圖9 HKQT站多模多頻GNSS聯(lián)合潮位反演以及誤差分析Fig.9 Multi-GNSS sea level retrieval series and error analysis of HKQT
針對GNSS-MR 技術(shù)在實際用于潮位監(jiān)測中存在精度不高和時間分辨率低兩大問題,本文提出一種基于IGGIII 模型的穩(wěn)健回歸方法,在綜合考慮海面動態(tài)變化和對流層延遲的情況下對多模多頻GNSS-MR 反演結(jié)果進(jìn)行融合,有效改善了潮位反演結(jié)果的精度和時間分辨率。本文主要得出如下結(jié)論:
(1)在潮位反演中多模多頻數(shù)據(jù)相較于單系統(tǒng)單頻數(shù)據(jù)可以大大提高反演結(jié)果的時間分辨率,但即使經(jīng)過海面動態(tài)改正和對流層延遲改正,反演結(jié)果的精度依舊沒有得到有效提高。在穩(wěn)健回歸模型下考慮海面動態(tài)變化誤差和對流延遲,利用GNSS-MR 技術(shù)進(jìn)行潮位反演,BRST站和HKQT站反演結(jié)果RMSE分別為12.43 cm和7.09 cm,與傳統(tǒng)方法相比精度提高約40%—70%。
(2)傳統(tǒng)方法獲取潮位反演值由于SNR 數(shù)據(jù)采樣時間不規(guī)律,潮位值也無法實現(xiàn)均勻采樣。在基于穩(wěn)健回歸模型的數(shù)據(jù)融合算法中可以在特定的時間窗口內(nèi)建立方程求解潮位值,BRST 站和HKQT站均可實現(xiàn)10 min的均勻采樣,這大大提高了GNSS-MR 技術(shù)用于潮位監(jiān)測的時間分辨率,滿足了潮位監(jiān)測需求。
利用這種基于IGGIII 模型的穩(wěn)健回歸方法對潮位反演值進(jìn)行融合目前也存在一定局限性:(1)IGGIII模型也無法完全去除異常值對結(jié)果的影響,系統(tǒng)內(nèi)存在的低質(zhì)量信號會影響整體反演精度。(2)某一個時間點可能存在衛(wèi)星數(shù)量少或者信號質(zhì)量差的情況,這會導(dǎo)致這一時間段狀態(tài)方程無法求解從而導(dǎo)致潮位數(shù)據(jù)缺失。(3)反演潮位時獲得的潮位值是根據(jù)一段時間的SNR 序列獲得,一般將其與SNR 序列的平均時刻對應(yīng);而驗潮站測量結(jié)果也為短時間的平均值,因此反演值和觀測值之間會存在時間尺度上的表示性誤差。對于此類誤差,需要在未來進(jìn)一步研究。