范小猛，胡川，張重陽，李成洪

(重慶交通大學 土木工程學院,重慶 400074)

0 引言

在全球衛(wèi)星導航系統(tǒng)(GNSS)坐標時間序列中，不僅包含真實信號，還含有白噪聲(WN)、閃爍噪聲(FN)、隨機漫步噪聲(RWN)等多種有色噪聲.合理的分離出信號和噪聲，有助于獲取更準確的測站運動參數(shù)估值，對于精化速度場、建立地球框架以及揭示構造變形運動等有著重要的意義[1-2].

經(jīng)驗模態(tài)分解(EMD)、整體經(jīng)驗模態(tài)分解(EEMD)和小波降噪是三種較常用的噪聲分離方法.HUANG等[3]在1998 年討論了EMD 在時間序列分析中的應用問題.張雙成等[4]將基于EMD 的時間序列分析法引入GNSS 坐標序列分析，對中國境內(nèi)及周邊9 個國際GNSS 服務(IGS)站三個坐標分量上的序列數(shù)據(jù)進行處理，有效地分離出線性趨勢項、周期項和噪聲項，證明了EMD 法降噪的有效性.隨著研究的不斷深入，研究人員發(fā)現(xiàn)EMD 存在模式混疊的現(xiàn)象[5-6]，為了改善這個問題，WU 等[7]在EMD 的基礎上提出了EEMD.張恒璟等[8]采用EMD 和EEMD分別對BJFS 站和URUM 站10 a 左右的高程數(shù)據(jù)進行處理，分離序列趨勢項和周期項運動，發(fā)現(xiàn)EEMD可以明顯減弱模式混疊現(xiàn)象，但不能完全消除模式混疊，同時EEMD分解時所加入的WN 會導致其分解結果不唯一.小波分析具有多分辨率、低熵性、去相關性等特點，得到了GNSS 時間序列分析研究人員的重視[9-11].魏冠軍等[12]采用小波變換閾值降噪對高程分量的坐標時序進行處理，有效減弱了噪聲信息，證明了該方法降噪的可行性，但樣本序列數(shù)據(jù)時間長度在3 a 左右，沒有涉及到更長時間跨度時序的降噪分析.馬俊等[13]提出了小波包系數(shù)信息熵法，對KMIN 站近5 a 的三個坐標分量上的時序進行了降噪處理，得到了僅受WN 影響的時間序列，消除了有色噪聲的影響，但同樣沒有對更長時間跨度的坐標時間序列進行分析.

目前，已經(jīng)有眾多學者對EMD、EEMD 和小波降噪有了深入研究，但少有文獻考慮到數(shù)據(jù)長度可能會對降噪效果產(chǎn)生的影響.隨著GNSS基準站的不間斷觀測，相關機構已積累了20 多年的坐標時間序列數(shù)據(jù)，為研究上述方法在不同時序長度下的降噪效果提供了可能.針對上述問題，本文通過對比分析EMD、EEMD 和小波降噪在不同GNSS 高程坐標時序長度下的降噪效果，采用信噪比(SNR)、相關系數(shù)和均方根誤差(RMSE)定量評價它們的降噪性能，為相關研究提供一定參考.

1 降噪方法

1.1 經(jīng)驗模態(tài)分解

EMD 是一種基于非線性、非平穩(wěn)數(shù)據(jù)的自適應時頻信號分解方法，在不同尺度上分析數(shù)據(jù)，提取有用信息[3].其基本思想是通過篩選，把原始信號X(t)分解成m個不同時間尺度的本征模態(tài)分量(IMF)和殘差項rn，分解流程如下：

1)尋找原始信號中所有的局部極大、極小值點，利用三次樣條函數(shù)分別擬合上述點，得到關于時間的上、下包絡線，計算兩條包絡線的平均包絡線m1(t)，用原始信號X(t)減去平均包絡得到新的待分解信號h1(t)

2)判斷h1(t)是否滿足以下條件：一是待分解的信號中的極值點的個數(shù)與過零點的個數(shù)應相等或最多相差1；二是在任意點上，上下包絡線的均值為0，即上下包絡線應對稱.如果滿足上述條件，將h1(t)作為分解得到的第一個分量 I MF1(t)，然后進行第3)步操作，否則將h1(t)作為新的原始信號重復第1)步操作.

3)計算剩余信號r1=X(t)-IMF1(t)，然后將r1作為新的待分解信號.

4)判斷r1是否滿足以下條件：r1為單調(diào)變化的殘差序列，且不能繼續(xù)分解得到IMF 分量.如果滿足，EMD 分解過程結束，殘差項rm=r1，否則重復流程1)～3)操作.原始信號經(jīng)過m次分解可以表示為

為了獲取不含噪聲的GNSS 坐標時間序列，需要對上式中分解出的IMF 分量進行重構.首先，計算各分量與原始序列間的相關系數(shù)；然后識別第一個取局部極小值的相關系數(shù)對應的分量IMFk，將其作為信號與噪聲之間的分界層；最后，將分界層與之前的IMF 分量層相加作為高頻噪聲，然后用原始信號減去重構高頻噪聲得到降噪后信號[3].

1.2 整體經(jīng)驗模態(tài)分解

EMD 的缺點之一是模式混疊：信號被EMD 分解為m個不同的IMF 分量后，不同頻率的信號會存在于同一個IMF 分量中，即一個IMF 分量包含著不同尺度的信號.EEMD 通過在分解前向原始信號中添加白噪聲w(t)，生成待分解信號，改善模式混疊問題.

式中，下標i表示第i次分解.其余分解步驟與EMD 相同，此處不再贅述.將多次分解得到的尺度相同的IMF 分量整體取平均，作為當前尺度的IMF 分量，從而得到最終的一系列不同尺度的IMF 分量數(shù)據(jù)，按照相關系數(shù)準則重構信號即可得到降噪之后的GNSS坐標時間序列.本文根據(jù)文獻[14]實驗結果，將加入的WN 與原始信號的標準差之比設為0.1，整體平均次數(shù)設置為100 次.

1.3 小波降噪

多分辨率分析是小波分析的重要概念，它是利用特定的小波基函數(shù)的伸縮變換，將GNSS 坐標時間序列在不同分辨率上進行分解，得到能夠反映時間序列概貌的低頻部分以及能夠反映時間序列細節(jié)的高頻部分，進而對不同分辨率上的信號進行分析[12].小波閾值降噪是在多分辨率分析的基礎上，通過設置一定的閾值對分解之后的各個尺度上的信號進行處理，達到分離噪聲的目的.小波閾值降噪的具體步驟如下：

1)小波分解.選取合適的小波基函數(shù)，對GNSS坐標時間序列進行小波分解，得到不同尺度上的小波系數(shù).

2)閾值降噪.對分解得到的各個尺度上的信號進行非線性閾值降噪處理，剔除小波系數(shù)小于閾值的小波項.

3)重構信號.將剔除小波項之后的各尺度上的信號重構得到降噪后的GNSS 坐標時間序列.

經(jīng)過實驗對比，本文采用Coif3 小波對原始數(shù)據(jù)進行分解，根據(jù)軟閾值規(guī)則對數(shù)據(jù)進行降噪處理.

2 實驗數(shù)據(jù)和方案

2.1 實驗數(shù)據(jù)

由于我國境內(nèi)的IGS 站時間序列在水平方向上主要以線性為主，而在高程方向上表現(xiàn)出復雜的非線性運動[15-16]，因此本文以高程方向上的坐標時序為例進行分析.為保證實驗所需時序長度，選取我國境內(nèi)BJFS、KMIN、LHAS、SHAO、URUM、WUHN 共6 個IGS 站點所對應的高程方向單日解坐標時序，各站點完整數(shù)據(jù)的起止歷元均為1999.17—2020.39，截取其中2000—2019 年的觀測數(shù)據(jù)作為本次實驗的原始數(shù)據(jù).

采用線性擬合去除原始數(shù)據(jù)中的趨勢項，并且根據(jù)3σ 準則剔除孤立值.另外，為了避免由于數(shù)據(jù)缺失造成模式混疊現(xiàn)象，以及滿足小波分析要求數(shù)據(jù)具有均勻取樣和零均值的特性，利用線性插值對數(shù)據(jù)缺失的部分進行插補，得到適用于降噪分析的連續(xù)數(shù)據(jù).

2.2 實驗方案

實驗方案共分為三種：

方案一：將每個站點共20 a 的坐標時間序列數(shù)據(jù)依次分為4 段，分別為2000—2004、2005—2009、2010—2014、2015—2019 時間段，得到24 個時序長度為5 a 的樣本序列；

方案二：將每個站點共20 a 的坐標時間序列數(shù)據(jù)平均分為2 段，分別為2000—2009、2010—2019，得到12 個時序長度為10 a 的樣本序列；

方案三：將每個站點2000—2019 共20 a 的坐標時間序列數(shù)據(jù)分別整體作為一段，得到6 個時序長度為20 a 的樣本序列.

降噪效果通過降噪后信號與原始信號之間的SNR、相關系數(shù)以及RMSE 來判斷.其計算公式分別為：

式中：X(t)為原始序列；為降噪后的真實信號；N為采樣點個數(shù).SNR 體現(xiàn)的是噪聲信號在信號中所占的比重，其值越大，表示噪聲含量越低；相關系數(shù)越接近1，表示兩者的相似度越高，降噪效果越好；RMSE 表示的是兩者之間的偏差程度，其值越小代表降噪效果越好.

3 實驗結果與分析

3.1 降噪結果

經(jīng)過計算，三種方法在不同樣本序列長度下所分解出的分量個數(shù)相同，從5 a 到20 a 分別為10、11、12 個.對于EMD 和EEMD 而言，根據(jù)各IMF 分量與原樣本序列的相關系數(shù)得到的信號與噪聲間分界層不同：EMD 法得到的不同樣本序列分界層數(shù)分別為第5、4、5 層，而EEMD 法得到的不同樣本序列分界層數(shù)分別為第2、5、5 層，根據(jù)上述分界層對樣本序列進行降噪處理.

圖1 和圖2 描述了不同IGS 站點高程方向上經(jīng)過降噪的坐標時間序列.從圖中可以發(fā)現(xiàn)三種方法都可以有效的對坐標時間序列數(shù)據(jù)進行降噪，降噪后信號表現(xiàn)出明顯的周期性變化，曲線變化也較為平滑.其中，圖1 中BJFS 站的5 a-2 樣本序列以及KMIN 站的5 a-4 樣本序列的光滑度明顯優(yōu)于其余相同時間跨度的樣本序列，其原因可能是在這兩個樣本數(shù)據(jù)中缺失歷元數(shù)較少，數(shù)據(jù)較為完整.例如BJFS 站的5 a-2樣本序列的缺失歷元數(shù)為54，僅占樣本總數(shù)的3%，遠小于其它曲線光滑度較低的樣本序列缺失歷元數(shù)，說明時間序列中缺失數(shù)據(jù)的存在會對降噪效果產(chǎn)生一定影響.

圖1 BJFS 站和KMIN 站降噪結果對比

圖2 SHAO 站和URUM 站降噪結果對比

圖2 中SHAO 和URUM 兩個站點，EMD 降噪后的曲線在個別時間段呈現(xiàn)出明顯的異常波動.其中SHAO 站的5 a-2 樣本序列數(shù)據(jù)在2007 年前后波動最大，最大達到了80 mm，5 a-4 樣本序列數(shù)據(jù)中異常波動較多，但波動幅度較小，其余站點也存在相似的異常波動情況.這種異常波動現(xiàn)象可能是原始時間序列數(shù)據(jù)中存在長時間的連續(xù)中斷，如SHAO 站連續(xù)中斷的歷元數(shù)為411，缺值率達到了22%，即便是后續(xù)對缺失處進行了插補，依然會對EMD 降噪結果造成偏差.

此外，在上述出現(xiàn)異常波動的時間點，EEMD 和小波降噪得到的序列趨勢相似，振幅為3～5 mm，趨勢較平緩，說明即使數(shù)據(jù)存在長時間的連續(xù)缺失，這兩種方法也可以有效減弱插值的影響，削弱信號中存在的噪聲.

3.2 降噪結果對比分析

為了進一步說明EMD、EEMD 和小波降噪的效果，分別計算不同樣本長度下各方法降噪結果的SNR、相關系數(shù)和RMSE，其中樣本長度為5 a 和10 a 時，將各指標的平均值作為最終結果，如表1 所示.

由表1 可知：

表1 各站點評價參數(shù)統(tǒng)計 m

1) EMD 受初始數(shù)據(jù)質(zhì)量的影響較大.在數(shù)據(jù)較完整的站點，EMD 降噪序列各評價指標略差或持平于EEMD 和小波降噪結果；而在缺失比例大的站點，各評價指標值明顯較差.

2)相比于EMD，EEMD 不僅改善了模式混疊問題，還在一定程度上提高了EMD 分解的精度，即便是在缺值數(shù)較多的SHAO 站，EEMD 降噪序列的相關系數(shù)值大部分在0.7 以上，降噪效果穩(wěn)定.

3)在5 a 和10 a 的樣本時間序列數(shù)據(jù)中，小波降噪序列的評價參數(shù)值明顯優(yōu)于EMD 和EEMD；當時間跨度為20 a 時，EEMD 降噪序列的評價參數(shù)與小波相近.這說明在這三種方法中，小波降噪的適用性更好，而EEMD 更適用于對時間跨度大的坐標時間序列進行降噪分析.

4)大部分站點的EEMD 降噪序列和小波降噪序列隨著時間跨度的增加，三個評價參數(shù)值均有所增加，難以根據(jù)某一個參數(shù)準確說明其降噪性能的強弱.

為了進一步探究EMD、EEMD 和小波降噪對有色噪聲的剔除能力，分別計算三種方法降噪序列的譜指數(shù)，譜指數(shù)越接近零，在一定程度上說明有色噪聲的振幅越小，對真實信號的影響也越小[13]，結果如圖3 所示.

圖3 三種方法得到的降噪序列的譜指數(shù)

由圖3 可以得出以下結論：1)EEMD 降噪后時序的譜指數(shù)大于EMD 降噪序列的譜指數(shù)，剔除有色噪聲的能力更好，但是否因EEMD 方法所加入WN 產(chǎn)生的影響還有待進一步研究；2)小波分析降噪在不同時間跨度的時序數(shù)據(jù)下都可以更好的抑制有色噪聲，使得降噪后的信號更加可靠.

4 結論

本文通過設置不同樣本長度的降噪實驗，對EMD、EEMD 以及小波分析三種方法的降噪效果進行對比，并利用SNR、RMSE 和相關系數(shù)定量分析了其降噪性能.實驗表明：

1)EMD 易受坐標時間序列自身質(zhì)量的影響，降噪效果相對較差；EEMD 和小波降噪在時間序列存在長時間中斷，質(zhì)量較差時，依然可以得到較好的降噪結果.

2)隨著時序樣本長度的增加，三種方法的降噪性能均有提升.其中，小波降噪的適用性最廣泛，在不同的樣本數(shù)據(jù)中，降噪性能均優(yōu)于EMD 和EEMD；EEMD 在時間跨度較長的情況下降噪性能更好，與小波降噪性能接近；EMD 降噪效果相對最差.

3)小波降噪相比于EMD 和EEMD，所得到的降噪序列中含有的有色噪聲更少，更有利于對數(shù)據(jù)進行后續(xù)的分析和使用.