數(shù)學形態(tài)學分形維數(shù)廣義估算方法及其應用

王 冰， 李洪儒， 許葆華

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

引 言

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

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

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

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

1 基于數(shù)學形態(tài)學的分形維數(shù)計算方法

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

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

(1)

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

(2)

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

(3)

2 數(shù)學形態(tài)學的分形維數(shù)廣義估算方法

2.1 基于數(shù)學形態(tài)學的分形維數(shù)廣義估算方法

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

(4)

(5)

(6)

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

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

2.2 仿真分析

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

(7)

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

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

圖1 不同分形維數(shù)下的FBM時域波形

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

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

圖2 尺度為64時的不同定義下的形態(tài)學覆蓋

圖3描述了形態(tài)學覆蓋面積與結構元素尺度的變化關系?？梢钥闯觯叨仍酱?，形態(tài)學覆蓋面積越大。4種方法相比，基于形態(tài)梯度的原定義覆蓋面積最大，分析其原因，基于形態(tài)膨脹和形態(tài)腐蝕的定義分別用上覆蓋和下覆蓋代替原有覆蓋，導致覆蓋面積減小。基于形態(tài)差值的定義采用形態(tài)閉和開運算，運算量加倍，覆蓋面積最小。

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

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

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

表1 FBM信號的分形維數(shù)估計

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

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

3 基于分形維數(shù)和灰色關聯(lián)分析的軸承退化狀態(tài)識別

3.1 基于數(shù)學形態(tài)學分形維數(shù)的特征提取

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

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

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

(8)

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

圖4 仿真信號時域波形

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

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

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

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

3.2 性能退化狀態(tài)識別

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

圖6 基于分形維數(shù)和灰色關聯(lián)分析的軸承退化狀態(tài)識別流程圖

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

4 實驗數(shù)據(jù)驗證

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

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

圖7 滾動軸承內(nèi)圈在4種工作狀態(tài)下振動信號時域波形

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

4.1 特征提取

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

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

圖8 軸承內(nèi)外圈標準性能退化實驗數(shù)據(jù)特征參數(shù)

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

4.2 灰色關聯(lián)分析與退化狀態(tài)識別

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

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

圖9 軸承內(nèi)圈測試樣本Xij退化狀態(tài)識別結果

圖10 軸承外圈測試樣本Yij退化狀態(tài)識別結果

5 結 論

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

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

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

參考文獻：

[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設備退化狀態(tài)識別與故障預測方法[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的機械故障演化規(guī)律分析建模與預測[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] 潘玉娜,陳進. 小波包-支持向量數(shù)據(jù)描述在軸承性能退化評估中的應用研究[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] 潘玉娜,陳進. 結合循環(huán)平穩(wěn)和支持向量數(shù)據(jù)描述的軸承性能退化評估研究[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] 訾艷陽,胥永剛,何正嘉.離散振動信號分形盒維數(shù)的改進算法和應用[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] 李兵,張培林,任國全，等. 基于數(shù)學形態(tài)學的分形維數(shù)計算及在軸承故障診斷中的應用[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] 李兵,張培林,米雙山，等. 齒輪故障信號多重分形維數(shù)的形態(tài)學計算方法[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] 沈路,周曉軍,張文斌，等.基于形態(tài)濾波與灰色關聯(lián)度的滾動軸承故障診斷[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.