基于奇異譜分析法的GPS時間序列周期項探測

湯文娟

(1.廣州市房地產測繪院，廣東 廣州 510030； 2.廣州市測繪產品質量檢驗中心，廣東 廣州 510030)

1 引 言

高精度、實時的非線性運動特性研究和監測是大地測量學科的熱門研究課題。只有在采用更完善的非線性運動模型的基礎上，才能構建更高精度參考框架。

GPS參考框架點坐標變化具有復雜的非線性特征，很多經典的時間序列分析方法并不能有效地分離出不同尺度的周期運動特性，從而無法對所得結果給出合理的地球物理解釋。奇異譜分析已經被證明為分解時間序列的有力工具，能彌補常用譜分析的不足。奇異譜分析方法的優越性主要在于①不需要預先給定濾波周期，只需根據資料自身確定，具有較強的自適應性。②對原始序列要求比較寬松，不需要對統計分布和平穩性做假設[1]。奇異譜分析方法根據序列自身的時間相關特性可以對序列進行動力重構，進行不同振蕩頻率的信號分離，在測繪領域中廣泛用于序列插值、濾波去噪、趨勢識別、周期項提取以及預報模型的建立。按照時間序列分析理論，每一個時間序列經過合理的變換后都可以分解為趨勢項、周期項和隨機噪聲三個部分[2]。本文在時域和頻域內，應用奇異譜分析對GPS連續觀測站的位置變化情況進行分解，提取不同尺度的周期項，并與經典GPS時間序列模型最小二乘擬合結果進行對比。

2 奇異譜分析原理及周期項探測方法

2.1 奇異譜分析的基本原理

奇異譜分析(Singular-Spectrum Analysis，SSA)是時間序列中常用的分析與預測技術，組合了經典時間序列分析、多元幾何、多元統計、動態系統和信號處理等多種元素。奇異譜分析將原始序列分解為緩慢變化趨勢、周期項、噪聲序列，基于時間序列的動力重構出發、與經驗正交函數相聯系，從事先未知物理本質的、包含噪聲的有限長觀測序列中，濾去非周期性的異常現象，對方差譜信號有強化和放大，將頻域信號分解為時頻信號加以識別和估計，提取盡可能多的可靠信息[3]。

應用奇異譜分析法研究分析框架非線性變化的特征，從框架點坐標序列中提取其周年運動、半周年運動等多時間、多尺度的非線性運動，從而構建非線性運動參考框架。對于長度N(N>2)的時間序列X=XN=(x1，…，xN)，取窗體長度為L(1

(1)嵌入(Embedding)

嵌入過程是將一個一維的時間序列X=(x1，…，xN)轉換為多維序列X1，…，XK，即將原時間序列映射為K=N-L+1個L滯后向量，時間序列X的L滯后軌跡矩陣X為：

(1)

其中Xi=(xi，xi+1，…，xi+L-1)。

(2)軌跡矩陣的奇異值分解(SVD)

定義矩陣S=XXT，XT為X的轉置矩陣。計算矩陣S的特征值λi和特征向量Ui，并將特征值按從大到小排列，其特征值依次為λ1≥…≥λL≥0，對應的特征向量為U1，…，UL。記Vi=XTUi(i=1，…，L)。其中Ui是矩陣的左特征向量，又稱為時間經驗正交函數(Temporal Empirical Orthogonal Function，T-EOF)；Vi是矩陣的右特征向量，又稱為時間主成分(Temporal Principal Components，T-PC)，軌跡矩陣X可以由初等矩陣合成X=X1+…+Xd。初等矩陣為：

(2)

(3)特征值三要素分組

將初等矩陣Xi的下標{1，2，…，d}分成p個不相交的子集I1，I2，…，Ip，設I={i1，i2，…，im}，那么合成矩陣XI=Xi1+Xi2+…+Xim，計算集合I=I1，I2，…，Ip的每個合成矩陣，那么序列的分解式可寫為：

X=XI1+XI2+…+XIp

(3)

其中，選取集合I1，I2，…，Ip的過程稱為分組。

(4)通過對角平均重構時間序列

奇異譜分析將序列分離成各種單一頻率的振蕩信號以及噪聲項，利用時間主成分和時間經驗正交函數展開各分量，可以獲得對應于各頻率振蕩信號的重構分量序列，從而達到提取有用信息、過濾噪聲的目的。將分組后的矩陣XIj重構為長度為N的新序列。設Y是一個L×K的矩陣，其元素為yij，其中1≤i≤L，1≤j≤K，L*=min(L，K)，K*=max(L，K)且N=L+K-1，根據對角平均公式將矩陣Y轉換為y1，…，yN的時間序列，對角平均公式為：

(4)

(5)

2.2 基于SSA的周期項探測

(1)對應特征值接近相等；

3 GPS坐標序列中周期項探測

以美國西北部阿拉斯加州觀測質量好、精度高且穩定性強的GPS連續觀測臺站AV06為例，對其周期性運動特性進行分析，測站位置如圖1所示。

圖1 AV06觀測臺站位置分布圖

根據實際分析處理經驗所知，在周期項探測的過程中，用于分析的時間序列跨度和窗口長度L的選取對于分析結果有著很大的影響。如果需要分析序列中的年周期項或者更大尺度的周期項，則需要較長的時間序列，并且嵌入維數最好設置為周期的整數倍，這樣可以獲得較為精確的周期項信號。

對AV06觀測臺站的日觀測坐標序列進行預處理，采用奇異譜插值法得到2005年～2016年12年間連續的3 902組日觀測坐標序列，如圖2所示。

圖2 AV06臺站日觀測坐標序列

3.1 AV06測站N方向的周期項探測

基于分析年周期的需要，在構造軌跡矩陣時，將嵌入維數設置為 1 500。分析特征值的分布情況，發現有幾組明顯成對的特征值，分別對前三對分量進行周期分析。周期探測結果如表1所示。

AV06-N主要成對RC的檢驗參數和周期 表1

綜合三對RC分量的奇異譜分析結果，可以得出測站AV06在北方向上存在顯著的年周期、半年周期信號。其他的RC成分在分析過程中并不滿足周期波動成分檢驗的三條判別標準，暫時無法探測出其他周期成分。

(1)將第3、4主分量進行序列重構得到年周期項，并與最小二成擬合所得年周期進行對比，如圖3所示。

圖3 AV06-N年周期項對比圖

其中，藍色曲線為SSA分析所得年周期，綠色曲線為最小二乘擬合所得年周期信號。從圖3中可以看出，經SSA分析，GPS時間序列中存在370天的周期項，即年周期項(由于頻率設置的原因，如果頻率之間的顆粒度更小，可以獲得更精確的周期)。經SSA探測出的周期項與最小二乘擬合的周期項大體一致，但是最小二乘擬合結果為振幅恒定的周期振幅，而SSA所探測到的周期項振幅有隨時間而增大的傾向。這對于框架的非線性運動特性研究有很好的指導意義，并且振幅逐年增大的趨勢也可以反過來配合一些地球物理因素的解釋，本文就不做進一步探討。

(2)將第5、6主分量進行序列重構得到半年周期信號，并與最小二成擬合所得半年周期進行對比，如圖4所示。

圖4 AV06-N半年周期項對比圖

通過周期圖可以看出，GPS時間序列中存在175天的周期信號，即半年周期信號。半年周期對比圖與年周期對比圖相比，區別要更明顯，這是因為SSA探測出的周期與半年有一定差距，隨著時間的推移，兩者半年周期信號開始出現峰值偏移。另外，奇異譜分析所得的半年周期信號的振幅除了有著明顯隨時間增大的趨勢，其中也暗含著大約2.5年的周期項。這個現象有待進一步討論。

3.2 AV06測站E、U方向的周期項探測

對于AV06測站E、U方向2005年～2016年12年間的日觀測坐標序列(共3 902組數據)，構造嵌入維數為1 500的軌跡矩陣。分析特征值的分布情況，分別對前四對、前三對特征值明顯成對的分量進行周期分析，分析結果如表2、表3所示。

AV06-E主要成對RC的檢驗參數和周期 表2

AV06-U主要成對RC的檢驗參數和周期 表3

從表2、表3中結果可以看出，測站AV06在E、U方向的分量坐標序列，除了常見的半年周期和年周期，還有1.5年周期、3年周期等等。由此可見，同一測站不同方向上坐標序列的周期性也不盡相同，但都含有近似年周期、半年周期項。

綜合測站AV06三個坐標分量上的SSA周期探測結果，有以下兩條結論：

①各分量序列中都含有半年周期、年周期項，除此之外，各坐標分量的周期性還有細小差別，即不同方向上的坐標序列可能還包含1.5年周期、3年周期，具體情形又略有不同。

②跟最小二乘擬合結果有差異。奇異譜探測到的周期項并非振幅恒定，周期信號的振幅有隨時間增大的趨勢，并且在總體增大趨勢中還暗含周期變化。

4 周期項影響因素分析

通過奇異譜分析方法提取出的周期項可知，不同測站的不同方向上的周期項也不盡相同，也不是嚴格的周年、半周年周期項，這些相近的周期項可能由同一個物理因素所引起，由于受分析方法的分辨率和定位精度所限，而表現出略微的差異[5]。下面對于周期性變化的幾個主要影響因素進行分析。

(1)非模型化的地殼周期運動。奇異譜分析提取出的近似年周期、半年周期就屬于地殼周期運動的結果。實驗結果表明，GPS三個方向上時間序列并不存在嚴格的周年、半周年運動，而且不同測站之間的周期項具體周期值也互不相同。這種情況是由地殼運動的復雜性所致，也是嚴格的地殼周期運動與局部人類活動疊加的結果。

(2)測站周期變化。各種地球物理參數會對GPS測站時間序列中坐標有所影響，未模型化的固體潮、海洋潮汐以及大氣質量負荷所引起的荷載會使測站表現出周期性變化。Amiri-Simkooei等[6]通過諧波估計，提取出很多介于170天～180天之間的周期項，劉偉、李昭等[5]也通過功率譜分析方法獲得相似的結論，這恰好與Penna等給出的未模型化S2海洋潮汐荷載效應相吻合。

(3)多路徑效應。由于GPS星座幾何形狀呈周期性變化，這導致多路徑效應也呈周期性變化，而GPS星座周期性變化與計算GPS測站坐標所基于的太陽日之間存在約 247 s的差異，因此會產生350天的周期項[7]。因為本文中計算頻率的精度在±0.000 1，因此提取出的周期項應該在338天～363天之間。

(4)周期性變化的潛在因素。SSA所探測出的長周期項沒有物理背景所對應，可能是由其他很多潛在的影響因素所導致的，比如臺站下基巖的熱膨脹、天線相位中心模型誤差、對流程濕分量影響、軌道模型誤差等等[8]。

關于所探測出的周期結果，本文并未做進一步探討，但是可為參考框架點的非線性運動研究提供一定的數據支持。同時，這樣的周期變化特性又該如何結合實際的地球物理因素給出模型化的改正，這也是有待解決的問題。