利用奇異譜分析探測(cè)GPS 坐標(biāo)時(shí)間序列時(shí)變周期信號(hào)

劉磊

（安徽省基礎(chǔ)測(cè)繪信息中心 安徽合肥 230031）

1 引言

GPS 基準(zhǔn)站坐標(biāo)時(shí)間序列廣泛應(yīng)用于參考框架的建立與維持、地殼形變監(jiān)測(cè)等高精度大地測(cè)量和地球動(dòng)力學(xué)研究領(lǐng)域。幾乎所有的GPS 基準(zhǔn)站都呈現(xiàn)出顯著的非線性運(yùn)動(dòng)變化（尤其是高程方向的周年、半周年特征）。隨著地學(xué)等研究領(lǐng)域?qū)Υ蟮販y(cè)量成果所要求的精度越來(lái)越高，GPS 坐標(biāo)時(shí)間序列中的非線性變化越來(lái)越受到關(guān)注[1]。

針對(duì)坐標(biāo)時(shí)間序列非線性運(yùn)動(dòng)分析，眾多學(xué)者利用最小二乘擬合（Least Squares Fitting, LSF）方法研究了基準(zhǔn)站的非線性運(yùn)動(dòng)特征。然而，基于LSF 獲取的基準(zhǔn)站周期振幅和相位信息反映的是恒定的周期特征，無(wú)法顧及季節(jié)性信號(hào)的年際變化特征。為此，有關(guān)學(xué)者提出了利用卡爾曼濾波、小波分析等方法探測(cè)GPS 坐標(biāo)時(shí)間序列中的年際非線性信號(hào)，取得了一定的成果。

奇異譜分析（Singular Spectrum Analysis, SSA）是用于研究非線性時(shí)間序列的有效方法，在大地測(cè)量學(xué)、氣象學(xué)，海洋學(xué)等領(lǐng)域有著廣泛的應(yīng)用，可以直接從短而有噪聲的時(shí)間序列中提取有效信息，不會(huì)被正弦波假定約束，更適合于時(shí)變周期的提取[2]。針對(duì)利用SSA 探測(cè)GPS 序列中的非線性信號(hào)，國(guó)內(nèi)外學(xué)者也進(jìn)行了一定的研究。王解先等利用SSA 對(duì)BJSH 站坐標(biāo)時(shí)間序列缺失的數(shù)據(jù)進(jìn)行插值，并分析了序列中的趨勢(shì)項(xiàng)和周期項(xiàng)[3]。Chen 等利用SSA 提取GPS 坐標(biāo)時(shí)間序列中的周期信號(hào)，通過(guò)多種方法的比較分析，驗(yàn)證了SSA 在提取時(shí)變季節(jié)信號(hào)中的價(jià)值[4]。周茂盛等提出了利用多通道奇異譜分析（MSSA）提取共模誤差的新思路，該方法提高了坐標(biāo)時(shí)間序列的精度[5]。

本文詳細(xì)闡述了SSA 方法的基本原理，結(jié)合GPS 時(shí)間序列分析的實(shí)際，分析了SSA 中的周期信號(hào)分組與重構(gòu)方法。進(jìn)而以BJSH 站為例，將其應(yīng)用于GPS 基準(zhǔn)站坐標(biāo)時(shí)間序列的分析，并與最小二乘擬合結(jié)果比較，探討基準(zhǔn)站的三維坐標(biāo)周期變化特征。

2 原理與方法

2.1 奇異譜分析

奇異譜分析是對(duì)維度為N 的GPS 坐標(biāo)時(shí)間序列C=（C1，C2，C3…CN） 進(jìn)行分析并轉(zhuǎn)化成多維坐標(biāo)時(shí)間序列，從而將原始GPS 坐標(biāo)時(shí)間序列進(jìn)行分解并提取其中可靠信息。SSA 具體步驟主要分為四個(gè)部分：

1）嵌入

首先，將原始GPS 坐標(biāo)時(shí)間序列C 轉(zhuǎn)化成L×K的軌跡矩陣，其中L 為窗口長(zhǎng)度，K=N-L+1（1

其中Cij=Ci+j-1 即反對(duì)角線上的元素都相等。

2）奇異值分解

3）特征三要素分組

將矩陣下標(biāo){1,2,3…d}分成m 個(gè)子集I1，I2，…Im，設(shè)I={i1，i2，…ip}，得到合成矩陣CI=Ci1+Ci2+…+Cip。計(jì)算集合I=I1，I2，…，Im的每個(gè)合成矩陣，則分解后的序列為[6]：

4）重構(gòu)GPS 坐標(biāo)時(shí)間序列

將分組得到的矩陣轉(zhuǎn)換成長(zhǎng)度為N 的新序列，令X 為L(zhǎng)×K 的矩陣，其中K=N-L+1，設(shè)L*=min{L，K}，K*=max{L，K}，根據(jù)對(duì)角平均化公式將矩陣X 轉(zhuǎn)化為X1，X2…XN的時(shí)間序列，Xi,j為矩陣中的各元素（1≤i≤L，1≤j≤K）。對(duì)角平均公式為：

2.2 基于SSA 的周期信號(hào)分離方法

當(dāng)原始GPS 坐標(biāo)時(shí)間序列C1，C2，C3…CN中存在一個(gè)周期信號(hào)時(shí)，可以得到一對(duì)近似相等的特征值[6]。令這對(duì)特征值在序列中的位置為k 和k+1，那么滿足所述條件的第k 個(gè)和第k+1 個(gè)RC（重構(gòu)成分）之和構(gòu)成了原始GPS 坐標(biāo)時(shí)間序列的一個(gè)周期信號(hào)。

由于時(shí)間序列往往是離散的，根據(jù)上述條件確定周期信號(hào)具有一定缺陷。因此，根據(jù)Vautard 等的研究，補(bǔ)充三個(gè)準(zhǔn)則來(lái)識(shí)別周期信號(hào)[7]：

（1） 兩個(gè)連續(xù)特征值近似相等，且在序列中足夠大；

（2） 對(duì)應(yīng)的兩個(gè)T-EOF 必須頻率相近，對(duì)TEOFk和T-EOFk+1傅里葉變換后得到Ek（f）和Ek+1（f），找出達(dá)到最大值的|Ek（f）|2和|Ek+1（f）|2所對(duì)應(yīng)的頻率fk和fk+1，其中，兩頻率的差值δfk=|fk+1-fk|應(yīng)很小。

（3） |Ek（f）|2和|Ek+1（f）|2足夠大，也就是說(shuō)介于fk和fk+1之間的一個(gè)頻率f*可以通過(guò)這一對(duì)特征成分反映出來(lái)。

說(shuō)明原始序列的中間頻率f^*的周期性變化至少有2/3 可通過(guò)這對(duì)特征成分所表示。

2.3 最小二乘擬合

GPS 基準(zhǔn)站坐標(biāo)時(shí)間序列可以表示為：

式中，ti是以一年為單位的觀測(cè)值i 的時(shí)間，Y0是一個(gè)恒定的初始偏移量，v 是一個(gè)恒定的速度，ak和bk是周期項(xiàng)的系數(shù)，fk是周期頻率，（ti） 是噪音的含量。如果將線性趨勢(shì)與年周期（f1=1 cpy）和半年周期（f2=2 cpy）信號(hào)一起分解，那么矩陣A 變?yōu)椋?/p>

未知向量x 是

最小二乘的主要優(yōu)點(diǎn)是易于實(shí)現(xiàn)，并能直觀地解釋線性趨勢(shì)和季節(jié)振幅的估計(jì)，但是最小二乘只能獲得恒定的振幅或相位，其缺點(diǎn)在于長(zhǎng)周期的變化可能被誤認(rèn)為是一種線性趨勢(shì)。

3 結(jié)果與分析

本文選擇陸態(tài)網(wǎng)提供的基準(zhǔn)站坐標(biāo)時(shí)間序列（以BJSH 站為例），利用SSA 對(duì)其進(jìn)行周期性運(yùn)動(dòng)分析。考慮到本文研究的是時(shí)間序列的周期變化特征，選擇去趨勢(shì)的BJSH 站結(jié)果進(jìn)行分析。此外，由于數(shù)據(jù)存在缺失，本文采用三次多項(xiàng)式插值得到BJSH 站2000 年-2017 年連續(xù)的單天坐標(biāo)時(shí)間序列。

3.1 SSA 窗口長(zhǎng)度L 的選擇

應(yīng)用SSA 分析GPS 坐標(biāo)時(shí)間序列時(shí)，窗口長(zhǎng)度L 的選擇對(duì)于獲取正確的結(jié)果至關(guān)重要。L 越大，原始序列的分解越精細(xì)，但過(guò)大的L 不利于奇異值分解的運(yùn)算；如果 L 太小，則奇異值分解將導(dǎo)致不同的成分相互混雜，不利于信號(hào)的提取。Vautard 等研究發(fā)現(xiàn)[8]，SSA 能夠分解出周期介于L/5 和L 之間的信號(hào)。Ghil 等提出，在選擇窗口滯后時(shí)間時(shí)，應(yīng)兼顧兩個(gè)因素[9]：所提取的信息量與對(duì)該信息的統(tǒng)計(jì)置信度。前者要求窗口滯后應(yīng)盡可能寬，即一個(gè)較大的L，而后者則要求盡可能多地重復(fù)感興趣的特征，即盡可能大的N/L 比。文獻(xiàn)[3]通過(guò)模擬實(shí)驗(yàn)也證明了從GPS 坐標(biāo)時(shí)間序列中提取年度和半年周期時(shí)，選擇2 年或3 年的滯后窗口是合適的。因此，本文在利用SSA 研究GPS 坐標(biāo)時(shí)間序列的年變化和半年變化時(shí)，選擇兩年的窗口長(zhǎng)度將GPS 時(shí)間序列分解成730 種時(shí)間模式，進(jìn)而獲取其年周期和半年周期信號(hào)特征。

3.2 周期項(xiàng)探測(cè)

3.2.1 SSA 分析結(jié)果

圖1（左）為BJSH 站高程方向經(jīng)SSA 分析時(shí)數(shù)據(jù)方差各部分的歸一化特征值按降序排列結(jié)果。由圖可知，前兩個(gè)最大的特征值分別可以解釋數(shù)據(jù)方差的34.21%和6.13%。圖1（右）顯示了前50 個(gè)相應(yīng)的特征值EOFs 對(duì)應(yīng)的主周期頻率，從圖中可以清楚地看出年周期的特征值和半年周期的特征值。取前兩對(duì)數(shù)值相近的特征值，將其對(duì)應(yīng)的重構(gòu)分量分為兩對(duì)，如圖2 圖3 顯示了特征值對(duì)應(yīng)的兩對(duì)經(jīng)驗(yàn)正交函數(shù)（EOFs）和相應(yīng)的重構(gòu)分量（RCs），分別代表年周期信號(hào)和半年周期信號(hào)。圖4 描述了原始GPS基準(zhǔn)站高程時(shí)間序列、利用SSA 提取的年周期和半年周期。由圖可知，BJSH 站的周期振幅隨時(shí)間變化而不同，具有顯著的時(shí)變特征，而SSA 能夠有效地提取GPS 坐標(biāo)時(shí)間序列中的時(shí)變周期。

圖1 （左）GPS 高程數(shù)據(jù)協(xié)方差矩陣得到的歸一化特征值；（右）特征值EOFs 對(duì)應(yīng)的主周期頻率

圖2 特征值對(duì)應(yīng)的兩對(duì)EOFs 年周期和半年周期

圖3 特征值對(duì)應(yīng)的兩對(duì)RCS 年周期和半年周期

圖4 GPS 高程坐標(biāo)時(shí)間序列、重建年信號(hào)、重建半年信號(hào)

3.2.2 測(cè)站N，E，U 三個(gè)方向上的周期項(xiàng)探測(cè)

根據(jù)2.2，當(dāng)存在一對(duì)近似相等的特征值時(shí)，符合補(bǔ)充準(zhǔn)則的一對(duì)特征成分的兩個(gè)RC 之和表示原始坐標(biāo)時(shí)間序列中的一個(gè)周期項(xiàng)。表1～表3 分別給出了U、N、E 三個(gè)方向上主要成對(duì)RC 的周期項(xiàng)分析結(jié)果，其中K 為選擇的重構(gòu)分量，λk為對(duì)應(yīng)的特征值，f 為|E（f）|2最大值對(duì)應(yīng)的頻率，|E（f）|2為TEOFk經(jīng)過(guò)傅里葉變換后的值，f*為一對(duì)特征值對(duì)應(yīng)的頻率的均值。

表1 U 方向上主要成對(duì)RC 的周期項(xiàng)分析

根據(jù)表中的三對(duì)RC 的分析結(jié)果，發(fā)現(xiàn)U 方向上前兩個(gè)特征值較大的重構(gòu)分量RC 對(duì)應(yīng)年周期，第二對(duì)RC 對(duì)應(yīng)半年周期，而其余的RC 則不能滿足上述的檢驗(yàn)標(biāo)準(zhǔn)，無(wú)法探測(cè)出周期項(xiàng)信號(hào)。

采用同樣的方法對(duì)BJSH 站N、E 分量的時(shí)間序列分別進(jìn)行周期項(xiàng)探測(cè)（表2、表3）。結(jié)果表明N 方向上無(wú)法探測(cè)出1/3 年周期，而且不同方向上選取的RC 也不同，但前幾個(gè)較大的特征值對(duì)應(yīng)的重構(gòu)分量都對(duì)應(yīng)著年周期、半年周期。

表2 N 方向上主要成對(duì)RC 的周期項(xiàng)分析

表3 E 方向上主要成對(duì)RC 的周期項(xiàng)分析

綜合BJSH 站三個(gè)方向上的周期項(xiàng)探測(cè)，可以得出以下結(jié)論：BJSH 站在不同方向都近似存在年周期、半年周期項(xiàng)，但不同方向上存在的周期項(xiàng)也有略微不同，比如1/3 年周期等。

3.3 SSA 與LSF 比較分析

本文進(jìn)一步分析了SSA 與LSF 提取GPS 坐標(biāo)時(shí)間序列周期特征的差異，如圖5 所示。

圖5 U 方向上LSF 和SSA 對(duì)比圖，（a）SSA 與LSF 年周期對(duì)比圖；（b）SSA 與LSF 半年周期對(duì)比圖；（c）SSA 與LSF 主周期對(duì)比圖

（1）SSA 得到的年周期與LSF 得到的年周期作對(duì)比，如圖5（a）所示。從圖中可以看出，SSA 與LSF都具有提取周期項(xiàng)的能力。并且可以清楚地觀察到，2006-2011 年間SSA 擬合的波峰比LSF 擬合的結(jié)果更接近原始序列的波峰，2012-2016 年間SSA 擬合的波谷比LSF 擬合的結(jié)果更接近原始序列的波谷。

（2）SSA 得到的半年周期與LSF 作對(duì)比，如圖5（b）所示。圖5（b）比圖5（a）的區(qū)別要更明顯，隨著時(shí)間的推移，SSA 與LSF 得到的峰值產(chǎn)生明顯的偏移。

（3）比較SSA 得到的主周期信號(hào)（年周期加半年周期）與最小二乘擬合得到的主周期信號(hào)（年周期加半年周期），如圖5（c）所示。從圖中可以看出，在最小二乘擬合中，季節(jié)性信號(hào)在每年的5 月中旬達(dá)到峰值；在SSA 中，獲得的季節(jié)性信號(hào)隨時(shí)間而變化。

綜合SSA 與LSF 比較的結(jié)果發(fā)現(xiàn)，SSA 得到的周期信號(hào)能隨原始坐標(biāo)時(shí)間序列的變化而變化，而最小二乘擬合的結(jié)果只能得到振幅和相位恒定的周期信號(hào)，因此SSA 在提取GPS 坐標(biāo)時(shí)間序列非線性信號(hào)上有明顯的優(yōu)勢(shì)。

4 結(jié)束語(yǔ)

本文詳細(xì)研究了SSA 方法應(yīng)用于GPS 坐標(biāo)時(shí)間序列分析的基本原理，通過(guò)SSA 對(duì)BJSH 站GPS坐標(biāo)時(shí)間序列的處理和分析，研究了GPS 基準(zhǔn)站的周期項(xiàng)特征。與LSF 結(jié)果比較發(fā)現(xiàn)，SSA 提取的周期項(xiàng)與原始信號(hào)的周期項(xiàng)吻合程度優(yōu)于LSF 估計(jì)的結(jié)果，尤其是對(duì)于年際變化特征而言，SSA 顯著優(yōu)于LSF。根據(jù)本文SSA 對(duì)BJSH 站的周期項(xiàng)探測(cè)發(fā)現(xiàn)，BJSH 站在U、N、E 三個(gè)方向上存在顯著的周期變化特征，但各方向上探測(cè)的周期項(xiàng)存在差異。BJSH 站的周期項(xiàng)并不是嚴(yán)格的年周期、半年周期。本文沒(méi)有深入分析影響周期項(xiàng)的主要因素，如何根據(jù)影響周期項(xiàng)的因素給出解決問(wèn)題的辦法有待進(jìn)一步研究。