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

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

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

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

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

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

1 原理理論與方法

1.1 EMD降噪原理

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

原始信號(hào)被分解后的表達(dá)式為:

(1)

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

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

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

(2)

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

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

1.2 IVMD降噪原理

1.2.1 VMD算法

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

(3)

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

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

(4)

1.2.2 IVMD算法

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

(5)

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

(6)

(7)

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

利用樣本熵(SE)能夠反映信號(hào)復(fù)雜程度的特性，在該改進(jìn)的VMD方法中采用樣本熵作為判斷噪聲部分與水位序列中真實(shí)信號(hào)的標(biāo)準(zhǔn).

2 EMD-IVMD算法

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

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

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

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

3)分別將低頻分量與高頻分量進(jìn)行重構(gòu)。

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

5)最后將EMD重構(gòu)的低頻信號(hào)與IVMD處理后的重構(gòu)信號(hào)疊加作為最終濾波后信號(hào)。EMD-IVMD的具體實(shí)施步驟見圖1。

圖1 EMD-IVMD組合方法流程

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

(8)

式中：x(t)為原始時(shí)間序列；s(t)為降噪后信號(hào)；N為時(shí)序信號(hào)長(zhǎng)度。

(9)

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

3 實(shí)驗(yàn)數(shù)據(jù)分析

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)

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

(10)

圖2 構(gòu)成模擬信號(hào)的分量信號(hào)波形

表1 模擬信號(hào)在不同方法下的精度指標(biāo)對(duì)比

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

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

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

圖4 模擬信號(hào)與各方法重構(gòu)信號(hào)對(duì)比

3.2 實(shí)測(cè)站水位序列實(shí)驗(yàn)分析

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

圖5 76站重構(gòu)信號(hào)與原始信號(hào)對(duì)比

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

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

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

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

表3 4個(gè)實(shí)測(cè)站降噪指標(biāo)值對(duì)比

4 結(jié) 論

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