基于X-11-ARIMA 模型在GNSS 定位數(shù)據(jù)后處理的應(yīng)用

匡宇龍，雷孟飛

( 湖南聯(lián)智科技股份有限公司, 長沙 410200 )

0 引 言

導(dǎo)航定位在測繪和監(jiān)測領(lǐng)域受到廣泛應(yīng)用，其全天候提供位置服務(wù)的特性讓測繪和監(jiān)測工作能夠?qū)崿F(xiàn)24 h 實時監(jiān)控位置信息，獲取長時間的位置變化信息或者通過長時間靜態(tài)定位獲取高精度的測繪點位坐標. 在理想情況下，所選取的靜態(tài)監(jiān)測點位不存在位移或者位移可以忽略. 但是在地災(zāi)監(jiān)測的點位，比如邊坡及橋梁等區(qū)域，因其通行車輛的特點使得震動不可避免，如何區(qū)分正常波動和災(zāi)害發(fā)生導(dǎo)致的位移將是導(dǎo)航定位在地災(zāi)監(jiān)測領(lǐng)域中的關(guān)鍵部分[1].

在導(dǎo)航定位后處理算法領(lǐng)域中，可行的數(shù)據(jù)處理步驟是：首先通過差分定位獲得基線向量改正值；然后通過卡爾曼濾波或其改進算法對坐標參數(shù)進行濾波平滑處理，能夠?qū)崿F(xiàn)高精度實時數(shù)據(jù)輸出，在監(jiān)測中能夠滿足探測位移的發(fā)生，并且能夠剔除粗差等干擾[2]. 但是以年為單位的數(shù)據(jù)記錄表明，單純通過一次濾波，或單純使用基于原始數(shù)據(jù)的濾波處理雖然能夠剔除粗差和異常位移的干擾，但切實發(fā)生的長期性波動如車輛通行不同帶來的震動幅度變化、潮汐變化、電離層變化、對流層變化、衛(wèi)星星座分布變化及可能存在的地殼變化帶來的干擾則不屬于“誤差”，而是切實發(fā)生的位移，這些變化難以通過濾波算法進行剔除. 而對該數(shù)據(jù)進行時間序列分析則能夠針對這些長期性變化的波動干擾進行分離和處理，并獲得在長時間尺度上提高定位精度的效果.

時間序列分析被廣泛應(yīng)用于數(shù)據(jù)流的分析和預(yù)測[3]. 在各種時間序列分析方法中，差分自回歸移動平均模型(ARIMA)模型則是由Box 與Jenkins 提出的一種著名時間序列預(yù)測模型. 游賢菲等[4]. 利用X-11-ARIMA 模型對手足口疫病發(fā)病數(shù)預(yù)測和張立欣等[5]利用X-11-ARIMA 模型鐵路貨運量進行分析的應(yīng)用. 而在GNSS 領(lǐng)域的運用，文獻[6-8]利用ARMA 模型改善最小二乘解算觀測量數(shù)據(jù)，文獻[8-9]利用ARMA 模型分析預(yù)測大氣電離層電子總量(TEC). 關(guān)于ARIMA 和季節(jié)性拆分方法的運用，有HILLMER S C[9]和JUNAIDI[10]以及PENG Y 等[11]運用ARIMA 及數(shù)據(jù)拆分進行時間序列分析. Xin J Z 等[12]針對橋梁結(jié)構(gòu)監(jiān)測中衛(wèi)星定位的應(yīng)用，采用了卡爾曼濾波結(jié)合ARIMA 模型的方法預(yù)測橋梁形變，Li Q S 等[13]利用ARIMA-GRACH 模型探測周跳，JANPAULE I 等[14]分析了十年跨度以日為間隔的定位解算數(shù)據(jù)，指出地磁暴對全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)數(shù)據(jù)影響的時間跨度為一天.

文中綜合分析了國內(nèi)外研究情況，結(jié)合衛(wèi)星定位監(jiān)測數(shù)據(jù)的特點，采用時間序列分析的方式拆分監(jiān)測數(shù)據(jù). 通過ARIMA 算法預(yù)測補償平滑算法數(shù)據(jù)缺失，得到原始的完整趨勢項數(shù)據(jù). 通過實驗對比分離周期性波動和誤差后的數(shù)據(jù)改善情況.

1 數(shù)據(jù)處理方案

1.1 基于季節(jié)結(jié)構(gòu)的分解方法

在分離出趨勢因素Tt后，將乘法模型轉(zhuǎn)換可得Xt/Tt=St·It，即將原序列除去趨勢項后便得到季節(jié)因素和不規(guī)則因素項. 顯然繼續(xù)對St·It進行移動平均，便能夠剔除不規(guī)則因素項而分離出季節(jié)因素項.不同的平均項數(shù)的選擇將會對剔除結(jié)果產(chǎn)生顯著的影響，當(dāng)項數(shù)越大則剔除不規(guī)則因素的效果越好，但是越大的項數(shù)意味著損失更多的信息，尤其以尾項的數(shù)據(jù)損失最為嚴重. 所以移動項數(shù)既不能過大，也不宜過小. 在X-11 方案中，一般選擇(3×3)、(3×5)或者(3×9)項移動平均法，視原始序列的不規(guī)則因素大小和分析目的而定. 由于在剔除不規(guī)則因素之前，已經(jīng)通過十二項移動平均法分離趨勢項，不規(guī)則因素經(jīng)過兩次移動平均處理幾乎趨近于零，所以一般不用擔(dān)心存在不規(guī)則因素消除得不夠徹底的問題. 對于一組時間序列的原始數(shù)據(jù)進行分析需要循環(huán)多次進行分解運算，即重復(fù)上述的分離趨勢和剔除不規(guī)則項運算，在多次迭代后得到被平滑和拆分的數(shù)據(jù)，以進行后續(xù)的分析預(yù)測工作.

1.2 季節(jié)性ARIMA 預(yù)測模型

時間序列分析模型ARIMA 在GNSS 領(lǐng)域的應(yīng)用按照其參與的數(shù)據(jù)運算環(huán)節(jié)，參考GNSS 于INS融合的類型，將其分為松耦合與緊耦合. 文中提出的利用X-11 方案的ARIMA 模型進行季節(jié)性分析后處理平滑，因其不參與原始數(shù)據(jù)定位解算而定義為松耦合. 在文獻[6-7]的研究中， 利用ARMA 模型進行隨機模型參數(shù)估計，實時修正GNSS 觀測數(shù)據(jù). ANSARI K 等[8]利用ARMA 模型分析TEC，研究電離層延遲與衛(wèi)星角度關(guān)系，這類運用將ARMA 模型用于協(xié)助原始數(shù)據(jù)解算獲取坐標參數(shù)，可認為是緊耦合的模式. 從國內(nèi)外研究可知，GNSS 定位的確存在一定的周期性波動，無論該波動是來自于電離層周期性變化還是其他因素. 于是可以推論基于一種時間序列分析的數(shù)據(jù)平滑后處理在原理上是可行有效的.

ARIMA 模型建模需要確定兩組參數(shù)，一組是自回歸移動參數(shù) (p,d,q) ，另一組是季節(jié)性自回歸移動參數(shù) (P,D,Q)S.

通常季節(jié)性ARIMA 模型被寫作SARIMA(p,d,q)(P,D,Q)S，S表示Seasonal，即季節(jié)項數(shù)據(jù)序列，其數(shù)學(xué)算式可寫作如式(10)形式：

同位素示蹤法則是一種以同位素作為示蹤物質(zhì)，對研究對象的特征和行為進行示蹤觀察的信息獲取方法，分為兩類： 放射性同位素示蹤法和穩(wěn)定同位素示蹤法。

式中：p為自回歸階數(shù)；d為差分階數(shù)；q為移動平均階數(shù)；P為季節(jié)性自回歸階數(shù)；D為季節(jié)性差分階數(shù)；Q為季節(jié)性移動平均階數(shù).

式(10)的符號含義如式(11)～(16)所示：

1.3 數(shù)據(jù)處理

對北斗定位數(shù)據(jù)進行時間序列分析，將定位數(shù)據(jù)中的趨勢項和周期項以及不規(guī)則項通過X-11 方法進行拆分. 其中趨勢項數(shù)據(jù)可以視為剝離干擾和周期性波動的穩(wěn)定平滑數(shù)據(jù). 但是由于X-11 固有特點導(dǎo)致首尾各缺失6 歷元數(shù)據(jù)，無法滿足監(jiān)測的實時性需求.

本算法包括三個部分，首先對原始數(shù)據(jù)進行拆分，得到原始數(shù)據(jù)的趨勢項數(shù)據(jù)、季節(jié)項數(shù)據(jù)和不規(guī)則項數(shù)據(jù)；接下來通過ARIMA 算法對趨勢項和不規(guī)則項數(shù)據(jù)進行建模，并給出12 個歷元的預(yù)測值；最后將預(yù)測得到的趨勢項和不規(guī)則項與季節(jié)項數(shù)據(jù)進行結(jié)合，季節(jié)項相應(yīng)的12 個歷元數(shù)值通過季節(jié)項自身周期性推導(dǎo). 趨勢項、不規(guī)則項以及季節(jié)項重組的數(shù)據(jù)再次進行拆分運算，得到缺失的6 個歷元趨勢項.對首部缺失的6 個歷元同樣可以通過該方法進行補充. 由此可以得到原始數(shù)據(jù)完整歷元的趨勢項數(shù)據(jù)，達到剝離干擾并獲得貼近真實位移的完整數(shù)據(jù). 具體流程如下.

1.3.1 數(shù)據(jù)拆分

針對北斗基線向量殘差觀測數(shù)據(jù)進行X-11 拆分，由于數(shù)據(jù)過零點且有負值，不適合使用乘法模型，故而采用Xt=Tt+St+It加法模型.

X-11 拆分算法采用(12×2)的模型. 進行兩輪運算，首輪運算窗口寬度為12 個歷元，從頭到尾遍歷原始數(shù)據(jù). 按式(3)獲得新時間序列，排列可記為6.5、7.5、···、[(n–6)+0.5].

接著對新數(shù)列進行第二輪運算，按式(5)得到序列T：7、8、···、n–6，即一次拆分的趨勢項數(shù)據(jù)序列，可以看到首尾各有6 個歷元的數(shù)據(jù)缺失.

1.3.2 數(shù)據(jù)預(yù)測

在前一步驟中通過拆分運算依次獲得的T趨勢項數(shù)據(jù)序列、S季節(jié)項數(shù)據(jù)序列以及I不規(guī)則項數(shù)據(jù)序列中，均存在由于拆分算法特點而導(dǎo)致的數(shù)據(jù)缺失問題. 由于北斗觀測數(shù)據(jù)并非真正意義上的周期性時間序列，所以傳統(tǒng)時間序列分析中選擇滯后歷元替代缺失歷元的方法并不適用. 為了補充該部分數(shù)據(jù)，同時保證數(shù)據(jù)準確性，采用數(shù)據(jù)預(yù)測的方法對缺失數(shù)據(jù)進行填充.

以趨勢項T為例，其數(shù)據(jù)序列按對應(yīng)歷元記為：7、8、···、n–6. 對序列T進行ACF 和PACF 分析，得到ARMA(p,d,q)模型的對應(yīng)參數(shù)(3,1,3). 結(jié)合數(shù)據(jù)T，通過模型ARMA(3,1,3)對前12 歷元和后12 歷元進行預(yù)測，補足的T數(shù)據(jù)序列可記為：–5、···、0、1、···、6、7、8、···、n–6、n–5、···、n、n+1、···、n+6. 其中–5、···、0、1、···、6 為前向補足數(shù)據(jù)，n–5、···、n、n+1、···、n+6 為后向補足數(shù)據(jù).

選取補足數(shù)據(jù)T中序列為：1～n的數(shù)據(jù)序列與原始序列X進行差分運算，便能獲得序列為1～n的混合序列Y. 隨后Y序列拆分得到季節(jié)項序列S同樣首尾缺失6 歷元數(shù)據(jù)，由于季節(jié)項序列S具有周期性，便能推導(dǎo)出相應(yīng)歷元的數(shù)值，從而得到完整序列，進而與序列Y差分得到完整長度的不規(guī)則項序列I. 由于I僅通過差分剝離季節(jié)項序列S得到，所以沒有數(shù)據(jù)缺失.

1.3.3 數(shù)據(jù)重組與再拆分

經(jīng)過上述步驟得到了與原始數(shù)據(jù)序列長度相同的趨勢項序列T、季節(jié)項序列S和不規(guī)則序列I，但由于趨勢項T僅經(jīng)過一次預(yù)測得到缺失的6 歷元數(shù)據(jù)，其可靠性相對而言不夠充分. 所以需要對三組序列進一步延伸，得到未來6 歷元的數(shù)據(jù). 由于在上一步驟中序列T已經(jīng)獲得了后向12 歷元的數(shù)據(jù)，分為6 歷元補足數(shù)據(jù)，6 歷元預(yù)測數(shù)據(jù). 季節(jié)項序列S按其周期性推導(dǎo)6 歷元預(yù)測數(shù)據(jù)，而不規(guī)則序列I則通過ARMA 預(yù)測得到相應(yīng)的6 歷元預(yù)測數(shù)據(jù).

將得到的三組數(shù)據(jù)進行疊加還原為原始數(shù)據(jù)類型，隨后對“還原的”原始數(shù)據(jù)進行拆分運算，得到缺失的6 個預(yù)測歷元的趨勢項序列，易見該數(shù)據(jù)長度貼合真實數(shù)據(jù)長度，從而達到了獲取當(dāng)期歷元趨勢項序列數(shù)據(jù)的目的. 由于趨勢項數(shù)據(jù)剝離了其他波動，可認為該數(shù)據(jù)為貼近真實位移的平滑數(shù)據(jù).

經(jīng)實驗分析，拆分后預(yù)測能夠?qū)?shù)據(jù)補充到滿足拆分后處理得到當(dāng)下歷元數(shù)據(jù). 與直接ARIMA 預(yù)測對比，精度有一定程度提升.

2 實驗分析

2.1 數(shù)據(jù)獲取

本文實驗通過基于ublox-f9p 芯片研制的接收機所接收的GNSS 數(shù)據(jù)進行測試，接收機接收北斗B1、B2、B3；GPS L1、L2 信號，采用單基站實時動態(tài)(RTK)工作模式，數(shù)據(jù)經(jīng)過自適應(yīng)卡爾曼濾波處理. 選取3 個項目上的歷史數(shù)據(jù)進行對比測試，如表1 所示，均是經(jīng)過自適應(yīng)卡爾曼濾波[14]后的長時間跨度數(shù)據(jù)樣本. 實驗分為兩大部分，第一部分為X-11 季節(jié)性拆分，其原理是加權(quán)滑動平均濾波. 展示拆分季節(jié)項和不規(guī)則項后的趨勢曲線相較于僅經(jīng)過卡爾曼濾波的數(shù)據(jù). 第二部分為ARIMA 預(yù)測數(shù)據(jù)，包括直接對原始數(shù)據(jù)進行ARIMA 預(yù)測，以及將原始數(shù)據(jù)進行季節(jié)性拆分后，分別對趨勢項和不規(guī)則項進行ARIMA 預(yù)測，并將預(yù)測結(jié)果結(jié)合季節(jié)項還原為原類型數(shù)據(jù).

表1 數(shù)據(jù)樣本分類

2.2 X-11 拆分

如圖1 所示，Y軸為位移刻度，反映的是基線向量在X方向的坐標殘差，單位為mm，X軸為歷元，歷元間隔為1 h. 圖中藍色曲線為只經(jīng)過卡爾曼濾波處理的衛(wèi)星監(jiān)測數(shù)據(jù)，紅色曲線為在卡爾曼濾波基礎(chǔ)之上利用X-11 拆分后得到的趨勢項. 由于長時間跨度，圖1 中藍色曲線所代表的卡爾曼濾波本身每個點間隔為1 h. 而藍色曲線的“粗差”實際上跨度為1～2 個點，即1～2 h，這種情況可以視為動態(tài)定位數(shù)據(jù)本身反映的“位移”而非誤差所造成的. 無論引起該誤差的原因何在，如此長時間跨度的位置誤差無法用短窗口的實時平滑算法處理，如果并非真實位移所引起，則數(shù)值會在該類誤差源影響消失后回到正常水平. 利用這點便能統(tǒng)計出誤差源影響周期，并能夠剔除該類長跨度周期性的誤差干擾.

圖1 數(shù)據(jù)“波動”的時間跨度

經(jīng)過卡爾曼濾波數(shù)據(jù)的方差為4.733，X-11 濾去卡爾曼濾波的季節(jié)項和不規(guī)則項后數(shù)據(jù)的方差為2.718. 整體精度提升了42.6%，應(yīng)對存在的“粗差”也能進行較為明顯的消除. 但是由于自回歸及滑動平均算法的固有特點，即生成一個平滑數(shù)據(jù)點需要該點之后6 個數(shù)據(jù)點的支持(以X-11(12×2)方法為例)，如圖2 所示(假設(shè)平滑1 000 歷元)，反映到實際情況就是6 h 的延遲，這在地災(zāi)監(jiān)控中是無法使用的. 而為了嘗試解決該算法固有的滯后性問題，運用ARIMA算法的預(yù)測能力，基于現(xiàn)有數(shù)據(jù)給出需要平滑的歷元之后6 歷元的數(shù)據(jù)，以供數(shù)據(jù)平滑算法保證實時性.

圖2 算法原理

2.3 ARIMA 預(yù)測分析

按數(shù)據(jù)不同構(gòu)建不同的擬合模型，原始數(shù)據(jù)均是SARIMA 模型，即季節(jié)性ARIMA 模型，季節(jié)周期設(shè)為13，其根據(jù)是對原始數(shù)據(jù)進行自相關(guān)函數(shù)ACF分析，由圖3 中所示，當(dāng)滯后歷元設(shè)置為13 時自相關(guān)到達一個高峰，隨后自相關(guān)程度下降，即該時間序列約存在13 歷元的周期，判斷周期項數(shù)值設(shè)為13 較為合適.

圖3 坐標殘差數(shù)據(jù)ACF 分析

對趨勢項和不規(guī)則項的建模均為ARIMA 模型，由于原始數(shù)據(jù)存在零點與負值數(shù)據(jù)，所以分解采用的方法為相加模型(Additive Decompose). 顯然，由于剔除了季節(jié)項，所以對趨勢項和不規(guī)則項的建模不需要引入季節(jié)性因素，只需要一般ARIMA 算法進行模型的構(gòu)建. 構(gòu)建的參數(shù)如表2 所示. 由于季節(jié)項數(shù)據(jù)屬于周期性數(shù)據(jù)，所以不進行預(yù)測，直接按歷元所處周期進行數(shù)據(jù)獲取.

表2 預(yù)測模型

趨勢項以及不規(guī)則項的預(yù)測效果如圖4(a)所示.藍色曲線為原始數(shù)據(jù)拆分而來的各項數(shù)據(jù)，包括Trend 趨勢項和Resid 不規(guī)則項，紅色曲線為模型擬合后進行預(yù)測給出的數(shù)據(jù).

將上述兩項數(shù)據(jù)結(jié)合季節(jié)項數(shù)據(jù)后，便能還原為原始數(shù)據(jù). 如圖4(c)所示，藍色曲線為原始數(shù)據(jù)，紅色曲線表示對原始數(shù)據(jù)直接進行SARIMA 建模所輸出的預(yù)測數(shù)據(jù)，綠色曲線則表示拆分項分別預(yù)測的數(shù)據(jù)加以整合后還原的數(shù)據(jù).

圖4 數(shù)據(jù)預(yù)測結(jié)果

圖4(c)中兩種預(yù)測方式的性能經(jīng)量化后，得到的直觀數(shù)據(jù)如表3 所示. 三組數(shù)據(jù)樣本的預(yù)測結(jié)果分別進行四項誤差分析，并對比二者性能差距.

表3 直接預(yù)測與拆分后預(yù)測效果對比

直接預(yù)測是直接對原始數(shù)據(jù)進行ARIMA 建模并預(yù)測，拆分后預(yù)測是通過X-11 拆分原始數(shù)據(jù)為趨勢項、季節(jié)項和不規(guī)則項，然后分別進行ARIMA 建模并預(yù)測，將各自預(yù)測結(jié)果按拆分原理進行還原，得到經(jīng)拆分后預(yù)測的結(jié)果.

數(shù)據(jù)對比方法采用四類方法，分別是平均絕對誤差(MAE)、均方誤差(MSE)、平均絕對百分比誤差(MAPE)和均方根誤差(RMSE). 對比直接預(yù)測和拆分預(yù)測的改進程度采用如下公式計算：

當(dāng) Δ 數(shù)值為正表明拆分預(yù)測較于直接預(yù)測沒有改善，當(dāng) Δ 數(shù)值為負時則表示拆分預(yù)測對比直接預(yù)測誤差程度得到了改進. 可以看出三組數(shù)據(jù)在擬合上均是拆分更優(yōu)，而拆分預(yù)測相較于直接預(yù)測都有著一定程度的改善. 在串絲斜坡數(shù)據(jù)上拆分預(yù)測提升顯著，整體提升約在50%. 可見在數(shù)據(jù)存量在300 歷元左右時，拆分預(yù)測較直接預(yù)測有著顯著的改善，而2 600歷元以上的數(shù)據(jù)中，拆分預(yù)測改進程度較小. 據(jù)此可以大致推斷用作平滑的數(shù)據(jù)范圍可以選定當(dāng)前歷元往前300 個歷元的數(shù)據(jù)為最佳.

3 結(jié) 論

時間序列分析法在其他領(lǐng)域有著成熟的運用[15]，基于該算法的特性，結(jié)合GNSS 行業(yè)中通過卡爾曼濾波等算法進行數(shù)據(jù)處理[16]的運用經(jīng)驗和總結(jié)出的一些問題，嘗試使用本文對時間序列分析法的運用方式解決實時性不足的問題. 理論上經(jīng)過預(yù)測獲得所需數(shù)據(jù)后，能夠?qū)崿F(xiàn)無數(shù)據(jù)缺失的完整平滑處理. 但是在地災(zāi)監(jiān)測中，重要的是真實和迅速，所以不能單純依賴預(yù)測，本文尚未解決預(yù)測數(shù)據(jù)與實時數(shù)據(jù)結(jié)合并互相修正的問題. 不過根據(jù)GNSS 行業(yè)其他方向的經(jīng)驗，可參考慣性導(dǎo)航系統(tǒng)(INS)/GNSS耦合方式，將ARIMA 和GNSS 經(jīng)卡爾曼濾波濾波后的數(shù)據(jù)進行解算層面的耦合，構(gòu)建一套以預(yù)測數(shù)據(jù)輔助后處理的運行系統(tǒng).