一種IGS站高程時間序列分析方法

張恒璟，陸帝，文漢江，程鵬飛，崔東東

(1.遼寧工程技術大學 測繪與地理科學學院，遼寧 阜新 123000；2.中國測繪科學研究院，北京 100036)

0 引言

對IGS站坐標時間序列的分析表明，在水平方向上時間序列主要表現出線性運動的特征，在高程方向則存在著復雜的周期性運動和非線性運動，這種運動特征在不同的IGS站上有著不同的表現[1-2]。對于IGS站高程時間序列的非平穩非線性信號，在使用周期擬合模型去除序列的非平穩性時，難免會使序列的分析結果出現錯誤。

在處理與分析非線性非平穩的信號時，希爾伯特-黃變換[3]通過使用自適應的方法避免了以往方法的缺點。經驗模態分解方法(empirical mode decomposition，EMD)分解，會因為部分信號中斷而導致混疊模態分量問題[4-5]；整體經驗模態分解方法(ensemble empirical mode decomposition，EEMD)通過引入大量的白噪聲來協助分析時間序列，利用整體平均的思想很好地解決了混疊模態的問題[6-7]。然而，引入大量白噪聲協助分析，使得求整體均值的過程繁瑣，導致了計算困難和費時的問題；同時，引入大量噪聲對原始序列存在一定的破壞，且引入的噪聲經處理后會有殘留。本文引入互補集合經驗模態分解方法(complementary ensemble empirical mode decomposition，CEEMD)，該方法通過引入獨立分布的N對正、負輔助白噪聲[8]，減少EEMD方法在信號重構時導致的噪聲殘留量，提高信號重構的效果。

1 EEMD和CEEMD分解

1.1 EEMD分解

EEMD算法解決了EMD方法中的因信號中斷導致的混疊模態分量問題。EEMD分解算法是通過將高斯白噪聲均勻地添加到時間序列信號中去，在有限次添加白噪聲之后進行EMD分解。參考文獻[6-7]的分解過程如下；

(1)

②假設高斯白噪聲的添加次數為M，對經過M次噪聲添加并經EMD分解過后的IMF分量進行求整體均值，作為EEMD分解的IMF分量

(2)

③整體均值計算后得到各階IMF分量及剩余量，原信號即可寫為：

(3)

EEMD在進行整體平均時，默認每次經噪聲添加再分解所得到的都是等幅度的IMF分量，但由于實際信號中存在著一定的間歇性高頻分量，導致混疊高頻分量信號，這樣得到的EEMD分解結果中存在著未消除的之前添加的白噪聲信號，噪聲在分解結果中殘留量較大，會降低序列信號的重構效果[9]。

1.2 CEEMD分解

CEEMD的分解過程與EEMD相比，改進之處在于在序列信號進行EMD分解前，向待處理的序列信號中加入正負相對的白噪聲序列信號[10]。過程如下：

①在原始序列信號xt中添加n組幅值相同、相位角相差180°的白噪聲，得到新的2組信號，即

(4)

式中：xt為原始信號；xn為白噪聲；X1，X2為得到的新的信號。由此可得信號集合中的數量為2n組。

元代聲律學家周德清在《中原音韻》中，對宮調的聲情特色作出歸納，其中有“仙呂調清新綿邈，黃鐘宮富貴纏綿”，與姜夔“清空峭拔”詞風相偕，本詞的宮調應合仙呂調。不論是否入樂，“燕”作去聲時似乎更能將音韻延長，以達“綿邈”，且“燕子”輕靈之感亦增清新面貌。姜夔頗曉音律，對平仄音節應較為敏感。“燕”作平聲時，與“雁”字連讀，則似有不順，略為拗口。盡管“燕地”渺遠，有似開拓意境之嫌，但既思之意境，何不用“幽雁”“北雁”呢？

②對新得到的信號X1和X2分別使用EMD方法進行分解，將分解后得到的每組信號記為IMF1和IMF2。

③將CEEMD分解的最終結果記為Xi，Xi的值為每組IMF1和IMF2分量的均值，即：

(5)

CEEMD分解充分保留了EEMD的優點，解決了混疊模態的現象，相對于EEMD算法，CEEMD方法節省了信號處理時間，并且隨著白噪聲添加的數量增加，最后重構出的信號數據中噪聲的殘留量減小，最終的噪聲殘留量可以忽略不計，提高了信號重構的效果。

2 信號分析結果評價標準

IGS站高程時間序列信號分析效果可從以下兩方面進行考慮：

一是重構后信號的信噪比(signal-to-noise ratio，SNR)，這里采用的信噪比計算方法為：

(6)

二是重構信號的均方根誤差(root mean square error，RMSE)。利用降噪后的時間序列Xi和原始時間序列xt定義RMSE為：

(7)

RMSE表現出了重構信號和原始信號的區別，RMSE值越小，則表明重構后的信號質量越好，精度越高。反之，RMSE值越大，重構效果越不理想。

3 實驗分析

本文實驗數據采用BJFS站2000—2018年和SHAO站2000—2013年的高程時間序列作為原始信號，數據來源于SOPAC網站，利用3倍中誤差剔除數據中的粗差，同時剔除序列均值、同震與非地震跳變的影響。時間序列如圖1、圖2所示。

圖1 BJFS站高程時間序列

圖2 SHAO站高程時間序列

3.1 分界點n的判定

對得到的BJFS站和SHAO站的數據進行CEEMD分解，得到時間序列的IMF分量圖，如圖3、圖4所示。

圖3 BJFS站分解的IMF分量

圖4 SHAO站分解的IMF分量

IGS站高程時間序列信號經分解后為一系列從高頻到低頻排列的IMF分量[11]，為使重構后的信號達到更好的效果，如何正確識別噪聲與信號的分界點是關鍵[12]。Wu和Huang等提出了一種使用平均周期與能量密度乘積來判定噪聲與信號分界點的方法，當某一IM-F分量的平均周期與能量密度的乘積發生顯著變化時即為噪聲與信號的分界點n[13-15]，然后從該分界點的下一IMF分量進行信號重構[16]。

定義IMF分量平均周期與能量密度乘積為：

(8)

表1 BJFS站部分IMF分量平均周期與能量密度乘積

表2 SHAO站IMF分量平均周期與能量密度乘積

BJFS站：IMF1～IMF4的CONST均值為0.046 9，IMF1～IMF4與其均值偏差不超過84%，而IMF5的CONST值為0.128 5，大約是IMF1～IMF5均值的2倍，故可認為IMF1～IMF4為噪聲模態分量，分界點為5。可認為IMF1～IMF4均為噪聲分量，其余為信號分量。

SHAO站：IMF1～IMF3的CONST均值為0.070 3，IMF1～IMF3與其均值偏差不超過44%，而IMF4的CONST值為0.227 2，其與IMF1～IMF4的均值0.109 5偏差2倍以上，故可認為IMF1～IMF3為噪聲模態分量，分界點為4。舍棄噪聲分量IMF1～IMF3，對剩余的信號分量進行重構。

3.2 信號重構分析

對經過分界點n識別并去除噪聲分量的余下信號分量進行重構，重構效果如圖5和圖6所示。

圖5 BJFS站信號重構圖

圖6 SHAO站信號重構圖

根據信噪比和均方根誤差的重構評價指標如表3、表4所示。

表3 BJFS站信號重構評價指標

表4 SHAO站信號重構評價指標

分析表3、表4可見：

BJFS站：對使用CEEMD方法分解后的信號進行重構，其信噪比相比較于EEMD分解和EMD分解分別提高了17%和27%，均方根誤差指標分別降低了41%、55%。

SHAO站：對使用CEEMD方法處理后的信號進行信號重構的信噪比比EEMD、EMD分別提高了19%、32%，均方根誤差相比于EEMD、EMD分別降低了20%、37%。

4 結束語

采用CEEMD方法分解高程時間序列信號避免了EMD分解的混疊模態分量問題。相較于EEMD分解，CEEMD通過添加獨立同分布的、正負相對的互補白噪聲序列，提高了計算效率，減弱了EEMD由于添加大量白噪聲對原始序列的破壞并減少了噪聲的殘留量。

使用CEEMD方法對時間序列進行處理，并進行噪聲的識別和舍去，對剩余的信號分量進行重構，其重構效果相比較于EEMD、EMD方法在信噪比上有明顯的提高，在均方根誤差上有明顯的降低。可以看出，CEEMD分解和重構效果優于EEMD和EMD，為后續對CORS站高程時間序列的研究提供了思路。但對于序列含有色噪聲影響的問題還需作進一步的研究。