許煥盛

(廣東省地質(zhì)局第二地質(zhì)大隊，廣東 汕頭 515000)

山東省于2007年啟動了山東省衛(wèi)星定位連續(xù)運行綜合應用服務系統(tǒng)(Shan Dong Continuously Operating Reference Station,SDCORS)項目，積累的觀測數(shù)據(jù)不斷豐富，為研究區(qū)域形變、地球板塊運動和地球自轉(zhuǎn)等提供重要的數(shù)據(jù)支撐[1]。研究GPS坐標時間序列的非線性變化，可獲得基準站的位置變化及三維速度場，為工程應用基準站的選擇、基準站的速度及數(shù)據(jù)分析提供參考[2-3]。相關研究表明：GPS坐標時間序列不僅具有線性變化趨勢，還具有年、半年周期等的非線性變化特征，其中以U方向表征的最為明顯[4-9]。本文以12個SDCORS站連續(xù)4年的觀測序列為研究對象，計算坐標時間序列的功率譜，通過功率譜圖分析其周期項，并將周期項特征進行統(tǒng)計，這對于分析基準站的穩(wěn)定性及位置變化具有重要意義。

1 數(shù)據(jù)處理策略

1.1 數(shù)據(jù)預處理

選取的各基準站的基本信息如表1所示。連續(xù)均勻的坐標時間序列是數(shù)據(jù)分析的基礎，因此在數(shù)據(jù)分析之前需對數(shù)據(jù)進行預處理，包括剔除孤立值與插值、階躍項改正、趨勢項擬合及去除趨勢項等步驟。

表1 站點基本信息

1.1.1 孤立值的剔除與數(shù)據(jù)補全

由于外界觀測條件、多路徑效應及傳輸信號干擾等因素的影響，坐標時間序列中不可避免含有孤立值，也稱粗差，這些粗差的存在會造成數(shù)據(jù)分析結(jié)果的偏差。因此，在數(shù)據(jù)進行周期分析和求站點速度之前，應將其剔除。GPS坐標時間序列由于基站的改造、觀測條件不充分等干擾因素，導致在某天或某段時間出現(xiàn)部分數(shù)據(jù)缺失或不完整的現(xiàn)象，且在剔除粗差的步驟中，部分點被剔除，需要將這兩部分的數(shù)據(jù)進行擬合補齊。

本文采用3倍中誤差為判別原則，簡稱“3S”方法，探測坐標時間序列中的粗差，并將其剔除[9]?！?S”基本原理如下：

(1)

則認為Sa為粗差，需將其從觀測序列中剔除。通常采用迭代的方法剔除坐標時間序列含有的粗差。

對于缺失和剔除粗差的數(shù)據(jù)，經(jīng)過多次實驗與分析，采用三次多項式插值的方法對數(shù)據(jù)進行補全，保證數(shù)據(jù)的連續(xù)性和完整性，為后續(xù)的研究和分析奠定基礎。

1.1.2 階躍項的修正

由于接收機故障、天線更換、地震等原因，GPS坐標時間序列在某時段出現(xiàn)階躍或稱中斷，會導致測站速度估計有偏，需要對其進行修正。具體做法是從坐標突變的位置開始，往后的坐標時間序列統(tǒng)一加上階躍值進行修正。由于sdrc基準站U方向的坐標時間序列中含有明顯的坐標階躍，本文就以sdrc為例對階躍項的修正作以說明。圖1所示為sdrc站U方向的原始坐標序列，從圖可以看出存在兩處坐標突變的地方，分別在第307天和第775天。圖2所示為sdrc階躍修正后的坐標時間序列，由圖可知，經(jīng)過階躍修正后的坐標時間序列在66.42 m附近波動，說明具有良好的連續(xù)性與一致性。

圖1 sdrc原始坐標時間序列

圖2 sdrc階躍修正后的坐標時間序列

1.1.3 趨勢項擬合

在進行譜分析時要求坐標時間序列為零均值，即沒有線性趨勢項。為此，需要對坐標時間序列進行線性擬合，其原理如下：

yt=β0+β1t+xt

(2)

式中，β0、β1為線性趨勢擬合的常數(shù)項和一次項系數(shù)。利用該模型計算得出各個時間點的擬合值，然后將對應時間點的坐標時間序列值減去擬合值，得到殘差時間序列，即

(3)

圖3所示為sdrc基準站三個方向經(jīng)過粗差剔除、數(shù)據(jù)插值和階躍修正后的坐標時間序列及線性趨勢擬合。由圖可知：sdrc站的N、E方向具有明顯的線性變化，U方向不僅有線性變化，還表現(xiàn)出周期性變化，因此U方向的運動是十分復雜的；U方向數(shù)據(jù)的離散程度要比N、E方向大，這主要是因為U方向的精度要低于其它兩個方向。在進行功率譜分析時要求數(shù)據(jù)零均值，因此需要將線性趨勢項從坐標時間序列中去除，以滿足功率譜分析的條件。圖4所示為sdrc去除線性趨勢項的坐標時間序列，由圖可知：坐標時間序列經(jīng)過上述步驟處理后，變化量在零的上下進行波動，說明具備了零均值的特性。同時圖4的坐標時間序列表現(xiàn)出明顯的季節(jié)性變化，說明存在一定的周期項，這就需要對坐標時間序列進行功率譜估計，探測潛在的周期。

圖3 sdrc坐標時間序列及擬合直線

圖4 sdrc殘差序列

1.2 坐標時間序列的功率譜估計

本文采用傳統(tǒng)功率譜估計——直接法，探測坐標時間序列中存在的潛在周期，基本思想為：直接對坐標時間序列作離散傅里葉變換，計算求得坐標時間序列的幅值，然后將幅值的平方除以時間長度作為功率譜密度，得到功率譜密度與周期之間的關系圖，也稱之為周期圖法。周期圖的峰值則為坐標時間序列可能存在的周期，且能量越高，對應的周期越明顯[10]。其優(yōu)點是計算簡單，而且不依賴于相關函數(shù)?；驹砣缦拢?/p>

設x(n)為坐標時間序列，其長度為N的離散傅里葉變換為：

(4)

則功率譜密度定義為式(5)：

(5)

2 坐標時間序列周期性分析

坐標時間序列經(jīng)過數(shù)據(jù)預處理得到部分基準站3個方向的殘差序列及各方向的功率譜圖，如圖5所示。

圖5 殘差序列與功率譜圖

圖5所示(a)(c)(e)(j)為各站的殘差序列圖，橫坐標為年份，縱坐標為各方向的振幅，單位為mm；(b)(d)(f)(h)為各站的功率譜圖，橫坐標為周期，單位為a，縱坐標為單位周期內(nèi)的能量，單位為mm2/a。由圖5可知：各站點在N方向的能量較小，均小于5×107，說明周期項不明顯。其中sdrc站半年周期項能量最大，為4.963×107，說明sdrc站半年周期運動較明顯，其他幾個站半年周期項的能量差別不大。另外jimo站、sdrc站、rizh站還存在2年的周期項，cash站、jimo站、rizh站可能還存在更長的周期項。E方向上，各站半年周期項對應的能量突出，功率譜處于8～10×107之間，能量值相差不大，說明各站的半年周期項運動在E方向明顯。U方向上，各基準站2年周期項的能量最高，其中sdrc站能量最小，為8.565×107，其他三個站的量級為108。

對其它基準站的功率譜圖進行分析，統(tǒng)計各基準站的周期項及功率譜如表2所示，表中功率譜采用科學計數(shù)法，量級為107。

表2 基準站周期項及功率譜

由表2可知：

(1)基準站在N、E、U方向均表現(xiàn)出一定的周期性。N方向的周期運動不明顯；E方向主要以0.5年周期項運動最為突出，且能量值差別不大；U方向存在2年周期項，絕大部分基準站存在1年周期項，且2年周期項的能量高于1年周期項的能量，說明U方向的2年周期運動較為明顯。

(2)坐標時間序列中并不是嚴格存在半年周期項和年周期項。從功率譜可以看出，能量最強對應的周期并不是嚴格的半年周期項、年周期項等周期，而是在周期項的附近，稱為類半年周期項、類年周期項等。

(3)坐標時間序列功率譜在E方向的0.5年周期項的能量普遍高于N方向的，說明E方向的0.5年周期性運動比N方向的要明顯。

(4)N方向整體周期項不明顯，且大多數(shù)周期項不明顯的基準站具有更長周期項的趨勢，具體原因有待分析。

3 結(jié) 論

本文對SDCORS站坐標時間序列的周期特征進行了分析與統(tǒng)計，得到：SDCORS站坐標時間序列存在周期項，且N、E方向年、半年周期項明顯，U方向2年周期項明顯，功率譜能量明顯高于其它兩個方向，這是由于U方向的影響因素較多。在后續(xù)的研究中，將對更長時間跨度的坐標時間序列進行周期特征的研究與分析，探討周期項的影響因素，為得到精確的噪聲模型提供數(shù)據(jù)基礎。