基于經(jīng)驗(yàn)?zāi)B(tài)與奇異值分解的振動(dòng)源數(shù)估計(jì)方法

劉維新，葉超

(中國工程物理研究院，機(jī)械制造工藝研究所，四川綿陽 621900)

0 前言

轉(zhuǎn)輪模塊是某大型光學(xué)實(shí)驗(yàn)裝置空間濾波器的核心部件，對(duì)裝置實(shí)現(xiàn)預(yù)定功能具有重要作用，其可靠性問題備受關(guān)注。轉(zhuǎn)輪模塊由電機(jī)、減速器、聯(lián)軸器、圓光柵、轉(zhuǎn)輪等部件組成，運(yùn)行在真空環(huán)境中，對(duì)可靠性的要求極高。

本文作者在某大型光學(xué)實(shí)驗(yàn)裝置部件可靠性實(shí)驗(yàn)時(shí)發(fā)現(xiàn)，由于機(jī)械結(jié)構(gòu)緊湊、低速重載、工況轉(zhuǎn)換頻繁、沖擊頻率高，轉(zhuǎn)輪模塊中的行星減速器是其主要故障源，是影響裝置可靠性的關(guān)鍵因素。

目前，行星減速器故障診斷領(lǐng)域已有較多學(xué)者取得良好的研究成果，但大多針對(duì)均勻載荷和伺服電機(jī)驅(qū)動(dòng)下的減速器，對(duì)非均勻載荷和步進(jìn)電機(jī)驅(qū)動(dòng)下的行星減速器的研究較少。由于振動(dòng)傳感器只能布置在行星減速器箱體外，行星減速箱內(nèi)一系列激勵(lì)源產(chǎn)生的振動(dòng)響應(yīng)與干擾、噪聲等耦合在一起，振動(dòng)傳感器測(cè)量得到的是統(tǒng)計(jì)意義上的耦合源信號(hào)，很難直接測(cè)量得到單一振動(dòng)激勵(lì)源。非均勻載荷和步進(jìn)電機(jī)推力輸出特性造成的干擾與轉(zhuǎn)輪模塊行星減速器故障特征耦合在一起難以剝離，實(shí)現(xiàn)振動(dòng)監(jiān)測(cè)信號(hào)內(nèi)激勵(lì)源(振源)與外激勵(lì)源(干擾)的盲源分離，是實(shí)現(xiàn)行星減速器故障檢測(cè)與辨識(shí)的重要前提。

傳統(tǒng)盲源分離中，為計(jì)算簡單，一般假設(shè)源的數(shù)量等于觀測(cè)傳感器數(shù)，則混合/解耦矩陣是方陣，容易估計(jì)。然而，在文中研究對(duì)象轉(zhuǎn)輪模塊行星減速器中，這種假設(shè)無法成立，因?yàn)槭馨惭b位置和空間限制，傳感器數(shù)量只能是有限的，而源的數(shù)量很容易超過傳感器數(shù)量。

為實(shí)現(xiàn)轉(zhuǎn)輪模塊行星減速器振動(dòng)監(jiān)測(cè)信號(hào)盲源分離，需準(zhǔn)確估計(jì)內(nèi)部激勵(lì)振源數(shù)。針對(duì)振動(dòng)傳感器數(shù)小于本底振源數(shù)的源數(shù)估計(jì)問題，本文作者提出一種基于經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)的虛擬通道擴(kuò)展方法。通過經(jīng)驗(yàn)?zāi)B(tài)分解得到的固有模態(tài)函數(shù)構(gòu)建振動(dòng)信號(hào)觀測(cè)矩陣，以擴(kuò)充振動(dòng)傳感器觀測(cè)通道數(shù)量。針對(duì)奇異值分解(Singular Value Decomposition，SVD)特征矩陣中噪聲和數(shù)據(jù)觀測(cè)誤差以及不確定性導(dǎo)致的源數(shù)估計(jì)不準(zhǔn)確問題，提出一種基于類內(nèi)散度與類間距離比值優(yōu)化的聚類分析方法。通過對(duì)奇異值分解后特征值矩陣中對(duì)角線特征值的聚類分析，獲得盲源數(shù)估計(jì)結(jié)果。通過仿真實(shí)驗(yàn)和模擬實(shí)驗(yàn)，驗(yàn)證文中所提振源數(shù)估計(jì)方法的有效性。

1 基于經(jīng)驗(yàn)?zāi)B(tài)分解的虛擬通道擴(kuò)展方法

經(jīng)驗(yàn)?zāi)B(tài)分解能夠把復(fù)雜信號(hào)分解為一系列內(nèi)稟函數(shù)(Intrinsic Mode Function,IMF)之和，每個(gè)IMF分量代表原信號(hào)的一個(gè)固有信號(hào)模態(tài)，且包含了不同的特征時(shí)間尺度，讓原信號(hào)特征以不同的分辨率表現(xiàn)出來。由于每個(gè)IMF所包含的頻率成分不僅與采樣頻率有關(guān)，還隨著信號(hào)的變化而變化，所以EMD方法是一種自適應(yīng)的信號(hào)處理方法，適用于時(shí)變、非平穩(wěn)信號(hào)分解。

本文作者分析轉(zhuǎn)輪模塊振動(dòng)監(jiān)測(cè)系統(tǒng)的特點(diǎn)和振動(dòng)機(jī)制，分別選取行星減速器2個(gè)振動(dòng)傳感器監(jiān)測(cè)信號(hào)的前5個(gè)IMF，將原有的2個(gè)振動(dòng)觀測(cè)通道擴(kuò)展為10個(gè)振動(dòng)觀測(cè)通道。

虛擬通道擴(kuò)展方法具體實(shí)現(xiàn)過程如下：

(1)采集轉(zhuǎn)輪模塊行星減速器箱體表面振動(dòng)信號(hào)()；

(2)辨識(shí)信號(hào)()上所有局部極大值和極小值點(diǎn)，用三次樣條曲線分別連接所有局部極大值點(diǎn)和極小值點(diǎn)，形成上下包絡(luò)線；

(3)計(jì)算上下包絡(luò)線的均值，然后計(jì)算=()-，并按條件1、條件2判斷是否滿足IMF的條件，若滿足則為振動(dòng)信號(hào)()的1個(gè)IMF分量；

條件1：具有相同數(shù)量的極值點(diǎn)和過零點(diǎn)，或最多相差一個(gè)；

條件2：任意時(shí)刻極大值和極小值的包絡(luò)線關(guān)于時(shí)間軸局部對(duì)稱，平均值為0；

(4)若不滿足IMF的條件，則將作為原始數(shù)據(jù)，重復(fù)步驟(1)—(3)次，直到1滿足IMF的條件，記為=1，為振動(dòng)信號(hào)()的1個(gè)IMF分量；

(5)將IMF分量從振動(dòng)信號(hào)()中分離出來，得到=()-，將作為原始數(shù)據(jù)重復(fù)步驟(1)—(4)，重復(fù)循環(huán)次，得到個(gè)IMF分量，即得到=-;…;=-1-；

(7)選取1路振動(dòng)傳感器經(jīng)驗(yàn)?zāi)B(tài)分解后的前5個(gè)IMF，構(gòu)建該路傳感器新的5路虛擬通道觀測(cè)矩陣；

(8)循環(huán)結(jié)束。

在經(jīng)驗(yàn)?zāi)B(tài)分解過程中，原振動(dòng)信號(hào)()相鄰2個(gè)極大值點(diǎn)或極小值點(diǎn)之間的時(shí)間跨度被稱為時(shí)間特征尺度，反映的是振動(dòng)信號(hào)不同模態(tài)的局部振動(dòng)特性。經(jīng)驗(yàn)?zāi)B(tài)分解過程就是按照特征時(shí)間尺度從小到大的順序依次分離各IMF，包含原振動(dòng)信號(hào)從高到低不同頻率段的振動(dòng)成分，也可以把經(jīng)驗(yàn)?zāi)B(tài)分解方法理解為是一組濾波器。

2 基于奇異值分解的盲源數(shù)估計(jì)方法

轉(zhuǎn)輪模塊行星減速器可視為一個(gè)多輸入多輸出系統(tǒng)(MIMO)。行星減速器中個(gè)振源(=1,…,)作為MIMO系統(tǒng)的輸入，經(jīng)過從源到觀測(cè)的頻率響應(yīng)函數(shù)(FRF)矩陣={}(=1,…,;=1,…,)的作用，在振動(dòng)傳感器處混合形成系統(tǒng)輸出(=1,…,)。各源分量之間可能相關(guān)也可能不相關(guān)。

由譜分析理論，個(gè)振源的自譜和互譜可構(gòu)成一個(gè)×的譜矩陣，譜矩陣的線性無關(guān)行(列)的數(shù)目，即矩陣的秩，即為振動(dòng)源中的不相關(guān)源數(shù)。譜矩陣通常不可能直接觀測(cè)，它與傳感器觀測(cè)矩陣的關(guān)系為=。由線性代數(shù)相關(guān)定理可知，如果的秩rank()=,則rank()≥rank()，工程應(yīng)用中可近似認(rèn)為rank()=rank()，其中：和分別對(duì)應(yīng)和的頻率響應(yīng)矩陣。因此,工程中可以利用傳感器觀測(cè)信號(hào)的譜矩陣的秩估計(jì)不相關(guān)振源數(shù)。

奇異值分解是一種對(duì)矩陣的正交化分解方法，是一種數(shù)值穩(wěn)定的矩陣求秩方法。一個(gè)×矩陣的SVD表達(dá)為=，其中：是一個(gè)由矩陣的正交特征向量組成的×對(duì)角矩陣；是一個(gè)由矩陣的奇異值組成的×對(duì)角矩陣;是一個(gè)由矩陣的正交特征向量組成的×矩陣。

理論上，譜矩陣經(jīng)過奇異值分解后，非零的奇異值數(shù)目即為待求的秩。然而，對(duì)于實(shí)際采集的傳感器觀測(cè)數(shù)據(jù)，譜矩陣的奇異值大多不會(huì)出現(xiàn)為零的情況，原因是工程應(yīng)用中采集的傳感器觀測(cè)信號(hào)，由于有限數(shù)據(jù)采樣精度及奇異值分解計(jì)算的誤差累積效應(yīng)，觀測(cè)信號(hào)譜矩陣的獲得涉及對(duì)自譜和互譜密度的估計(jì)，導(dǎo)致奇異值分解的非零誤差。一般來說，由數(shù)模轉(zhuǎn)換引入的誤差占比較大。

通常，工程應(yīng)用中采集到的傳感器觀測(cè)信號(hào)奇異值分解后得到的奇異值呈現(xiàn)大小不等的兩端分布，其中較大的奇異值與待估計(jì)的振源有對(duì)應(yīng)關(guān)系。文獻(xiàn)[14]中建議采用閾值方法對(duì)不同奇異值進(jìn)行分類，即設(shè)置一固定閾值，大于該閾值的值認(rèn)為是有效的奇異值，對(duì)應(yīng)某一振源；小于該閾值的值認(rèn)為是無效的奇異值，對(duì)應(yīng)誤差或噪聲。然而，實(shí)際應(yīng)用中，閾值的設(shè)定必須綜合考慮信號(hào)處理與參數(shù)估計(jì)誤差對(duì)譜矩陣奇異值分解的影響，而影響譜估計(jì)精度的參數(shù)和因素較多且相互間并無明確的函數(shù)關(guān)系，導(dǎo)致閾值的合理設(shè)定較為困難。

針對(duì)上述奇異值分解特征矩陣中噪聲和數(shù)據(jù)觀測(cè)誤差以及不確定性導(dǎo)致的源數(shù)估計(jì)不準(zhǔn)確問題，提出一種基于類內(nèi)散度與類間距離比值優(yōu)化的聚類分析方法。通過對(duì)奇異值分解后特征值矩陣中對(duì)角線特征值的聚類分析，獲得振動(dòng)源數(shù)估計(jì)結(jié)果。所提方法的實(shí)施流程如圖1所示。

圖1 所提振源數(shù)估計(jì)方法流程

所提振源數(shù)估計(jì)方法的具體實(shí)現(xiàn)過程如下：

(1)采集轉(zhuǎn)輪模塊減速箱箱體表面90°夾角布置的2個(gè)振動(dòng)傳感器的振動(dòng)信號(hào)()、()；

(2)采用第1節(jié)所述的基于經(jīng)驗(yàn)?zāi)B(tài)分解的虛擬通道擴(kuò)展方法，得到振動(dòng)信號(hào)的固有模態(tài)；采用每個(gè)單獨(dú)振動(dòng)信號(hào)的前5個(gè)固有模態(tài)擴(kuò)展振動(dòng)觀測(cè)信號(hào)矩陣；

(3)計(jì)算振動(dòng)觀測(cè)信號(hào)矩陣中各固有模態(tài)的自/互相關(guān)系數(shù)，構(gòu)建自/互相關(guān)系數(shù)矩陣:

其中：IMF，IMF，代表原始振動(dòng)信號(hào)的第個(gè)固有模態(tài)函數(shù)與第個(gè)固有模態(tài)函數(shù)的互相關(guān)系數(shù)，、取值范圍均為[1,2,3,4,5]。

(4)根據(jù)所構(gòu)建的自/互相關(guān)系數(shù)矩陣，分別計(jì)算矩陣和矩陣的特征值和特征向量;計(jì)算得到的矩陣的特征向量為矩陣奇異值分解的左奇異向量，計(jì)算得到的矩陣的特征向量為矩陣奇異值分解的右奇異向量，特征值組成的對(duì)角線矩陣為矩陣的奇異值矩陣，即=，rank()=10；

(5)以對(duì)角矩陣的10個(gè)對(duì)角線特征元素作為均值聚類的對(duì)象，基于類內(nèi)散度與類間距離的比值，即用DB指標(biāo)來評(píng)價(jià)聚類效果，確定最優(yōu)聚類數(shù)；

(6)最優(yōu)聚類數(shù)減1(噪聲干擾形成的奇異值矩陣特征元素聚集為一類)即為振動(dòng)信號(hào)盲源數(shù)估計(jì)結(jié)果。

3 仿真實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證文中所提振動(dòng)信號(hào)盲源數(shù)估計(jì)方法的有效性，采用仿真實(shí)驗(yàn)驗(yàn)證盲源數(shù)估計(jì)結(jié)果是否正確。仿真信號(hào)的振源已知，便于評(píng)價(jià)算法的有效性。構(gòu)造仿真信號(hào)如下：

=sin(2π)+sin(2π)+sin(2π)+0.1

(1)

式中：=250；=100；=50;為隨機(jī)干擾，范圍為[0,1)。由式(1)可看出該仿真信號(hào)具有3個(gè)振源，其時(shí)域波形如圖2所示。

圖2 仿真信號(hào)時(shí)域波形

原仿真信號(hào)觀測(cè)傳感器數(shù)量為1，由盲源分離原理可知，最多只能估計(jì)出一個(gè)振動(dòng)源。為準(zhǔn)確估計(jì)振動(dòng)源，采用文中第1節(jié)所述的基于經(jīng)驗(yàn)?zāi)B(tài)分解的虛擬通道擴(kuò)展方法，擴(kuò)展傳感器觀測(cè)矩陣，得到仿真信號(hào)固有模態(tài)分解結(jié)果如圖3所示。

圖3中IMF8為殘余分量，IMF1～I(xiàn)MF7為仿真信號(hào)固有模態(tài)。采用這7個(gè)固有模態(tài)函數(shù)擴(kuò)展原觀測(cè)傳感器通道，可擴(kuò)展為7通道觀測(cè)信號(hào)。采用文中第2節(jié)所述的基于奇異值分解的源數(shù)估計(jì)方法，構(gòu)建自/互相關(guān)系數(shù)矩陣，然后對(duì)矩陣進(jìn)行奇異值分解，得到奇異值矩陣。對(duì)奇異值矩陣的7個(gè)對(duì)角線特征值進(jìn)行聚類分析，基于類內(nèi)散度與類間距離的比值來評(píng)價(jià)聚類效果，結(jié)果如圖4所示。

圖3 仿真信號(hào)經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果

由DB指標(biāo)的定義可知，DB指標(biāo)越小，聚類效果越好，但DB指標(biāo)不應(yīng)為0。因此，由圖4可判斷仿真信號(hào)觀測(cè)矩陣7個(gè)奇異值特征最優(yōu)分類數(shù)為4，分類結(jié)果如圖5所示。

圖4 仿真信號(hào)DB指標(biāo)聚類效果 圖5 仿真奇異值矩陣特征值聚類結(jié)果

由圖5可以看出：仿真信號(hào)奇異值矩陣對(duì)角線特征值聚類為4類；有3個(gè)聚類結(jié)果分別對(duì)應(yīng)不同的振動(dòng)源，第4個(gè)聚類結(jié)果數(shù)值接近為0，對(duì)應(yīng)噪聲和數(shù)據(jù)觀測(cè)誤差等因素；文中方法振源辨識(shí)結(jié)果與仿真信號(hào)真實(shí)情況一致。結(jié)果表明：對(duì)于仿真信號(hào)，文中所提方法可有效實(shí)現(xiàn)振動(dòng)源數(shù)盲估計(jì)。

4 故障模擬實(shí)驗(yàn)對(duì)比驗(yàn)證

為進(jìn)一步驗(yàn)證文中所提振源數(shù)估計(jì)方法的有效性，采用某大型科學(xué)裝置上的轉(zhuǎn)輪模塊原型樣機(jī)開展實(shí)驗(yàn)驗(yàn)證。轉(zhuǎn)輪模塊結(jié)構(gòu)與構(gòu)建的原型樣機(jī)狀態(tài)監(jiān)測(cè)系統(tǒng)分別如圖6、圖7所示。

圖6 轉(zhuǎn)輪模塊結(jié)構(gòu)

圖7 轉(zhuǎn)輪模塊狀態(tài)監(jiān)測(cè)系統(tǒng)

該監(jiān)測(cè)系統(tǒng)由轉(zhuǎn)輪模塊、8個(gè)振動(dòng)傳感器、1臺(tái)DASP-V11振動(dòng)信號(hào)采集分析儀組成。8個(gè)振動(dòng)傳感器分別布置在前支撐軸承、后支撐軸承、行星減速器、步進(jìn)電機(jī)上，其中，本文作者重點(diǎn)研究的行星減速器上呈90°夾角布置有2個(gè)振動(dòng)傳感器。為獲取減速器故障狀態(tài)下的振動(dòng)信號(hào)，通過在減速器太陽輪、行星輪表面進(jìn)行加工的方式模擬減速器齒輪故障，并獲取對(duì)應(yīng)工況下的減速器振動(dòng)監(jiān)測(cè)信號(hào)。實(shí)驗(yàn)過程：轉(zhuǎn)輪旋轉(zhuǎn)360°后反向往復(fù)運(yùn)動(dòng)，采樣頻率1 000 Hz。齒輪故障模擬效果如圖8所示。

圖8 一級(jí)行星輪崩齒故障模擬效果

以崩齒故障為例，驗(yàn)證文中方法的有效性。一級(jí)行星齒輪崩齒故障振動(dòng)信號(hào)監(jiān)測(cè)結(jié)果如圖9所示。

圖9 崩齒故障行星減速器振動(dòng)監(jiān)測(cè)信號(hào)

行星減速器具有2個(gè)振動(dòng)觀測(cè)傳感器，根據(jù)盲源分離原理分析，最多只能估計(jì)出2個(gè)振動(dòng)源。為準(zhǔn)確估計(jì)振動(dòng)源，采用第1節(jié)所述的基于經(jīng)驗(yàn)?zāi)B(tài)分解的虛擬通道擴(kuò)展方法，擴(kuò)展傳感器觀測(cè)矩陣。考慮到算法的復(fù)雜度和運(yùn)算速度，選取每個(gè)振動(dòng)傳感器監(jiān)測(cè)信號(hào)的前5個(gè)固有模態(tài)函數(shù)擴(kuò)展觀測(cè)信號(hào)，結(jié)果如圖10所示。

圖10 崩齒故障振動(dòng)信號(hào)經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果

采用垂向和橫向振動(dòng)信號(hào)的前5個(gè)固有模態(tài)IMF1～I(xiàn)MF5擴(kuò)展原觀測(cè)傳感器通道，可擴(kuò)展為10通道觀測(cè)信號(hào)。采用文中第2節(jié)所述的基于奇異值分解的源數(shù)估計(jì)方法，得到奇異值矩陣。對(duì)奇異值矩陣的10個(gè)對(duì)角線特征值進(jìn)行聚類分析，基于類內(nèi)散度與類間距離的比值評(píng)價(jià)聚類效果，結(jié)果如圖11所示。

由圖11可判斷仿真信號(hào)觀測(cè)矩陣10個(gè)奇異值特征最優(yōu)分類數(shù)為5，分類結(jié)果如圖12所示。

圖11 振動(dòng)信號(hào)DB指標(biāo)聚類效果 圖12 振動(dòng)信號(hào)奇異值矩陣特征值聚類結(jié)果

由圖12可以看出：奇異值矩陣對(duì)角線特征值聚類為5類，其中5個(gè)聚類結(jié)果分別對(duì)應(yīng)不同的振動(dòng)源，第5個(gè)聚類結(jié)果數(shù)值接近0，對(duì)應(yīng)噪聲和數(shù)據(jù)觀測(cè)誤差等因素，即文中方法估計(jì)轉(zhuǎn)輪模塊崩齒故障振動(dòng)源有4個(gè)。而采用傳統(tǒng)的閾值方法處理奇異值矩陣對(duì)角線特征時(shí)，特征閾值如設(shè)為0.4，則小于0.4的對(duì)角線特征值都被歸類于噪聲和觀測(cè)誤差而不被認(rèn)為是振動(dòng)源。傳統(tǒng)方法振源估計(jì)結(jié)果為1，顯著不符合轉(zhuǎn)輪模塊的物理特性，并且該特征閾值只能人為估計(jì)，很難準(zhǔn)確設(shè)定。結(jié)果表明，文中所提振動(dòng)信號(hào)盲源數(shù)估計(jì)方法效果優(yōu)于傳統(tǒng)閾值方法。

5 結(jié)論

本文作者針對(duì)轉(zhuǎn)輪模塊行星減速器振動(dòng)源數(shù)盲估計(jì)問題，提出一種基于經(jīng)驗(yàn)?zāi)B(tài)分解與奇異值分解相結(jié)合的振動(dòng)信號(hào)盲源數(shù)估計(jì)方法。針對(duì)振動(dòng)傳感器數(shù)小于本底振源數(shù)的源數(shù)估計(jì)問題，通過經(jīng)驗(yàn)?zāi)B(tài)分解得到的固有模態(tài)函數(shù)構(gòu)建振動(dòng)信號(hào)觀測(cè)矩陣，以擴(kuò)充振動(dòng)傳感器觀測(cè)通道數(shù)量。針對(duì)奇異值分解特征矩陣中噪聲和數(shù)據(jù)觀測(cè)誤差以及不確定性導(dǎo)致的源數(shù)估計(jì)不準(zhǔn)確問題，通過對(duì)奇異值分解后特征值矩陣中對(duì)角線特征值的聚類分析，獲得盲源數(shù)估計(jì)結(jié)果。通過仿真實(shí)驗(yàn)和故障模擬實(shí)驗(yàn)驗(yàn)證所提方法的有效性。結(jié)果表明：利用傳統(tǒng)直接采用原始振動(dòng)信號(hào)進(jìn)行源數(shù)估計(jì)的方法只能計(jì)算得到2個(gè)振動(dòng)源；傳統(tǒng)的基于閾值的奇異值分解方法難以準(zhǔn)確估計(jì)振動(dòng)源數(shù)；所提的基于虛擬通道擴(kuò)展與奇異值分解相結(jié)合的方法可準(zhǔn)確估計(jì)振動(dòng)源數(shù)，效果優(yōu)于傳統(tǒng)方法。