EMD-IVMD方法在驗潮站水位序列降噪中的應用

郇常敏,周世健,魯鐵定,賀小星,徐華卿

(1.東華理工大學 測繪工程學院，南昌 330013；2.南昌航空大學，南昌 330063；3.江西理工大學 土木與測繪工程學院，江西 贛州 341000)

當前海平面的急劇變化使人們對驗潮站數(shù)據(jù)的研究逐步加深，驗潮站水位序列受板塊邊界的構造運動、冰川均衡調(diào)整以及其它誤差的影響變得十分復雜[1]，且其有非平穩(wěn)、非線性等特點。因此，如何對驗潮站中的水位序列進行有效地降噪，使其發(fā)揮積極作用，具有重要的現(xiàn)實意義。驗潮站坐標時間序列在海平面變化、潮汐模型建立及海洋垂直基準建立等科學研究和工程建設等方面應用廣泛，其作用對于氣候變化、海洋自然災害、海岸帶工程的建設具有重大意義[2-4]。

在信號分析處理領域，許多專家學者進行了深入的研究，比如Huang等[5]提出了一種用于信號分析的方法—經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)，在處理非平穩(wěn)、非線性信號數(shù)據(jù)的降噪方面應用廣泛。在EMD方法的應用方面，張雙成等[6-7]針對GPS時間序列具有的一系列問題，采用EMD的方法對中國區(qū)域內(nèi)的IGS站時間序列的N,E,U方向進行處理，最終實現(xiàn)GNSS時間序列信號的有效降噪；同樣魯鐵定等[8]針對采用相關系數(shù)準則確定臨界分量K值不確定的問題，提出一種利用綜合評價指標T的新型EMD降噪方法。最終在EMD方法的基礎之上,Wu等[9]在原始信號中添加不同強度的高斯白噪聲，提出整體經(jīng)驗模態(tài)分解的降噪方法。

在Dragomiretskiy等[10]提出一種新型自適應復雜信號分解方法—變分模態(tài)分解方法之后；Xu等[11]針對傳統(tǒng)的VMD(variational mode decomposition)方法在降噪過程中存在無法確定分解模態(tài)函數(shù)(IMF)個數(shù)以及懲罰因子的取值的問題，提出了一種以能量熵互信息(eneergy entropy mutual information,EEMI)作為目標函數(shù)的改進變分模態(tài)分解(VMD)與小波包結合的降噪算法；魯鐵定等[12]針對變形監(jiān)測數(shù)據(jù)中分解含有的噪聲難以去除等問題，提出一種變分模態(tài)分解與樣本熵相結合的方法；羅亦泳等[13]提出改進變分模態(tài)分解算法(improved variational mode decomposition，IVMD)構建變形特征提取新方法，這些方法模型的結合與改進使得VMD方法的性能有了比較大的提升。為了識別軸承早期故障，Zhou等[14]提出了一種基于自適應變分模態(tài)分解(VMD)的軸承早期特征提取方法,該方法不僅改善了蝗蟲優(yōu)化算法(grasshopper optimiztion algorithm,GOA)的局部最優(yōu)問題，而且可以自適應地確定VMD的模態(tài)數(shù)和懲罰參數(shù)。徐新等[15]針對變壓器故障診斷中傳統(tǒng)的BP神經(jīng)網(wǎng)絡算法效率低、收斂速度慢等缺點，提出一種基于蝗蟲優(yōu)化算法的神經(jīng)網(wǎng)絡的故障判斷方法。

基于以上研究，綜合考慮變分模態(tài)分解的參數(shù)優(yōu)化和驗潮站水位序列的特性，以及經(jīng)驗模態(tài)分解的適用性，采用一種EMD與IVMD相結合的降噪方法，并結合實驗比較本文方法與傳統(tǒng)的EMD與EEMD(ensemble empirical mode decomposition)方法的降噪效果。

1 原理理論與方法

1.1 EMD降噪原理

經(jīng)驗模態(tài)分解(EMD)方法可將時間序列分解為若干分量，在處理非平穩(wěn)、非線性信號方面非常有效。該方法可將原始信號分解為若干個頻率由高到低的分量(本征模態(tài)函數(shù))和一個殘余項，最終將殘余項與低頻的IMF分量進行重構可實現(xiàn)信號的降噪。分量必須滿足的兩個條件是[5]：①在整個數(shù)據(jù)段內(nèi)，極值點的個數(shù)和過零點的個數(shù)必須相等或相差最多不能超過1個。②在任意時刻，由局部極大值點形成的上包絡線和由局部極小值點形成的下包絡線的平均值為0，即上、下包絡線相對于時間軸局部對稱。

原始信號被分解后的表達式為:

(1)

式中：x(t)為原始時間序列；t=1,2,…,N；m為EMD分解所得本征模態(tài)函數(shù)的個數(shù)；K為本征模態(tài)函數(shù)分量的序號；r(t)為殘余項;IMFk為分解后得到的本征模態(tài)函數(shù)分量。

原始信號進行EMD分解的過程為：①找出原始信號x(t)的極大值和極小值點，計算出上、下包絡線的均值T1,將原始信號減去該包絡均值，得到1個新的序列C1(t)；②重復進行步驟①)，直到滿足IMF閾值條件為止，則可得到IMF分量；③將原始時間序列x(t)減去第1個IMF分量形成1個新的數(shù)據(jù)序列x′(t),重復步驟①和②，最終得到m個IMF分量，殘余項滿足單調(diào)條件才可行。

原始數(shù)據(jù)序列經(jīng)EMD分解之后，需要確定臨界的IMF即低頻信號與高頻信號的分界點，可采用相關系數(shù)準則來判斷，互相關系數(shù)的首個極小值點時對應的IMF分量即為分界IMF的最佳預估項。判斷公式如下：

(2)

式中:Pk為第k個IMF分量與原始序列的相關系數(shù);IMFk(t)為分解后得到的本征模態(tài)函數(shù)分量;x(t)為原始數(shù)據(jù)序列。

針對EMD模態(tài)混疊等問題，利用高斯白噪聲具有頻率均勻分布的統(tǒng)計特性，通過每次加入同等幅值的不同白噪聲來改變信號的極值點特性，然后對多次EMD得到的相應IMF進行總體平均來抵消加入的白噪聲，從而有效抑制模態(tài)混疊的產(chǎn)生，即為集合經(jīng)驗模態(tài)分解 (EEMD)方法，其本質(zhì)則是一種疊加高斯白噪聲的多次經(jīng)驗模式分解[9]。

1.2 IVMD降噪原理

1.2.1 VMD算法

VMD(變分模態(tài)分解)分解是一種自適應分解方法，本征模態(tài)函數(shù)(Intrinsic mode function，IMF)被定義為一個有帶寬限制的調(diào)幅-調(diào)頻函數(shù),VMD算法是通過設置模態(tài)數(shù)、懲罰參數(shù)和上升步長等參數(shù)將信號分解成K個中心頻率為ωk的模態(tài)函數(shù)，其實質(zhì)為變分問題的構造和求解。VMD算法構造約束變分問題可表示為[6]：

(3)

式中:x(t)為原始數(shù)據(jù)序列；μk(t)為模態(tài)函數(shù)；ωk為各模態(tài)函數(shù)的實際中心頻率；e-jωkt為各解析信號的預估中心頻率。為求解該約束變分問題，引入拉格朗日乘子λ(t)和二次懲罰因子α，將該約束變分問題轉換為無約束變分問題，因此增廣的拉格朗日表達式為：

L(μk(t),ωk,λ(t))=

(4)

1.2.2 IVMD算法

鑒于驗潮站水位序列具有非線性、非平穩(wěn)的特點，且其包含的噪聲復雜，用單一指標來獲取信號的特征相對困難，多個單一指標的融合可提供更強的魯棒性。因此對于VMD的參數(shù)確定問題，文中采用由能量熵和互信息組成的能量熵互信息(EEMI)指標，采用兩個模態(tài)函數(shù)的能量熵互信息之和作為目標函數(shù)，利用GOA算法對VMD參數(shù)進行優(yōu)化?；认x優(yōu)化算法(GOA)是由Shahrzad Saremi等[15]提出的一種啟發(fā)式仿生優(yōu)化算法，其優(yōu)點是搜索效率高和收斂速度快，該算法模擬的是自然界中蝗蟲群體的捕食行為。該算法主要分為探索和開發(fā)兩個步驟：在探索過程中蝗蟲群體長距離快速跳躍，這一階段有利于全局搜索；在開發(fā)中蝗蟲群體只在小范圍跳躍，有利于局部搜索。

(5)

式中，fitness為目標函數(shù)；γ=(K,α)為VMD參數(shù)的取值范圍；EEMI=EEMI,EE、mi分別為能量熵和互信息，其計算式為：

(6)

(7)

式中：(X;Y)～p(x,y)，其中p(x,y)是X和Y的聯(lián)合概率分布函數(shù)，而p(x)和p(y)分別是X和Y的邊緣概率分布函數(shù)；Ei是頻域中振動信號的能量分布；IMFi(t)是不同頻率帶寬的模態(tài)分量。

利用樣本熵(SE)能夠反映信號復雜程度的特性，在該改進的VMD方法中采用樣本熵作為判斷噪聲部分與水位序列中真實信號的標準.

2 EMD-IVMD算法

EMD是一種自適應性信號處理方法，在處理非線性、非平穩(wěn)信號方面有其獨特的優(yōu)勢，且VMD同樣是一種新型信號處理方法，針對傳統(tǒng)VMD無法確定固有本征模態(tài)函數(shù)的個數(shù)和懲罰因子的數(shù)值，采用蝗蟲優(yōu)化算法(GOA)優(yōu)化目標函數(shù)的一種新型改進的VMD方法，并將EMD與改進的VMD進行結合，分別處理不同頻率階段的信號，以展示出其優(yōu)勢。

經(jīng)驗模態(tài)分解(EMD)與基于蝗蟲優(yōu)化算法優(yōu)化目標函數(shù)的變分模態(tài)分解(VMD)在驗潮站水位序列的降噪步驟如下:

1)對原始序列進行EMD分解，得到一定數(shù)量的本征模態(tài)函數(shù)(IMF)。

2)采用相關系數(shù)準則確定低頻分量與高頻分量的分界點，即相關系數(shù)首次出現(xiàn)局部極小值則認為是信噪分界點的最佳預估項。

3)分別將低頻分量與高頻分量進行重構。

4)設置VMD算法參數(shù)的范圍，初始化GOA算法參數(shù)；在文獻[16]中模態(tài)數(shù)K∈[2,8],并且懲罰因子α∈[1 000,10 000]，通過文獻[17]可得K和α的取值范圍應為[2,8]和[1 000,10 000],GOA算法的種群數(shù)N=30,最大循環(huán)數(shù)L=10。利用最優(yōu)參數(shù)的VMD將重構后的高頻IMF分量進行再次分解，計算各IMF分量的樣本熵，設定SE的閾值，將大于該閾值的視為噪聲部分，小于該閾值的分量視為信號將其進行重構。

5)最后將EMD重構的低頻信號與IVMD處理后的重構信號疊加作為最終濾波后信號。EMD-IVMD的具體實施步驟見圖1。

圖1 EMD-IVMD組合方法流程

選取信噪比(SNR)與均方根誤差(RMSE)來評價該組合方法的降噪效果，算式為：

(8)

式中：x(t)為原始時間序列；s(t)為降噪后信號；N為時序信號長度。

(9)

式中:s(t)為經(jīng)EMD-IVMD分解后重構的信號，特別說明仿真實驗部分式(15)中的x(t)為未加噪的真實信號，實測數(shù)據(jù)實驗部分x(t)則表示含噪數(shù)據(jù)序列。對于仿真實驗，信噪比越高，RMSE越小，表明濾噪效果越好，提取的信號越接近于真實信號；對于實測數(shù)據(jù)中，若信噪比越高，RMSE越小，表明濾噪效果越好且提取的有用信號更多。

3 實驗數(shù)據(jù)分析

3.1 模擬數(shù)據(jù)實驗

模擬數(shù)據(jù)主要有3個恒定振幅的周期項和高斯白噪聲來組成，其設置采樣頻率為1 HZ，采樣點數(shù)為1 024個，加入信噪比為6 dB的白噪聲(Noise)，該模擬信號構成分量波形如圖2所示。表1為模擬信號在不同方法下的精度指標對比，模擬數(shù)據(jù)表達式為：

(10)

圖2 構成模擬信號的分量信號波形

表1 模擬信號在不同方法下的精度指標對比

圖3為模擬數(shù)據(jù)的加噪信號波形，圖4為模擬信號與采用3種方法降噪后的信號波形對比圖。從圖4可以看出由EMD-IVMD方法處理后的重構信號與加噪模擬信號相比其波形非常相近且其信號曲線的光滑程度要更高，而通過傳統(tǒng)的EMD和EEMD方法處理后的信號波形與加噪模擬信號的波形存在一定的差異且其信號曲線的光滑程度過高，分析其原因可能是過分濾除真實信號，使其丟失部分真實信號。由上述分析得出，EMD-IVMD方法比傳統(tǒng)的EMD和EEMD方法降噪效果更明顯，提取有效信號的效果更優(yōu)。

模擬數(shù)據(jù)實驗中，在IVMD進行分解時，GOA算法優(yōu)化后的K=3,α=1 051。為了定量說明這3種方法的去噪效果，采用信噪比(SNR)和均方根誤差(RMSE)兩個指標，表1為模擬數(shù)據(jù)經(jīng)3種方法去噪后的兩種降噪效果評價指標值。由表1可知，在模擬數(shù)據(jù)中，EMD-IVMD方法與EMD、EEMD方法相比，EMD-IVMD方法重構后的信號的SNR相對分別提高了11.40%和10.53%，且RMSE相對分別減小了22.45%和20.83%。綜上所述，模擬實驗結果表明EMD-IVMD方法相比于EMD和EEMD方法去噪效果更優(yōu)。

圖3 模擬數(shù)據(jù)的加噪信號波形

圖4 模擬信號與各方法重構信號對比

3.2 實測站水位序列實驗分析

為驗證本方法對實際驗潮站水位變化序列的降噪效果，采用OLANDS NORRA UDDE(69站),CUXHAVEN 2(7站),AARHUS(76站),ESBJERG(80站)4個驗潮站的實測數(shù)據(jù)進行分析，數(shù)據(jù)來源于PSMSL網(wǎng)站，(http://www.psmsl.org/data/obtaining/stations/52.php)。OLANDS NORRA UDDE站所選數(shù)據(jù)為1887—2020年跨度為133年的時間序列，CUXHAVEN 兩站所選數(shù)據(jù)為1843—2018年跨度為175年的時間序列數(shù)據(jù)，AARHUS站所選數(shù)據(jù)為1888—2017年跨度為129年的時間序列數(shù)據(jù)，ESBJERG站所選數(shù)據(jù)為1889—2017年跨度為128年的時間序列數(shù)據(jù)。在進行實驗之前，去除序列中較大的粗差，將剔除粗差的時間序列作為本實驗的原始時間序列，經(jīng)觀察分析剔除粗差后的信號振幅明顯要小于未剔除前，說明達到了較好的剔除粗差的效果。限于篇幅有限，僅以76站(AARHUS)為例進行展示分析，圖5為76站原始序列信號與經(jīng)3種方法降噪后的重構信號對比圖，特別說明圖5中展示的為水位變化序列(扣除序列平均值)。通過圖5觀察分析可知，EMD，EEMD方法對水位序列處理后其形成的水位序列曲線雖平滑，但其過量濾去了其中的真實有用信號，其振幅波動范圍與原始信號有著較大的差異。而EMD-IVMD方法所形成的序列曲線，在濾去部分噪聲后，其信號振幅波動范圍與原始信號有著相近的趨勢且一致性較好。

圖5 76站重構信號與原始信號對比

圖6 76站經(jīng)EMD-IVMD分解后IMF分量圖

如圖6所示為76站經(jīng)EMD-IVMD方法處理后的IMF分量圖，且第一次經(jīng)EMD分解后產(chǎn)生了11個分量，通過相關系數(shù)法確定臨界的IMF分量為第6個，將前6個IMF分量進行重構，而且通過觀察分析圖6，前6個IMF表現(xiàn)出明顯的隨機特性，所以再結合相關系數(shù)法可以推斷出前6個IMF為噪聲部分。第二次IVMD對重構后的時間序列進行分解，且分解后產(chǎn)生3個IMF分量。圖6(a)為76站原始時間序列經(jīng)EMD分解后的IMF分量圖，圖6(b)為經(jīng)IVMD分解后的IMF分量圖。因為EMD-IVMD方法要進行兩次分解，在對第76站進行第二次分解時GOA算法優(yōu)化后的K=3,α=9 688，在第二次分解后，計算各模態(tài)分量的樣本熵，如表2所示為IVMD分解后各模態(tài)分量樣本熵，且分解后得到的3個模態(tài)分量的樣本熵分別為0.69、0.70和0.74，所以通過大量實驗分析最終確定閾值為0.71，將IMF1、IMF2重構后與經(jīng)過第一次分解后的重構信號再次重構得到最終降噪信號。

表2 IVMD分解后各模態(tài)分量樣本熵

表3給出了4個實測驗潮站經(jīng)3種方法處理后的精度評價指標，通過表3和圖5可以看出對于4個驗潮站，EMD-IVMD與EMD和EEMD兩種方法相比，均方根誤差較小且信噪比較大，表明EMD-IVMD方法能夠更好地濾除噪聲并提取更多的有用信號。此外發(fā)現(xiàn)EMD與EEMD相比，EMD評價指標略優(yōu)于EEMD，兩者較為接近，可能源于驗潮站數(shù)據(jù)序列中噪聲以有色噪聲為主，需進一步分析。76站(AARHUS)水位變化序列通過文中提出的方法處理后，SNR降噪指標相比于其余兩種方法提高了0.46%和0.62%，EMD-IVMD濾噪效果較好，RMSE相對分別減少了8.16%和14.41%。

表3 4個實測站降噪指標值對比

4 結 論

針對驗潮站水位變化序列的非線性非平穩(wěn)特性，提出一種融合EMD與蝗蟲優(yōu)化算法的VMD方法的濾波算法。該方法需要對數(shù)據(jù)序列進行兩次分解，首次利用EMD分解提取出高頻含噪部分和低頻信號部分，并重構得到以高頻噪聲為主的分量以及以信號為主的低頻分量。對重構的高頻噪聲分量進行IVMD分解，將其含有的部分有效信號提取出來。最后將EMD首次提取的低頻重構分量與IVMD提取的有效信號疊加得到濾波后最終信號。文中通過模擬實驗以及實測驗潮站水位序列數(shù)據(jù)進行分析，并評估EMD、EEMD和EMD-IVMD 3種方法的濾波效果。仿真結果顯示，文中提出的EMD-IVMD方法更能有效地提取出高頻噪聲部分，降噪效果要明顯優(yōu)于傳統(tǒng)的EMD和EEMD方法。實測數(shù)據(jù)結果表明，與傳統(tǒng)的EMD和EEMD方法相比，4個驗潮站EMD-IVMD方法的信噪比平均分別提高了1.67%和1.52%，均方根誤差平均分別減小了9.59%和13.51%。