張明敏,劉 盼,周海龍,從建鋒

(山東科技大學(xué) 測繪科學(xué)與工程學(xué)院，山東 青島 266590)

目前，IGS測站在全球分布越來越多，積累近二十年的坐標(biāo)時間序列，準(zhǔn)確分析其變化特征，從而獲得精確的測站坐標(biāo)，不僅對于研究大陸板塊運(yùn)動具有重要的意義，而且可以對各種坐標(biāo)參考框架進(jìn)行修正，進(jìn)一步提高IGS測站坐標(biāo)的精度[1-2]。鑒于IGS測站的廣泛性及重要性，對IGS站點(diǎn)坐標(biāo)的準(zhǔn)確預(yù)測不僅在氣象預(yù)報，高精度快速測繪等領(lǐng)域具有重要作用，而且可以為地震分析預(yù)報、地殼變形監(jiān)測、研究地球物理現(xiàn)象提供一定的參考信息。GPS坐標(biāo)時間序列的研究始于上世紀(jì)九十年代，隨后越來越多的學(xué)者通過研究GPS坐標(biāo)時間序列獲得更多測繪學(xué)相關(guān)方面的成果[3]。國內(nèi)外學(xué)者研究發(fā)現(xiàn)GPS 坐標(biāo)時間序列中存在明顯的周期性(周年、半周年)運(yùn)動，且在坐標(biāo)垂向分量最為明顯[4-5]。根據(jù)不同研究角度，相關(guān)學(xué)者對GPS 坐標(biāo)時間序列時頻分析提出諸多方法，如主成分空間濾波法、最大似然估計(jì)法、小波相干分析法及功率譜分析法等[6-9]。此外，在時間序列模型研究方面，傳統(tǒng)時間序列模型存在建模過程復(fù)雜、無法動態(tài)調(diào)整模型參數(shù)、穩(wěn)定性弱等缺陷[10]，因此，學(xué)者們又提出自回歸綜合移動平均模型、奇異普分析模型和基于BP神經(jīng)網(wǎng)絡(luò)的時間序列融合模型等改進(jìn)模型[11-14]。為達(dá)到對站點(diǎn)高程精準(zhǔn)預(yù)報及穩(wěn)定性監(jiān)測目的，選擇高精度的時間序列模型至關(guān)重要?；诖?，本文根據(jù)SOPAC提供的時間序列數(shù)據(jù)，分別確立多項(xiàng)式周期模型參數(shù)與自回歸移動平均(Autoregressive and Moving Average，ARMA)模型參數(shù)，并對這兩種模型的預(yù)測精度進(jìn)行實(shí)驗(yàn)驗(yàn)證及對比分析，從而選出最優(yōu)模型。

1 數(shù)據(jù)獲取及預(yù)處理

1.1 數(shù)據(jù)來源

SOPAC提供IGS基準(zhǔn)站單天解坐標(biāo)時間序列數(shù)據(jù)產(chǎn)品，主要分為RAW-TRE、CL-TRE、CL-RES及CL-DETRE等不同類型，為了實(shí)驗(yàn)驗(yàn)證及對比分析兩種模型的預(yù)測精度，在區(qū)域上，本文選取我國及周邊7個IGS站CL-TRE產(chǎn)品(經(jīng)數(shù)據(jù)分析中心對原始觀測數(shù)據(jù)降噪后的坐標(biāo)時間序列)的高程坐標(biāo)。所選站點(diǎn)如圖1所示。由于2014—2017年IGS站時間序列數(shù)據(jù)缺失嚴(yán)重，所以選用于擬合模型的時間序列范圍為2009-01-01—2012-12-31，并對2013全年數(shù)據(jù)進(jìn)行預(yù)測及精度檢驗(yàn)。

圖1 所選IGS站點(diǎn)分布

1.2 數(shù)據(jù)預(yù)處理

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

圖2 各測站預(yù)處理前后時間序列對比

2 模型參數(shù)建立及數(shù)據(jù)測試分析

2.1 多項(xiàng)式周期模型

2.1.1 功率譜分析

功率譜分析是一種通過傅里葉變換將時域內(nèi)的信號轉(zhuǎn)換到頻域內(nèi)進(jìn)行研究分析的方法，主要用來分析離散數(shù)據(jù)的周期特性。若一組離散數(shù)據(jù)具有周期性，則其周期運(yùn)動對應(yīng)的功率在全部功率中占很大比重，在功率譜圖中對應(yīng)著功率的峰值[14]，從而在時域內(nèi)不易反映出的特性，可以在頻域內(nèi)很容易觀察出來。

對IGS站點(diǎn)消除趨勢項(xiàng)并進(jìn)行傅里葉變換后，得到各個測站U分量的功率譜，見圖3。通過觀察各個測站高程方向坐標(biāo)時間序列的功率譜圖可以明顯地看到時間序列周期特性趨于一致，除了具有明顯的年周期項(xiàng)外，還有相對年周期能量較小的半年周期項(xiàng)。雖然不同站點(diǎn)的周期并非完全相同，即并不是嚴(yán)格意義上周年項(xiàng)周期為365.25 d，半周年項(xiàng)周期為182.62 d，但可從圖3中看出3個測站的年周期在363～368 d內(nèi)，半年周期在181～184 d內(nèi)，即可以進(jìn)行周期擬合。

2.1.2 模型構(gòu)建

各站點(diǎn)高程方向坐標(biāo)以年周期運(yùn)動為主，同時存在半年周期運(yùn)動。對于趨勢項(xiàng)通常采用多項(xiàng)式函數(shù)進(jìn)行擬合，而周期項(xiàng)可采用周期函數(shù)進(jìn)行擬合，故本文采用的擬合函數(shù)模型為：

圖3 4個IGS站高程方向上的功率譜

yi=a+b·ti+A1sin(2πti+φ1)+

A2sin(4πti+φ2).

(1)

式中：a為常數(shù)項(xiàng)，b為線性運(yùn)動速率，A1，φ1分別為年周期項(xiàng)的振幅和相位，A2，φ2分別為半年周期項(xiàng)的振幅和相位。對式(1)線性化后變?yōu)椋?/p>

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測站高程坐標(biāo)代入上式，可求得各測站的擬合函數(shù)模型的參數(shù)，見表1。由表1可知每個站點(diǎn)的年周期項(xiàng)振幅比半年周期項(xiàng)振幅大，說明年周期運(yùn)動是站點(diǎn)的主要運(yùn)動規(guī)律。

表1 各測站擬合參數(shù)

2.2 ARMA模型

ARMA模型為

Xt=φ1Xt-1+φ2Xt-2+…+φpXt-p-

θ1∈t-1-θ2∈t-2-…-θq∈t-q.

(5)

式中:p和q是模型的自回歸階數(shù)和移動平均階數(shù)；φ，θ是不為零的待定系數(shù)；Xt是平穩(wěn)、正態(tài)、零均值的時間序列；∈t是獨(dú)立的誤差項(xiàng)。建模的基本步驟：①消除IGS站點(diǎn)高程坐標(biāo)時間序列的趨勢項(xiàng)；②采用AIC準(zhǔn)則確定p，q階數(shù)的值；③利用Matlab軟件確定模型參數(shù)，同時，查看殘差∈t是否為白噪聲序列判斷模型的適用性；④增加模型趨勢項(xiàng)。

2.3 模型驗(yàn)證及對比分析

利用以上兩種模型分別預(yù)報所選測站2013-01-01—2013-12-31的高程坐標(biāo)，并與SOPAC提供的時間序列數(shù)據(jù)進(jìn)行對比分析。預(yù)報結(jié)果如表2所示。

表2 兩種預(yù)測模型的預(yù)報中誤差對比

由表2可知，多項(xiàng)式周期模型的預(yù)報中誤差絕對值在2～4.5 mm之內(nèi)，其平均值為±3.63 mm，預(yù)報中誤差絕對值的最大值為4.10 mm，最小值為2.93 mm，其中，只有BJFS站的預(yù)報中誤差絕對值小于3 mm，而CHAN站與KUNM站的預(yù)報中誤差絕對值則大于4 mm。在ARMA模型中，預(yù)報中誤差絕對值均小于4 mm，平均值為±3.33 mm，此外，該模型預(yù)報中誤差絕對值的最大值為3.75 mm，最小值為2.71 mm，除BJFS站的預(yù)報中誤差絕對值小于3 mm外，LHAZ站的預(yù)報中誤差絕對值同樣小于3 mm。通過對比兩種模型的預(yù)報中誤差，發(fā)現(xiàn)ARMA模型的預(yù)報精度整體高于多項(xiàng)式周期模型，且精度提升在5%以上，其中，BJFS站、CHAN站及LHAZ站的精度提升在6%～7%之間。KUNM站與SHAO站的精度提升甚至高于10%，這可能是由于預(yù)處理前的KUNM站與SHAO站時間序列完整性強(qiáng)所造成的，同時，在多項(xiàng)式周期模型中，所選IGS測站預(yù)報中誤差的標(biāo)準(zhǔn)差為0.51 mm，而在ARMA模型中，所選IGS測站預(yù)報中誤差的標(biāo)準(zhǔn)差為0.33 mm，可見，ARMA模型的穩(wěn)定性強(qiáng)于多項(xiàng)式周期模型。綜上可知ARMA模型更加適合預(yù)報測站高程坐標(biāo)。

3 結(jié)束語

本文通過對多項(xiàng)式周期模型與ARMA模型進(jìn)行實(shí)驗(yàn)驗(yàn)證及對比分析。發(fā)現(xiàn)兩種模型的坐標(biāo)預(yù)報中誤差絕對值均在4.5 mm以內(nèi)，即符合預(yù)報精度，同時，ARMA模型的預(yù)報精度較多項(xiàng)式周期模型的預(yù)報精度提升5%以上，且穩(wěn)定性高于多項(xiàng)式周期模型，由此可見ARMA模型更加適合作為站點(diǎn)高程預(yù)報模型。本文用7個IGS站4年的數(shù)據(jù)擬合各測站高程方向的時間序列模型，有待用更長的時間序列和更加精確的模型加以檢驗(yàn)及水平方向的分析。