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

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

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

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

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

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

1 基礎(chǔ)理論

1.1 EMD

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

對于待分解信號x(t),令i=1,則EMD 分解過程如下:

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

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

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

4 ) 得到第i個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).

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

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

1.2 EEMD 與CEEMD

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

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

2 ) 對添加白噪聲的信號進行EMD 分解得到IMFs;

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

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

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

1 )成對地添加互補(符號相反)白噪聲至待分解信號;

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

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

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

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

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

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

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

圖1 EMD 分解噪聲輔助信號得到的前5 個互補IMFsFig.1 The first five complementary IMFs obtained from the noise-assisted signal through EMD

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

1 )pi(f): 互補IMFs 中第i個分量內(nèi)存在原信號的概率,i=2,3,4.對于每個頻率f,當互補IMFs中第i個分量內(nèi)存在任意窗口的指標高于噪聲閾值,認為該分量包含原信號并計數(shù).進行1 000 次獨立實驗后,設(shè)mi為包含原信號的分量個數(shù),通過計算其所占百分比可以估計得到pi(f),即:pi(f)=mi/1 000.其中,互補IMFs 中第i個分量的噪聲閾值為

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

2 )ri(f): 互補IMFs 中第i個分量內(nèi)存在原信號x(t) 的時間占比,即

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

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

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

進一步定義Pi(f) 為互補IMFs 中第i個分量內(nèi)任意半波窗口中存在原信號x(t) 的概率,表達式為

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

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

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

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

其中,MSD(f) 值越大,x(t) 分裂的幅度和程度越大.CEEMD 算法使用平均算子處理N組互補IMFs,得到最終IMFs,其MSD(f) 曲線如圖4 所示.對照圖3與圖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)明顯的峰值,即該頻帶內(nèi)MS 現(xiàn)象最為嚴重;其次,在重疊頻帶以外,MSD(f) 曲線呈現(xiàn)明顯的下降趨勢,尤其在2.9 Hz≤f ≤3.1 Hz 頻帶內(nèi),曲線接近于0.但是,無論頻率如何改變,使用平均算子都將出現(xiàn)MS 現(xiàn)象.

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

3 中值CEEMD

3.1 中值CEEMD 算法

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

1 )初始化,令n=1.

2 )添加成對互補白噪聲至原始信號中

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

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

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

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

圖5 為該算法框圖.

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

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

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

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

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

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

3.3 MCEEMD 分解完備性

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

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

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

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

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

3.4 仿真分析

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

其中,設(shè)定原始信號為: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).信號總采樣點512 個,采樣頻率fs=100 Hz.圖8(a)～ 8(d)依次對應(yīng)EEMD、CEEMD、MEEMD 和MCEEMD 分解該仿真信號所得結(jié)果 (由于篇幅原因只顯示到d5分量).由圖8 可以發(fā)現(xiàn),4 種方法的分解結(jié)果中d3分量為主尺度.然而,EEMD 和CEEMD單獨使用平均算子,導(dǎo)致d2分量中存在部分原始信號x(t),表現(xiàn)出明顯的MS 問題.而采用中值算子的MEEMD 和MCEEMD 兩種方法,所得到的d3分量與原始信號x(t) 高度吻合,既精確地提取出周期信號,又實現(xiàn)了周期信號與間歇信號的尺度分離.

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

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

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

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

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

4 實際案例分析

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

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

4.1 機械故障信號

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

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

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

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

此外,僅MEEMD 將分解得到的主要信號功率(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 包含了功率相當?shù)某煞?當IMF4 部分分裂至IMF5 后,將導(dǎo)致20～ 40 Hz頻率的成分主要分布在IMF5,使IMF5 占主導(dǎo)地位.相應(yīng)地,使用中值算子將會在集合過程中把IMF4歸于IMF5 中.該現(xiàn)象可能對后續(xù)的特征提取產(chǎn)生不利的影響.

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

4.2 超聲多普勒血流信號

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

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

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

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

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

圖12 原始信號的頻率歸一化功率譜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

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

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

表1 4 種方法的性能指標Table 1 Performance indicators of the four methods

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

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

表2 4 種方法的計算時間Table 2 Calculation time of the four methods

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

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

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

5 結(jié)束語

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