擬穩平差法在GNSS時間序列修復中的研究與應用

齊 聰 黨亞民 楊 強

1 山東科技大學測繪與空間信息學院,青島市前灣港路579號,266590 2 中國測繪科學研究院,北京市蓮花池西路28號,100036

隨著全球導航衛星系統GNSS的快速發展,許多國家和組織在全球建立了大量GNSS觀測站[1],GNSS觀測站的長時間坐標時間序列可為板塊運動、地震以及速度場、形變場等地球動力學研究提供準確數據。由于觀測站在觀測過程中會出現各種不確定因素,使得時間序列出現數據缺失[2],因此需要對缺失的時間序列進行修復。

目前比較常見的插值方法有拉格朗日法、三次樣條法和正交多項式法[3]:武艷強等[4]基于三次樣條插值法提出多點三次樣條法,此方法能夠在連續缺失數據時提高插值精度;Dong等[5]基于坐標序列的缺失情況,使用不同的插值方法,采用周圍站點的坐標平均值作內插處理;田慧等[6]采用二階正交多項式擬合最小二乘曲線法進行插值,能夠較好地反映原始時序的運動趨勢。但上述方法大多基于數學模型[7],沒有考慮到實際的變化因素。近年來,國內外許多學者提出了更符合實際的方法:Carrizosa[8]使用優化后的變鄰域搜索算法(variable neighborhood search, VNS)對時間序列缺失值進行修復,保留了原始數據的統計學特性;王方超等[9]將調整最大似然法(regularized EM algorithm, RegEM)引入到GNSS時間序列插值中,該方法不依靠數學模型和先驗信息,而是依賴數據特性進行插值。

通過同步觀測衛星的方式在GNSS衛星定位中選取任意2個測站,可以得到其三維坐標和一條基線向量,通過建立基線向量網對待求測站的坐標進行平差計算。考慮到待求點與周圍測站具有相同的變化趨勢,聯合待求點及其周圍區域已知測站構建GNSS網進行擬穩平差計算,得到待求點坐標,并與拉格朗日法和三次樣條法進行對比,驗證其可行性。

1 基本原理和精度評定指標

1.1 擬穩平差法基本原理

由于存在各種誤差,測量過程中數據可能有一定程度的波動,會對平差結果產生影響。本文提出一種基于間接平差法,根據相對穩定的測站坐標時間序列變化預測待求點坐標變化的時間序列分析方法。其平差準則為:

(1)

式中,V為觀測誤差,P為權向量,X為平差網中擬穩點的未知參數。

在公共歷元數據中選出與待求點總體運動趨勢相近的測站(即擬穩點),將上述測站作為基準站,以基準站第n天的三維坐標和待求點的近似三維坐標作為初始值,設待求點第(n+1)天的坐標為參數,可以表示為:

(2)

(3)

得到基線的誤差方程為:

(4)

當GNSS觀測網中有m個待求點、n條基線向量時,GNSS基線網的誤差方程為:

(5)

(6)

式中,D為基線向量協方差陣。令NBB=BTPB,W=BTPl,使用間接平差法組建法方程:

(7)

(8)

式中,Q為NBB的逆矩陣。

1.2 精度指標

通過人為打斷時間序列的方式模擬缺失的時間序列,并使用平均絕對誤差MAE、均方根誤差RMSE以及Pearson相關系數3種指標對插值結果進行分析。MAE與RMSE用來衡量插值與實際觀測值之間的偏差,MAE為L1范數,RMSE為L2范數,次數越高對于異常值越敏感。Pearson相關系數用來衡量2個數據間的相關程度,系數取值為-1～1。

平均絕對誤差MAE的計算公式為:

(9)

均方根誤差RMSE的計算公式為:

(10)

Pearson相關系數的計算公式為:

(11)

2 模擬數據實驗分析

為體現插值方法的客觀性,采用中國地殼運動觀測網絡的真實數據進行模擬實驗。選擇我國東南沿海地區的FJDH、FJYX、FJZP、FJLC四個測站為研究對象,與周圍相關性較強的測站組成GNSS相對觀測網。測區內各測站2021年的數據完整情況如圖1所示。

圖1 數據完整情況Fig.1 Data integrity

由圖1可見,各測站在2021年的觀測數據完整率均在70%以上,大部分測站數據缺失率約為10%,可以保證每個歷元至少有一半以上的測站用于組建GNSS網進行平差計算。

為模擬缺失數據,將2021年剔除粗差后的FJDH、FJYX、FJZP、FJLC觀測站的時間序列人為打斷成5 d、30 d、90 d的數據。其中FJDH站刪除doy16～20、doy31～60、doy91～180的數據,FJYX站刪除doy16～20、doy31～60、doy71～160的數據,FJZP站刪除doy91～180、doy251～255、doy301～330的數據,FJLC站刪除doy31～35、doy91～120、doy151～240的數據。分別使用拉格朗日法、三次樣條法、擬穩平差法進行插值處理,基于一致性和隨機性進行對比分析。在使用拉格朗日法與三次樣條法時,根據缺失天數選擇缺失數據前后相同區域的數據,平均分成3份后選出6～8 d的數據作為樣本構建函數進行插值。采取擬穩平差法計算時,使用含有待求點缺失序列的測站進行計算,可以分段進行插值。由于篇幅有限,模擬實驗僅在高程方向進行。原始時間序列與3種插值方法的結果如圖2所示。

圖2 4個觀測站插值結果Fig.2 Interpolation results of four observation stations

由圖2可見,在缺失5 d數據的情況下,拉格朗日法、三次樣條法和擬穩平差法的插值結果與原始序列之間不存在偏移,表現穩定;在缺失30 d數據的情況下,拉格朗日法和三次樣條法無法反映原始序列的趨勢,會產生凸起并存在較大偏移,而擬穩平差法的插值結果與原始序列一致;在缺失90 d數據的情況下,3種方法的插值結果與缺失30 d數據的結果相似,但擬穩平差法的一致性更明顯。由此可見,擬穩平差法的優勢在于能夠利用周圍測站的信息估算待求點的運動狀態,從而減小數據缺失對插值結果的影響。而基于多項式的線性插值方法只能保證時間序列的連續性和光滑性,無法保證與原始序列的關聯性。因此在數據缺失較多時,可以使用擬穩平差法進行插值修復。

為了更直觀地體現3種插值方法的效果,將插值部分進行放大處理,如圖3所示。在數據缺失較少的情況下,3種方法的插值結果與原始數據之間的誤差較小,其中拉格朗日法與三次樣條法的結果基本一致。在缺失30 d、90 d數據的情況下,拉格朗日法與三次樣條法的插值結果與原始序列偏差較大,最高達23 mm和17 mm,而擬穩平差法單日最大誤差為13 mm,說明擬穩平差法的一致性更好。

表1(單位mm)為3種插值方法的精度指標。從插值方法的精度看,在缺失5 d數據的情況下,3種方法的MAE與RMSE差距不大,說明3種方法在缺失少量數據時插值效果相當;在缺失大量數據的情況下,3種方法的修復效果存在顯著差異,擬穩平差法的MAE與RMSE指標誤差最小,在6 mm以內,相比其他2種方法精度分別提高10%和60%。從相關系數來看,在缺失5 d數據的情況下,3種方法在原始序列波動較大的FJLC、FJZP站的相關系數較小,拉格朗日法與三次樣條法在時間序列較為平緩的FJDH、FJYX站的相關性略高于擬穩平差法;在缺失30 d和90 d數據的情況下,擬穩平差法的相關系數要遠遠大于其他2種插值方法,說明擬穩平差法與原始序列的運動趨勢具有較高的相似性。

3 擬穩平差算法仿真實驗

首先使用擬穩平差法補齊測區內2021年測站的時間序列,得到完整的時間序列;然后選擇FJZP站作為待求點,使用擬穩平差法對補齊后的觀測站進行仿真實驗;最后將仿真值與真實值進行比較,評估擬穩平差法在時間序列仿真中的可信度。

由圖4可知,通過擬穩平差法仿真出的序列與真實序列基本重疊,且3個方向上的運動趨勢基本一致。將每個歷元的仿真值與真實值作差得到FJZP站2021年的時間序列殘差,如圖5所示。由圖可見,3個方向上仿真值與真實值的殘差均在0 mm附近波動,E、N方向上的時間序列殘差在10 mm以內,U方向上的時間序列殘差在20 mm以內。

圖4 仿真序列與真實序列對比Fig.4 Comparison between simulated sequence and real sequence

圖5 仿真序列與真實序列殘差Fig.5 Residual of simulated sequence and real sequence

由表2(單位mm)可見,N方向上的誤差最小,約為2 mm;E方向次之,誤差約為3.5 mm;U方向誤差最大,約為7～8 mm。3個方向上的皮爾遜相關系數達0.7～0.8,表明仿真值與真實值之間的相關性較高。從出現異常值的情況看,3個方向上的RMSE與MAE相差不大,說明出現異常值的概率非常小。由此可知,通過擬穩平差法得到的仿真數據與真實數據基本同步,在精度上能滿足測量要求,可以用作更進一步的研究。

表2 仿真值與真實值精度分析

4 結 語

1)當數據缺失較少時,3種方法的插值結果都很理想,與真實值的偏差較小,且對于較為平穩的時間序列,拉格朗日法和三次樣條法與真實值的相關性更高;在缺失30 d和90 d數據的情況下,拉格朗日法與三次樣條法出現凸起,與真實值出現較大幅度的偏移,而擬穩平差法可以根據周圍測站的運動趨勢較好地插值出待求點的時間序列。

2)使用擬穩平差法得到待求點完整的時間序列,計算得到仿真值與真實值的殘差,E和N方向上的殘差約為3 mm,U方向上的殘差約為7～8 mm,且仿真值與真實值的相關性較高。

擬穩平差法在時間序列空缺插值中具有較高的準確性和可靠性,能夠有效補齊缺失數據,并能在一定程度上模擬真實觀測值,提供高精度的數據基礎。