EEMD-四線性平行因子的機(jī)械多故障欠定盲分離方法

EEMD-四線性平行因子的機(jī)械多故障欠定盲分離方法

朱亞靜，李澤東，李志農(nóng),谷士鵬，馬亞平（1.南昌航空大學(xué)無損檢測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，南昌330063；2.中國飛行試驗(yàn)研究院，西安710089）平行因子（ParallelFactor,PARAFAC）在寬松條件下具有分解唯一性的優(yōu)勢(shì)，在機(jī)械故障診斷中獲得了成功應(yīng)用。李曉明等[1]將傳統(tǒng)奇異譜分解結(jié)合標(biāo)準(zhǔn)PARAFAC張量分解方法進(jìn)行擴(kuò)展，對(duì)機(jī)械故障特征進(jìn)行提取。苗育茁等[2]利用連續(xù)小波變換構(gòu)造3階張量，通過PARAFAC對(duì)其進(jìn)行分解，提出了一種多尺度平行因子分析算法，可以有效地從機(jī)械非線性多故障狀態(tài)中提取特征，進(jìn)行故障診斷。楊誠等[3]結(jié)合Volterra模型與PARAFAC構(gòu)建了故障預(yù)測(cè)模型，主要解決了PARAFAC在分解機(jī)械非線性系統(tǒng)信號(hào)時(shí)估計(jì)參數(shù)較多的問題。李志農(nóng)等[4]結(jié)合變分模態(tài)分解(VariationalModeDecomposition，VMD)和PARAFAC各自優(yōu)點(diǎn)，提出了基于VMDPARAFAC欠定盲分離方法，并將其應(yīng)用到滾動(dòng)軸承的復(fù)合故障特征提取中。孫宜權(quán)等[5]采用自適應(yīng)PARAFAC和BSS算法對(duì)柴油機(jī)噴油多故障進(jìn)行診斷。然而，上述方法中的PARAFAC分解都是通過構(gòu)造3階張量完成的，采用三線性PARAFAC進(jìn)行機(jī)械故障診斷時(shí)，只能對(duì)振動(dòng)信號(hào)的部分信號(hào)建模，存在信息不完整的問題。在平行因子的發(fā)展過程中，四線性PARAFAC的提出具有重要的意義。四線性PARAFAC具有三線性平行因子所不具有的優(yōu)勢(shì)，且采用四線性平行因子盲分離算法可有效地對(duì)機(jī)械多故障進(jìn)行源數(shù)估計(jì)，更好地分解信號(hào)，提取故障特征[6]。然而，在研究機(jī)械多故障欠定盲分離時(shí)，通常情況下對(duì)四線性PARAFAC模型采用QALS迭代擬合，而在欠定條件下，收斂過程中會(huì)出現(xiàn)陷入局部收斂的問題，得到的估計(jì)源信號(hào)往往不是很理想。集成經(jīng)驗(yàn)?zāi)B(tài)分解(EnsembleEmpiricalModeDecomposition，EEMD)是針對(duì)經(jīng)驗(yàn)?zāi)B(tài)分解的不足而提出的1種噪聲輔助數(shù)據(jù)分析方法，該方法充分利用白噪聲均值為零的特點(diǎn)，使得由于干擾導(dǎo)致的任何尺度的信號(hào)所出現(xiàn)間斷，都能得到白噪聲的補(bǔ)充，使間斷信號(hào)在相應(yīng)尺度形成連續(xù)信號(hào)，該方法很好抑制了由于EMD（EmpiricalModeDecomposition，EMD）分離所導(dǎo)致各分量之間出現(xiàn)混頻現(xiàn)象[7]。文獻(xiàn)[8]結(jié)合EEMD和具有自適應(yīng)增長特性的輸出隱藏反饋Elman網(wǎng)絡(luò)對(duì)滾動(dòng)軸承進(jìn)行診斷。文獻(xiàn)[9]將EEMD應(yīng)用于機(jī)械故障中的復(fù)雜信號(hào)處理，提取故障特征。文獻(xiàn)[10]提出1種改進(jìn)的EEMD方法，以消除在EEMD方法中只添加噪聲而導(dǎo)致的局限性，結(jié)果表明所提方法可以抑制機(jī)械信號(hào)模態(tài)混淆并且在故障診斷方面效果明顯。文獻(xiàn)[11]將EEMD和散布熵結(jié)合，構(gòu)建故障信號(hào)的高維特征，在結(jié)合LPP-KNN進(jìn)行故障診斷。文獻(xiàn)[12]將集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)和低相干K-SVD相結(jié)合進(jìn)行齒輪故障特征提取。上述研究均表明EEMD在進(jìn)行機(jī)械振動(dòng)信號(hào)分解時(shí)具有獨(dú)特的優(yōu)勢(shì)。在此，結(jié)合四線性平行因子盲分離算法和集合經(jīng)驗(yàn)?zāi)B(tài)分解算法(EEMD)的的獨(dú)特優(yōu)勢(shì)，提出了一種基于EEMD分解的四線性PARAFAC欠定盲源分離方法。該方法首先利用EEMD分解傳感器所采集的觀測(cè)信號(hào)，接著從分解得到的IMF分量中選出恰當(dāng)相關(guān)IMF分量信號(hào)。根據(jù)選出的相關(guān)IMF分量與原觀測(cè)信號(hào)重新構(gòu)造信號(hào)，使得經(jīng)四線性分解處理的重組信號(hào)數(shù)大于或等于源信號(hào)數(shù)，這樣盲分離中的欠定就被轉(zhuǎn)換為正定或超定。最后利用四線性交替最小二乘法迭代(QuadrilinearAlternatingLeastSquares，QALS)進(jìn)行擬合迭代，得到混合矩陣估計(jì)，用最短路徑法得到源信號(hào)估計(jì)。然后，通過仿真和實(shí)驗(yàn)研究以驗(yàn)證該方法的有效性。1EEMD-四線性平行因子盲分離方法信號(hào)線性混疊模式可以表示為如下形式：式中：S(t)=[s1(t),…,sM(t)]T為源信號(hào)矩陣，每個(gè)源信號(hào)包含L個(gè)數(shù)據(jù)點(diǎn)，A(t)是N×M維的混合矩陣(N<p>在觀測(cè)信號(hào)X(t)中加入白噪聲信號(hào)，利用EEMD方法將含噪的觀測(cè)信號(hào)分解成多個(gè)IMF信號(hào)分量。具體步驟如下：第一步：首先在處理觀測(cè)信號(hào)X(t)=[x1(t),…,xN(t)]T前設(shè)定平均處理次數(shù)M，其中i=1,2…,M；第二步：接著在觀測(cè)信號(hào)X(t)添加白噪聲ni(t)，將加入白噪聲的信號(hào)作為一個(gè)整體，然后與原觀測(cè)信號(hào)組成新的觀測(cè)信號(hào)Xi(t)=X(t)+ni(t)；第三步：用EMD分解Xi(t)得到IMF，得到子分量ci,n(t)；其中：ri,n(t)是得到的殘余分量；第四步：重復(fù)前面3個(gè)步驟，通過分解得到新的一系列的IMF；第五步：對(duì)得到的IMF分量進(jìn)行M次平均獲得cn(t)；其中：i=1,2…,M；n=1,2…M。第六步：利用相似系數(shù)公式求通過分解得到的IMF分量相關(guān)系數(shù)，然后選取相關(guān)性較大的有效IMF。選取相關(guān)性較強(qiáng)的IMF信號(hào)分量，并將其與原來的傳感器數(shù)據(jù)組合，構(gòu)建新的觀測(cè)信號(hào)X(t)。這樣，實(shí)現(xiàn)了將欠定條件下的盲源分離轉(zhuǎn)化為正定或超定條件下盲源分離。然后將經(jīng)中心化處理后的觀測(cè)信號(hào)進(jìn)行分段處理，其被均勻劃分成J個(gè)不會(huì)重疊的數(shù)據(jù)段，經(jīng)計(jì)算可知每段的數(shù)據(jù)點(diǎn)數(shù)為L/J個(gè)，記其為K，即K=L/J，具體段數(shù)可以任意選擇，若最后一段的點(diǎn)數(shù)不夠，可以補(bǔ)零。每一個(gè)數(shù)據(jù)段對(duì)應(yīng)的采集時(shí)間為h，則其共被分為t/h個(gè)時(shí)間段，記其為I，此時(shí)：信號(hào)塊的時(shí)滯協(xié)方差矩陣為：其中：k=1,2,…,K，序號(hào)tk表示t內(nèi)第k個(gè)數(shù)據(jù)點(diǎn)。傳感器在h段采集到的各源信號(hào)數(shù)據(jù)時(shí)滯協(xié)方差Rn可以通過計(jì)算得到，有：式(6)中：tk表示t內(nèi)第k個(gè)源信號(hào)。將塊時(shí)滯協(xié)方差矩陣轉(zhuǎn)換成為對(duì)角陣，并且使該矩陣等于矩陣B，將源時(shí)滯協(xié)方差矩陣轉(zhuǎn)換成為對(duì)角矩陣并且使該矩陣等于矩陣C，則式(5)、式(6)可以分別改寫為式(7)和式(8)。將源時(shí)滯協(xié)方差矩陣以及塊時(shí)滯協(xié)方差矩陣疊加，疊加后所得的3階矩陣計(jì)為RJ,K其中：n=1,2,…,N,k=1,2,…,K。相鄰對(duì)角陣相乘可以換位，則式(9)可表示為：將式中N×K個(gè)切片累積成N×J×K×I維的R（四維數(shù)據(jù)集）。其中：n=1,2,…,N,k=1,2…,K，j=1,2,…,J，i=1,2,…I，式(11)即為所構(gòu)造的EEMD-四線性盲分離模型。如上實(shí)現(xiàn)了集合經(jīng)驗(yàn)?zāi)B(tài)四線性平行因子盲源分離算法的構(gòu)建，其分解過程中采用QALS得到混合矩陣估計(jì)A′。分解算法的主要思想為讓上一次估計(jì)更新的矩陣參與下一次的矩陣求解，同時(shí)使根據(jù)最終估計(jì)得到的4個(gè)載荷矩陣與被分解矩陣的差值盡量最小，迭代過程步驟分為4步。第二步：將A重新代入式(12)，接著最小化式(12)得到估計(jì)的，令其等于B，更新矩陣B。第三步：將求得的矩陣B重新代入式(12)，接著最小化式(12)得到矩陣估計(jì)，令其等于C，更新矩陣C。第四步：將按照以上循環(huán)步驟得到的載荷矩陣代入式(17)，若滿足式(17)，則退出，反之繼續(xù)。一般認(rèn)為ε為1×10-6，??3是殘余誤差的平方和，表示含噪的切片，‖?‖F(xiàn)為Frobenius范數(shù)，i=1,2,…,I，j=1,2,…,J，上標(biāo)+為Moore-Penrose偽逆，和表示A、B和C的上一次估計(jì)，m為迭代次數(shù)。當(dāng)式（17）所示的停止準(zhǔn)則得到滿足即表示矩陣更新結(jié)束算法收斂，得到最終的混合矩陣估計(jì)A′，令A(yù)′=A。令ai為A第i個(gè)列向量，si(t)表示源信號(hào)的第i個(gè)源向量。線性混疊模型X(t)=AS(t)可以表示為：其中：X(t)是觀測(cè)信號(hào)，它可以表示成a1s1(t)、a2s2(t)、…、ansn(t)的累加。ai表示矢量方向，si(t)表示向量長度。對(duì)于欠定混合的情況，采用最短路徑法[9]恢復(fù)源信號(hào)。推斷源的問題可以用式（19）表示：由式(18)以及式(19)即可得到源信號(hào)估計(jì)。其主要思想是求每個(gè)時(shí)刻最優(yōu)解，然后得到恢復(fù)的源信號(hào)。2仿真研究構(gòu)造的兩個(gè)仿真信號(hào)s1(t)、s2(t)如式(20)所示，式中：f1=20Hz，f2=30Hz，f3=50Hz，f4=60Hz。選取采樣頻率fs=5000Hz，采樣點(diǎn)數(shù)N=20000。繪制兩個(gè)源信號(hào)s1(t)、s2(t)的時(shí)域波形如圖1所示。相應(yīng)的頻譜如圖2所示。圖1源信號(hào)時(shí)域波形圖圖2源信號(hào)頻譜圖用1×2維的隨機(jī)混合矩陣A，將兩個(gè)源信號(hào)s1(t)、s2(t)進(jìn)行隨機(jī)混合，混合后得到1路的混合信號(hào)如圖3所示。已知源信號(hào)的個(gè)數(shù)為2，混合后觀測(cè)信號(hào)的個(gè)數(shù)為1，因此，該盲源分離問題為欠定求解問題。圖3混疊信號(hào)時(shí)域波形圖利用四線性平行因子盲分離算法將該混合信號(hào)分離，得到恢復(fù)源信號(hào)如圖4所示。相應(yīng)的頻譜如圖5所示。圖4恢復(fù)信號(hào)時(shí)域波形圖圖5恢復(fù)信號(hào)的頻譜對(duì)比圖4和圖1可知，采用四平行因子盲分離算法得到的恢復(fù)信號(hào)與源信號(hào)相似度不是特別高，分離效果不理想。由圖5可知，分離得到的恢復(fù)信號(hào)的特征頻率分別為20Hz、30Hz、50Hz、60Hz，與源信號(hào)特征頻率不一致，分離得到的頻譜圖效果不是很理想，因此，在欠定條件下直接采用四線性平行因子盲源分離算法分解效果不理想，需對(duì)其進(jìn)一步改進(jìn)。為了改進(jìn)四線性平行因子盲源分離算法的這種不足，用集合經(jīng)驗(yàn)?zāi)B(tài)平行因子盲源分離算法分解混合信號(hào)。該算法首先用EEMD算法分解混合信號(hào)，通過分離得到14個(gè)子分量，然后選擇相關(guān)子分量。選擇相關(guān)分量的具體做法是利用相關(guān)系數(shù)求各子分量與觀測(cè)信號(hào)的相關(guān)度，一般認(rèn)為信號(hào)與觀測(cè)信號(hào)越相似，則相關(guān)度的值越接近為1。根據(jù)相關(guān)度的值定義其在0到0.3之間為微相關(guān)，0.3到0.5之間為實(shí)相關(guān)，0.5到0.8之間為顯著相關(guān)，0.8到1之間為高度相關(guān)。根據(jù)EEMD算法分解得到的14個(gè)子分量的相關(guān)系數(shù)如表1所示。由表1可知，子分量IMF1、IMF6、IMF7、IMF8的相關(guān)系數(shù)值相對(duì)較高，其他分量均為不相關(guān)，所以選擇IMF1、IMF6、IMF7、IMF8為相關(guān)分量并將其與原觀測(cè)信號(hào)組合得到新觀測(cè)信號(hào)，這樣就將欠定轉(zhuǎn)化為超定。繪制相關(guān)分量IMF1、IMF6、IMF7、IMF8的時(shí)域波形圖和頻譜圖，得到的圖形分別如圖6以及圖7所示。圖6根據(jù)EEMD分解得到相關(guān)子分量時(shí)域波形圖圖7根據(jù)EEMD分解得到相關(guān)子分量的幅值譜圖表1各子分量相關(guān)系數(shù)將通過分解得到相關(guān)子分量和混疊信號(hào)組合，得到新觀測(cè)信號(hào)，用核一致算法估計(jì)振源數(shù)，如圖9所示。圖9組件數(shù)由圖9可知，振源數(shù)為2。將得到的觀測(cè)信號(hào)用四線性平行因子盲分離算法處理，用QALS算法擬合迭代，當(dāng)達(dá)到迭代終止條件時(shí)得到4個(gè)載荷矩陣A、B、C、D。由圖8可知，4個(gè)載荷矩陣中，載荷矩陣B、D分解形式相同，表明建立的模型符合式(11)所示的四線性模型，并未偏離四線性模型。圖8載荷矩陣用集合經(jīng)驗(yàn)?zāi)B(tài)四線性平行因子盲源分離算法處理通過組合得到的觀測(cè)信號(hào)，依據(jù)模型求得混合矩陣估計(jì)，以及根據(jù)式(19)求得估計(jì)源信號(hào)，繪制恢復(fù)信號(hào)的時(shí)域波形圖如10圖所示。相應(yīng)的恢復(fù)信號(hào)頻譜如圖11所示。由圖11可知，圖11(a)中的特征頻率分別為50Hz、60Hz，圖11(b)中的特征頻率分別為20Hz、30Hz，除排列順序有所不同外，其分別與兩源信號(hào)的特征頻率相一致，這種排列順序改變并不影響分離結(jié)果的正確性。對(duì)比圖10與源信號(hào)時(shí)域波形，恢復(fù)信號(hào)的時(shí)域波形圖和源信號(hào)時(shí)域波形圖非常相似。由此可知，提出的方法能夠有效分離混合信號(hào)，分離效果非常令人滿意。對(duì)比提出的方法和傳統(tǒng)的四線性平行因子盲源分離算法，發(fā)現(xiàn)欠定條件下采用四線性平行因子盲源分離算法并不能得到很好的分離效果，而采用集合經(jīng)驗(yàn)?zāi)B(tài)四線性平行因子盲源分離算法可以很好分離混合信號(hào)并且分離效果比四線性平行因子盲源分離算法好很多，證明了該算法的有效性。圖10恢復(fù)信號(hào)時(shí)域波形圖11恢復(fù)信號(hào)的頻譜圖3實(shí)驗(yàn)研究為了驗(yàn)證所提算法的有效性，將集合經(jīng)驗(yàn)?zāi)B(tài)四線性平行因子盲分離算法應(yīng)用于滾動(dòng)軸承復(fù)合故障盲源分離中。實(shí)驗(yàn)所需要的數(shù)據(jù)由滾動(dòng)軸承診斷試驗(yàn)臺(tái)采集，試驗(yàn)臺(tái)由旋轉(zhuǎn)電機(jī)、滾子軸承和加速度傳感器組成。軸承的代號(hào)是N205EM。電機(jī)轉(zhuǎn)速約為1300r/min，轉(zhuǎn)動(dòng)頻率約為21.7Hz，采樣頻率fs=100kHz。數(shù)據(jù)中包含一組內(nèi)滾道和滾子的復(fù)合故障信號(hào)。使用線切割機(jī)技術(shù)人為制造缺陷以產(chǎn)生故障源，內(nèi)圈故障大小為0.3mm×0.15mm，滾子故障大小為0.3mm×0.15mm，利用軸承故障頻率機(jī)理公式可計(jì)算得到滾子的故障特征頻率fb=102.3Hz，內(nèi)圈故障特征頻率fi=145.7Hz。所用信號(hào)的采樣點(diǎn)數(shù)為1001000，繪制采集的混合觀測(cè)信號(hào)時(shí)域波形圖，如圖12所示。圖12觀測(cè)信號(hào)時(shí)域波形圖將觀測(cè)信號(hào)用EEMD算法分解，得到16個(gè)IMF分量。利用相關(guān)系數(shù)法選擇相關(guān)性較大的子分量，各子分量的相關(guān)系數(shù)如表2所示。由表2可知，子分量IMF1、IMF2、IMF3、IMF4、IMF5的相關(guān)系數(shù)均大于0.3，為有效子分量，其他子分量信號(hào)的相似系數(shù)值均在0.3以下，所以選擇IMF1、IMF2、IMF3、IMF4、IMF5為相關(guān)分量。將選擇的相關(guān)子分量與原觀測(cè)信號(hào)組合得到新的觀測(cè)信號(hào)。繪制IMF1、IMF2、IMF3、IMF4、MF5的時(shí)域圖，得到的圖形如圖13所示。表2各子分量相關(guān)系數(shù)圖13采用EEMD分解觀測(cè)信號(hào)得到的有效子分量時(shí)域波形圖將分解得到的這5個(gè)相關(guān)的子分量信號(hào)和混疊信號(hào)重新構(gòu)造成新的觀測(cè)信號(hào)，計(jì)算觀測(cè)信號(hào)段的時(shí)滯協(xié)方差矩陣，并將所有數(shù)據(jù)段的時(shí)滯協(xié)方差疊加成4階張量形式，得到四線性模型，用QALS擬合，當(dāng)滿足式(17)時(shí)迭代結(jié)束，求得估計(jì)矩陣。利用所求的估計(jì)矩陣和最短路徑法得到恢復(fù)信號(hào)，繪制恢復(fù)信號(hào)的時(shí)域圖和頻譜如圖13、圖14所示。圖14恢復(fù)信號(hào)時(shí)域波形