李 兵, 張培林, 任國全, 劉東升, 米雙山
(1.石家莊軍械工程學院 自行火炮教研室,石家莊 050003;2.石家莊軍械工程學院 導彈機電工程教研室,石家莊 050003)
振動是發動機在運轉過程中產生的必然現象,振動信號包含了發動機運行狀態的豐富信息,發動機工作性能的變化可以通過振動信號表現出來,而其采集和分析不影響發動機的正常運轉。因此,振動信號分析法成為對發動機進行結構、故障分析和狀態監測最為廣泛、也是最行之有效的方法[1]。
傳統的故障診斷技術是針對振動信號的時域或頻域特征來提取特征量進行故障識別的。但對于多部件、多激勵的往復式發動機而言,由于其運動件多而復雜,工作條件惡劣,激勵力眾多且其頻率范圍寬廣,各種激勵經相應的傳遞及耦合均被反映在發動機表面的振動中,使得所測得的表面振動信號是極其復雜的混合信號[2,3],用傳統的理論方法去刻畫這種復雜性具有很大的局限性。
分形幾何為復雜信號提供了一種幾何結構分析方法,它在很多領域已經有了成功的應用。在機械設備故障診斷領域中,人們也開始用分形幾何方法對振動信號進行分析,并取得了一定的成果[4,5]。然而,在振動信號分析中,人們一般只考慮信號的單重分形特征,它只能從整體上反映信號的不規則性,缺乏對局部奇異性的刻畫。多重分形維數能更精細的刻畫信號的局部尺度行為,可以更加全面反映信號的分形特性。描述多重分形的主要方法有多重分形譜和廣義分形維數兩套參數體系,研究證明這兩套參數體系之間存在著勒讓德變換的關系,即已知一套參數即可計算出另一套參數,本文主要研究基于Renyi熵的廣義分形維數的計算方法。
計算廣義分形維數最常用的方法是覆蓋法即盒計數法(Box-counting method),這種方法有其不可避免的缺點,這主要是因為盒計數法采用了規則劃分網格的方法,因此對分形維數的估計在某些情況下非常不穩定[6]。
針對盒計數法計算廣義分形維數存在的缺點,本文提出了基于數學形態學操作的廣義分形維數的計算方法,并對發動機正常、失火和氣門間隙過大故障信號進行了分析,結果表明,基于數學形態學的廣義分形維數可以非常有效的表征不同狀態下發動機故障信號的特性。
自從曼德爾布羅特在70年代提出分形概念以來,分形在物理、天文、地理、數學、計算機等科學領域取得了廣泛的應用,并且取得了大量富有新意的成果。但是隨著理論研究和應用的深入,研究者們越來越清楚地意識到,對于大多數客觀存在的分形物體而言,僅用一個分形維數并不能完全刻畫其結構。多重分形是定義在分形上的,由多個標量指數的奇異測度(即不存在測度密度的測度)所組成的集合。它刻畫的是分形測度在支集上的分布情況,即用一個譜函數來描述分形不同層次的特征。這就是從形體的部分(小尺度)出發,根據自相似性質,研究其最終整體(大尺度)特征的理論基礎。
描述多重分形的一種方法是廣義分形維數,覆蓋法是分形研究中最通用的方法,它既用于簡單分形,也用于復雜分形。覆蓋法就是用尺度為ε的相同大小的盒子對整個集合進行覆蓋,所需盒子總數是N(ε),設點落入第i個盒子的概率為Pi(ε),對給定參數q,可計算出廣義信息熵Kq(ε)的表達式為[7]

從而廣義分形維數定義為

由于具有不同標度指數的子集可通過q的改變進行區分,當q=0時:

為容量維數(盒維數)。
當q=1時:

即為信息維數。
當q=2時,式(2)即為關聯維數:

由此說明廣義分形維數Dq包含了自相似分形理論涉及的大部分分形維數。
然而,由于盒維數的計算本質是對信號進行規則的網格劃分,文獻[8] 指出,這種規則的劃分網格的方法必將會導致對所需盒子數估計的誤差,從而影響分形維數估計的精度,而且誤差相當可觀。
為克服盒維數計算方法本質上的缺陷,必須尋找新的切入點來實現對分維數的精確估計。因此,本文提出了基于數學形態學的廣義分形維數估計方法。
數學形態學是以集合論、積分幾何和拓撲學為基礎發展起來的一種有別于時空域分析和頻域分析的數學方法。它的基本思想是用具有一定形態的結構元素去度量和提取信號中的對應形態,以達到對信號進行分析和識別的目的。膨脹和腐蝕是數學形態學的兩種基本運算,它們可以在保持信號的基本特征的前提下,去除信號中尺度上小于結構元的細節,從而得到結構簡化的信號數據。這就如同在不同的尺度下觀察信號,當使用尺度較小的結構元時,可以看到較多的細節如果換用尺度較大的結構元,看到的就只有粗的輪廓了。而各種分形維數估計算法的基本思想都是在不同尺度下對分形集合進行度量,而數學形態學正好為我們提供了一種在不同尺度下度量信號的數學工具[9]。所以,采用數學形態學來計算信號的分形維數應該是一個非常自然而恰當的選擇。
數學形態學的基本運算包括腐蝕和膨脹兩種算子。設f(n)和 g(m)分別為定義在 F={0,1,2,…,N -1}和 G={0,1,2,…,M -1}上的離散函數,且 N?M。這里為f(n)輸入信號,g(m)為結構元素。
f(n)關于g(m)的腐蝕定義為:

f(n)關于g(m)的膨脹定義為:

腐蝕和膨脹運算等價于離散函數在滑動濾波窗(相當于結構元素)內的最小值和最大值濾波。一般情況下,腐蝕運算減小了信號的峰值,加寬了谷域;而膨脹運算增大了信號的谷值,擴展了峰頂,關于數學形態學更加詳細的理論描述可參考文獻[10,11] 。
文獻[8] 提出了基于數學形態學的單一分形維數估計方法,其主要計算過程如下。
假設信號為f,g單位結構元素,定義尺度ε下的結構元為:

則在不同尺度ε下對信號f的膨脹和腐蝕分別為f⊕εg(n)和 fΘεg(n)。其中 ε =1,2,…,εmax,且 εmax≤N/2。
根據文獻[8] 的方法,選擇g為長度為3的扁平結構元素,即g={0,0,0},選擇扁平結構元素的好處是可以減少一部分計算量。
定義尺度ε對信號的覆蓋面積為:

文獻[8] 證明各尺度下Ag(ε)滿足:

則DM可由下式求得:

實際計算中對 log(Ag(ε)/ε2)和 log(1/ε)進行最小二乘線性擬合得到對信號分形維數DM的估計。
由前述的在不同尺度ε下對信號f的膨脹和腐蝕分別為f⊕εg(n)和 fΘεg(n),我們可以定義一個反映局部度量的分布函數ui(ε):

很顯然,上式中的 f⊕εg(n)-fΘεg(n)表示了信號膨脹結果與腐蝕結果之間的差異,其作用就像是單個網格上所需的盒子數,分布函數ui(ε)則描述了這種差異的分布,信號在尺度上的不均勻性可以通過分布函數的高階矩所表現出的奇異性來刻畫。
對給定參數q,可計算出形態學廣義信息熵Kq(ε)的表達式為:

由式(12)可以看出,基于數學形態學的網格劃分是固定的,即所有尺度下的網格數均為信號的點數N(ε)=N。如果采用計盒方法計算信號的覆蓋面積,Ag(ε)相當于盒子數 N(ε)與單位盒子面積 ε2的乘積,則:

從而保證了與廣義信息熵式(1)描述的一致性。
作為一個多重分形度量Kq(ε)與尺度之間必定滿足如下所示的指數關系:

則基于數學形態學的廣義分形維數可由下式求得:

實際計算中對Kq(ε)和log(ε)進行最小二乘線性擬合得到對廣義分形維數Dq的估計。
可以證明當 q=0 時,K0(ε)= - log(Ag(ε)/ε2),則多重分形維數Dq退化為基于數學形態學的單一分形維數。
實驗在F3L912型三缸四沖程柴油機上進行,加速度傳感器粘貼在第一缸缸蓋處,實驗在空間較大的車間內進行,計算機、電荷放大器等裝置在車間里的控制間內放置。發動機故障設置為一缸失火和進氣門間隙過大兩種故障。實驗時發動機轉速為1 200 r/min,采樣頻率為40 kHz,采樣點數為發動機一個工作周期。

圖1 發動機三種狀態下振動加速度信號時域波形圖Fig.1 Waveform of vibration signal from engine with three states
圖1為發動機在正常、失火和進氣門間隙過大三種狀態下的時域波形圖。從圖中可以看出,三種狀態下缸蓋振動加速度信號的時域波形具有較明顯的差異。失火時,燃爆信號消失,進氣門間隙過大時,進氣門關閉響應明顯變大。
以進氣門間隙過大時的信號為例說明基于數學形態學的分形維數計算方法,圖2為采用不同尺度的結構元素對信號進行形態學變換的結果。其中實線為采用尺度為20的扁平結構元、虛線為采用尺度為100的扁平結構元運算結果。

圖2 發動機氣門間隙故障信號及其腐蝕、膨脹圖Fig.2 Dilation and erosion waveform of the vibration signal from engine with valve fault
本文選擇分析信號的單位結構元素g為長度為3的扁平結構元素,即 g={0,0,0},最大尺度為100,即εmax=100,參數q取值范圍為[-20:2:20] ,然后根據2.3所描述方法計算信號的廣義分形維數。圖3給出了采用數學形態學操作所計算的進氣門間隙過大故障信號的廣義分形維數圖。從圖中可以看出,信號的廣義分形維數Dq隨q呈嚴格的遞減,當q=0時Dq就退化為基于數學形態學的單一分形維數。

圖3 發動機氣門間隙故障信號的形態學廣義分形維數Fig.3 Morphological generalized fractal spectrum of signal from engine with valve fault
圖4 為采用數學形態學方法計算的發動機三種狀態下振動信號的廣義分形維數譜,每種狀態取樣本數為4。

圖4 發動機三種狀態下信號廣義分形維數(盒計數法)Fig.4 Generalized fractal spectrums of signals from three engine states using box-counting method
從圖中可以看出,在q=0時即采用單一的分形維數很難將三種狀態進行區分;在q<0時三種狀態的分形維數存在著很大的混疊,完全不可分;但是當q>0,三種狀態的分形維數可以進行明顯的區分。這說明了采用單一的分形維數包含的信號奇異性信息十分有限,在部分情況下很難對不同狀態的信號進行區分,而廣義分形維數可以在不同的層次上對信號的奇異性進行分析,因此采用廣義分形維數對復雜信號進行分析具有很大的優越性和必要性。
同時,本文還采用傳統的盒計數法計算了發動機在三種狀態下的廣義分形維數。在計算中,盒子的尺度選擇為[2,4,8,16,32,64,128,256,512] 個采樣點。圖5給出了相應的廣義分形維數譜,與圖4相比可以發現,盒計數法計算的廣義分形維數對發動機的三種狀態區分性較差,不利于對發動機故障狀態進行分類。

圖5 發動機三種狀態下信號的形態學廣義分形維數Fig.5 Morphological generalized fractal spectrums of signals from three engine states
數學形態學為數字信號處理提供了一種新的快速分析手段,本文提出了一種基于數學形態學的廣義分形維數計算方法,并證明了與盒計數法計算廣義分形維數的一致性。通過對實際的發動機正常、失火故障和氣門間隙過大故障振動加速度信號的分析結果表明,與傳統的盒計數方法相比,基于數學形態學計算方法的廣義分維數估計能夠更加準確的將不同狀態下的發動機振動加速度信號進行區分,可以有效的用于發動機故障狀態的判別。此外,數學形態學只涉及簡單的加減和取大取小運算,因此與傳統的盒計數方法相比,基于數學形態學的廣義分形維數計算效率更高,為發動機故障特征提取和故障診斷提供了一種更為快速、有效的分析方法。
[1] Wu J D,Chen J C.Continuous wavelet transform technique for fault signal diagnosis of internal combustion engines[J] .NDT and E International,2006,39(4):304 -311.
[2] 金 萍,陳怡然,白 燁.內燃機表面振動信號的性質[J] .天津大學學報(自然科學與工程技術版),2000,33(01):63-68.
[3] 夏 勇,張振仁,陳衛昌,等.分形維數在內燃機振動診斷中的應用[J] .振動、測試與診斷,2001,21(03):209-212.
[4] 李 娜,方彥軍.利用關聯維數分析機械系統故障信號[J] .振動與沖擊,2007,26(4):136-139.
[5] 胥永剛,何正嘉.分形維數和近似熵用于度量信號復雜性的比較研究[J] .振動與沖擊,2003,22(03):25-29.
[6] Chaudhuri BB,Sarkar N.Texture segmentation using fractal dimension[J] .IEEE Transactions on Pattern Analysis and Machine Intelligence,1995,17(1):72-77.
[7] Grassberger P.Generalized dimensions of strange attractors[J] .Physics Letters A,1983,97(6):227-230.
[8] Maragos P,Sun F K.Measuring the fractal dimension of signals:morphological covers and iterative optimization[J] .IEEE Transactions on Signal Processing,1993,41(1):108-121.
[9] Radhakrishnan P,Lian T L,Daya Sagar B S.Estimation of fractal dimension through morphological decomposition[J] .Chaos,Solitons and Fractals,2004,21(3):563 -572.
[10] Maragos P,Schafer R W.Morphological filters-PartⅠ:Their set-theoretic analysis and relations to linear shiftinvariant filters[J] . IEEE Transactions on Acoustics,Speech,and Signal Processing,1987,ASSP - 35(8):1153-1169.
[11] Maragos P,Schafer R W.Morphological filters-PartⅡ:their relationsto median,order-statistic,and stack filters[J] .IEEE Transactions on Acoustics, Speech, and Signal Processing,1987,ASSP-35(8):1170-1184.