數學形態學分形維數廣義估算方法及其應用

王 冰， 李洪儒， 許葆華

(軍械工程學院導彈工程系， 河北 石家莊 050003)

引 言

滾動軸承是旋轉機械中的關鍵性元件，由于長期受到載荷、安裝、潤滑狀態等因素的影響，軸承極易受到損傷，進而影響整個機械設備的性能。因此，軸承一直都是設備狀態監測與故障診斷領域的熱門研究對象[1]。

隨著維修理論的發展，基于狀態的維修越來越得到人們關注。故障預測技術是實現基于狀態維修的關鍵。主要包括三個方面：運行狀態的特征提取；退化狀態識別；故障預測[2，3]。其中，運行狀態特征提取是實現退化狀態識別和故障預測的基礎，直接影響退化狀態識別的準確性和故障預測的可信性[4]。其著眼點在于特征信息對設備性能退化程度的量化反應能力，而非對不同故障模式的區分能力。文獻[4，5]提取并選取小波相關特征尺度熵作為預測特征指標；文獻[6，7]以小波包分解后的節點能量作為特征向量，定量評估樣本的退化程度；文獻[8]將通過循環平穩分析得到的組合切片累積能量作為預測特征值，均取得了較好的效果。

滾動軸承故障振動信號是一種典型的非平穩、非線性信號。分形維數則可以定量描述其復雜性和非線性。常見的分形維數有關聯維數、盒維數、信息維數等，其中，盒維數應用最為廣泛[9]。但由于現有的盒維數估計方法不準確[10]，文獻[11]提出一種基于數學形態學的分形維數計算方法，該方法計算量小，計算結果穩定且準確，已在軸承故障診斷領域得到了有效的應用[12，13]。

本文在基于數學形態學的分形維數計算方法基礎上，提出三種形態學分形維數廣義估算方法，以FBM仿真信號為例對其準確性和適用性進行分析。以此為基礎，提出基于數學形態分形維數與灰色關聯分析的退化狀態識別方法。并采用實測數據對該方法進行驗證。

1 基于數學形態學的分形維數計算方法

分形維數可以定量描述分形對象的復雜程度。在各種分形維數中，盒維數的應用最為廣泛。然而盒維數的計算過程由于對信號進行規則的網格劃分，因此影響了估算的精度[10]。針對此問題，文獻[11]提出一種基于數學形態學的分形維數估計方法。具體過程如下。

假設離散時間信號為f(n)，n=1,2,…,N，單位結構元素為g，則在尺度ε下所使用的結構元素定義為[14]

(1)

定義尺度ε對信號的覆蓋面積為

(2)

根據文獻[6]證明，覆蓋面積滿足如下條件

(3)

2 數學形態學的分形維數廣義估算方法

2.1 基于數學形態學的分形維數廣義估算方法

從幾何意義上分析，在上述方法中，由式(2)所定義的形態學覆蓋面積是為信號f(n)以結構元素εg進行膨脹和腐蝕運算后的差分，即形態梯度，本文稱該方法為基于形態梯度的形態覆蓋定義方法。為了分析其他形態學算子在形態學覆蓋定義中的有效性，本文提出3種形態學覆蓋廣義定義，其表達式如下：

(4)

(5)

(6)

根據覆蓋面積的定義方式，式(4)～(6)分別定義為基于形態膨脹、形態腐蝕和形態差值的形態學覆蓋。

以形態學覆蓋面積的廣義定義為基礎，通過最小二乘線性擬合即可得到對信號分形維數的廣義估計。

2.2 仿真分析

為了驗證三種形態學覆蓋廣義定義在計算中的有效性，采用FBM信號進行仿真分析。分數布朗運動(FBM)是描述時間(或空間)長程相關、功率譜滿足指數規律、增量服從正態分布的非平穩隨機過程的數學模型之一，是用來描述自然界隨機分形一種常用的分形模型，滿足如下定義[15]

(7)

式中k為常量，ξ為標準高斯隨機過程。H(0

本文采用分形工具箱Fraclab中的函數fbmlevinson來生成BH(t)[16]，采樣頻率為f=2 048 Hz，采樣點數為N=2 048，圖1顯示了3種不同維數(D=1.1，1.5，1.9)下的FBM分形信號。

圖1 不同分形維數下的FBM時域波形

分別根據原定義和3種廣義估算方法計算FBM信號分形維數。選擇扁平結構元素g=[0,0,0]，尺度進行離散化設置，取值為[2,4,8,16,32,64,128,256]。

圖2對比了在4種形態學覆蓋的定義下，對D=1.1的FBM信號進行尺度為64時的形態學覆蓋結果。在圖2(a)所示的原定義中，形態學覆蓋為膨脹和腐蝕之間所包含的區域，而其上覆蓋和下覆蓋分別為圖2(b)和(c)所描述的基于形態膨脹和形態腐蝕的形態學覆蓋。對比分析同一尺度下的形態學覆蓋面積大小，圖2(a)覆蓋最大，圖2(b)和(c)次之，圖2(d)中所定義的形態差值覆蓋最小。

圖2 尺度為64時的不同定義下的形態學覆蓋

圖3描述了形態學覆蓋面積與結構元素尺度的變化關系。可以看出，尺度越大，形態學覆蓋面積越大。4種方法相比，基于形態梯度的原定義覆蓋面積最大，分析其原因，基于形態膨脹和形態腐蝕的定義分別用上覆蓋和下覆蓋代替原有覆蓋，導致覆蓋面積減小。基于形態差值的定義采用形態閉和開運算，運算量加倍，覆蓋面積最小。

圖3 不同離散尺度下形態學覆蓋面積變化趨勢圖

由于分形維數是通過對覆蓋面積和尺度的最小二乘線性擬合而計算出的，因此形態覆蓋面積的變化趨勢影響著算法的有效性。從圖3可以看出，4種曲線具有近似的變化趨勢，從定性上分析，3種廣義估算方法是可行的。下面將進行定量對比分析。

表1為4種定義方法對FBM信號分形維數的定量估計值和相對誤差。實驗所采用計算機配置為Intel(R) Core i5-2400 CPU @3.10 GHz，內存為4 GB，軟件為Matlab R2010b。

表1 FBM信號的分形維數估計

可以看出，基于形態差值的定義方法計算結果明顯偏低，誤差較大，耗時較長。而基于形態膨脹和形態腐蝕的方法與原算法相比，計算精度相當，誤差均較低。并且由于單次運算只進行一次腐蝕或膨脹操作，計算量減半，運算速度更快。

綜上分析，基于形態膨脹和形態腐蝕的分形維數廣義估算方法在運算精度方面與原算法相當，運算速度更快，可以在實際中替代原算法。基于形態差值的分形維數廣義估算方法效率精度較低，運算速度較慢，不適合實際運用。

3 基于分形維數和灰色關聯分析的軸承退化狀態識別

3.1 基于數學形態學分形維數的特征提取

如何選取能夠表征軸承性能退化程度的特征指標，是對軸承進行性能退化狀態識別的基礎。所選取的特征指標應具有良好的穩定性和關聯性。穩定性是指該指標應對工況變化以及環境噪聲具有很強的抗干擾能力；關聯性是指該指標的變化趨勢應與設備的性能退化過程具有相關聯，隨性能的退化表現出單調性的變化趨勢。

從軸承振動信號的本質分析，其信號在全壽命周期中的變化過程，就是其內部隨機成分所占比例不斷變化的過程。因此，提出基于數學形態學分形維數的特征提取方法，用可以表征分形對象復雜度的分形維數表征軸承的性能退化程度。

為了該方法的有效性，采用仿真信號進行分析[17]

(8)

式中 0.1t2cos(2π×10t+2)為故障模擬信號，cos(2π×50t)為常規振動信號，n(t)代表高斯噪聲，幅值0.1t2用來近似反映故障隨時間的變化過程。信號采樣點數為N=10 240，采樣頻率為1 024 Hz。為了研究白噪聲的影響，分別設定3種噪聲強度1，3，6，圖4顯示了噪聲強度為1時的仿真信號時域波形。

圖4 仿真信號時域波形

將信號等分為10段并順序標記，用10組數據近似描述仿真信號故障程度不斷加深的性能退化過程。采用基于形態膨脹的數學形態學分形維數計算方法計算每一組數據的分形維數。選擇扁平型單位結構元素g=[0,0,0]，尺度為[2,4,8,16,32,64,128,256]，不同噪聲強度下的仿真信號分形維數變化趨勢如圖5所示。

圖5 不同噪聲強度下分形維數變化趨勢圖

可以看出，隨著故障程度的不斷加深，基于數學形態學的分形維數數值不斷減小，其值在區間[1.2,1.58]內變化，呈現出遞減的趨勢。且在故障初期變化緩慢，分析認為這與仿真信號以0.1t2進行幅值變化緊密相關。另外，隨著噪聲強度的不斷加強，其中，當噪聲強度為1時，分形維數的數值在區間[1.48,1.72]內變化，噪聲強度為6時，其數值在區間[1.6,1.78]范圍內變化，可見，噪聲強度越大，形態分形維數的計算結果越大，但從整體而言，依然保持遞減趨勢。綜上分析，基于數學形態學的分形維數與性能退化過程具有較好的關聯性，并且對噪聲具有一定的抗干擾能力。

為了增強預測特征向量的完備性，本文另選取已在文獻[18]中驗證的峰值、峰峰值、有效值3個故障預測特征參數，以此構成對電機軸承性能退化狀態識別的四維預測特征向量。

3.2 性能退化狀態識別

性能退化狀態識別是一種多類分類問題，與故障模式識別側重于對多故障類型分類不同，性能退化狀態識別更強調對性能退化狀態的分類。灰色關聯分析是灰色系統理論進行系統分析的重要方法，它是根據系統各個因素之間的內部聯系或發展態勢的相似程度來度量因素之間關聯度的方法[14]。本節在前文基礎上提出基于分形維數和灰色關聯分析的軸承性能退化狀態識別方法。其步驟如圖6所示。

圖6 基于分形維數和灰色關聯分析的軸承退化狀態識別流程圖

首先確定樣本對象，將性能退化過程振動信號分為標準退化實驗樣本數據和待識別樣本實驗數據。之后按照性能退化的程度將性能退化過程分為n個退化狀態，分別提取每個退化狀態下的特征向量，并對結果進行歸一化處理。特征向量記為F={F1,F2,…,Fi}，其中i為向量中特征參數的個數。本文取i=4。在此基礎上，建立標準退化模式矩陣和待識別樣本特征向量矩陣。計算待識別樣本特征向量與標準模式退化模式矩陣的灰色關聯度，關聯度最大值即為待識別樣本的性能退化狀態。

4 實驗數據驗證

為驗證本文所提方法在電機滾動軸承性能退化狀態識別中的有效性，采用美國凱斯西儲大學電氣工程實驗室的滾動軸承振動加速度實驗數據進行實例分析[19]。測試軸承為SKF6205-2RS深溝球軸承，采樣頻率為48 kHz。轉速為1 730 r/min，載荷為2.25 kW。以4種不同故障程度的軸承內圈和外圈振動信號近似模擬內外圈由正常狀態逐漸經歷一系列不同的退化狀態直至完全失效的這一性能退化過程。故障程度通過軸承內外圈點蝕的直徑來描述。分別以直徑0，0.18，0.36，0.54 mm的點蝕故障模擬軸承內外圈正常狀態、輕度故障、中度故障、重度故障4種工作狀態。

選取軸承內外圈在4種工作狀態下的振動加速度數據，數據長度為6 000，分別記為Si(i=1,…，4)，Ti(i=1,…，4)，其中i為軸承工作狀態序號，以此作為軸承內外圈標準性能退化實驗數據；另外選取內外圈在4種工作狀態下長度為1 000的振動加速度數據，分別記為Xij(i=1,…，4;j=1,…，10)，Yij(i=1,…，4;j=1,…，10)，其中i為軸承工作狀態序號，j為待識別樣本序號，以此作為軸承內圈和外圈的待識別樣本數據。

圖7 滾動軸承內圈在4種工作狀態下振動信號時域波形

圖7為軸承內圈在4種工作狀態下的振動加速度數據Si的原始時域波形。由圖中可以明顯看出，4種狀態下的軸承振動信號在時域結構上具有較為明顯的差異。正常信號呈現出幾乎隨機性的分布，隨著內圈故障程度的不斷加深，信號的規律性明顯增強，幅值也不斷增大。

4.1 特征提取

采用3.1節提出的特征提取方法，分別計算標準性能退化實驗數據(Si，Ti)和待識別樣本數據(Xij，Yij)的特征向量，建立標準退化模式矩陣和待識別樣本特征向量矩陣。圖8(a)和(b)分別為經歸一化處理后的軸承內外圈標準性能退化實驗數據特征參數。其中，橫坐標代表不同的性能退化程度，序號1，2，3，4分別代表軸承內外圈正常狀態、輕度故障、中度故障、重度故障4種工作狀態。

可以看出，隨著軸承內外圈的故障程度的不斷加深，4個特征參量均發生明顯的關聯性變化，分形

圖8 軸承內外圈標準性能退化實驗數據特征參數

維數隨著故障退化程度的加深而不斷減小，峰值、峰峰值和有效值隨著故障退化程度的加深而不斷增加。且分形維數的分布更加平均，區分效果更好，有效值的效果次之，峰值和峰峰值對輕度和重度故障的區分效果不明顯。

4.2 灰色關聯分析與退化狀態識別

根據特征提取結果，分別建立軸承內外圈的標準退化模式矩陣Fin和Fout，軸承內外圈待識別樣本特征向量矩陣Aij(i=1,…,4;j=1,…,10),Bij(i=1,…,4;j=1,…,10)，其中i為實際軸承工作狀態序號，j為待識別樣本序號。

分別計算Aij與Fin,Bij與Fout的灰色關聯度，分辨系數ρ=0.5。根據灰色關聯度的大小識別待識別樣本的性能退化狀態。圖9(a)～(d)顯示了不同內圈退化狀態樣本的識別結果。可以看出，除中度和重度故障狀態測試樣本中各有2個未被正確識別外，其余各個樣本均可被正確識別。其總體狀態識別成功率達到90%；在圖10(a)～(d)所顯示的軸承外圈故障性能退化狀態識別結果中，除重度故障狀態測試樣本有2個樣本未被正確識別外，其余各個樣本均被正確的識別。其總體狀態識別成功率達到95%。由此可見，在軸承內外圈性能退化狀態識別中，基于數學形態學分形維數與灰色關聯分析的方法可以達到較高的識別率，效果理想。

圖9 軸承內圈測試樣本Xij退化狀態識別結果

圖10 軸承外圈測試樣本Yij退化狀態識別結果

5 結 論

(1)基于形態膨脹、形態腐蝕、形態差值的分形維數廣義估算方法是從形態覆蓋的角度對原有算法的廣義擴展。通過仿真實驗可以看出，基于形態膨脹、形態腐蝕的估算方法在精確度和運算速度方面均具有一定的優勢，而基于形態差值的估算方法則運算量大，精確度較之其他方法也略低。

(2)基于形態分形維數與灰色關聯分析的性能退化狀態識別方法將形態學分形維數引入到預測特征提取中，在仿真和實例數據的分析處理中取得了較為理想的效果。

(3)機械設備在實際運行中的性能退化過程要復雜得多，需要考慮轉速、載荷及環境因素時變性帶來的影響，且不易采集。采用不同損傷程度實驗數據是一種用離散近似模擬連續的方法，必然有一定的局限性。

參考文獻：

[1] Patil M S,Mathew J,Rajendra Kumar P K. Bearing signature analysis as a medium for fault detection: a review[J].Journal of Tribology,2008,130(1) : 1—7.

[2] QIU H, LEE J, LIN J, et al. Robust performance degradation assessment methods for enhanced rolling element bearing prognostics[J]. Advanced Engineering Informatics, 2003, 17: 127—140.

[3] ZHANG X D, XU R, KWAN C, et al. An integrated approach to bearing fault diagnostics and prognostics[A].2005 American Control Conference[C].Portland, 2005: 2 750—2 755.

[4] 曾慶虎,邱靜,劉冠軍. 基于小波相關特征尺度熵的HSMM設備退化狀態識別與故障預測方法[J]. 儀器儀表學報,2008,12(29):2 559—2 564.ZENG Qing-hu, QIU Jing, LIU Guan-jun. Equipment degradation state recognition and fault prognosis method based on wavelet correlation feature scale entropy and HSMM[J]. Chinese Journal of Scientific Instrument,2008,12(29):2 559—2 564.

[5] 曾慶虎,邱靜,劉冠軍. 基于HSMM的機械故障演化規律分析建模與預測[J]. 機械強度,2010,32(5):695—701.ZENG Qing-hu, QIU Jing, LIU Guan-jun. Modeling and prognostics of mechanical fault based on HSMM[J]. Journal of Mechanical Strength,2010,32(5):695—701.

[6] 肖文斌,陳進，周宇，等. 小波包變換和隱馬爾可夫模型在軸承性能退化評估中的應用[J]. 振動與沖擊，2011,30(8):32—35.XIAO Wen-bin, CHEN Jin， ZHOU Yu,et al. Wavelet packet transform and hidden Markov model based bearing performance degradation assessment[J]. Journal of Vibration and Shock，2011,30(8):32—35.

[7] 潘玉娜,陳進. 小波包-支持向量數據描述在軸承性能退化評估中的應用研究[J]. 振動與沖擊，2009,28(4):164—167.PAN Yu-na, CHEN Jin. Wavelet package-support vector data description applied in bearing performance degradation assessment[J]. Journal of Vibration and Shock， 2009,28(4):164—167.

[8] 潘玉娜,陳進. 結合循環平穩和支持向量數據描述的軸承性能退化評估研究[J]. 機械科學與技術,2009,28(4):442—445.PAN Yu-na, CHEN Jin. Assessment of bearing performance degradation by cyclostationarity analysis and support vector data description[J]. Mechanical Science and Technology for Aerospace Engineering, 2009,28(4):442—445.

[9] 訾艷陽,胥永剛,何正嘉.離散振動信號分形盒維數的改進算法和應用[J].機械科學與技術, 2001, (03): 373—376.ZI Yan-yang, XU Yong-gang, HE Zheng-jia, Application and improved algorithm of fractal box dimension for disperse vibration signal[J]. Mechanical Science and Technology, 2001,(03): 373—376.

[10] Chaudhuri B B, Sarkar N. Texture segmentation using fractal dimension[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995. 17(1):72—77.

[11] 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.

[12] 李兵,張培林,任國全，等. 基于數學形態學的分形維數計算及在軸承故障診斷中的應用[J]. 振動與沖擊,2010,29(5):191—194.LI Bing, ZHANG Pei-lin, REN Guo-quan, et al. Mathematic morphology-based fractal dimension calculation and its application in fault diagnosis of roller bearings[J]. Journal of Vibration and Shock, 2010，29(5):191—194.

[13] 李兵,張培林,米雙山，等. 齒輪故障信號多重分形維數的形態學計算方法[J]. 振動、測試與診斷, 2011,31(4)：450—453.LI Bing, ZHANG Pei-lin, MI Shuang-shan,et al. Mathematical morphology based on multifractal dimensions for gear fault diagnosis[J]. Journal of Vibration, Measurement & Diagnosis，2011,4(31):450—453.

[14] 沈路,周曉軍,張文斌，等.基于形態濾波與灰色關聯度的滾動軸承故障診斷[J].振動與沖擊. 2009,28(11):17—20.SHEN Lu, ZHOU Xiao-jun, ZHANG Wen-bin, et al. Fault diagnosis of rolling elem ent bearing based on morphological filter and grey incidence[J]. Journal of Vibration and Shock. 2009,28(11):17—20.

[15] 朱波華. 基于分形理論的信號處理與分析的研究[D]. 武漢：武漢理工大學,2008.ZHU Bo-hua. Research on signal process and analysis based on fractal theory [D].Wuhan : Wuhan University of Technology, 2008.

[16] Fraclab: A fractal analysis toolbox for signal and image processing[EB/OL]. http://fraclab.saclay.inria.fr/, 2009.

[17] 徐東,徐永成，陳循，等. 基于EMD的灰色模型的疲勞剩余壽命預測方法研究[J]. 振動工程學報，2011,24(1):104—110.XU Dong, Xu Yong-cheng, CHEN Xun, et al. Residual fatigue life prediction based on grey model and EMD[J]. Journal of Vibration Engineering, 2011,24(1):104—110.

[18] 潘玉娜.滾動軸承的性能退化特征提取及評估方法研究[D]. 上海：上海交通大學，2011.PAN Yu-na. Study on feature extraction and assessment method of rolling element bearing performance degradation[D]. Shanghai: Shanghai Jiao Tong University,2011.

[19] Kenneth A Loparo.Bearing data center [EB/OL].Case western reserve university , 2000[2010-03-26] . http://www.eecs.case.edu/laboratory/bearing/welcome_overview.htm.