許 煒 呂志平 董 行 鄺英才 楊凱淳 黃 輝

1 信息工程大學(xué)地理空間信息學(xué)院，鄭州市科學(xué)大道62號，450001 2 中國人民解放軍96712部隊(duì)，江西省景德鎮(zhèn)市，333300 3 哈爾濱工業(yè)大學(xué)(深圳)空間科學(xué)與應(yīng)用技術(shù)研究院，深圳市平山一路6號，518055 4 中國人民解放軍61206部隊(duì)，北京市，100042 5 中國人民解放軍96711部隊(duì)，安徽省池州市，247100

針對GNSS坐標(biāo)時間序列粗差探測與剔除的研究很多，Pernetti[1]利用DIA(detection identification adaptation)方法對GPS坐標(biāo)時間序列的粗差和中斷進(jìn)行檢測，結(jié)論顯示，早期的GPS坐標(biāo)時間序列中大約有6%的數(shù)據(jù)為粗差；Mao等[2]及李曉光等[3]采用3倍中誤差(3σ)準(zhǔn)則進(jìn)行粗差探測，但3σ準(zhǔn)則易受粗差的影響；為增強(qiáng)粗差探測的穩(wěn)健性，四分位距(interquartile range, IQR)準(zhǔn)則得到了更加廣泛的應(yīng)用，其對粗差有一定的抵御能力。IQR準(zhǔn)則也衍生出不少新的算法，如經(jīng)典的最小二乘結(jié)合IQR準(zhǔn)則的方法(LS-IQR)[4]、抗差1范數(shù)結(jié)合IQR準(zhǔn)則的方法(L1-ModIQR)[5]、數(shù)據(jù)驅(qū)動的奇異譜分析技術(shù)結(jié)合IQR準(zhǔn)則的方法(SSA-IQR)[6]及在頻域上與小波分析結(jié)合的方法(WA-IQR)[7]等；崔太岷等[8]還提出基于假設(shè)檢驗(yàn)的方法對觀測數(shù)據(jù)進(jìn)行粗差探測。近年來，隨著人工智能的廣泛應(yīng)用，甚至發(fā)展出了利用機(jī)器學(xué)習(xí)和深度學(xué)習(xí)對各種時間序列進(jìn)行異常探測的算法[9]。

本文通過檢驗(yàn)常用的3σ準(zhǔn)則和IQR準(zhǔn)則的探測效果，對三維GNSS大地坐標(biāo)時間序列進(jìn)行粗差探測，并分析三維方向的探測效果。

1 GNSS坐標(biāo)時間序列粗差的特點(diǎn)

三維GNSS坐標(biāo)時間序列的粗差是空間位置的粗差，如圖1(單位mm)所示，假設(shè)對某點(diǎn)進(jìn)行多次GNSS定位，并模擬其中4次定位結(jié)果的粗差。

圖1 三維GNSS定位

為方便說明，以二維平面坐標(biāo)為例。假設(shè)圖2(單位mm)為GNSS多次測量某點(diǎn)的定位結(jié)果，A、B、C、D為4個模擬粗差。

圖2 二維GNSS定位

以3σ準(zhǔn)則為例，如圖2所示，假設(shè)通過平面定位某點(diǎn)100次，N方向和E方向的中誤差均為1 mm，并服從正態(tài)分布，且分布于3σ之外的數(shù)據(jù)不大于1%，其中1σ區(qū)域?yàn)橹虚g淡藍(lán)色區(qū)域，3σ為外圈淺黃色區(qū)域，淺紅色正方形區(qū)域?yàn)橥ㄟ^N、E方向判定不是粗差而分布于3σ之外的區(qū)域。A、B、C、D為4個人為加入的偏離超過3σ的粗差點(diǎn)，如僅考慮N方向，A與B可被識別為粗差，但C與D將不能被正確識別；同樣，如果僅考慮E方向，B和D無法被正確探測出來。雖然2個方向的探測均不能識別出在3σ之外的D點(diǎn)，但這樣的區(qū)域占比很小，僅同時探測N、E方向可大大簡化計(jì)算量，且不損失過多的探測效率。

類似于平面位置的3σ準(zhǔn)則時間序列粗差探測，三維GNSS坐標(biāo)的粗差探測應(yīng)探測各個方向的偏差是否大于某個統(tǒng)計(jì)量指標(biāo)，需要求出模型化后三維各方向的統(tǒng)計(jì)量信息。但這種處理方法計(jì)算復(fù)雜，若簡化為僅考慮N、E、U三個相互正交方向的偏差，計(jì)算量可被大大簡化，且能覆蓋大多數(shù)情形，所以本文采用各個分量分別探測出來的時間并集作為粗差發(fā)生歷元。

為更有效地檢測GNSS坐標(biāo)時間序列中的粗差，可在時域或頻域上對時間序列進(jìn)行建模，然后針對得到的殘差序列進(jìn)行探測。由于在頻域上對時間序列進(jìn)行分析要求采樣均勻[10]，但GNSS時間序列中存在不少缺失數(shù)據(jù)，破壞了其均勻性，進(jìn)行插值則可能會引入新的粗差。所以，本文分別依據(jù)3σ準(zhǔn)則和IQR準(zhǔn)則，利用迭代最小二乘法求取參數(shù)和剔除粗差。

2 綜合N、E、U三方向的迭代最小二乘粗差探測方法

在實(shí)際應(yīng)用中，時間序列的自相關(guān)性隨觀測歷元間隔的增大而降低[11]，坐標(biāo)時間序列中存在與時間相關(guān)的有色噪聲成分，而且對于跨度較長的坐標(biāo)時間序列，其前后的年周期和半年周期振幅往往不一樣，所以一般采用移動開窗法對粗差進(jìn)行探測?？紤]到GNSS坐標(biāo)時間序列中的年周期信號最明顯[12]，所以窗口及步長均設(shè)為1 a，不足1 a的數(shù)據(jù)以離目標(biāo)數(shù)據(jù)最近的1 a數(shù)據(jù)作為探測窗口。

僅靠1次探測在很多情況下不能成功探測出所有粗差，所以在單方向上采取迭代求取參數(shù)和剔除粗差的方式。迭代的優(yōu)勢在于可以逐步剔除粗差，并消除粗差對最小二乘法估計(jì)參數(shù)的影響，最終得到正確的殘差序列。而三維粗差探測為N、E、U方向探測出的粗差并集，可大幅降低單方向漏判的概率，具體探測步驟如下：

1)輸入觀測方程和權(quán)。

2)利用最小二乘法估計(jì)出GNSS時間序列模型，求取參數(shù)和殘差。

3)按3σ準(zhǔn)則或IQR準(zhǔn)則探測和剔除N、E、U方向的粗差。

3σ準(zhǔn)則是指當(dāng)大量GNSS坐標(biāo)時間序列數(shù)據(jù)模型化后的殘差服從正態(tài)分布時，誤差分布于±3σ之內(nèi)的數(shù)據(jù)占99.73%，所以將分布于3σ之外的數(shù)據(jù)認(rèn)為是粗差造成的，其誤判比例很小。具體判別式為：

(1)

(2)

4)重復(fù)步驟1)～3)，直至探測出的粗差個數(shù)為0。

5)求取N、E、U三方向探測所得的粗差并集。

3 算例分析

3.1 模擬數(shù)據(jù)分析

對于單站單方向GNSS坐標(biāo)時間序列，一般采取式(3)進(jìn)行建模[10]：

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

(3)

式中，ti為序列的歷元，單位a；a對應(yīng)于截距，其值與預(yù)設(shè)的初始坐標(biāo)位置有關(guān)；b為站點(diǎn)坐標(biāo)線速度；c、d和e、f分別為年周期項(xiàng)和半年周期項(xiàng)系數(shù)，能夠完整描述周期函數(shù)的振幅和初始相位；gi為由各種原因引起的階躍式坐標(biāo)突變；Thj為階躍發(fā)生的開始時刻；H為階躍函數(shù)，發(fā)生突變前H值為0，發(fā)生突變后H值為1；vti為觀測噪聲，可表示成不同噪聲模型的組合。

為說明粗差探測方法的有效性并驗(yàn)證三維GNSS坐標(biāo)時間序列的粗差剔除效果，對模擬數(shù)據(jù)進(jìn)行處理。為使問題不過于復(fù)雜，又具有一定的代表性，實(shí)驗(yàn)根據(jù)北京房山(BJFS)站的線速度、周年和半周年振幅模擬了2009～2019年共10 a的三維GNSS坐標(biāo)時間序列，截距(a)在初始位置N、E、U方向坐標(biāo)均設(shè)為0，加入中誤差均為3 mm的白噪聲，且不添加階躍信號，依據(jù)式(3)的參數(shù)模擬結(jié)果如表1(單位mm)所示。

表1 模擬時間序列參數(shù)值

粗差模擬過程如下：1)將三維方向的噪聲按照勾股定理求取空間位置誤差的模長，將得到的模長放大6倍，并從中隨機(jī)選取200個大于3倍空間位置誤差模長的時間點(diǎn)作為粗差，粗差占比為5.5%。粗差的天頂距φ在0～π中均勻分布，水平方向θ在0～2π中均勻分布。2)將粗差的模長分別乘以cosφ、sinφ·cosθ及sinφ·sinθ，并加入U(xiǎn)、N、E方向選取的200個粗差時間點(diǎn)，得到各方向被粗差污染的時間序列。原序列和被粗差污染的序列如圖3所示。

圖3 模擬數(shù)據(jù)及粗差

為檢驗(yàn)2種準(zhǔn)則的檢測性能和GNSS三維位置粗差的探測效果，設(shè)計(jì)模擬實(shí)驗(yàn)對2種方案進(jìn)行對比：1)利用最小二乘結(jié)合3σ準(zhǔn)則方法(LS-3σ)進(jìn)行單方向及綜合三方向探測；2)利用最小二乘結(jié)合IQR準(zhǔn)則方法(LS-IQR)進(jìn)行單方向及綜合三方向探測。對比各方法探測的粗差個數(shù)和誤判個數(shù)。

各方案粗差探測結(jié)果如表2所示，3σ準(zhǔn)則和IQR準(zhǔn)則在各方向的探測結(jié)果如圖4所示，綜合三方向的探測結(jié)果如圖5所示。

表2 模擬GNSS坐標(biāo)時間序列粗差探測效果

圖4 3σ準(zhǔn)則與IQR準(zhǔn)則各方向粗差探測效果

圖5 3σ準(zhǔn)則與IQR準(zhǔn)則綜合三方向粗差探測誤判情況

結(jié)合表2和圖4可以看出，若僅考慮單方向，2種準(zhǔn)則均能較好地識別出偏離較大的粗差，但無法識別偏離較小的粗差，所以單方向的探測率均不高。其中，3σ準(zhǔn)則水平方向的探測率在55%左右，存在5%左右的誤判，而垂直方向的探測率在80%左右，誤判率也在5%左右；IQR準(zhǔn)則探測率較低，水平方向在45%左右，垂直方向約為75%，不過IQR準(zhǔn)則在本次實(shí)驗(yàn)中沒有出現(xiàn)誤判的情況，可認(rèn)為其誤判率顯著低于3σ準(zhǔn)則。低誤判率的主要原因?yàn)镮QR統(tǒng)計(jì)量并未進(jìn)行標(biāo)準(zhǔn)化，相當(dāng)于對粗差判定的置信區(qū)間進(jìn)行了壓縮。2種方法在垂直方向的探測效率均好于水平方向，主要原因是在粗差的模擬過程中，垂直分量粗差的占比要大于水平分量。將各方向探測得到的粗差求取并集后，3σ準(zhǔn)則可探測所有粗差的發(fā)生時刻，但誤判率也達(dá)到了13.5%；IQR準(zhǔn)則的探測率超過98%，沒有出現(xiàn)誤判。圖5中的粗差誤判和漏判歷元均是在3個方向上的，粗差與模擬正常值沒有明顯差別。綜合來看，IQR準(zhǔn)則要好于3σ準(zhǔn)則，其以稍低的探測率換取了低誤判率。

3.2 實(shí)測數(shù)據(jù)分析

根據(jù)模擬實(shí)驗(yàn)結(jié)果，實(shí)測數(shù)據(jù)采用IQR準(zhǔn)則進(jìn)行粗差探測。實(shí)測數(shù)據(jù)為SOPAC(Scrips Orbit and Permanent Array Center)利用GAMIT軟件處理的全球參考站坐標(biāo)時間序列，選取其中9個近年來數(shù)據(jù)相對完整且穩(wěn)定的站點(diǎn)作為分析對象，時間跨度為2009.001 4～2018.998 6(2009-01-01～2018-12-31)，具體歷元數(shù)因各站點(diǎn)的數(shù)據(jù)缺失程度不一而不同。

SOPAC提供4類坐標(biāo)時間序列：1)僅去除較大粗差的原始序列(raw)，即平差前對水平方向偏離大于1 000 mm、垂直方向大于3 000 mm的觀測值進(jìn)行剔除；2)移除了均值、粗差、階躍的干凈時間序列(cleaned trended)；3)在干凈時間序列的基礎(chǔ)上移除趨勢的時間序列(cleaned detrended)；4)對干凈時間序列進(jìn)行主成分分析濾波得到的過濾時間序列(filtered)。其中，cleaned trended數(shù)據(jù)的粗差剔除策略是觀測值減去計(jì)算值(O-C)，以檢測模型化后的殘差是否大于所設(shè)置的閾值，其閾值在N、E、U方向分別為25.0 mm、25.0 mm和35.0 mm。

本文對SOPAC提供的raw數(shù)據(jù)進(jìn)行粗差探測，并與提供的cleaned trended數(shù)據(jù)進(jìn)行對比來進(jìn)行驗(yàn)證。為保證模型的一致性，在進(jìn)行粗差探測前，通過查詢cleaned trended數(shù)據(jù)中的階躍信息進(jìn)行模型改正，再利用IQR準(zhǔn)則對時間序列的三方向進(jìn)行粗差探測，并求得并集，探測結(jié)果見表3。由表3可知，每個測站都存在數(shù)據(jù)缺失的情況，缺失比例在2%～32%，不滿足均勻采樣的條件，表明實(shí)測的IGS站原始序列數(shù)據(jù)不適合直接在頻域上進(jìn)行分析。利用本文方法對9個IGS站數(shù)據(jù)進(jìn)行粗差探測得出，原始坐標(biāo)時間序列中含有的粗差占比在0.1%～4.7%之間。

表3 IGS參考站粗差探測效果

為詳細(xì)說明探測效果，以BUCU站為例，對其三維位置粗差探測效果進(jìn)行分析，檢驗(yàn)探測結(jié)果的有效性，如圖6所示?？梢钥闯觯?/p>

圖6 BUCU站粗差探測效果

1)由于水平方向存在非常明顯的線性趨勢，局部特征容易被線性趨勢掩蓋，僅靠目測不易判斷其中存在的粗差。在僅考慮垂直方向的情況下，2個水平方向的明顯粗差未能被識別出來。綜合探測3個方向后，識別的粗差個數(shù)明顯增多，包括很多從圖上看起來不像粗差的值。

2)根據(jù)模擬實(shí)驗(yàn)的結(jié)果，在多個方向上同時探測出某個歷元存在粗差的比例大概在70%左右，而BUCU站的實(shí)測數(shù)據(jù)中僅有10%左右，說明存在不小比例的誤判。其原因可能是存在一些未探明的階躍信號，如2010.053 4～2010.072 6連續(xù)8 d被標(biāo)記為粗差，通過檢查N方向數(shù)據(jù)發(fā)現(xiàn)，其很有可能是某種類型的階躍，表明在cleaned trend數(shù)據(jù)中還存在沒有被探明的階躍。

3)坐標(biāo)時間序列中的噪聲并不是單純的白噪聲，其先驗(yàn)隨機(jī)模型不準(zhǔn)確，導(dǎo)致其與模擬實(shí)驗(yàn)結(jié)果存在差異。為更明顯地對粗差剔除效果進(jìn)行對比，實(shí)驗(yàn)通過最小二乘法估計(jì)出利用IQR準(zhǔn)則剔除粗差后序列的趨勢及截距，將得到的去趨勢數(shù)據(jù)與SOPAC提供的cleaned detrend數(shù)據(jù)進(jìn)行對比，結(jié)果如圖7所示。從圖7可以看出，相比于cleaned detrend數(shù)據(jù)，IQR準(zhǔn)則綜合三方向剔除粗差后的時間序列不僅剔除了孤立的異常點(diǎn)，細(xì)節(jié)特征也更加明顯。

圖7 IQR準(zhǔn)則綜合三方向剔除粗差后時間序列與cleaned detrend時間序列對比

4 結(jié) 語

本文對比了最小二乘分別結(jié)合3σ準(zhǔn)則和IQR準(zhǔn)則對三維GNSS坐標(biāo)時間序列進(jìn)行粗差探測的方法。由模擬實(shí)驗(yàn)結(jié)果可見，綜合三方向的探測結(jié)果明顯好于僅考慮單方向的結(jié)果；3σ準(zhǔn)則的探測效率高，但存在較多的誤判；IQR準(zhǔn)則探測效率略低，但沒有出現(xiàn)誤判的情況?？傮w來看，IQR準(zhǔn)則要好于3σ準(zhǔn)則。

然而實(shí)測數(shù)據(jù)情況要更為復(fù)雜，其中含有大量的階躍信號，如果未準(zhǔn)確建模，將會對參數(shù)估計(jì)產(chǎn)生很大影響，進(jìn)而影響探測效果。而在SOPAC提供的干凈時間序列中依然存在很多沒有被探明的階躍，造成了探測中粗差與階躍的混淆，因此將階躍和粗差進(jìn)行區(qū)分十分重要。

另外，GNSS坐標(biāo)時間序列中的噪聲成分復(fù)雜，其隨機(jī)模型不準(zhǔn)確會使參數(shù)估值和殘差的計(jì)算受到很大影響，進(jìn)而影響粗差探測的效果。獲取更加準(zhǔn)確的函數(shù)模型和隨機(jī)模型也有助于提高粗差探測的準(zhǔn)確性。

GNSS坐標(biāo)時間序列中粗差產(chǎn)生的原因常被籠統(tǒng)地歸結(jié)為儀器故障或測量條件的限制，對粗差產(chǎn)生的原因進(jìn)行定性和定量分析，有助于在利用GNSS觀測數(shù)據(jù)獲取坐標(biāo)時間序列階段避免或減弱粗差的影響。