羅 勇 匡翠林 盧辰龍 曾凡河

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

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

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

1 SSA基本原理

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

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

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

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

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

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

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

2.1 小波和EMD 基本原理

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

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

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

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

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

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

其中，N為信號的長度；yt是去噪后的信號；y是原始不含噪聲的信號；ut為原始含噪信號；SRMSE為濾波后信號的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越大，表示去噪后的信號和原始未含噪的參考信號越相似。

3）信噪比（SNR）為：

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

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

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

從圖1對比可看出，SSA 去噪后的序列與原始未含噪序列的殘差序列比較平穩(wěn)且趨近于零，而小波和EMD 去噪后的殘差序列不平穩(wěn)，且EMD 在端部的殘差比較大。為了進(jìn)一步證明SSA 方法去噪的有效性，采用3種評價指標(biāo)進(jìn)行定量比較分析。從表1中可看出，3 種方法去噪后的信號與真實(shí)信號的相關(guān)系數(shù)都達(dá)到了0.99以上，說明3種方法均是有效濾波法。從濾波后信號的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é)果對比Fig.1 Comparison between the simulation denoise results of SSA，db10wavelet and EMD

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

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

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

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

4 實(shí)例分析

圖3是采用SSA、db10小波和EMD 三種方法分別對我國IGS站BJFS站的高程序列進(jìn)行濾波去噪的結(jié)果。3種方法去噪后的序列都比較平滑，均能有效剔除序列中每個時間點(diǎn)的高頻噪聲，但SSA 去噪的序列和db10小波去噪后的序列相關(guān)性更高，而EMD 方法由于分解不穩(wěn)定，存在過度去噪的現(xiàn)象，如圖3中紅色橢圓標(biāo)記處。表3是3種方法去噪后的信號互相關(guān)系數(shù)比較，SSA與小波去噪后的相關(guān)系數(shù)達(dá)到0.995，高于SSA與EMD、db10 小波與EMD 的相關(guān)系數(shù)，說明SSA 與db10小波去噪后的序列更吻合。SSA 與原始時間序列、db10小波與原始時間序列的互相關(guān)系數(shù)均為0.868，稍高于EMD 的0.857，說明SSA、db10小波去噪后的信號與原始信號更接近。SSA、db10小波和EMD 去噪后的信號與噪聲的信噪比分別為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 方法分解的信號從上往下分別為0.25a、0.5a、1a和似2a周期信號，結(jié)果與文獻(xiàn)［1］一致；db10 小波分解的信號從上往下分別為：0.5 a、1a、似2a周期信號和周期更長的信號；EMD分解的的信號從上往下分別為：0.5a、1a、似2a周期信號和周期更長的信號。從圖可見，SSA 方法分解得到的0.25a、0.5a和1a周期信號在其對應(yīng)的頻率存在明顯且單一的峰值，而另兩種方法僅1a周期信號比較明顯，其余信號則存在周期混疊現(xiàn)象?？傊?種方法均能分解出季節(jié)性信號，從貢獻(xiàn)率大小看，1a周期信號所占比重最大，0.5a周期信號次之；從分解效果看，SSA 方法提取的季節(jié)信號更加穩(wěn)定平滑，且信號混疊現(xiàn)象不明顯。

表4是3種方法對分解的年周期信號進(jìn)行的相關(guān)性分析，SSA 提取的年周期信號與db10 小波提取的年周期信號之間的互相關(guān)系數(shù)達(dá)到了0.96，而SSA 和db10 小波與EMD 之間的互相關(guān)系數(shù)分別為0.89和0.86，說明SSA 提取的年周期信號和db10小波提取的年周期項(xiàng)信號更接近，且SSA 提取的年周期信號與原始時間序列的相關(guān)性要高于db10小波與EMD 的相關(guān)性，這表明SSA 提取的年周期項(xiàng)信號更接近于原始時間序列，也間接說明了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去噪后信號的互相關(guān)系數(shù)RTab.3 Correlation coefficients of three kinds of method denoising for BJFS station

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

5 結(jié) 語

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

［1］王解先，連麗珍，沈云中.奇異譜分析在GPS站坐標(biāo)監(jiān)測序列分析中的應(yīng)用［J］.同濟(jì)大學(xué)學(xué)報：自然科學(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ù)時間序列的奇異譜分析和小波分析［J］.測繪科學(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］江志紅，丁裕國.奇異譜分析的廣義性及其應(yīng)用特色［J］.氣象學(xué)報，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ù)報方法［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 振動監(jiān)測數(shù)據(jù)中的應(yīng)用［J］.大地測量與地球動力學(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)測應(yīng)用中的比較［J］.大地測量與地球動力學(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］.測繪學(xué)報，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）