湯文娟

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

1 引 言

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

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

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

2.1 奇異譜分析的基本原理

奇異譜分析(Singular-Spectrum Analysis，SSA)是時(shí)間序列中常用的分析與預(yù)測技術(shù)，組合了經(jīng)典時(shí)間序列分析、多元幾何、多元統(tǒng)計(jì)、動(dòng)態(tài)系統(tǒng)和信號處理等多種元素。奇異譜分析將原始序列分解為緩慢變化趨勢、周期項(xiàng)、噪聲序列，基于時(shí)間序列的動(dòng)力重構(gòu)出發(fā)、與經(jīng)驗(yàn)正交函數(shù)相聯(lián)系，從事先未知物理本質(zhì)的、包含噪聲的有限長觀測序列中，濾去非周期性的異?，F(xiàn)象，對方差譜信號有強(qiáng)化和放大，將頻域信號分解為時(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)參考框架。對于長度N(N>2)的時(shí)間序列X=XN=(x1，…，xN)，取窗體長度為L(1

(1)嵌入(Embedding)

嵌入過程是將一個(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，對應(yīng)的特征向量為U1，…，UL。記Vi=XTUi(i=1，…，L)。其中Ui是矩陣的左特征向量，又稱為時(shí)間經(jīng)驗(yàn)正交函數(shù)(Temporal Empirical Orthogonal Function，T-EOF)；Vi是矩陣的右特征向量，又稱為時(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è)合成矩陣，那么序列的分解式可寫為：

X=XI1+XI2+…+XIp

(3)

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

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

奇異譜分析將序列分離成各種單一頻率的振蕩信號以及噪聲項(xiàng)，利用時(shí)間主成分和時(shí)間經(jīng)驗(yàn)正交函數(shù)展開各分量，可以獲得對應(yīng)于各頻率振蕩信號的重構(gòu)分量序列，從而達(dá)到提取有用信息、過濾噪聲的目的。將分組后的矩陣XIj重構(gòu)為長度為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ù)對角平均公式將矩陣Y轉(zhuǎn)換為y1，…，yN的時(shí)間序列，對角平均公式為：

(4)

(5)

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

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

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

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

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

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

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

圖2 AV06臺站日觀測坐標(biāo)序列

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

通過奇異譜分析方法提取出的周期項(xiàng)可知，不同測站的不同方向上的周期項(xiàng)也不盡相同，也不是嚴(yán)格的周年、半周年周期項(xiàng)，這些相近的周期項(xiàng)可能由同一個(gè)物理因素所引起，由于受分析方法的分辨率和定位精度所限，而表現(xiàn)出略微的差異[5]。下面對于周期性變化的幾個(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)，而且不同測站之間的周期項(xiàng)具體周期值也互不相同。這種情況是由地殼運(yùn)動(dòng)的復(fù)雜性所致，也是嚴(yán)格的地殼周期運(yùn)動(dòng)與局部人類活動(dòng)疊加的結(jié)果。

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

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

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

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