羅 勇 匡翠林 盧辰龍 曾凡河

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙市麓山南路932號(hào)，410083

2 鄭州市市政工程勘測(cè)設(shè)計(jì)研究院，鄭州市民生路1號(hào)，450052

SSA是在Karhumen－Loeve分解理論基礎(chǔ)上逐漸發(fā)展起來的，是對(duì)一維時(shí)間序列進(jìn)行周期震蕩的主成分分析方法，它從時(shí)間序列的動(dòng)力重構(gòu)出發(fā)，并與經(jīng)驗(yàn)正交函數(shù)相關(guān)，不受正弦波假定的約束，從含噪的數(shù)據(jù)中盡可能多地提取可靠信息，將具有顯著震蕩的分量提取出來，并且選擇其中若干有意義的分量進(jìn)行重建，從而降低噪聲［1－2］。本文基于模擬數(shù)據(jù)和IGS跟蹤站BJFS站的實(shí)測(cè)數(shù)據(jù)，采用SSA、db10小波和EMD 方法進(jìn)行對(duì)比分析，并通過精度統(tǒng)計(jì)指標(biāo)比較3種方法在GPS坐標(biāo)序列去噪和周期信號(hào)提取方面的優(yōu)劣勢(shì)。

1 SSA基本原理

SSA 的分析對(duì)象是經(jīng)過中心化后的一維時(shí)間序列x1，x2，x3.....xN，對(duì)其進(jìn)行滯后排列，選擇合適的窗口大小M（M＜N／2），建立一個(gè)滯后矩陣［3］：

求出滯后矩陣Y的協(xié)方差矩陣C（它是一個(gè)類似于Toeplizt的矩陣，其對(duì)角線上的元素都相同）：

再求出C的特征值λ1≥λ2≥…≥λM≥0和特征向量。式就是時(shí)間序列的奇異譜，對(duì)其奇異值進(jìn)行的運(yùn)算稱為奇異譜分析。λk對(duì)應(yīng)的特征向量Ek稱為時(shí)間經(jīng)驗(yàn)正交函數(shù)T－EOF（temporal empirical orthogonal function）［4］。之后計(jì)算滯后序列Y在Ek上的投影：

式中，0≤i≤N－M，1≤j≤M，aik為Ekj所反映的時(shí)間演變權(quán)重，稱為時(shí)間主成分，記為TPC。

SSA 最重要的功能是重建（reconstruction component，RC）［5］。由第k個(gè)時(shí)間主成分（TEOF）和時(shí)間主分量（TPC）重建的成分記為：

根據(jù)SSA 原理，當(dāng)原序列中存在一個(gè)周期振蕩成分時(shí)，SSA 得到一對(duì)接近相等的特征值，對(duì)應(yīng)的一對(duì)TEOF 正交，一對(duì)TPC 正交。為不失一般性，設(shè)這對(duì)特征成分序號(hào)為k和k＋1。滿足上述條件的第k和k＋1的RC 之和是一個(gè)周期振蕩成分。然而，用離散的有限長(zhǎng)度序列作SSA時(shí)，即使是周期信號(hào)，也得不到λk＝λk＋1，而TEOF及TPC之間總是正交的。所以，實(shí)際應(yīng)用中依據(jù)這些條件識(shí)別周期振蕩成分較困難。因此，Vautard［3］和Ghil［6］根據(jù)譜性質(zhì)提出了相應(yīng)的補(bǔ)充判據(jù)。由于SSA 具有穩(wěn)定識(shí)別周期信號(hào)的優(yōu)點(diǎn)，SSA 及MSSA 在趨勢(shì)提取與周期信號(hào)提取中有著廣泛的應(yīng)用［6－7］，在大地測(cè)量領(lǐng)域也有初步的應(yīng)用［8－9］。

2 比較分析方法與精度評(píng)定指標(biāo)

2.1 小波和EMD 基本原理

對(duì)于任何的信號(hào)或者可積函數(shù)f（t），它的連續(xù)小波變換為［10］：

Wf（x，y）是基本小波函數(shù)，x和y分別是伸縮和平移因子，ψ＊（t）是ψ（t）的共軛函數(shù)。y確定小波的中心位置，x確定小波的時(shí)域廣度。因此，在分析時(shí)間序列時(shí)，只要恰當(dāng)?shù)剡x取小波函數(shù)，便能得到較好的頻域和時(shí)域局部性［11］。

經(jīng)驗(yàn)?zāi)Ｊ椒纸猓‥MD）是對(duì)一組時(shí)間序列進(jìn)行平穩(wěn)化處理，將原始信號(hào)分解成一組從高頻到低頻的基本信號(hào)，每個(gè)信號(hào)都代表了一個(gè)本征式的分量，即固有模態(tài)函數(shù)（intrinsic mode function，IMF），IMF要滿足兩個(gè)條件［12］：1）待分析信號(hào)的極值數(shù)與過零值數(shù)相等或最多相差1；2）在任意點(diǎn)上，由局部極值確定的上下包絡(luò)線均值為零。其中最低頻的IMF普遍為原信號(hào)的均值或者趨勢(shì)項(xiàng)，高頻的IMF分量一般為噪聲。計(jì)算式為：

式中，y（t）為原始序列，ai（t）表示獲得的IMF本征分量，n為分量總數(shù)，r（t）為趨勢(shì)分量。

2.2 統(tǒng)計(jì)評(píng)定指標(biāo)

1）均方根誤差（RMSE）為：

其中，N為信號(hào)的長(zhǎng)度；yt是去噪后的信號(hào)；y是原始不含噪聲的信號(hào)；ut為原始含噪信號(hào)；SRMSE為濾波后信號(hào)的RMS值，SRMSE越小說明yt和y越相似；NRMSE 為濾波后噪聲部分的RMS 值，NRMSE越接近噪聲水平，去噪效果越好。

2）互相關(guān)系數(shù)（R）為：

式中，cov（yt，y）是yt和y的協(xié)方差，σyt和σy分別是yt和y的標(biāo)準(zhǔn)差。R越大，表示去噪后的信號(hào)和原始未含噪的參考信號(hào)越相似。

3）信噪比（SNR）為：

3 模擬數(shù)據(jù)分析

構(gòu)造模擬數(shù)據(jù)模型：

模擬數(shù)據(jù)的采樣頻率為1s，樣本大小為3 000，主要由兩個(gè)周期信號(hào)疊加而成，e（t）為服從正態(tài)分布N（0，1.52）的高斯白噪聲，并且加入一個(gè)低頻的趨勢(shì)項(xiàng)。3種方法去噪結(jié)果如圖1所示，從上至下分別是：模擬的原始時(shí)間序列；未加噪聲的模擬時(shí)間序列；去噪后的序列；去噪后的序列與未含噪聲序列的殘差序列；原始時(shí)間序列與去噪后序列的殘差序列。

從圖1對(duì)比可看出，SSA 去噪后的序列與原始未含噪序列的殘差序列比較平穩(wěn)且趨近于零，而小波和EMD 去噪后的殘差序列不平穩(wěn)，且EMD 在端部的殘差比較大。為了進(jìn)一步證明SSA 方法去噪的有效性，采用3種評(píng)價(jià)指標(biāo)進(jìn)行定量比較分析。從表1中可看出，3 種方法去噪后的信號(hào)與真實(shí)信號(hào)的相關(guān)系數(shù)都達(dá)到了0.99以上，說明3種方法均是有效濾波法。從濾波后信號(hào)的SRMSE值可以看出，SSA 方法最小，且其濾波后噪聲的NRMSE 更加接近噪聲水平，SSA方法去噪后序列的信噪比也略高于小波和EMD，以上都說明SSA 方法去噪性能的優(yōu)越性。

表1 3種方法模擬去噪效果的定量分析Tab.1 Quantitative analysis of the simulation denoise results of three kinds of method

圖1 SSA、db10小波和EMD 模擬去噪結(jié)果對(duì)比Fig.1 Comparison between the simulation denoise results of SSA，db10wavelet and EMD

模擬數(shù)據(jù)有兩個(gè)周期函數(shù)，分別是：

圖2是3種方法分解出來的周期信號(hào)分別與原始真實(shí)信號(hào)s1和s2進(jìn)行的對(duì)比?？傮w上，3種方法均能有效地從原始含噪信號(hào)中分離出s1和s2信號(hào)。但SSA 分解出來的s1和s2信號(hào)與原始真實(shí)信號(hào)的殘差都比較平穩(wěn)且趨于零；而通過小波和EMD 方法得到的殘差序列不平穩(wěn)且存在端部效應(yīng)，說明SSA 方法提取周期項(xiàng)的性能要優(yōu)于小波分析與EMD 方法。從表2 也可看出，用SSA 方法提取的信號(hào)s1和s2與對(duì)應(yīng)參考信號(hào)的RMSE要遠(yuǎn)低于db10小波和EMD，而互相關(guān)系數(shù)R大于db10小波和EMD，以上都表明SSA 方法對(duì)周期性信號(hào)提取的優(yōu)越性。

表2 3種方法分解的信號(hào)與對(duì)應(yīng)參考信號(hào)的RMSE（mm）和互相關(guān)系數(shù)RTab.2 RMSE（mm）and correlation coefficients of decomposed signal and the reference signal

圖2 SSA、db10小波和EMD 分解的s1信號(hào)和s2信號(hào)Fig.2 Signal s1and s2derived from decomposition of SSA，db10wavelet and EMD

4 實(shí)例分析

圖3是采用SSA、db10小波和EMD 三種方法分別對(duì)我國(guó)IGS站BJFS站的高程序列進(jìn)行濾波去噪的結(jié)果。3種方法去噪后的序列都比較平滑，均能有效剔除序列中每個(gè)時(shí)間點(diǎn)的高頻噪聲，但SSA 去噪的序列和db10小波去噪后的序列相關(guān)性更高，而EMD 方法由于分解不穩(wěn)定，存在過度去噪的現(xiàn)象，如圖3中紅色橢圓標(biāo)記處。表3是3種方法去噪后的信號(hào)互相關(guān)系數(shù)比較，SSA與小波去噪后的相關(guān)系數(shù)達(dá)到0.995，高于SSA與EMD、db10 小波與EMD 的相關(guān)系數(shù)，說明SSA 與db10小波去噪后的序列更吻合。SSA 與原始時(shí)間序列、db10小波與原始時(shí)間序列的互相關(guān)系數(shù)均為0.868，稍高于EMD 的0.857，說明SSA、db10小波去噪后的信號(hào)與原始信號(hào)更接近。SSA、db10小波和EMD 去噪后的信號(hào)與噪聲的信噪比分別為11.118、11.133 和9.863，EMD 的信噪比最低。從以上分析可看出，前兩種方法的去噪效果相當(dāng)，且略優(yōu)于EMD。

圖4是3種方法進(jìn)行多尺度分解后的結(jié)果，由于篇幅限制，只給出db10 小波分解得到的d7～d 1 0分量和EMD的IMF 7～I(xiàn)MF 1 0分量。SSA 方法分解的信號(hào)從上往下分別為0.25a、0.5a、1a和似2a周期信號(hào)，結(jié)果與文獻(xiàn)［1］一致；db10 小波分解的信號(hào)從上往下分別為：0.5 a、1a、似2a周期信號(hào)和周期更長(zhǎng)的信號(hào)；EMD分解的的信號(hào)從上往下分別為：0.5a、1a、似2a周期信號(hào)和周期更長(zhǎng)的信號(hào)。從圖可見，SSA 方法分解得到的0.25a、0.5a和1a周期信號(hào)在其對(duì)應(yīng)的頻率存在明顯且單一的峰值，而另兩種方法僅1a周期信號(hào)比較明顯，其余信號(hào)則存在周期混疊現(xiàn)象。總之，3種方法均能分解出季節(jié)性信號(hào)，從貢獻(xiàn)率大小看，1a周期信號(hào)所占比重最大，0.5a周期信號(hào)次之；從分解效果看，SSA 方法提取的季節(jié)信號(hào)更加穩(wěn)定平滑，且信號(hào)混疊現(xiàn)象不明顯。

表4是3種方法對(duì)分解的年周期信號(hào)進(jìn)行的相關(guān)性分析，SSA 提取的年周期信號(hào)與db10 小波提取的年周期信號(hào)之間的互相關(guān)系數(shù)達(dá)到了0.96，而SSA 和db10 小波與EMD 之間的互相關(guān)系數(shù)分別為0.89和0.86，說明SSA 提取的年周期信號(hào)和db10小波提取的年周期項(xiàng)信號(hào)更接近，且SSA 提取的年周期信號(hào)與原始時(shí)間序列的相關(guān)性要高于db10小波與EMD 的相關(guān)性，這表明SSA 提取的年周期項(xiàng)信號(hào)更接近于原始時(shí)間序列，也間接說明了SSA 方法提取年周期項(xiàng)的優(yōu)越性。

圖3 3種方法的BJFS站高程序列去噪結(jié)果Fig.3 The denoise results of BJFS station of three kinds of method

圖4 3種方法的多尺度分解結(jié)果Fig.4 Multi－scale decomposition results of three kinds of method

表3 SSA、db10小波和EMD去噪后信號(hào)的互相關(guān)系數(shù)RTab.3 Correlation coefficients of three kinds of method denoising for BJFS station

表4 SSA、db10小波和EMD分解的1a周期信號(hào)的互相關(guān)系數(shù)RTab.4 Correlation coefficients of three kinds of method decomposing annual cycle

5 結(jié) 語

奇異譜分析是主成分分析方法的一種，主要對(duì)一維時(shí)間序列進(jìn)行分析，適用于從短時(shí)間序列中提取信息?；诖朔椒ǎ疚膶?duì)模擬數(shù)據(jù)和GPS實(shí)測(cè)數(shù)據(jù)進(jìn)行去噪和周期項(xiàng)提取的研究，并與小波和EMD 方法對(duì)比分析。結(jié)果表明，3種方法都是有效的去噪方法，且均能較好地分解出季節(jié)性信號(hào)，但SSA 方法去噪效果更好，分解得到的周期信號(hào)更加穩(wěn)定和平滑，且不存在信號(hào)混疊現(xiàn)象。

［1］王解先，連麗珍，沈云中.奇異譜分析在GPS站坐標(biāo)監(jiān)測(cè)序列分析中的應(yīng)用［J］.同濟(jì)大學(xué)學(xué)報(bào)：自然科學(xué)版，2013（2）：282－288（Wang Jiexian，Lian Lizhen，Shen Yunzhong.Application of Singular Spectral Analysis to GPS Station Coordinate Monitoring Series［J］.Journal of Tongji University（Natural Science），2013（2）：282－288）

［2］徐克紅，程鵬飛，文漢江.太陽黑子數(shù)時(shí)間序列的奇異譜分析和小波分析［J］.測(cè)繪科學(xué)，2007，32（6）：35－38（Xu Kehong，Cheng Pengfei，Wen Hanjiang.Singular Spectrum Analysis and Wavelet Analysis on Time Series of Sunspot［J］.Science of Surveying and Mapping，2007，32（6）：35－38）

［3］Vautard R，Yiou P，Ghil M.Singular－Spectrum Analysis：A Toolkit for Short，Noisy Chaotic Signals［J］.Physica D：Nonlinear Phenomena，1992，58（1）：95－126

［4］江志紅，丁裕國(guó).奇異譜分析的廣義性及其應(yīng)用特色［J］.氣象學(xué)報(bào)，1998，56（6）：736－745（Jiang Zhihong，Ding Yuguo.Generality and Applied Features for Singular Spectrum Analysis［J］.Acta Meteorologica Sinica，1998，56（6）：736－745）

［5］吳洪寶.奇異譜分析——最大熵預(yù)報(bào)方法［J］.甘肅氣象，2000，18（1）：1－5（Wu Hongbao.Singular Spectral Analysis－Maximum Entropy Forecast Method［J］.Gansu Meteorology，2000，18（1）：1－5）

［6］Ghil M，Vautard R.Interdecadal Oscillations and the Warming Trend in Global Temperature Time Series［J］.Nature，1991，350（6316）：324－327

［7］Jevrejeva S，Moore J C.Singular Spectrum Analysis of Baltic Sea ice Conditions and Large‐Scale Atmospheric Patterns Since 1708［J］.Geophysical Research Letters，2001，28（23）：4 503－4 506

［8］Yoo J C，Dorico P.Trends and Fluctuations in the Dates of Ice Break－up of Lakes and Rivers in Northern Europe：the Effect of the North Atlantic Oscillation［J］.Journal of Hydrology，2002，268（1）：100－112

［9］Rangelova E，Wal W，Sideris M G，et al.Spatiotemporal Analysis of the GRACE－derived Mass Variations in North America by Means of Multi－Channel Singular Spectrum Analysis［M］.Springer Berlin Heidelberg，2010

［10］黨星海，趙麗潔，孔令杰，等.小波分析在GPS 振動(dòng)監(jiān)測(cè)數(shù)據(jù)中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2013，33（2）：147－150（Dang Xinghai，Zhao Lijie，Kong Lingjie，et al.Application of Wavelet Analysis in GPS Dynamic Deformation Data Processing［J］.Journal of Geodesy and Geodynamics，2013，33（2）：147－150）

［11］羅飛雪，戴吾蛟.小波分解與EMD 在變形監(jiān)測(cè)應(yīng)用中的比較［J］.大地測(cè)量與地球動(dòng)力學(xué)，2010，30（3）：137－141（Luo Feixue，Dai Wujiao.Comparison of EMD with Wavelet Decomposition for Dynamic Deformation Monitoring［J］.Journal of Geodesy and Geodynamics，2010，30（3）：137－141）

［12］戴吾蛟，丁曉利，朱建軍，等.基于經(jīng)驗(yàn)?zāi)Ｊ椒纸獾臑V波去噪法及其在GPS 多路徑效應(yīng)中的應(yīng)用［J］.測(cè)繪學(xué)報(bào)，2007，35（4）：321－327（Dai Wujiao，Ding Xiaoli，Zhu Jianjun，et al.EMD Filter Method and Its Application in GPS Multipath［J］.Acta Geodaetica et Cartographica Sinica，2007，35（4）：321－327）