GPS 坐標(biāo)時間序列粗差剔除方法比較分析

舒 穎

（中鐵第四勘察設(shè)計院集團有限公司，武漢 430063）

0 引言

全球定位系統(tǒng)（global positioning system,GPS）測量過程中會產(chǎn)生多種誤差，如電離層延遲誤差、對流層延遲誤差、多路徑效應(yīng)[1-3]等。這些誤差都已經(jīng)建立大量誤差模型來消除和減弱其對GPS 數(shù)據(jù)的影響。然而在定位測量過程中受多種因素的影響仍然會存在一些粗差，其存在會使得坐標(biāo)、速度估值存在偏差[4]，給測量結(jié)果的精度帶來很大的影響；并且GPS 速度不確定性受GPS 噪聲模型的影響比較大，如果不對噪聲模型進行準(zhǔn)確估計，會引起速度的有偏估計。因此有必要對GPS坐標(biāo)時間序列進行粗差探測并選擇準(zhǔn)確的噪聲模型，來保證獲取的GPS 坐標(biāo)時間序列的可靠性。GPS 坐標(biāo)時間序列中常用的粗差探測方法如5 倍中誤差法（5σ）、3 倍中誤差法[5]（3σ）。由于GPS 坐標(biāo)時間序列主要表現(xiàn)為有色噪聲特性，并不符合正態(tài)分布的高斯白噪聲，上述方法在粗差探測中，應(yīng)用受到一定限制[6]。為了更加有效地探測和剔除GPS 臺站位移序列的粗差，本文采用常見的幾種粗差剔除方法對GPS 坐標(biāo)時間序列進行處理并比較各種方法的優(yōu)劣性，接著對剔除粗差后的GPS 坐標(biāo)序列進行噪聲模型分析，總結(jié)粗差剔除前后時間序列的最佳噪聲模型，以獲得更加可靠的結(jié)果。

1 粗差剔除方法基本原理

1.1 3 倍中誤差法

3 倍中誤差法基本原理是假設(shè)測得的觀測數(shù)據(jù)為L（ii=1,2,???,n），相應(yīng)的殘差為v（ii=1,2,???,n），根據(jù)誤差理論，隨機誤差σ服從正態(tài)分布，σ為標(biāo)準(zhǔn)差，一般是未知的，通常用貝塞爾公式算得S代替σ，以L代替真值，其計算公式[7]為

對某個時刻的觀測值Li，若其殘差滿足大于3 倍標(biāo)準(zhǔn)差的條件，則認為觀測值Li為粗差。

1.2 中位數(shù)絕對偏差法

中位數(shù)絕對偏差法簡稱MAD 法。MAD 法最突出的特性是抵抗粗差，它不要求對原始時間序列進行插值[8]，其計算公式為

式中：Xi、Xj（i,j=1,2,???,n）為某個時刻的觀測值；MMAD為觀測值序列的中位數(shù)絕對偏差值。對正態(tài)分布的數(shù)據(jù)使用MAD 法得出的值與標(biāo)準(zhǔn)差相似，本文使用的MAD 值實際上為3×1.482 6×MAD，這樣算出的中位數(shù)絕對偏差與3 倍中誤差接近，但不等于3 倍中誤差。

1.3 四分位距法

為了更加有效地探測GPS 臺站位移序列的粗差，本文嘗試采用四分位距（interquartile range,IQR）粗差剔除方法對其進行處理[9]，首先采用最小二乘方法對位移時間序列進行線性擬合求出趨勢項，對趨勢項進行去除得到位移序列的殘差序列，接著將殘差序列按照從小到大的順序進行排列，計算殘差序列的四分位數(shù)間距，P25 稱下四分位數(shù)（Q1），P50 稱中位數(shù)（Q2），P75 稱上四分位數(shù)（Q3），下四分位數(shù)、中位數(shù)和上四分位數(shù)分別為排序后序列第25%、50%、75%的數(shù)值，將P75與P25 之差定義為四分位數(shù)間距（IQR）。具體計算方法如下：

1）上四分位數(shù)Q3的位置為(n+1) ×0.75，n表示項數(shù)；

2）下四分位數(shù)Q1的位置為(n+1) ×0.25；

3）四分位距IQR的計算，IIQR=Q3?Q1，IIQR為殘差序列的四分位數(shù)間距值；

4）計算a、b，a=(Q1?1.5IIQR)，b=(Q1+1.5IIQR)；

5）原始序列中位于（a，b）區(qū)間之外的值，則為粗差[10]。

1.4 GPS 坐標(biāo)序列精度評價指標(biāo)

本文通過計算 GPS 坐標(biāo)時間序列的均方根（root mean square,RMS）來分析共模誤差對GPS坐標(biāo)時間序列的影響程度[11]。RMS的計算公式為

式中：RRMS為序列的均方根值；vi為時間序列第i個值的殘差；N為公共歷元的數(shù)量。

1.5 GPS 最優(yōu)噪聲模型估計準(zhǔn)則

為了對噪聲模型進行更加準(zhǔn)確的估計，本文在傳統(tǒng)極大似然估計法（maximum likelihood estimation,MLE）的基礎(chǔ)上，采用赤池信息量（Akaike information criteria,AIC）和貝葉斯信息量（Bayesian information criteria,BIC）的方法進行模型判定，以獲得更加穩(wěn)健的噪聲模型估計結(jié)果[12]。其基本原理為：

式中：AAIC為樣本的赤池信息量；BBIC為樣本的貝葉斯信息量；L為對應(yīng)模型下的似然函數(shù)；n為樣本觀測的個數(shù)；k為變量個數(shù)。

當(dāng)似然函數(shù)L越大，對應(yīng)的模型越趨近于真實模型，此時AIC/BIC 值最小。當(dāng)AIC、BIC的結(jié)果不一致時，此時需要BIC 準(zhǔn)則才能獲得較為準(zhǔn)確的估計結(jié)果，即BIC值最小的模型是更加可靠的模型?；诖耍疚膶V波后坐標(biāo)時間序列的噪聲模型進行估計分析，考慮多種噪聲模型組合，探討不同條件下GPS 坐標(biāo)時間序列的最佳噪聲模型，為GPS臺站速度估計以及相關(guān)參數(shù)估計提供有力的根據(jù)。

2 實驗與結(jié)果分析

為了比較GPS 坐標(biāo)時間序列采用5 倍中誤差法、3 倍中誤差法、中位數(shù)絕對偏差法、四分位距法這4 種粗差剔除的效果，選取了國際全球衛(wèi)星導(dǎo)航系統(tǒng)（global navigation satellite system,GNSS）服務(wù)組織（International GNSS Service,IGS）分析中心、美國噴氣推進實驗室（Jet Propulsion Laboratory,JPL）通過GIPSY 軟件解算的40 個GPS基準(zhǔn)站坐標(biāo)時間序列進行分析，其坐標(biāo)序列可以通過文獻[12]給出的網(wǎng)站下載。對GPS 坐標(biāo)時間序列中的粗差采用上述方法進行剔除，通過比較粗差的剔除率和去除粗差后序列的RMS 值以獲得效果最優(yōu)的方法。

2.1 數(shù)據(jù)選取及預(yù)處理

經(jīng)過計算和挑選，本文最終在全球區(qū)域內(nèi)，選取1996—2016 年這20 a的40 個GPS基準(zhǔn)站的數(shù)據(jù)來進行粗差的探測與剔除，這些基準(zhǔn)站的名稱及所在位置如圖1 所示，站點坐標(biāo)信息如表1 所示。

圖1 GPS基準(zhǔn)站名稱及其所在位置

表1 本節(jié)實驗所用GPS基準(zhǔn)站的經(jīng)緯度

2.2 4 種粗差探測方法的比較和分析

為了較為直觀地描述粗差剔除前后的GPS 坐標(biāo)時間序列，本文列出了部分站點時間序列垂向分量上的粗差剔除情況，如圖2 所示。從圖中可以看出：GPS 時間序列中或多或少存在大幅度偏離正常值的觀測值，即我們要剔除的粗差，4 種粗差剔除方法在一定程度上能探測出并且去除這些“不正?！钡挠^測值；且由于4 種方法選取閾值的方式和大小不同，造成了粗差剔除效果上的不同。

圖2 部分測站粗差剔除前后的序列

為了定量分析4 種濾波方法去除粗差的效果，本文計算并對比4 種剔除方法的粗差剔除率（百分比），在無錯誤剔除數(shù)據(jù)的情況下，該值越大，表明剔除效果越好。表2 為用4 種粗差剔除方法所處理站點剔除率的結(jié)果。

表2 4 種粗差探測方法的剔除率比較 %

（續(xù)）

由表可以看出4 種方法效果的優(yōu)劣性：四分位距（IQR）法的效果優(yōu)于MAD 法；MAD 法優(yōu)于3 倍中誤差法；3 倍中誤差法優(yōu)于5 倍中誤差法。為了使統(tǒng)計結(jié)果更加可靠，我們進一步對40 個站去除粗差前后序列的RMS 結(jié)果進行統(tǒng)計分析。表3 給出了使用不同方法去除粗差前后序列RMS的最大值（Max）、最小值（Min）、均值（Mean）。

表3 4 種方法去除粗差前后站點序列的RMS 統(tǒng)計分析結(jié)果 mm

從表可知，就RMS 均值而言，用4 種方法去除粗差后序列的 RMS的均值在東（E）方向、北（N）方向上相對于原始序列有小幅度降低，在垂直（U）方向分量上有大幅度減小，即去除粗差后序列的均方根明顯變小。

用4 種方法去除粗差后序列的RMS的最大值、最小值在每個方向上相對于原始序列均有大幅度降低，且以四分位距法的效果最佳。

對40 個GPS基準(zhǔn)站序列的RMS 最大值、最小值進行統(tǒng)計：東方向濾波前最大值為87.52 mm，最小值為 10.89 mm；北方向濾波前最大值為49.30 mm，最小值為2.34 mm；垂向分量濾波前最大值為215.81 mm，最小值為6.43 mm。

四分位距法剔除粗差后序列的RMS 在東方向最大值為37.12 mm，最小值為2.95 mm；北方向最大值為42.53 mm，最小值為2.25 mm；垂向分量最大值為16.57 mm，最小值為4.75 mm。

剔除粗差后的序列RMS 平均值比原始序列的RMS 平均值在E 方向分別下降了56.1%、62.3%、61.1%、63.7%，在N 方向上分別下降了15.6%、20.3%、18.0%、23.7%，在U 方向分別下降了57.9%、69.0%、76.7%、82.4%。

綜上所述，GPS基準(zhǔn)站的時間序列經(jīng)粗差去除后序列的觀測值同真值之間的偏差明顯減小，即殘余坐標(biāo)時間序列的不確定性和可靠性得到了較大的提高。

2.3 粗差對噪聲模型估計的影響

為了分析粗差對噪聲模型的影響，實驗數(shù)據(jù)依然采用表1 所述GPS 站點，對粗差剔除后的坐標(biāo)序列進行相應(yīng)的噪聲模型估計分析。

噪聲數(shù)據(jù)處理依然采用極大似然估計進行，采用白噪聲（white noise,WN）、冪律噪聲（power law noise,PL）、閃爍噪聲/白噪聲（flicker noise/white noise,FN/WN）、閃爍噪聲/隨機游走噪聲/白噪聲（flicker noise/random walk noise/white noise,FN/RW/WN）4 種常用噪聲模型進行統(tǒng)計，建立其最優(yōu)噪聲模型。表4 為圖1 所示中區(qū)域內(nèi)站點采用四分位距法去除粗差前后3 個方向坐標(biāo)分量的最優(yōu)噪聲模型結(jié)果。

表4 4 種方法去除粗差前后最佳噪聲模型站點數(shù)統(tǒng)計結(jié)果 mm

由表可以看出，GPS 坐標(biāo)時間序列并非都是單一的噪聲模型，而是呈現(xiàn)出多樣性。GPS 坐標(biāo)時間序列在東方向主要呈現(xiàn)出白噪聲模型，占總噪聲模型個數(shù)40.0%，冪律噪聲模型占20%，閃爍噪聲/白噪聲模型占35%，閃爍噪聲/隨機游走噪聲/白噪聲模型約占5%。在北方向主要表現(xiàn)為冪律噪聲模型占20%，閃爍噪聲/白噪聲模型占52.5%，閃爍噪聲/隨機游走噪聲/白噪聲模型約占27.5%。在垂向分量上主要表現(xiàn)為白噪聲模型占2.5%，冪律噪聲模型占32.5%，閃爍噪聲/白噪聲模型占50%，閃爍噪聲/隨機游走噪聲/白噪聲模型約占15%?？梢钥闯鯣PS 坐標(biāo)時間序列在北方向與垂向分量上的最優(yōu)噪聲模型符合度較高。

除此之外，從表4 可知，經(jīng)粗差去除后，部分站點（約占總數(shù)的50%）3 個坐標(biāo)分量的噪聲模型發(fā)生了改變；且在粗差去除之前，東方向上的噪聲模型主要表現(xiàn)為白噪聲模型和閃爍噪聲/白噪聲模型，改正后主要呈現(xiàn)出閃爍噪聲/白噪聲型，70%以上的站點噪聲模型發(fā)生了改變；且改正后少部分站點（約占總數(shù)的12.5%）的最佳噪聲模型為“閃爍噪聲/隨機游走噪聲/白噪聲”模型。而根據(jù)現(xiàn)有研究表明，“閃爍噪聲/隨機游走噪聲/白噪聲”模型中的隨機游走噪聲主要是來自于觀測墩的不穩(wěn)定性。這表明經(jīng)過粗差改正之后，GPS 坐標(biāo)序列中噪聲的長周期分量（如隨機游走噪聲）變得顯著，且大跨度的時間序列為探測低頻噪聲的存在提供了條件，尤其是隨機游走噪聲振幅較小時不能被準(zhǔn)確地探測出來，從而被忽略，導(dǎo)致不能獲得最優(yōu)的噪聲模型。在垂直方向，去除粗差后主要表現(xiàn)為閃爍噪聲/白噪聲型模型，部分站最佳模型為冪律噪聲模型，改正后主要呈現(xiàn)出閃爍噪聲/白噪聲模型。上述站點經(jīng)粗差剔除之后，時間序列的最佳噪聲模型發(fā)生了變化，表明粗差對噪聲模型的影響較大，因此有必要對觀測序列中粗差對噪聲模型的影響進行進一步的研究。

3 結(jié)束語

本文以時間序列跨度為20 a的GPS基準(zhǔn)站數(shù)據(jù)為基礎(chǔ)，利用5 倍中誤差法、3 倍中誤差法、MAD 法、四分位距法4 種方法進行了GPS基準(zhǔn)站數(shù)據(jù)粗差的探測和剔除，并對結(jié)果精度進行了比較。結(jié)果表明：4 種粗差探測和剔除方法都能檢測出時間序列中的粗差，且以四分位距法的效果最為明顯。此外，對GPS基準(zhǔn)站噪聲模型的建立及粗差對其的影響進行探討，結(jié)果表明去除粗差之后的時間序列的最佳噪聲模型發(fā)生了改變，且 GPS 站坐標(biāo)序列噪聲模型呈現(xiàn)出多樣性并存在個體差異。