中值互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解

劉淞華 何冰冰 郎 恂 陳啟明 張榆鋒 蘇宏業(yè)

經(jīng)驗(yàn)?zāi)B(tài)分解 (Empirical mode decomposition,EMD)算法是由Huang 等[1]提出的一種數(shù)據(jù)驅(qū)動(dòng)分解的時(shí)-頻分析方法[2],適用于分析處理非平穩(wěn)非線(xiàn)性信號(hào).該方法能自適應(yīng)地將原信號(hào)分解為一組固有模態(tài)函數(shù)(Intrinsic mode functions,IMFs)和一個(gè)殘余量,再對(duì)各個(gè)IMF 進(jìn)行希爾伯特變換(Hilbert transform,HT),得到具有實(shí)際意義的瞬時(shí)頻率.EMD 方法自提出后在生物醫(yī)學(xué)工程,工業(yè)控制過(guò)程振蕩監(jiān)測(cè),機(jī)械設(shè)備故障診斷,圖像處理等眾多工程領(lǐng)域得到了廣泛的應(yīng)用[3-5].目前,EMD系列方法已與機(jī)器學(xué)習(xí)深度融合,在預(yù)測(cè)、分類(lèi)等問(wèn)題上實(shí)現(xiàn)了不小的突破[6-7].

文獻(xiàn)[8]指出EMD 存在模態(tài)混疊和分解局限性等問(wèn)題,其中模態(tài)混疊將導(dǎo)致時(shí)頻分布出現(xiàn)錯(cuò)位,使IMF 失去物理意義,并對(duì)后續(xù)的信號(hào)分解產(chǎn)生影響.Wu 等[8]基于白噪聲的統(tǒng)計(jì)特性研究結(jié)果[9],提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble EMD,EEMD).通過(guò)將獨(dú)立分布的高斯白噪聲加入至原信號(hào)中進(jìn)行EMD 分解,并對(duì)多次分解結(jié)果中相對(duì)應(yīng)的IMFs 進(jìn)行平均得到最終的IMFs.EEMD 利用高斯白噪聲的頻率均勻分布特性,改變信號(hào)極值點(diǎn)的分布,填充信號(hào)存在的間歇,可以有效抑制EMD的模態(tài)混疊效應(yīng).然而,有限的集合平均次數(shù)不能完全消除添加白噪聲的影響,導(dǎo)致重構(gòu)誤差等問(wèn)題.針對(duì)該問(wèn)題,Yeh 等[10]提出CEEMD (Complementary EEMD)方法,向待分解信號(hào)中加入若干對(duì)互補(bǔ)(符號(hào)相反)白噪聲,分別進(jìn)行EMD 分解,通過(guò)集合平均時(shí)互補(bǔ)噪聲相消,有效減少了白噪聲殘留,提高了分解的完備性.類(lèi)似地,利用輔助噪聲改進(jìn)EMD 的算法還有噪聲完備集合經(jīng)驗(yàn)?zāi)B(tài)分解(Complete EEMD with adaptive noise,CEEMDAN)[11]、噪聲輔助多元經(jīng)驗(yàn)?zāi)B(tài)分解(Multivariate EMD,MEMD)[12]和快速多元經(jīng)驗(yàn)?zāi)B(tài)分解(Fast MEMD,FMEMD)[13]等.

以上基于噪聲輔助的EMD 改進(jìn)算法有效抑制了模態(tài)混疊,但是它們?cè)诜椒ㄉ瞎餐嬖谀B(tài)分裂(Mode splitting,MS)問(wèn)題.其中模態(tài)混疊定義為一個(gè)IMF 包含不同的尺度,MS 定義為一個(gè)尺度分裂到兩個(gè)或更多的IMFs 上.文獻(xiàn)[14]提出中值集合經(jīng)驗(yàn)?zāi)B(tài)分解(Median EEMD,MEEMD),利用中值算子替代算術(shù)平均算子,有效抑制了MS 問(wèn)題.然而實(shí)驗(yàn)發(fā)現(xiàn),中值取決于數(shù)據(jù)集中間位置的值,這可能導(dǎo)致集合所得的IMFs 中存在毛刺,降低信號(hào)光滑性.其次,MEEMD 需要加入大量的噪聲組數(shù)以降低重構(gòu)誤差,并且達(dá)到與EEMD 一致的重建效果理論上需要添加后者1.25 倍的噪聲組數(shù),致使MEEMD 算法存在效率低、分解完備性差等問(wèn)題.因此,該算法在實(shí)際工程應(yīng)用中仍存在局限性.

事實(shí)上,CEEMD 方法提供了一種降低殘留白噪聲的思路,即: 通過(guò)集合平均中和互補(bǔ)白噪聲,提升分解完備性.基于此,本文提出中值互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解(Median CEEMD,MCEEMD)算法.首先向原信號(hào)中添加N對(duì)互補(bǔ)白噪聲,通過(guò)EMD分解得到 2N組IMFs;然后,分別對(duì)互補(bǔ)相關(guān)的兩組IMFs (添加一對(duì)互補(bǔ)白噪聲分解得到的兩組IMFs)進(jìn)行平均運(yùn)算得到N組IMFs;最后,對(duì)平均所得N組IMFs 進(jìn)行中值運(yùn)算,輸出最終的IMFs.實(shí)驗(yàn)結(jié)果驗(yàn)證了MCEEMD 方法的有效性.

本文的創(chuàng)新點(diǎn)歸納如下: 1)在EMD 系列方法的集合過(guò)程中,首次將不同算子結(jié)合使用,用以替代方法中的單一算子;2)所提方法結(jié)合了中值算子與平均算子的各自?xún)?yōu)點(diǎn),不但能夠在一定程度上抑制MS 問(wèn)題,而且克服了MEEMD 分解完備性差、IMFs 存在毛刺現(xiàn)象等缺陷;3)針對(duì)CEEMD 的模態(tài)分裂問(wèn)題,提出了一系列概率統(tǒng)計(jì)模型量化分析該現(xiàn)象,用以驗(yàn)證所提方法的合理性與有效性.

1 基礎(chǔ)理論

1.1 EMD

經(jīng)驗(yàn)?zāi)B(tài)分解根據(jù)原信號(hào)不同尺度的趨勢(shì)和波動(dòng)分解得到一組IMFs,其中,IMF 為滿(mǎn)足以下兩個(gè)條件的函數(shù): 1)信號(hào)中極值點(diǎn)個(gè)數(shù)和過(guò)零點(diǎn)個(gè)數(shù)必須相等或僅相差一個(gè);2)由信號(hào)局部極大值和極小值點(diǎn)分別擬合出上包絡(luò)線(xiàn)和下包絡(luò)線(xiàn),上、下包絡(luò)線(xiàn)的局部均值必須為零.

對(duì)于待分解信號(hào)x(t),令i=1,則EMD 分解過(guò)程如下:

1 )找到x(t) 的所有局部極大值點(diǎn)和極小值點(diǎn).

2 ) 利用三次樣條插值函數(shù)分別對(duì)信號(hào)x(t) 的極大值和極小值進(jìn)行擬合,構(gòu)成上、下包絡(luò)線(xiàn)e+(t),e-(t),然后求取上、下包絡(luò)均值m(t)=[e+(t)+e-(t)]/2,并計(jì)算x(t) 與包絡(luò)線(xiàn)均值的差,即h(t)=x(t)-m(t).

3 )若h(t) 滿(mǎn)足IMF 的兩個(gè)條件,執(zhí)行步驟4);不滿(mǎn)足則將x(t) 替換為h(t),然后執(zhí)行步驟1).

4 ) 得到第i個(gè)IMF:di(t)=h(t),以及剩余部分ri(t)=x(t)-di(t).

5 )判斷ri(t) 是否為常值或單調(diào)函數(shù),若是,則分解結(jié)束;否則,令x(t)=ri(t),i=i+1,返回步驟1).

信號(hào)x(t) 經(jīng)EMD 分解得到一組IMFs 和殘余分量rn(t),即

式(1)說(shuō)明原信號(hào)可被分解為n個(gè)IMF 分量和1 個(gè)殘余分量rn(t).由于EMD 方法具有完備性,上述分量能精確重構(gòu)出被分解信號(hào)x(t).

1.2 EEMD 與CEEMD

針對(duì)EMD 方法存在的模態(tài)混疊問(wèn)題,文獻(xiàn)[8]提出一種輔助噪聲改進(jìn)算法——EEMD.該方法利用白噪聲的頻率分布特性,使信號(hào)的極值點(diǎn)均勻分布,填充了信號(hào)中的頻率(尺度)間歇,能夠有效抑制EMD 的模態(tài)混疊問(wèn)題.其詳細(xì)算法流程參見(jiàn)文獻(xiàn)[8],簡(jiǎn)要步驟如下:

1 )添加一組白噪聲信號(hào)至待分解信號(hào)中;

2 ) 對(duì)添加白噪聲的信號(hào)進(jìn)行EMD 分解得到IMFs;

3 )循環(huán)迭代步驟1)和步驟2),直至N組獨(dú)立分布的白噪聲全部加入;

4 )將上述分解得到的N組IMFs 進(jìn)行集合平均,得到最終的分解結(jié)果.

該算法利用集合平均的中和作用,可以一定程度地抑制由添加白噪聲引起的重構(gòu)誤差.但完全消除該問(wèn)題需要大量的集合平均次數(shù),在實(shí)際工程應(yīng)用中具有局限性.文獻(xiàn)[9]在EEMD 算法基礎(chǔ)上,提出了CEEMD 方法,步驟如下:

1 )成對(duì)地添加互補(bǔ)(符號(hào)相反)白噪聲至待分解信號(hào);

2 )～ 4)同EEMD 算法中的步驟2)～ 4).

CEEMD 方法在向待分解信號(hào)中添加一組白噪聲的同時(shí)添加一組符號(hào)相反的白噪聲,并通過(guò)對(duì)2N組IMFs 進(jìn)行平均運(yùn)算,使任意互補(bǔ)相關(guān)的兩組IMFs 中白噪聲平均相消,有效減少了由白噪聲引起的重構(gòu)誤差,提升了分解完備性.然而,EEMD與CEEMD 均存在MS 問(wèn)題,且無(wú)論改變集合尺寸還是改變添加白噪聲的幅度都不能克服該問(wèn)題.

2 CEEMD 的模態(tài)分裂

在文獻(xiàn)[14]的基礎(chǔ)上,本文利用概率模型量化CEEMD 的MS 現(xiàn)象,用于驗(yàn)證中值算子抑制MS問(wèn)題的有效性.

根據(jù)CEEMD 算法流程: 1)向原信號(hào)中添加N對(duì)互補(bǔ)噪聲,分別經(jīng)過(guò)EMD 分解得到 2N組IMFs;2)分別對(duì)其中互補(bǔ)相關(guān)的兩組IMFs 進(jìn)行平均,可以得到N組相互獨(dú)立的IMFs;3)利用平均算子處理步驟2)中N組IMFs,得到最終的IMF 集合.為方便分析,定義步驟2) 中的任意一組IMFs 為互補(bǔ)IMFs.由于N組互補(bǔ)IMFs 之間相互獨(dú)立,所以?xún)H對(duì)其中一組互補(bǔ)IMFs 進(jìn)行量化分析即可準(zhǔn)確反映CEEMD 算法的MS 問(wèn)題.具體分析過(guò)程如下.

使用正弦信號(hào)x(t)=sin(2πft) 作為原信號(hào),其中頻率f=2 Hz,并添加一對(duì)互補(bǔ)白噪聲構(gòu)成噪聲輔助信號(hào)

其中,白噪聲w(t)～N(0,1),ε=0.2std{x(t)}表示添加白噪聲的幅度(信號(hào)x(t) 標(biāo)準(zhǔn)差的0.2 倍) (若無(wú)特殊說(shuō)明下同).信號(hào)(2)經(jīng)過(guò)EMD 分解得到兩組互補(bǔ)相關(guān)的IMFs,再進(jìn)行平均運(yùn)算可獲得一組互補(bǔ)IMFs,如圖1 所示.可以發(fā)現(xiàn): 1)d4分量與正弦信號(hào)x(t) 相似程度最高,視為主尺度.2)d3分量?jī)?nèi)包含少量與x(t) 相關(guān)的半波,這說(shuō)明出現(xiàn)了MS現(xiàn)象.3) EMD 算法由信號(hào)的極值點(diǎn)驅(qū)動(dòng)分解,并且信號(hào)x(t) 中一個(gè)半波內(nèi)僅存在一個(gè)極值點(diǎn).當(dāng)某一半波的極值點(diǎn)附近存在噪聲間歇時(shí),該極值點(diǎn)將參與到噪聲分量的篩選過(guò)程中,導(dǎo)致噪聲分量中出現(xiàn)該極值點(diǎn)所對(duì)應(yīng)的半波(如d3分量中虛線(xiàn)所圈部分).因此,分裂成分主要以半波的形式呈現(xiàn).4)分裂成分(虛線(xiàn)所圈部分) 與原信號(hào)的頻率、波形相似,進(jìn)而以原信號(hào)的半波為單位可以量化MS 現(xiàn)象.

圖1 EMD 分解噪聲輔助信號(hào)得到的前5 個(gè)互補(bǔ)IMFsFig.1 The first five complementary IMFs obtained from the noise-assisted signal through EMD

基于以上結(jié)論,對(duì)MS 現(xiàn)象進(jìn)行統(tǒng)計(jì)分析.假設(shè)x(t) 包含wx(f) 個(gè)半波,根據(jù)半波的位置從整個(gè)觀察區(qū)間中劃分出wx(f) 個(gè)窗口,如圖1 虛線(xiàn)劃分所示.如果某IMF 的半波窗口內(nèi)的信號(hào)絕對(duì)值之和大于噪聲窗口閾值,可以認(rèn)為該窗口中的信號(hào)與原信號(hào)相關(guān),進(jìn)而引出兩個(gè)關(guān)鍵的概率參數(shù).

1 )pi(f): 互補(bǔ)IMFs 中第i個(gè)分量?jī)?nèi)存在原信號(hào)的概率,i=2,3,4.對(duì)于每個(gè)頻率f,當(dāng)互補(bǔ)IMFs中第i個(gè)分量?jī)?nèi)存在任意窗口的指標(biāo)高于噪聲閾值,認(rèn)為該分量包含原信號(hào)并計(jì)數(shù).進(jìn)行1 000 次獨(dú)立實(shí)驗(yàn)后,設(shè)mi為包含原信號(hào)的分量個(gè)數(shù),通過(guò)計(jì)算其所占百分比可以估計(jì)得到pi(f),即:pi(f)=mi/1 000.其中,互補(bǔ)IMFs 中第i個(gè)分量的噪聲閾值為

式中,i=2,3,4,常數(shù)ρ,β分別為2.01,0.719.T1通過(guò)計(jì)算第1 個(gè)互補(bǔ)IMF 中所有半波窗口內(nèi)數(shù)據(jù)的絕對(duì)值之和,然后取平均獲得.

2 )ri(f): 互補(bǔ)IMFs 中第i個(gè)分量?jī)?nèi)存在原信號(hào)x(t) 的時(shí)間占比,即

其中,wi(f) 為第i個(gè)互補(bǔ)IMFs 中與原信號(hào)x(t) 相關(guān)的半波窗口個(gè)數(shù),wx(f) 為x(t) 全部半波窗口的數(shù)目.

上述互補(bǔ)IMFs 的pi(f) 和ri(f) 如圖2 所示,可以發(fā)現(xiàn): 1)低頻范圍內(nèi) (1.5 Hz≤f ≤2.1 Hz),d4對(duì)應(yīng)的p4(f) 和r4(f) 值最大,因此可以將d4分量視為主尺度;2)在頻率增加至3.8 Hz 的過(guò)程中,p3(f)和r3(f) 急劇增大,逐漸接近100%,即主尺度逐漸由d4分量轉(zhuǎn)移至d3分量;3) 在高頻范圍內(nèi),d2分量對(duì)應(yīng)的p2(f) 和r2(f) 增幅趨勢(shì)與之前d3分量相似,代表主尺度逐漸轉(zhuǎn)移至d2分量中.基于上述分析可以發(fā)現(xiàn),原信號(hào)頻率的改變使主尺度在不同分量之間進(jìn)行轉(zhuǎn)移,從而不同程度地導(dǎo)致了發(fā)生MS 現(xiàn)象.

圖2 互補(bǔ)IMFs (由噪聲輔助信號(hào)得到)的pi(f) 和 r i(f) 曲線(xiàn)Fig.2 The curves p i(f) and r i(f) corresponding to the complementary IMFs (obtained from the noise-assisted signal)

進(jìn)一步定義Pi(f) 為互補(bǔ)IMFs 中第i個(gè)分量?jī)?nèi)任意半波窗口中存在原信號(hào)x(t) 的概率,表達(dá)式為

通過(guò)pi(f) 和ri(f) 計(jì)算得到的Pi(f) 如圖3 所示.從圖3 中可以明顯看出,Pi(f) 與EMD 的濾波器組特性[9]相似.

在此基礎(chǔ)上,可以分析CEEMD 分解結(jié)果中各個(gè)分量存在原信號(hào)的程度.通過(guò)N次獨(dú)立實(shí)驗(yàn),對(duì)得到的N組互補(bǔ)IMFs 進(jìn)行平均運(yùn)算,即可得到CEEMD 輸出的IMFs.其中,對(duì)于第i個(gè)分量而言,N次獨(dú)立實(shí)驗(yàn)中包含x(t) 的互補(bǔ)IMFs 有Pi(f)×N個(gè),然后取平均得到其出現(xiàn)的比率PRi(f),表達(dá)式為

由于N組互補(bǔ)IMFs 之間相互獨(dú)立,所以通過(guò)平均運(yùn)算得到的PRi(f) 與Pi(f) 曲線(xiàn)一致.這說(shuō)明在N次獨(dú)立實(shí)驗(yàn)中,理論上有PRi(f) 占比的分量將包含部分x(t).此外,集合平均后得到的第i個(gè)IMF包含x(t) 的程度與PRi(f) 成正比關(guān)系,進(jìn)而可以使用PRi(f) 量化MS 現(xiàn)象.

通過(guò)上述分析發(fā)現(xiàn),MS 現(xiàn)象主要出現(xiàn)在主尺度前、后兩個(gè)分量中,因此通過(guò)計(jì)算主尺度及其前、后兩個(gè)分量的PRi(f) 即可量化MS 程度.針對(duì)本文研究的頻率范圍,設(shè)計(jì)MS 量化指標(biāo)MSD(f),即

其中,MSD(f) 值越大,x(t) 分裂的幅度和程度越大.CEEMD 算法使用平均算子處理N組互補(bǔ)IMFs,得到最終IMFs,其MSD(f) 曲線(xiàn)如圖4 所示.對(duì)照?qǐng)D3與圖4 容易發(fā)現(xiàn): 首先,在EMD 濾波器組的重疊頻帶內(nèi) (1.6 Hz≤f ≤1.8 Hz 或4.2 Hz≤f ≤4.9 Hz),使用平均算子獲得的MSD(f) 曲線(xiàn)出現(xiàn)明顯的峰值,即該頻帶內(nèi)MS 現(xiàn)象最為嚴(yán)重;其次,在重疊頻帶以外,MSD(f) 曲線(xiàn)呈現(xiàn)明顯的下降趨勢(shì),尤其在2.9 Hz≤f ≤3.1 Hz 頻帶內(nèi),曲線(xiàn)接近于0.但是,無(wú)論頻率如何改變,使用平均算子都將出現(xiàn)MS 現(xiàn)象.

然而,使用中值算子替換上述集合過(guò)程中的平均算子可以去除MSD(f) 的故障點(diǎn),如圖4 所示.其根本原因是: 對(duì)于第i個(gè)分量的任意半波窗口而言,經(jīng)過(guò)N次獨(dú)立實(shí)驗(yàn),只要包含x(t) 的分量數(shù)目不超過(guò)N/2,集合過(guò)程使用中值算子即可有效剔除半波窗口中的x(t) 成分.所以,當(dāng)原始的MSD(f)＜50%時(shí),使用中值算子得到的MSD(f) 曲線(xiàn)將被置零,可以有效抑制MS 現(xiàn)象.

3 中值CEEMD

3.1 中值CEEMD 算法

第2 節(jié)論證了中值算子抑制CEEMD 算法MS問(wèn)題的有效性.并且,鑒于使用集合平均可以中和互補(bǔ)白噪聲,彌補(bǔ)單一使用中值算子的缺陷,本文提出了中值CEEMD (MCEEMD)方法,步驟如下:

1 )初始化,令n=1.

2 )添加成對(duì)互補(bǔ)白噪聲至原始信號(hào)中

其中,w(t)～N(0,1),添加白噪聲的幅度為ε=0.2std{x(t)}.

式中,i=1,···,M,M是EMD 分解得到的IMF個(gè)數(shù).

4 ) 令n=n+1,如果n≤N,重復(fù)步驟2)和步驟3),其中,N是添加的互補(bǔ)噪聲組數(shù).

5 )對(duì)步驟4)輸出的N組互補(bǔ)IMFs 進(jìn)行中值運(yùn)算,得到最終的分解結(jié)果:.其中,第i個(gè)分量計(jì)算為

圖5 為該算法框圖.

3.2 MCEEMD 模態(tài)分裂率

第2 節(jié)利用概率模型從理論上分析了CEEMD的MS 程度.本節(jié)將基于仿真信號(hào)x(t) 進(jìn)一步對(duì)比MCEEMD、MEEMD 和CEEMD 的MS 程度,以驗(yàn)證中值算子抑制MS 問(wèn)題的有效性.

根據(jù)MS 的定義,主尺度所在IMF 的前、后兩個(gè)分量中可能存在原信號(hào)的分裂成分,因此實(shí)際MS 程度可由主尺度前、后兩個(gè)分量與原信號(hào)之間的比值來(lái)反映,定義模態(tài)分裂率[14]為

其中,i=1,···,M -1,M是分解得到的IMF 的個(gè)數(shù).

使用MCEEMD、MEEMD 和CEEMD 分解仿真信號(hào)x(t)=sin(2πft),初始化條件相同: 輔助噪聲幅度ε=0.2std{x(t)},總采樣點(diǎn)數(shù)512 個(gè),采樣頻率fs=100 Hz,頻率f=1.1,1.2,···,5.5 Hz,其中間隔為0.1 Hz.為測(cè)試集合尺寸N對(duì)MS 問(wèn)題造成的影響,采用N=10,20,50,100 分別進(jìn)行1 000 次獨(dú)立實(shí)驗(yàn)并獲得平均結(jié)果.圖6 顯示了3 種方法在不同集合尺寸下的MSR(f) 曲線(xiàn),可以看出: 1) MCEEMD實(shí)際模態(tài)分裂率與理論分析中圖4 的結(jié)果高度一致;2) MCEEMD 與MEEMD 相比,在濾波器組的重疊頻帶以外,MCEEMD 基本消除了噪聲的影響,使其MSR(f)曲線(xiàn)低于MEEMD 且趨近于零.但是在濾波器組的重疊頻帶內(nèi),MCEEMD 的MSR(f) 曲線(xiàn)略高于后者.其原因是,方法中的平均過(guò)程導(dǎo)致MS 程度略有增加;3) 與CEEMD方法相比,本文提出的方法能夠有效抑制不同集合尺寸下的MS 問(wèn)題.并且,本文方法的改善效果隨集合尺寸的增加而更加明顯.尤其當(dāng)N=100 時(shí),在濾波器組重疊頻帶以外,本文方法的MSR(f) 曲線(xiàn)以最快速率下降至零,這說(shuō)明MS 問(wèn)題得到了有效抑制.

圖6 不同集合尺寸下CEEMD、MEEMD 和MCEEMD 的 M SR(f) 曲線(xiàn)Fig.6 M SR(f) curves for different ensemble sizes within CEEMD,MEEMD and MCEEMD

3.3 MCEEMD 分解完備性

文獻(xiàn)[10]指出,成對(duì)地添加互補(bǔ)噪聲可以有效減小重建信號(hào)中的噪聲殘留,保證分解完備性.基于此,進(jìn)行數(shù)值實(shí)驗(yàn)對(duì)比分析MCEEMD、MEEMD 和CEEMD 的分解完備性,驗(yàn)證本文方法融合兩種算子(平均、中值算子)相對(duì)單一使用中值算子的有效性.

本次實(shí)驗(yàn)使用的仿真信號(hào)仍是x(t)=sin(2πft),其中,f=3.5 Hz,總樣本點(diǎn)數(shù)512 個(gè),輔助噪聲幅值為ε=0.2std{x(t)}.通過(guò)計(jì)算重建誤差的標(biāo)準(zhǔn)差與輔助噪聲幅值的比值,可以反映噪聲的殘留量.據(jù)此,定義噪聲殘留率SDR(N)[14]為

本次數(shù)值實(shí)驗(yàn)使用不同數(shù)值的集合尺寸檢驗(yàn)兩種方法的分解完備性,即N=10,···,100,其中間隔為5.最終的SDR(N) 結(jié)果由1 000 次獨(dú)立的實(shí)驗(yàn)取平均獲得,如圖7 所示.

圖7 MCEEMD、MEEMD 和CEEMD 在不同集合尺寸下的 S DR(N) 曲線(xiàn)Fig.7 S DR(N) curves for different ensemble sizes within MCEEMD,MEEMD and CEEMD

從圖7 中可以發(fā)現(xiàn),1) CEEMD 的SDR(N) 曲線(xiàn)始終最低,分解完備性最優(yōu);2) MCEEMD 由于結(jié)合了平均算子,相比MEEMD 的SDR(N) 曲線(xiàn)降低了約75%,分解完備性有明顯的提升.其根本原因是MEEMD 方法單一使用中值算子,導(dǎo)致噪聲殘留量較大,而本文所提方法利用平均算子中和了互補(bǔ)白噪聲,所以重建信號(hào)與原信號(hào)吻合程度更高,分解完備性?xún)?yōu)于MEEMD.

3.4 仿真分析

本節(jié)使用式(15)所示仿真信號(hào)進(jìn)一步驗(yàn)證本文方法的有效性.

其中,設(shè)定原始信號(hào)為:x(t)=cos(2πft),f=3.5 Hz,n(t)是兩段均值為0 且幅值為0.2 的高斯白噪聲,w(t)～N(0,0.2),t∈(100,200)∪(300,400).信號(hào)總采樣點(diǎn)512 個(gè),采樣頻率fs=100 Hz.圖8(a)～ 8(d)依次對(duì)應(yīng)EEMD、CEEMD、MEEMD 和MCEEMD 分解該仿真信號(hào)所得結(jié)果 (由于篇幅原因只顯示到d5分量).由圖8 可以發(fā)現(xiàn),4 種方法的分解結(jié)果中d3分量為主尺度.然而,EEMD 和CEEMD單獨(dú)使用平均算子,導(dǎo)致d2分量中存在部分原始信號(hào)x(t),表現(xiàn)出明顯的MS 問(wèn)題.而采用中值算子的MEEMD 和MCEEMD 兩種方法,所得到的d3分量與原始信號(hào)x(t) 高度吻合,既精確地提取出周期信號(hào),又實(shí)現(xiàn)了周期信號(hào)與間歇信號(hào)的尺度分離.

圖8 4 種方法分解仿真信號(hào)所得的前5 個(gè)IMFFig.8 The first five IMFs obtained by decomposing the simulated signal by four methods

進(jìn)一步地,上述方法的分解結(jié)果中,d3分量與x(t)的皮爾遜相關(guān)系數(shù)分別為0.9783,0.9066,0.9972,0.9978.顯然,采用中值算子的兩種方法能夠更好地從間歇噪聲中提取原始信號(hào),驗(yàn)證了中值算子抑制MS 問(wèn)題的有效性.

另外,本案例還對(duì)比了使用中值算子的兩種方法的性能.圖9(a)和圖9(b)分別顯示了MEEMD與MCEEMD 分解得到的d2分量.從圖9 中可以明顯發(fā)現(xiàn),MEEMD 的d2分量中存在大量毛刺,導(dǎo)致分量的精確性差.其根本原因是中值算子根據(jù)數(shù)據(jù)集內(nèi)中間位置的值來(lái)確定代表值,缺乏對(duì)數(shù)據(jù)集的概括能力,從而降低了信號(hào)光滑性.本文所提方法在中值算子的基礎(chǔ)上融合了平均算子,有效地改善了中值算子造成的毛刺現(xiàn)象.并且,兩種方法分解仿真信號(hào)的SDR(N) 曲線(xiàn)與圖7 一致,即MCEEMD的分解完備性?xún)?yōu)于MEEMD.

圖9 MEEMD、MCEEMD 分解結(jié)果中的 d2 模態(tài)Fig.9 The d2 mode in the decomposition results of MEEMD and MCEEMD

對(duì)上述仿真信號(hào)的分析初步表明,本文提出的MCEEMD 方法不僅可以有效抑制MS 現(xiàn)象,而且還能改善單一使用中值算子產(chǎn)生的毛刺現(xiàn)象和分解完備性差等問(wèn)題.此外,本文所提方法對(duì)“純凈信號(hào)加瞬態(tài)振蕩信號(hào)”的分解結(jié)果與上述結(jié)果高度一致.

4 實(shí)際案例分析

鑒于MCEEMD 在仿真實(shí)驗(yàn)中取得的良好性能,因此可期待更多的實(shí)際應(yīng)用.為了驗(yàn)證所提出方法在實(shí)際應(yīng)用中的有效性,本文研究了兩個(gè)不同領(lǐng)域的典型應(yīng)用,即: 1) 機(jī)械故障信號(hào)特征提取;2) 超聲多普勒血流信號(hào)的干擾抑制.由于機(jī)械系統(tǒng)中的故障涉及到部件的相對(duì)運(yùn)動(dòng)和磨損,同時(shí)監(jiān)測(cè)超聲血流信號(hào)依賴(lài)于多普勒效應(yīng),因此這兩種信號(hào)均呈現(xiàn)顯著的周期性,可用來(lái)測(cè)試分解算法的實(shí)際性能.

本文算法實(shí)驗(yàn)硬件環(huán)境是筆記本電腦,主要配置為Intel i5-8300H,2.3 GHz,以及8 GB 的內(nèi)存;軟件環(huán)境是32 位的Windows10 系統(tǒng),仿真運(yùn)行工具是MATLAB2020a.

4.1 機(jī)械故障信號(hào)

目前,國(guó)內(nèi)外對(duì)單向閥故障診斷研究的主要步驟包括: 信號(hào)預(yù)處理、特征提取和故障識(shí)別[14].其中,信號(hào)預(yù)處理階段通常使用EMD 系列方法分解原始信號(hào),對(duì)IMFs 進(jìn)行處理后得到故障的特征向量;然后,利用特征提取方法提取故障特征信息進(jìn)行訓(xùn)練;最后建立故障診斷模型,實(shí)現(xiàn)故障識(shí)別.在處理IMFs 過(guò)程中,為提升故障特征提取的準(zhǔn)確性,通常將EMD 系列方法分解結(jié)果中虛假分量、背景噪聲進(jìn)行分辨并去除[15-16].但是,EMD 系列方法存在MS問(wèn)題,可能導(dǎo)致部分故障特征分裂至被去除的分量中,降低信號(hào)預(yù)處理的效果.基于此,對(duì)比分析現(xiàn)有的噪聲輔助EMD 方法與本文方法在信號(hào)預(yù)處理階段的性能,驗(yàn)證本文方法在實(shí)際應(yīng)用中的有效性.

單向閥磨損擊穿狀態(tài)下的振動(dòng)信號(hào)采樣頻率為2 560 Hz,采樣數(shù)據(jù)長(zhǎng)度為2 048 個(gè)采樣點(diǎn)[17].使用EEMD、CEEMD、MEEMD 和MCEEMD 四種方法處理該信號(hào),并設(shè)置集合尺寸N=100,添加白噪聲幅度為振動(dòng)信號(hào)標(biāo)準(zhǔn)差的0.2 倍.由于分解結(jié)果中的MS 現(xiàn)象無(wú)法直觀判斷,所以通過(guò)計(jì)算各個(gè)分量的功率譜密度(Power spectral density,PSD)來(lái)反映該現(xiàn)象.

圖10(a)～ 10(d)分別對(duì)應(yīng)EEMD、CEEMD、MEEMD 和MCEEMD 分解結(jié)果中各個(gè)分量的PSD曲線(xiàn).從圖10 中可以看到: 1) EEMD 的MS 問(wèn)題最為嚴(yán)重.圖10(a)中,IMF4 分量在30 Hz 頻率處出現(xiàn)50%左右的功率泄漏,IMF2 分量在157 Hz 頻率處出現(xiàn)約37%的功率泄漏;2) CEEMD 在30 Hz頻率處表現(xiàn)出的模態(tài)分裂程度與EEMD 相同,并且類(lèi)似問(wèn)題同樣出現(xiàn)在80 Hz、157 Hz 頻率處;3) MEEMD 在 30 Hz、80 Hz 頻率處抑制了MS 問(wèn)題,但在157 Hz、194 Hz 頻率處仍表現(xiàn)出約23%和33%的功率泄漏;4) MCEEMD 在上述4 個(gè)頻率處均未出現(xiàn)明顯的MS 現(xiàn)象.

圖10 4 種方法分解結(jié)果的PSDFig.10 The PSD curves of the decomposition results from the four methods

此外,僅MEEMD 將分解得到的主要信號(hào)功率(20～ 40 Hz 頻率成分)集中在IMF5.這是由于該方法在IMF2 分量出現(xiàn)了一定程度的模態(tài)分裂現(xiàn)象,使得160～ 200 Hz 的成分部分被分裂到IMF3分量中,這將導(dǎo)致連鎖反應(yīng),使得之后所有的成分都分別往下一層IMF 分裂.圖10(a)和圖10(b)已經(jīng)表明,IMF4 和IMF5 包含了功率相當(dāng)?shù)某煞?當(dāng)IMF4 部分分裂至IMF5 后,將導(dǎo)致20～ 40 Hz頻率的成分主要分布在IMF5,使IMF5 占主導(dǎo)地位.相應(yīng)地,使用中值算子將會(huì)在集合過(guò)程中把IMF4歸于IMF5 中.該現(xiàn)象可能對(duì)后續(xù)的特征提取產(chǎn)生不利的影響.

文獻(xiàn)[18]指出,單向閥的故障特征主要出現(xiàn)在10～ 200 Hz.在該頻帶內(nèi),MCEEMD 方法有效抑制了MS 問(wèn)題,分解得到的IMFs 相對(duì)完整地保留了故障特征信息,使提取到的特征更具代表性,利于后續(xù)的故障識(shí)別研究.

4.2 超聲多普勒血流信號(hào)

超聲多普勒血流信號(hào)中經(jīng)常存在管壁信號(hào)的干擾,尤其在血流和管壁的臨界處,管壁信號(hào)會(huì)嚴(yán)重干擾血流信號(hào)中的低頻成分.因此,有效抑制管壁干擾對(duì)提高血流流速檢測(cè)的精度具有重要作用[19-20].

EMD 系列方法常用于超聲多普勒血流信號(hào)中血流成分的提取,主要包含三大步驟[21]: 1) 使用EMD 系列方法分解超聲多普勒血流信號(hào),得到一組IMFs.2)由于管壁信號(hào)功率通常比血流信號(hào)功率大約20 dB,所以通過(guò)管壁與血流功率比(Wall blood signal ratio,WBSR)[22]突變的位置,即可確定血流與管壁的IMF 分界點(diǎn).其中,分界點(diǎn)以前的高頻、低幅值IMFs 屬于血流信號(hào),而分界點(diǎn)以后的低頻、高幅值IMFs 屬于管壁信號(hào).3)將分界點(diǎn)以前的IMFs 相加構(gòu)成血流成分,并且,使用軟閾值降噪法對(duì)提取到的血流成分進(jìn)一步細(xì)分.然而,在步驟2)中,傳統(tǒng)EMD 方法的MS 現(xiàn)象可能導(dǎo)致分界點(diǎn)附近的IMF 中包含不同尺度的信息,即: 血流成分分裂至管壁成分中,或者,管壁成分分裂至血流成分中.這將影響提取血流成分的精度.基于此,對(duì)比分析現(xiàn)有噪聲輔助EMD 方法與本文方法提取血流成分的性能,驗(yàn)證本文方法提取血流成分的有效性.步驟詳述如下.

所分析的超聲多普勒人體頸動(dòng)脈血流信號(hào)從管壁和血流臨界點(diǎn)處測(cè)得,其中,采樣頻率為6 666 Hz,數(shù)據(jù)長(zhǎng)度為4 048 個(gè)樣本點(diǎn).使用MCEEMD 方法分解血流信號(hào)結(jié)果如圖11 所示.通過(guò)WBSR 方法判斷出d3分量為臨界點(diǎn),因此,將前三個(gè)IMFs 分量相加構(gòu)成血流成分.按上述步驟,分別獲得EEMD、CEEMD 和MEEMD 提取到的血流成分.

圖11 MCEEMD 分解血流信號(hào)所得的前8 個(gè)分量Fig.11 The first eight components of the blood flow signal decomposed by MCEEMD

由于血流成分與管壁干擾之間存在頻率分界點(diǎn),采用高通濾波器提取血流成分時(shí)的歸一化截止頻率通常設(shè)為f=0.09 Hz[23].因此,計(jì)算出原信號(hào)與各方法提取血流成分的歸一化PSD,如圖12 和圖13 所示.在此基礎(chǔ)上,可以通過(guò)兩個(gè)方面來(lái)評(píng)價(jià)方法的性能: 1) 圖13 中f=0.09 Hz 左側(cè)部分為管壁干擾,所以該部分的面積大小可以反映方法抑制干擾的性能;2) 右側(cè)部分可認(rèn)為是較為純凈的血流,通過(guò)與原信號(hào)該部分的PSD 曲線(xiàn)進(jìn)行互相關(guān)分析,可以評(píng)估方法提取血流的精度.

圖12 原始信號(hào)的頻率歸一化功率譜Fig.12 The frequency normalized PSD of the original signal

圖13 EEMD、CEEMD、MEEMD 和MCEEMD 提取的血流成分頻率歸一化功率譜Fig.13 Frequency normalized PSD of the blood flow component extracted by EEMD,CEEMD,MEEMD,and MCEEMD,respectively

對(duì)照?qǐng)D12 與圖13 容易發(fā)現(xiàn): 1) 在截止頻率f=0.09 Hz 左側(cè),MCEEMD、MEEMD 的功率譜面積較小,有效抑制了管壁干擾;2) 在其右側(cè),MEEMD 的功率譜與原信號(hào)相差較大,說(shuō)明該方法提取的血流成分嚴(yán)重失真.

進(jìn)一步地,使用皮爾遜相關(guān)系數(shù)(Pearson correlation coefficient,PCC)、均方根誤差(Root mean square error,RMSE)以及PSD 面積比(PSD area ratio)來(lái)定量分析方法性能.首先,對(duì)于PSD 中f=0.09 Hz右側(cè)部分,分別計(jì)算4 種方法提取的血流成分與原信號(hào)之間的PCC、RMSE.其次,對(duì)于f=0.09Hz 左側(cè)部分,分別計(jì)算4 種方法提取的血流成分與原信號(hào)之間的面積比,得到PSD 面積比.以上計(jì)算結(jié)果如表1 所示.

表1 4 種方法的性能指標(biāo)Table 1 Performance indicators of the four methods

從表1 中可以看出,MCEEMD 和CEEMD 方法的PCC 接近于1,說(shuō)明提取的血流成分與原信號(hào)吻合程度較高,并且兩種方法之間的差僅為萬(wàn)分位級(jí),可以忽略不計(jì).相比之下,MEEMD 由于單獨(dú)使用中值算子導(dǎo)致其相關(guān)系數(shù)最低.其次,MCEEMD方法的RMSE 相比MEEMD 降低了87%左右,方法的分解完備性得到大幅度提升,僅次于CEEMD.更重要的是,MCEEMD 的PSD 面積比在4 種方法中最小,分別比EEMD、CEEMD、MEEMD 減小了約50%,33%,42%,說(shuō)明該方法有效抑制了管壁干擾.

此外,表2 給出了4 種方法的計(jì)算時(shí)間.可以發(fā)現(xiàn),在方法初始參數(shù)(噪聲幅值、集合次數(shù)、篩選停止指標(biāo))相同的情況下,CEEMD 和MCEEMD的計(jì)算時(shí)間基本相同,約為EEMD 和MEEMD 的兩倍,其原因是: 前兩種方法之間的區(qū)別僅在于處理IMFs 集合的算子差異,調(diào)用EMD 子算法的次數(shù)相同,所以?xún)煞N方法計(jì)算時(shí)間高度相似;然而,由于它們加入了互補(bǔ)噪聲,使得調(diào)用EMD 子算法的次數(shù)是EEMD、MEEMD 的兩倍,計(jì)算時(shí)間也約為兩倍.但是,本文所提方法與MEEMD 相比,分解完備性有75%的提升,設(shè)置較少的集合次數(shù)即可獲得優(yōu)于MEEMD 的分解結(jié)果.當(dāng)要求兩種方法在達(dá)到相同的完備性標(biāo)準(zhǔn)時(shí),本文所提方法要求的獨(dú)立分解次數(shù)更少,運(yùn)行時(shí)間更短.

表2 4 種方法的計(jì)算時(shí)間Table 2 Calculation time of the four methods

綜上所述,本文方法既能精確地提取血流成分,又抑制了管壁干擾,對(duì)后續(xù)的噪細(xì)分、平均血流估計(jì)(血流測(cè)速)研究起到積極作用.

基于上述兩個(gè)實(shí)際案例可以看出,在實(shí)際應(yīng)用中EMD 系列方法常出現(xiàn)模態(tài)混疊、模態(tài)分裂現(xiàn)象,而MEEMD 又存在分解完備性差、IMFs 中出現(xiàn)毛刺等固有問(wèn)題.本文提出的MCEEMD 可以有效地抑制上述問(wèn)題,該方法對(duì)大部分存在模態(tài)分裂的應(yīng)用場(chǎng)合有效.

此外,我們還針對(duì)不同采樣頻率的仿真信號(hào)、工業(yè)控制過(guò)程振蕩信號(hào)(振蕩檢測(cè))[24]和心電信號(hào)(信號(hào)消噪)[25]進(jìn)行實(shí)驗(yàn),結(jié)果與本文結(jié)論高度一致,說(shuō)明本文所提方法對(duì)采樣頻率并不敏感,并且在多種領(lǐng)域中均具備實(shí)踐性,驗(yàn)證了方法的可推廣性及應(yīng)用價(jià)值.

5 結(jié)束語(yǔ)

針對(duì)MEEMD 單一使用中值算子造成的分解完備性差、IMFs 存在毛刺等問(wèn)題,本文提出MCEEMD 方法.該方法在使用中值算子的基礎(chǔ)上,利用集合平均中和了互補(bǔ)白噪聲.仿真實(shí)驗(yàn)結(jié)果表明,該方法既抑制了MS 問(wèn)題,又克服了單一使用中值算子造成的缺陷,提升了方法的分解完備性,使分解得到的IMF 更具物理意義.最后,使用兩個(gè)不同工程領(lǐng)域的案例,驗(yàn)證了本文所提方法相對(duì)EMMD、CEEMD 和MEEMD 的優(yōu)越性.因此,本文建議將MCEEMD 作為MEEMD 的標(biāo)準(zhǔn)形式.此外,針對(duì)特殊的實(shí)際信號(hào),可以嘗試更先進(jìn)的平均算子,例如: 模式平均、幾何平均和加權(quán)平均等.下一步的主要工作是研究其他平均算子的優(yōu)點(diǎn),以及如何完全抑制MS 問(wèn)題.