基于奇異譜分析法的GPS時(shí)間序列周期項(xiàng)探測(cè)

湯文娟

(1.廣州市房地產(chǎn)測(cè)繪院，廣東 廣州 510030； 2.廣州市測(cè)繪產(chǎn)品質(zhì)量檢驗(yàn)中心，廣東 廣州 510030)

1 引 言

高精度、實(shí)時(shí)的非線性運(yùn)動(dòng)特性研究和監(jiān)測(cè)是大地測(cè)量學(xué)科的熱門(mén)研究課題。只有在采用更完善的非線性運(yùn)動(dòng)模型的基礎(chǔ)上，才能構(gòu)建更高精度參考框架。

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

2 奇異譜分析原理及周期項(xiàng)探測(cè)方法

2.1 奇異譜分析的基本原理

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

應(yīng)用奇異譜分析法研究分析框架非線性變化的特征，從框架點(diǎn)坐標(biāo)序列中提取其周年運(yùn)動(dòng)、半周年運(yùn)動(dòng)等多時(shí)間、多尺度的非線性運(yùn)動(dòng)，從而構(gòu)建非線性運(yùn)動(dòng)參考框架。對(duì)于長(zhǎng)度N(N>2)的時(shí)間序列X=XN=(x1，…，xN)，取窗體長(zhǎng)度為L(zhǎng)(1

(1)嵌入(Embedding)

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

(1)

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

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

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

(2)

(3)特征值三要素分組

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

X=XI1+XI2+…+XIp

(3)

其中，選取集合I1，I2，…，Ip的過(guò)程稱(chēng)為分組。

(4)通過(guò)對(duì)角平均重構(gòu)時(shí)間序列

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

(4)

(5)

2.2 基于SSA的周期項(xiàng)探測(cè)

(1)對(duì)應(yīng)特征值接近相等；

3 GPS坐標(biāo)序列中周期項(xiàng)探測(cè)

以美國(guó)西北部阿拉斯加州觀測(cè)質(zhì)量好、精度高且穩(wěn)定性強(qiáng)的GPS連續(xù)觀測(cè)臺(tái)站AV06為例，對(duì)其周期性運(yùn)動(dòng)特性進(jìn)行分析，測(cè)站位置如圖1所示。

圖1 AV06觀測(cè)臺(tái)站位置分布圖

根據(jù)實(shí)際分析處理經(jīng)驗(yàn)所知，在周期項(xiàng)探測(cè)的過(guò)程中，用于分析的時(shí)間序列跨度和窗口長(zhǎng)度L的選取對(duì)于分析結(jié)果有著很大的影響。如果需要分析序列中的年周期項(xiàng)或者更大尺度的周期項(xiàng)，則需要較長(zhǎng)的時(shí)間序列，并且嵌入維數(shù)最好設(shè)置為周期的整數(shù)倍，這樣可以獲得較為精確的周期項(xiàng)信號(hào)。

對(duì)AV06觀測(cè)臺(tái)站的日觀測(cè)坐標(biāo)序列進(jìn)行預(yù)處理，采用奇異譜插值法得到2005年～2016年12年間連續(xù)的3 902組日觀測(cè)坐標(biāo)序列，如圖2所示。

圖2 AV06臺(tái)站日觀測(cè)坐標(biāo)序列

3.1 AV06測(cè)站N方向的周期項(xiàng)探測(cè)

基于分析年周期的需要，在構(gòu)造軌跡矩陣時(shí)，將嵌入維數(shù)設(shè)置為 1 500。分析特征值的分布情況，發(fā)現(xiàn)有幾組明顯成對(duì)的特征值，分別對(duì)前三對(duì)分量進(jìn)行周期分析。周期探測(cè)結(jié)果如表1所示。

AV06-N主要成對(duì)RC的檢驗(yàn)參數(shù)和周期 表1

綜合三對(duì)RC分量的奇異譜分析結(jié)果，可以得出測(cè)站AV06在北方向上存在顯著的年周期、半年周期信號(hào)。其他的RC成分在分析過(guò)程中并不滿足周期波動(dòng)成分檢驗(yàn)的三條判別標(biāo)準(zhǔn)，暫時(shí)無(wú)法探測(cè)出其他周期成分。

(1)將第3、4主分量進(jìn)行序列重構(gòu)得到年周期項(xiàng)，并與最小二成擬合所得年周期進(jìn)行對(duì)比，如圖3所示。

圖3 AV06-N年周期項(xiàng)對(duì)比圖

其中，藍(lán)色曲線為SSA分析所得年周期，綠色曲線為最小二乘擬合所得年周期信號(hào)。從圖3中可以看出，經(jīng)SSA分析，GPS時(shí)間序列中存在370天的周期項(xiàng)，即年周期項(xiàng)(由于頻率設(shè)置的原因，如果頻率之間的顆粒度更小，可以獲得更精確的周期)。經(jīng)SSA探測(cè)出的周期項(xiàng)與最小二乘擬合的周期項(xiàng)大體一致，但是最小二乘擬合結(jié)果為振幅恒定的周期振幅，而SSA所探測(cè)到的周期項(xiàng)振幅有隨時(shí)間而增大的傾向。這對(duì)于框架的非線性運(yùn)動(dòng)特性研究有很好的指導(dǎo)意義，并且振幅逐年增大的趨勢(shì)也可以反過(guò)來(lái)配合一些地球物理因素的解釋?zhuān)疚木筒蛔鲞M(jìn)一步探討。

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

圖4 AV06-N半年周期項(xiàng)對(duì)比圖

通過(guò)周期圖可以看出，GPS時(shí)間序列中存在175天的周期信號(hào)，即半年周期信號(hào)。半年周期對(duì)比圖與年周期對(duì)比圖相比，區(qū)別要更明顯，這是因?yàn)镾SA探測(cè)出的周期與半年有一定差距，隨著時(shí)間的推移，兩者半年周期信號(hào)開(kāi)始出現(xiàn)峰值偏移。另外，奇異譜分析所得的半年周期信號(hào)的振幅除了有著明顯隨時(shí)間增大的趨勢(shì)，其中也暗含著大約2.5年的周期項(xiàng)。這個(gè)現(xiàn)象有待進(jìn)一步討論。

3.2 AV06測(cè)站E、U方向的周期項(xiàng)探測(cè)

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

AV06-E主要成對(duì)RC的檢驗(yàn)參數(shù)和周期 表2

AV06-U主要成對(duì)RC的檢驗(yàn)參數(shù)和周期 表3

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

綜合測(cè)站AV06三個(gè)坐標(biāo)分量上的SSA周期探測(cè)結(jié)果，有以下兩條結(jié)論：

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

②跟最小二乘擬合結(jié)果有差異。奇異譜探測(cè)到的周期項(xiàng)并非振幅恒定，周期信號(hào)的振幅有隨時(shí)間增大的趨勢(shì)，并且在總體增大趨勢(shì)中還暗含周期變化。

4 周期項(xiàng)影響因素分析

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

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

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

(3)多路徑效應(yīng)。由于GPS星座幾何形狀呈周期性變化，這導(dǎo)致多路徑效應(yīng)也呈周期性變化，而GPS星座周期性變化與計(jì)算GPS測(cè)站坐標(biāo)所基于的太陽(yáng)日之間存在約 247 s的差異，因此會(huì)產(chǎn)生350天的周期項(xiàng)[7]。因?yàn)楸疚闹杏?jì)算頻率的精度在±0.000 1，因此提取出的周期項(xiàng)應(yīng)該在338天～363天之間。

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

關(guān)于所探測(cè)出的周期結(jié)果，本文并未做進(jìn)一步探討，但是可為參考框架點(diǎn)的非線性運(yùn)動(dòng)研究提供一定的數(shù)據(jù)支持。同時(shí)，這樣的周期變化特性又該如何結(jié)合實(shí)際的地球物理因素給出模型化的改正，這也是有待解決的問(wèn)題。