張明敏,劉 盼,周海龍,從建鋒
(山東科技大學 測繪科學與工程學院,山東 青島 266590)
目前,IGS測站在全球分布越來越多,積累近二十年的坐標時間序列,準確分析其變化特征,從而獲得精確的測站坐標,不僅對于研究大陸板塊運動具有重要的意義,而且可以對各種坐標參考框架進行修正,進一步提高IGS測站坐標的精度[1-2]。鑒于IGS測站的廣泛性及重要性,對IGS站點坐標的準確預測不僅在氣象預報,高精度快速測繪等領域具有重要作用,而且可以為地震分析預報、地殼變形監測、研究地球物理現象提供一定的參考信息。GPS坐標時間序列的研究始于上世紀九十年代,隨后越來越多的學者通過研究GPS坐標時間序列獲得更多測繪學相關方面的成果[3]。國內外學者研究發現GPS 坐標時間序列中存在明顯的周期性(周年、半周年)運動,且在坐標垂向分量最為明顯[4-5]。根據不同研究角度,相關學者對GPS 坐標時間序列時頻分析提出諸多方法,如主成分空間濾波法、最大似然估計法、小波相干分析法及功率譜分析法等[6-9]。此外,在時間序列模型研究方面,傳統時間序列模型存在建模過程復雜、無法動態調整模型參數、穩定性弱等缺陷[10],因此,學者們又提出自回歸綜合移動平均模型、奇異普分析模型和基于BP神經網絡的時間序列融合模型等改進模型[11-14]。為達到對站點高程精準預報及穩定性監測目的,選擇高精度的時間序列模型至關重要。基于此,本文根據SOPAC提供的時間序列數據,分別確立多項式周期模型參數與自回歸移動平均(Autoregressive and Moving Average,ARMA)模型參數,并對這兩種模型的預測精度進行實驗驗證及對比分析,從而選出最優模型。
SOPAC提供IGS基準站單天解坐標時間序列數據產品,主要分為RAW-TRE、CL-TRE、CL-RES及CL-DETRE等不同類型,為了實驗驗證及對比分析兩種模型的預測精度,在區域上,本文選取我國及周邊7個IGS站CL-TRE產品(經數據分析中心對原始觀測數據降噪后的坐標時間序列)的高程坐標。所選站點如圖1所示。由于2014—2017年IGS站時間序列數據缺失嚴重,所以選用于擬合模型的時間序列范圍為2009-01-01—2012-12-31,并對2013全年數據進行預測及精度檢驗。

圖1 所選IGS站點分布
通過觀察選取的GPS站點高程時間序列,發現時間序列觀測值中包含少量的突變點和間斷點。在對高程坐標時間序列進行周期特性分析時,都應剔除這些孤立值且將間斷點數據補充完整,以免降低分析精度。本文采用穩健孤立值探測法IQR準則剔除突變點[16],同時,利用正交多項式做最小二乘擬合完成間斷點插值。數據預處理前后對比如圖2所示,鑒于篇幅有限,文中列出4個測站預處理后高程坐標時間序列的對比圖,其中突變點明顯減少,間斷點也得到有效處理。






圖2 各測站預處理前后時間序列對比
2.1.1 功率譜分析
功率譜分析是一種通過傅里葉變換將時域內的信號轉換到頻域內進行研究分析的方法,主要用來分析離散數據的周期特性。若一組離散數據具有周期性,則其周期運動對應的功率在全部功率中占很大比重,在功率譜圖中對應著功率的峰值[14],從而在時域內不易反映出的特性,可以在頻域內很容易觀察出來。
對IGS站點消除趨勢項并進行傅里葉變換后,得到各個測站U分量的功率譜,見圖3。通過觀察各個測站高程方向坐標時間序列的功率譜圖可以明顯地看到時間序列周期特性趨于一致,除了具有明顯的年周期項外,還有相對年周期能量較小的半年周期項。雖然不同站點的周期并非完全相同,即并不是嚴格意義上周年項周期為365.25 d,半周年項周期為182.62 d,但可從圖3中看出3個測站的年周期在363~368 d內,半年周期在181~184 d內,即可以進行周期擬合。
2.1.2 模型構建
各站點高程方向坐標以年周期運動為主,同時存在半年周期運動。對于趨勢項通常采用多項式函數進行擬合,而周期項可采用周期函數進行擬合,故本文采用的擬合函數模型為:


圖3 4個IGS站高程方向上的功率譜
yi=a+b·ti+A1sin(2πti+φ1)+
A2sin(4πti+φ2).
(1)
式中:a為常數項,b為線性運動速率,A1,φ1分別為年周期項的振幅和相位,A2,φ2分別為半年周期項的振幅和相位。對式(1)線性化后變為:
yi=a+b·ti+csin(2πti)+dcos(2πti)+
esin(4πti)+fcos(4πti).
(2)


(3)

(4)
將2009-01-01—2012-12-31 IGS測站高程坐標代入上式,可求得各測站的擬合函數模型的參數,見表1。由表1可知每個站點的年周期項振幅比半年周期項振幅大,說明年周期運動是站點的主要運動規律。

表1 各測站擬合參數
ARMA模型為
Xt=φ1Xt-1+φ2Xt-2+…+φpXt-p-
θ1∈t-1-θ2∈t-2-…-θq∈t-q.
(5)
式中:p和q是模型的自回歸階數和移動平均階數;φ,θ是不為零的待定系數;Xt是平穩、正態、零均值的時間序列;∈t是獨立的誤差項。建模的基本步驟:①消除IGS站點高程坐標時間序列的趨勢項;②采用AIC準則確定p,q階數的值;③利用Matlab軟件確定模型參數,同時,查看殘差∈t是否為白噪聲序列判斷模型的適用性;④增加模型趨勢項。
利用以上兩種模型分別預報所選測站2013-01-01—2013-12-31的高程坐標,并與SOPAC提供的時間序列數據進行對比分析。預報結果如表2所示。

表2 兩種預測模型的預報中誤差對比
由表2可知,多項式周期模型的預報中誤差絕對值在2~4.5 mm之內,其平均值為±3.63 mm,預報中誤差絕對值的最大值為4.10 mm,最小值為2.93 mm,其中,只有BJFS站的預報中誤差絕對值小于3 mm,而CHAN站與KUNM站的預報中誤差絕對值則大于4 mm。在ARMA模型中,預報中誤差絕對值均小于4 mm,平均值為±3.33 mm,此外,該模型預報中誤差絕對值的最大值為3.75 mm,最小值為2.71 mm,除BJFS站的預報中誤差絕對值小于3 mm外,LHAZ站的預報中誤差絕對值同樣小于3 mm。通過對比兩種模型的預報中誤差,發現ARMA模型的預報精度整體高于多項式周期模型,且精度提升在5%以上,其中,BJFS站、CHAN站及LHAZ站的精度提升在6%~7%之間。KUNM站與SHAO站的精度提升甚至高于10%,這可能是由于預處理前的KUNM站與SHAO站時間序列完整性強所造成的,同時,在多項式周期模型中,所選IGS測站預報中誤差的標準差為0.51 mm,而在ARMA模型中,所選IGS測站預報中誤差的標準差為0.33 mm,可見,ARMA模型的穩定性強于多項式周期模型。綜上可知ARMA模型更加適合預報測站高程坐標。
本文通過對多項式周期模型與ARMA模型進行實驗驗證及對比分析。發現兩種模型的坐標預報中誤差絕對值均在4.5 mm以內,即符合預報精度,同時,ARMA模型的預報精度較多項式周期模型的預報精度提升5%以上,且穩定性高于多項式周期模型,由此可見ARMA模型更加適合作為站點高程預報模型。本文用7個IGS站4年的數據擬合各測站高程方向的時間序列模型,有待用更長的時間序列和更加精確的模型加以檢驗及水平方向的分析。