基于CEEMDAN和MPE的導(dǎo)墻振動(dòng)信號(hào)降噪方法及應(yīng)用

劉 昉， 張魯豐， 龐博慧， 梁 超， 姚 燁

(1.天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300354；2.華能瀾滄江水電股份有限公司，云南 昆明 650214)

0 引言

導(dǎo)墻作為輕型泄流結(jié)構(gòu)，長(zhǎng)期受水流、風(fēng)等環(huán)境因素的干擾，會(huì)由于疲勞和腐蝕發(fā)生開(kāi)裂損傷，而且導(dǎo)墻水下部分結(jié)構(gòu)損傷不易被直接發(fā)現(xiàn)[1]。Taxakana、Navajo等水利樞紐的導(dǎo)墻都因水流誘發(fā)振動(dòng)而破壞[2]。對(duì)導(dǎo)墻實(shí)施振動(dòng)監(jiān)測(cè)是保證水電站安全運(yùn)行的必要手段，而監(jiān)測(cè)數(shù)據(jù)受現(xiàn)場(chǎng)工作環(huán)境影響，會(huì)存在噪聲干擾，使得結(jié)構(gòu)振動(dòng)特征信息被淹沒(méi)，因此需要對(duì)信號(hào)進(jìn)行降噪以獲得真實(shí)信號(hào)。

國(guó)內(nèi)外學(xué)者對(duì)信號(hào)降噪進(jìn)行了很多研究，徐國(guó)賓等[3]利用小波變換有效地去除了水墊塘底板振動(dòng)信號(hào)的噪聲，同時(shí)保留原始信號(hào)的細(xì)節(jié)部分。李成業(yè)等[4]利用經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)[5]與小波變換相結(jié)合的濾波方法，對(duì)拱壩振動(dòng)信號(hào)進(jìn)行降噪，提取振動(dòng)特征。EMD無(wú)須像小波變換一樣設(shè)置基函數(shù)，克服了依賴主觀經(jīng)驗(yàn)的缺點(diǎn)，但其存在模態(tài)混疊和端點(diǎn)效應(yīng)問(wèn)題[6]。胡劍超等[7]對(duì)實(shí)測(cè)振動(dòng)信號(hào)進(jìn)行處理，采用完備總體經(jīng)驗(yàn)?zāi)B(tài)分解(complete ensemble empirical mode decomposition，CEEMD)與小波包閾值降噪結(jié)合的方法對(duì)信號(hào)進(jìn)行濾波降噪，CEEMD[8]通過(guò)向原信號(hào)中添加正負(fù)成對(duì)的高斯白噪聲以解決EMD存在的模態(tài)混疊現(xiàn)象，但加入的白噪聲難以徹底消除，重構(gòu)信號(hào)中含有殘余噪聲。自適應(yīng)噪聲的完備總體經(jīng)驗(yàn)?zāi)B(tài)分解(complete ensemble empirical mode decomposition with adaptive noise，CEEMDAN)[9]做出進(jìn)一步改進(jìn)，在EMD的每個(gè)階段自適應(yīng)地添加白噪聲，計(jì)算唯一余量信號(hào)獲得各IMF分量，分解過(guò)程完備。重構(gòu)誤差很小，消除了添加噪聲對(duì)信號(hào)的影響，提高了計(jì)算效率及準(zhǔn)確性。Bai等[10]使用CEEMDAN方法、排列熵(permutation entropy，PE)和峰值濾波聯(lián)合去噪算法，能有效去除齒輪箱振動(dòng)信號(hào)噪聲。

基于上述情況，本文使用CEEMDAN和多尺度排列熵(multi-scale permutation entropy，MPE)即CEEMDAN-MPE聯(lián)合降噪方法，對(duì)振動(dòng)數(shù)據(jù)進(jìn)行CEEMDAN處理；利用MPE檢驗(yàn)各IMF分量的隨機(jī)性，篩分為含噪聲的IMF和純凈IMF；然后采用小波閾值降噪處理含噪聲的IMF分量；最后將處理后的含噪聲的IMF分量與純凈IMF分量進(jìn)行重構(gòu)，達(dá)到降噪目的，并將其應(yīng)用于導(dǎo)墻的實(shí)測(cè)振動(dòng)位移信號(hào)，驗(yàn)證了該方法的有效性。

1 基本原理

1.1 CEEMDAN算法

CEEMDAN通過(guò)在EMD的每個(gè)階段添加自適應(yīng)白噪聲，在較少平均次數(shù)下其重構(gòu)誤差仍然很小，可以解決EEMD算法中的噪聲殘留問(wèn)題，運(yùn)算效率也有相應(yīng)提高，具體步驟總結(jié)如下。

(1)在原始信號(hào)x(t)中添加白噪聲ε0ni(t)，其中ε0為第一次添加的白噪聲幅值系數(shù)，則第i次信號(hào)可以表示為xi(t)=x(t)+ε0ni(t)，對(duì)各xi(t)通過(guò)EMD，得到各個(gè)信號(hào)第一階分量imfi1和殘差ri1(t)。

(1)

(2)定義算子Ei(·)為使用EMD后的第j個(gè)分量，對(duì)r1(t)中添加白噪聲ε1E1[ni(t)]，然后進(jìn)行分解，得到第二階分量：

(2)

以此類推，第k個(gè)殘余信號(hào)為

rk(t)=rk-1(t)-imfk。

(3)

(3)第(k+1)階分量為

(4)

(4)重復(fù)步驟(2)和(3)，直到信號(hào)無(wú)法繼續(xù)分解，運(yùn)算終止，共得到K個(gè)分量，最終殘余量為

(5)

此時(shí)，原始信號(hào)可以表示為

(6)

1.2 MPE算法

MPE能夠檢測(cè)信號(hào)的突變性和隨機(jī)性，是一種判斷信號(hào)復(fù)雜度的指標(biāo)[11]。與排列熵PE只能檢測(cè)單一尺度的隨機(jī)性相比，MPE能夠衡量信號(hào)不同尺度下的隨機(jī)性[12]，而泄流振動(dòng)信號(hào)在不同尺度均包含重要信息，因此MPE更適用于該信號(hào)分析。計(jì)算步驟如下。

對(duì)時(shí)間序列Z={z1,z2,…,zL}進(jìn)行多尺度粗?；幚恚?/p>

(7)

N=(L/s-(m-1)t)；

(8)

(9)

式中：m為嵌入維數(shù)；t為時(shí)間延遲。

(10)

式中：Pj為第j次符號(hào)出現(xiàn)的概率。

對(duì)式(10)計(jì)算排列熵進(jìn)行歸一化處理：

(11)

熵值大小代表時(shí)間序列的隨機(jī)和復(fù)雜程度，值越接近于1，隨機(jī)性越強(qiáng)，非平穩(wěn)程度越高。參考文獻(xiàn)[13]，將排列熵閾值設(shè)置為0.6，大于或等于閾值的IMF分量視為含噪聲的分量，需要進(jìn)一步處理。

1.3 小波閾值去噪

小波閾值去噪的原理為將信號(hào)經(jīng)過(guò)小波變換產(chǎn)生的小波系數(shù)中會(huì)含有信號(hào)的重要信息，真實(shí)信號(hào)的系數(shù)較大，噪聲信號(hào)的系數(shù)較小，通過(guò)對(duì)小波系數(shù)設(shè)置合適的閾值，小于閾值的小波系數(shù)認(rèn)為是由噪聲引起，將其置0從而達(dá)到去噪目的[14]。具體步驟主要包括：

(1)選擇合適的小波基及分解層數(shù)，對(duì)含噪信號(hào)進(jìn)行小波分解；

(2)選擇合適的閾值方法對(duì)各分解尺度下的小波系數(shù)進(jìn)行閾值處理；

(3)對(duì)小波系數(shù)進(jìn)行重構(gòu)，完成降噪。

小波閾值處理主要有軟閾值去噪和硬閾值去噪。軟閾值去噪與硬閾值去噪相比，克服了不連續(xù)的缺陷，去噪后結(jié)果更平滑，因此本文采用軟閾值去噪方法，其數(shù)學(xué)表達(dá)式為

(12)

2 CEEMDAN-MPE聯(lián)合降噪

CEEMDAN-MPE聯(lián)合降噪方法的步驟如下。

首先，利用CEEMDAN將導(dǎo)墻振動(dòng)信號(hào)進(jìn)行分解，獲得若干IMF分量，計(jì)算各IMF的MPE值ZMPE，檢驗(yàn)各IMF分量的隨機(jī)性，若MPE值大于等于設(shè)定的閾值0.6，為含噪聲的IMF分量；MPE值小于0.6則為純凈IMF。其次，對(duì)于含噪聲的IMF分量使用小波軟閾值方法進(jìn)行去噪。最后，將小波軟閾值處理后的結(jié)果與純凈IMF分量進(jìn)行重構(gòu)，得到去噪后的信號(hào)，流程見(jiàn)圖1。

圖1 CEEMDAN-MPE降噪方法流程Figure 1 Flow chart of CEEMDAN-MPE denoising method

利用仿真信號(hào)驗(yàn)證聯(lián)合降噪方法的效果。設(shè)置采樣時(shí)間為10 s，采樣頻率為100 Hz的仿真信號(hào)，其表達(dá)式為

y1=0.3sin πt+0.5sin 2πt+0.4sin 8πt。

(13)

加噪后信號(hào)表達(dá)式為

y2=y1+w。

(14)

式中:w=0.5randn(1,1 000)為白噪聲分量,由MATLAB軟件生成，白噪聲隨時(shí)間起伏變化快，其功率譜密度在頻帶內(nèi)是均勻的或者變換很小。

采用CEEMDAN-MPE方法對(duì)信號(hào)進(jìn)行處理，得到降噪后結(jié)果見(jiàn)圖2。由圖2可知，添加的白噪聲已被去除，信號(hào)波形與原始信號(hào)的相似程度較高。

圖2 降噪后信號(hào)波形對(duì)比Figure 2 Comparison of signal waveform after denoising

為驗(yàn)證降噪效果，分別使用CEEMDAN-相關(guān)系數(shù)方法及EEMD-MPE方法對(duì)信號(hào)進(jìn)行處理，并與本文方法進(jìn)行對(duì)比。含噪聲IMF分量與原信號(hào)相關(guān)性很差，因此可根據(jù)各IMF分量與原信號(hào)的相關(guān)系數(shù)來(lái)判定虛假分量[16]，相關(guān)系數(shù)計(jì)算公式為

(15)

式中：X為IMF分量；Y為原始信號(hào)；N為信號(hào)長(zhǎng)度。

選擇信噪比(SNR)、均方根誤差(RMSE)及降噪信號(hào)與原始信號(hào)的相關(guān)系數(shù)作為評(píng)價(jià)降噪效果的指標(biāo)。信噪比及相關(guān)系數(shù)越大，均方根誤差越小，則有較好的降噪效果。由表1可知本文采用方法降噪效果最優(yōu)。

表1 不同降噪方法的效果評(píng)價(jià)Table 1 Effect evaluation of different denoising methods

3 工程應(yīng)用

3.1 工程概況

某水電站泄水建筑物由主河床的4孔泄洪閘和右岸導(dǎo)流明渠內(nèi)3孔泄洪閘組成。以泄流導(dǎo)墻為研究對(duì)象，在其頂面布置10個(gè)測(cè)點(diǎn)，具體布置見(jiàn)圖3，各測(cè)點(diǎn)布置三向低頻振動(dòng)位移傳感器，導(dǎo)墻橫河向振動(dòng)量大，順河向及垂向振動(dòng)相對(duì)較小，所以主要對(duì)導(dǎo)墻橫河向振動(dòng)進(jìn)行研究。通過(guò)分析發(fā)現(xiàn),距離閘門(mén)較近的測(cè)點(diǎn)振動(dòng)量相對(duì)較大，且易受到噪聲干擾，淹沒(méi)有用信號(hào)，影響監(jiān)測(cè)結(jié)果。因此，選取位于導(dǎo)墻斜坡段的測(cè)點(diǎn)2進(jìn)行降噪處理。

圖3 導(dǎo)墻及測(cè)點(diǎn)布置Figure 3 Guide wall and measuring point layout

3.2 降噪處理與分析

測(cè)點(diǎn)2的振動(dòng)位移時(shí)程線及功率譜密度如圖4所示。由圖4可知，振動(dòng)信號(hào)有較多“毛刺”，存在噪聲，采用CEEMDAN-MPE聯(lián)合降噪方法對(duì)其進(jìn)行處理。首先進(jìn)行CEEMDAN分解，共得到13個(gè)IMF分量，IMF1～I(xiàn)MF13頻率依次降低。

圖4 測(cè)點(diǎn)2振動(dòng)信號(hào)時(shí)程圖與功率譜密度圖Figure 4 Time history and power spectral density of vibration signal at measuring point 2

通過(guò)MPE方法計(jì)算各IMF分量的MPE值(ZMPE)，需要設(shè)置合適的參數(shù)，本文使用互信息法和偽近鄰法[17]計(jì)算時(shí)間延遲τ，選擇嵌入維數(shù)m。信號(hào)的互信息與時(shí)間延遲τ的關(guān)系見(jiàn)圖5。τ=1時(shí)，互信息值取得的第一個(gè)極小值，為最優(yōu)延遲時(shí)間。

圖5 互信息隨時(shí)間延遲變化曲線Figure 5 Curve of mutual information with time delay

在時(shí)間延遲τ基礎(chǔ)上，嵌入維數(shù)通過(guò)分析偽近鄰率隨嵌入維數(shù)的變化規(guī)律來(lái)確定。判據(jù)1閾值設(shè)置為15，判據(jù)2閾值設(shè)置為2，最大嵌入維數(shù)設(shè)置為10。偽近鄰率計(jì)算結(jié)果如圖6所示。

圖6 偽近鄰率隨嵌入維數(shù)變化曲線Figure 6 Curve of FNNP with embedding dimension

圖6中，判據(jù)1的偽鄰近率逐漸降低，判據(jù)2的偽鄰近率在嵌入維數(shù)為6之后不再隨嵌入維數(shù)的增加而減少，因此，最佳嵌入維數(shù)m=6。選定嵌入維數(shù)及時(shí)間延遲后，經(jīng)過(guò)試算，選擇尺度因子s=12。各IMF分量的ZMPE如圖7所示。IMF1～I(xiàn)MF13的ZMPE逐漸降低，其中IMF1～I(xiàn)MF6的ZMPE大于0.6，視為含噪聲的IMF分量；IMF7～I(xiàn)MF13的ZMPE值小于0.6，為純凈IMF分量。

圖7 IMF分量的ZMPEFigure 7 ZMPE of the IMF component

采用小波軟閾值去噪方法對(duì)IMF1～I(xiàn)MF6分量進(jìn)行處理，選擇去噪效果較好的db4小波基，進(jìn)行4層分解，固定閾值為9.673 0。將降噪后的IMF1～I(xiàn)MF6與IMF7～I(xiàn)MF13進(jìn)行重構(gòu)，得到降噪后的導(dǎo)墻振動(dòng)信號(hào)如圖8所示。

圖8 降噪后信號(hào)時(shí)程圖與功率譜密度圖Figure 8 History diagram and power spectral density diagram of signal after denoising

對(duì)比圖4和圖8可知，降噪后信號(hào)噪聲減少，更好地反映了導(dǎo)墻振動(dòng)的波形特征。從功率譜密度圖中可知，大于6 Hz的高頻白噪聲成分已被去除，低頻成分信息保存完好，保留了信號(hào)的關(guān)鍵信息。

3.3 降噪效果評(píng)價(jià)

為進(jìn)一步驗(yàn)證本文算法的有效性，分別用EEMD-MPE算法及CEEMDAN-相關(guān)系數(shù)算法對(duì)信號(hào)進(jìn)行處理，降噪結(jié)果見(jiàn)圖9。

由圖9可知，CEEMDAN-MPE聯(lián)合降噪的效果優(yōu)于另外2種方法，EEMD-MPE降噪效果一般，仍存在有部分高頻噪聲，是因?yàn)镋EMD分解時(shí)加入的白噪聲無(wú)法徹底消除；CEEMDAN-相關(guān)系數(shù)降噪后信號(hào)峰值點(diǎn)處仍存在“毛刺”，而且部分低頻成分被濾除，是由于CEEMDAN分解能夠消除加入白噪聲的影響，但相關(guān)系數(shù)判斷各IMF分量屬性不準(zhǔn)確，使得效果不佳。在本文采用的方法中，通過(guò)計(jì)算各IMF分量的MPE值，能更好地判斷出含噪聲的IMF分量及純凈IMF分量。

圖9 降噪效果對(duì)比圖Figure 9 Comparison diagram of denoising effect

由于在實(shí)際工程中，實(shí)測(cè)信號(hào)中有效信號(hào)和噪聲信號(hào)的功率均是未知的，無(wú)法采用信噪比判斷降噪效果。因此本文引入降噪誤差比(dnSNR)判斷降噪質(zhì)量[18]，其計(jì)算公式如下：

dnSNR=10lg(Ps/Pg)。

(16)

式中：Ps為含噪信號(hào)的功率；Pg為去除噪聲信號(hào)的功率，dnSNR值越小，表明降噪效果越顯著。各方法的降噪誤差比計(jì)算結(jié)果見(jiàn)表2。

由表2可知，CEEMDAN-MPE算法的降噪誤差比最小，表明該算法降噪效果最好。

表2 降噪誤差比計(jì)算結(jié)果Table 2 Calculation results of denoising error ratio

4 結(jié)論

使用CEEMDAN和MPE聯(lián)合降噪的方法，對(duì)導(dǎo)墻振動(dòng)信號(hào)進(jìn)行CEEMDAN分解，利用MPE方法計(jì)算各IMF分量隨機(jī)程度，結(jié)合小波變換對(duì)含噪聲的IMF進(jìn)行降噪，并與余下純凈IMF分量進(jìn)行重構(gòu)，完成降噪。通過(guò)該方法對(duì)實(shí)測(cè)信號(hào)進(jìn)行分析，結(jié)果表明該方法在有效降噪的同時(shí)，更好地保留了結(jié)構(gòu)特征信息。比較于EEMD-MPE方法及CEEMDAN-相關(guān)系數(shù)方法，該方法更適用于導(dǎo)墻振動(dòng)信號(hào)降噪，提高了導(dǎo)墻振動(dòng)的監(jiān)測(cè)精度，為實(shí)現(xiàn)泄流導(dǎo)墻健康運(yùn)行及故障診斷提供依據(jù)。