張懷亮 鄒佰文 肖 雷

中南大學(xué)，長沙，410083

0 引言

機械設(shè)備在運轉(zhuǎn)過程中產(chǎn)生的磨粒群包含豐富的摩擦磨損信息，對磨粒群特征進行提取，可用于設(shè)備的狀態(tài)監(jiān)測和故障診斷。研究表明，磨粒群在一定尺度上表現(xiàn)出分形特性，傳統(tǒng)基于歐氏幾何的參數(shù)無法準確描述磨粒群特征［1－2］。分形維數(shù)是表征磨粒群分形特性的一個重要參數(shù)，因而提高磨粒群分形維數(shù)的計算精度，對于磨粒群特征的提取和設(shè)備磨損狀態(tài)監(jiān)測有重要意義。由于產(chǎn)生機制復(fù)雜，磨粒群的分形自相似特性是近似或統(tǒng)計意義上的，即在一定尺度范圍內(nèi)的自相似，這個尺度范圍即稱為無標度區(qū)，在無標度區(qū)外討論磨粒群的分形特性沒有意義，因此，無標度區(qū)的準確辨識是磨粒群分形維數(shù)計算中必須解決的關(guān)鍵問題之一。

辨識無標度區(qū)的常用方法有以下幾類：①基于r～N（r）（尺度～測度）雙對數(shù)曲線的人工識別法［3］；②基于r～N（r）點對相關(guān)性檢驗的排序逼近法［4－8］；③基于三折線、曲線－直線－曲線、反 S等模型的模型擬合法［9－11］。人工識別法主觀性強，排序逼近法在r～N（r）點對較多時計算量很大且較易出現(xiàn)局部最優(yōu)，模型擬合法在模型的選取方面缺乏理論依據(jù)，應(yīng)用時有一定局限性。因而，Yang等［12］利用無標度區(qū)一階導(dǎo)數(shù)近似不變的特點提出應(yīng)用聚類分析理論辨識無標度區(qū)的算法；Ji等［13］在此基礎(chǔ)上結(jié)合K－means和點斜率誤差法提高了算法的精度；王成棟等［14］根據(jù)無標度區(qū)二階導(dǎo)數(shù)在0附近波動的特點提出了一種自動辨識方法，這類利用r～N（r）雙對數(shù)曲線導(dǎo)數(shù)特點的辨識算法幾何意義明確，計算量較小，但大多對初始值敏感，容易出現(xiàn)局部最優(yōu)，可能導(dǎo)致辨識出來的無標度區(qū)雖然相關(guān)性檢驗參數(shù)符合要求，但并不是完整的無標度區(qū)，甚至不在真實無標度區(qū)內(nèi)，即無標度區(qū)的局部最優(yōu)。本文結(jié)合模擬退火和K－means算法對雙對數(shù)曲線的一階、二階導(dǎo)數(shù)進行全局聚類，可有效避免陷入局部最優(yōu)，從而提高磨粒群分形維數(shù)的計算精度。

1 磨粒群分形維數(shù)計算理論

磨粒群分散于設(shè)備的潤滑油中，通過鐵譜儀可將潤滑油樣中的磨粒群分離出來，制成鐵譜片，經(jīng)鐵譜顯微鏡拍照成像，得到磨粒群的圖像（圖1a），原始圖像背景復(fù)雜，通過背景去噪、閾值分割等圖像預(yù)處理過程，得到磨粒群的二值圖像（圖1b），磨粒群分形維數(shù)的計算過程基于其二值圖像。實際應(yīng)用中，為保證狀態(tài)監(jiān)測數(shù)據(jù)的可比性，獲取磨粒群圖片時應(yīng)遵守同樣的取樣倍數(shù)及視場選取規(guī)范。

分形維數(shù)的定義包括相似維、Hausdorff維、信息維、計盒維數(shù)等多種，其中，計盒維數(shù)（boxcounting dimension）由于其物理意義明確，便于程序化計算等優(yōu)點，廣泛應(yīng)用于圖像分形特征的表征。計盒維數(shù)的計算原理如下：將邊長為r的方形網(wǎng)格覆蓋所需測量的圖形，測算與圖形重疊的網(wǎng)格數(shù)目為測度N（r），不斷改變網(wǎng)格邊長，重復(fù)上述過程，得到一組不同尺度r下的測度N（r），那么計盒維數(shù)即為

圖1 磨粒群原始圖像及二值圖像

根據(jù)計盒維數(shù)的定義，磨粒群分形維數(shù)計算的一般過程如下：假設(shè)磨粒群圖像的分辨率為M×M，以邊長為r個像素的網(wǎng)格不斷覆蓋磨粒群圖像并統(tǒng)計重疊網(wǎng)格數(shù)目N（r），其中r的取值范圍為｛r｜1≤r＜ M／2，r＝1，2，4，…，2n｝［15］，根據(jù)r～N（r）點對序列繪出雙對數(shù)曲線，將曲線中的直線部分作為無標度區(qū)間，在無標度區(qū)間內(nèi)對r～N（r）點對進行線性回歸，取分形維數(shù)D＝－α，其中α為回歸直線的斜率。

在磨粒群分形維數(shù)的計算過程中，無標度區(qū)的辨識結(jié)果直接影響磨粒群分形維數(shù)的計算精度，辨識出來的無標度區(qū)若出現(xiàn)偏離、過短或過長等誤差，會導(dǎo)致無法正確評估磨粒群的分形維數(shù)。而應(yīng)用已有的無標度區(qū)辨識方法計算磨粒群的分形維數(shù)時，往往會出現(xiàn)局部最優(yōu)的現(xiàn)象，如圖2所示。圖2中的數(shù)據(jù)點為磨粒群圖像（圖1）應(yīng)用計盒維數(shù)算法得到的r～N（r）點對，A、B、C、D 4個區(qū)為不同的無標度區(qū)辨識結(jié)果，其中，A區(qū)、B區(qū)、C區(qū)數(shù)據(jù)點的線性相關(guān)系數(shù)均大于0.99，D區(qū)為0.97，將D區(qū)作為無標度區(qū)不嚴格，C區(qū)lnN（r）幾乎為定值，即分形維數(shù)為0，顯然不合理，B區(qū)雖然相關(guān)性符合要求，但長度太短，A區(qū)是較合理的無標度區(qū)。相對于A區(qū)、B區(qū)、C區(qū)為無標度區(qū)辨識過程中出現(xiàn)的兩種典型的局部最優(yōu)情況，即辨識結(jié)果雖在真實無標度區(qū)內(nèi)，但長度過短和偏離真實無標度區(qū)，4個區(qū)內(nèi)計算分形維數(shù)的結(jié)果為1.650、1.603、0.011、1.723，差異很大，所以，無標度區(qū)的準確辨識是磨粒群分形維數(shù)計算的關(guān)鍵環(huán)節(jié)。

圖2 無標度辨識中的局部最優(yōu)

2 分形無標度區(qū)辨識新算法

為避免無標度區(qū)辨識算法陷入局部最優(yōu)，將K－means聚類算法與全局搜索算法模擬退火算法相結(jié)合，從r～N（r）雙對數(shù)曲線中準確辨識磨粒群無標度區(qū)。

2.1 模擬退火K－means算法

K－means算法是一種基于距離的聚類算法，該算法以距離為相似性指標，將觀察對象分為k類［16］?；具^程如下：首先從研究對象中隨機選擇k個對象作為聚類中心，其余的對象根據(jù)它們到聚類中心的距離劃分到最近的聚類。然后，以各聚類中所有對象的平均值作為新的聚類中心，重復(fù)上述劃分過程，直到準則函數(shù)收斂。應(yīng)用K－means算法對雙對數(shù)曲線的導(dǎo)數(shù)進行聚類，可分離出無標度區(qū)，但是從K－means的聚類過程可以看出，算法對初始聚類中心的選擇依賴很大，初始聚類中心選擇不當，可能出現(xiàn)局部最優(yōu)，將它應(yīng)用到磨粒群無標度區(qū)的識別可能會出現(xiàn)無標度區(qū)過短等現(xiàn)象。

模擬退火算法是由Kirkpatrick提出的一種啟發(fā)式、隨機優(yōu)化算法［17］。該算法與初始值無關(guān)，算法求得的解與初始解狀態(tài)無關(guān)，具有漸近收斂性，已在理論上被證明是一種以概率1收斂于全局最優(yōu)解的全局優(yōu)化算法。因而，利用模擬退火算法能夠彌補K－means算法依賴初始聚類中心的缺點，跳出局部最優(yōu)。

模擬退火K－means算法以K－means聚類的結(jié)果作為模擬退火算法的初始解，通過隨機改變聚類中某個元素的類別產(chǎn)生新解，目標函數(shù)為當前聚類劃分的總類間離散度：

模擬退火K－means算法的實現(xiàn)步驟如下：

（1）將數(shù)據(jù)樣本的K－means聚類結(jié)果作為初始解w，根據(jù)式（2）計算目標函數(shù)值Jw。

（2）初始化溫度t0，t0＝Jw，初始化退火速度a和最大退火次數(shù)。

（3）在某一溫度t，迭代進行步驟（4）～ 步驟（6），直到最大迭代次數(shù)跳到步驟（7）。

（4）通過隨機改變一個聚類樣本的當前所屬類別產(chǎn)生新解w′，計算新的目標函數(shù)值J′w。

（5）判斷J′w是否為最優(yōu)目標函數(shù)值，若是，則保存聚類劃分w′為最優(yōu)聚類劃分、J′w為最優(yōu)目標函數(shù)值；否則跳到下一步。

（6）計算新的目標函數(shù)值與當前目標函數(shù)值的差ΔJ。判斷ΔJ是否小于0：若ΔJ＜0，則接受新解，即將新解作為當前解；若ΔJ≥0，則根據(jù)Metropolis準 則，以 概 率 p（p ＝ eΔJ／（kt））接 受新解。

（7）判斷是否達到最大退火次數(shù)，若是則結(jié)束算法，輸出最優(yōu)聚類劃分；否則對溫度t進行退火，返回步驟（3）繼續(xù)迭代。

2.2 無標度區(qū)辨識新算法

在計算磨粒群分形維數(shù)時，辨識無標度區(qū)可簡化為提取r～N（r）雙對數(shù)曲線的直線部分。由于直線的一階導(dǎo)數(shù)是常數(shù)，不隨自變量的改變而變化，而曲線的一階導(dǎo)數(shù)隨著自變量改變而變化，直線的二階導(dǎo)數(shù)為0，曲線的二階導(dǎo)數(shù)不為0，利用這個性質(zhì)，可以通過雙對數(shù)曲線的一階、二階導(dǎo)數(shù)特點識別出曲線的線性部分，作為無標度區(qū)。實際計算得到的并不是連續(xù)的曲線，而是離散的r～N（r）點對序列，因而定義雙對數(shù)曲線的一階、二階局部導(dǎo)數(shù)分別為

理論上，無標度區(qū)內(nèi)點的一階局部導(dǎo)數(shù)為常數(shù)，但磨粒群并不是完全理想的分形圖形，所以實際計算中，無標度區(qū)內(nèi)點的一階局部導(dǎo)數(shù)在一個常數(shù)附近做微小波動，而無標度區(qū)外點的波動幅度較大；無標度區(qū)內(nèi)點的二階局部導(dǎo)數(shù)在0附近作微小波動，而無標度區(qū)外的點偏離0較遠。因而，應(yīng)用模擬退火K－means算法對雙對數(shù)的一階局部導(dǎo)數(shù)聚類，分離出無標度區(qū)的大致范圍，對其二階局部導(dǎo)數(shù)聚類，辨識出無標度區(qū)的準確區(qū)間。磨粒群分形無標度區(qū)辨識的新算法如下：

首先，應(yīng)用模擬退火K－means算法對磨粒群分形維數(shù)計算中的雙對數(shù)曲線進行一階導(dǎo)數(shù)聚類，結(jié)果如圖3a所示。經(jīng)過聚類，雙對數(shù)曲線一階局部導(dǎo)數(shù)被分為兩類（圖中l(wèi)n·N（r）group1和ln·N（r）group2），兩類離散點交錯分布，將初始的雙對數(shù)曲線分成了許多段區(qū)間（圖3a中豎線分割），各區(qū)間長短不一。其中，當lnr較大時，出現(xiàn)了多段一階導(dǎo)數(shù)連續(xù)為0點區(qū)間，這是由于采用計盒法計算時，當盒子尺寸較大時，盒子尺寸的改變量相對盒子本身尺寸很小，覆蓋磨粒群的盒子個數(shù)并不隨盒子尺寸的變化而改變，這些區(qū)間不符合無標度區(qū)的定義，應(yīng)予以剔除。剩下的區(qū)間中，多為較短的區(qū)間，甚至只有一個點，由于磨粒群的分形無標度區(qū)一般較寬，因而選取包含連續(xù)點最多的區(qū)間作為無標度區(qū)，識別結(jié)果如圖3b所示。需要指出的是，雖然從圖3b上看，被剔除的區(qū)間較短，但由于雙對數(shù)曲線的壓縮特性，實際被剔除的點的數(shù)量是很大的，而這些點都明顯不屬于無標度區(qū)，說明計算過程中r的取值范圍可以進一步縮小以提高分形維數(shù)計算效率。在磨粒群分形維數(shù)計算過程中，如果對這些點不加以剔除，計算結(jié)果將出現(xiàn)很大的偏差。

通過上述的識別過程，無標度區(qū)的大概范圍已被確定，但區(qū)間還存在一些波動相對較大的部分（一階導(dǎo)數(shù)右端），為精確確定無標度區(qū)，還需進一步剔除這些點。由于一階導(dǎo)數(shù)在常數(shù)附近波動，其二階導(dǎo)數(shù)應(yīng)該在0附近作微小波動，若二階導(dǎo)數(shù)偏離0的幅度較大，則表明其一階導(dǎo)數(shù)波動劇烈，因而，對其二階導(dǎo)數(shù)進行模擬退火K－means聯(lián)合聚類，結(jié)果如圖3c所示，一階導(dǎo)數(shù)的波動在二階導(dǎo)數(shù)上表現(xiàn)得更明顯，與一階導(dǎo)數(shù)類似，通過聚類，區(qū)間分割，選取包含連續(xù)點數(shù)最多的區(qū)間，最終得到的無標度區(qū)如圖3d所示。磨粒群無標度區(qū)識別結(jié)果中，雙對數(shù)曲線二階導(dǎo)數(shù)在0附近作微小起伏，說明雙對數(shù)曲線在該區(qū)間內(nèi)接近直線，即磨粒群的無標度區(qū)。

對識別的無標度區(qū)內(nèi)的點進行線性回歸，相關(guān)系數(shù)R＝0.9994，說明無標度區(qū)內(nèi)r～N（r）雙對數(shù)曲線線性相關(guān)性很強，符合無標度區(qū)的定義。同時，整個辨識磨粒群分形無標度區(qū)的過程中，未出現(xiàn)辨識結(jié)果過短或落在無標度區(qū)外等局部最優(yōu)現(xiàn)象，在保證通過相關(guān)性檢驗的前提下盡可能選取了較長的區(qū)間，有效地避免了陷入局部最優(yōu)。對無標度區(qū)內(nèi)的點進行最小二乘擬合，即可得到直線的斜率α，根據(jù)盒維數(shù)的定義，取D＝－α，作為磨粒群的分形維數(shù)。

圖3 磨粒群分形無標度區(qū)辨識過程

2.3 方法檢驗及結(jié)果分析

為檢驗新方法對無標度區(qū)識別的準確性及精度，選取Koch曲線、Sierpinski三角形、Sierpinski方毯三個已知理論分形維數(shù)的圖形作為測試對象，三個圖形的大小一致，采用計盒法獲取三個圖形的r～N（r）雙對數(shù)曲線，利用新方法識別相應(yīng)的無標度區(qū)，對無標度區(qū)內(nèi)的點進行最小二乘擬合，確定圖形的分形維數(shù)。同時，為了對比研究本文無標度區(qū)識別方法的精度，同時用人工識別法和文獻［14］提出的二階導(dǎo)數(shù)法對三個圖形無標度區(qū)進行了識別，計算結(jié)果如表1所示，圖4示出了三種方法計算結(jié)果與理論值的誤差。

表1 測試圖形分形維數(shù)計算結(jié)果

圖4 三種無標度區(qū)識別方法相對誤差

應(yīng)用新方法辨識三個圖形無標度區(qū)的過程中均出現(xiàn)局部最優(yōu)現(xiàn)象，圖4表明，人工識別法在三個圖形的無標度區(qū)識別中都取得了較好的效果，平均誤差在3%～5%之間，但是人工識別法需與計算機進行交互操作，無法進行自動識別，且不同觀察者一般得到不同的結(jié)論，主觀性強；二階導(dǎo)數(shù)法對Koch曲線的識別效果很好，但對其他兩個圖形的無標度區(qū)識別效果較差；本文提出的識別方法雖對Koch曲線的識別效果與人工識別法效果相當，但對Sierpinski三角形和Sierpinski方毯的識別誤差均在1%以下，由于Sierpinski三角形和Sierpinski方毯的特點與磨粒群的特點相似，都是一種多孔結(jié)構(gòu)，因而應(yīng)用本文提出的方法對磨粒群分形無標度區(qū)的識別是合理的，分形維數(shù)的計算結(jié)果是可信的。

3 磨粒群的分形無標度區(qū)辨識

為分析新方法對磨粒群無標度區(qū)的辨識能力，應(yīng)用新方法計算了正常磨損、切削磨損、滾動疲勞磨損、滾動和滑動復(fù)合磨損及嚴重滑動磨損的磨粒群分形維數(shù)，磨粒群圖片來自文獻［18］，圖片大小為400像素×400像素，盒子尺寸增長方式為 ｛r｜1≤r＜M／2，r＝1，2，4，…，2n｝，計算程序基于MATLAB開發(fā)，計算機CPU主頻為3.00GHz，計算結(jié)果如表2所示，其中，無標度區(qū)中的數(shù)值范圍表示r～N（r）雙對數(shù)曲線中識別的無標度區(qū)取點范圍。

表2 應(yīng)用新方法計算磨粒群分形維數(shù)結(jié)果

計算結(jié)果表明，在圖片大小為400像素×400像素，盒子尺寸增長方式為1，2，4，…，2n時，磨粒群的分形無標度區(qū)大約為1～20個點之間，取此區(qū)間內(nèi)的點對進行擬合，可以得到較精確的磨粒群分形維數(shù)，這段區(qū)間內(nèi)，點對的線性相關(guān)系數(shù)大于0.99，同時區(qū)間長度占整個雙對數(shù)曲線的50%左右，未出現(xiàn)無標度區(qū)辨識的局部最優(yōu)。另外，試驗中不同磨損狀態(tài)的磨粒群分形維數(shù)差異較大，從1.44到1.71，這表明磨粒群的分形維數(shù)可作為摩擦副磨損狀態(tài)識別的一個判據(jù)。新算法對給定條件下磨粒群無標度區(qū)進行識別所耗用CPU時間約為0.24s，同樣條件下應(yīng)用文獻［8］、［13］和［14］提出的方法對磨粒群無標度區(qū)進行識別所用CPU 時間分別為0.3s、0.25s、0.8s左右，相對于已有算法，新算法在保持較高的無標度區(qū)識別精度的同時，識別效率也略有提高。

4 結(jié)論

（1）基于模擬退火K－means算法提出了一種辨識磨粒群分形無標度區(qū)的新算法，新算法能有效地避免無標度區(qū)辨識過程中陷入局部最優(yōu)，從而準確辨識磨粒群的分形無標度區(qū)，提高分形維數(shù)的計算精度。

（2）新算法對磨粒群圖像、Sierpinski方毯等多孔狀圖形無標度區(qū)辨識精度高，辨識得到的磨粒群分形無標度區(qū)較寬，區(qū)間內(nèi)點對的線性相關(guān)性好。

［1］Podsiadlo P，Stachowiak G W.Scale－invariant Analysis of Wear Particle Surface Morphology：II.Fractal Dimension［J］.Wear，2000，242（1／2）：180－188.

［2］Yuan C Q，Li J，Yan X P，et al.The Use of the Fractal Description to Characterize Engineering Surfaces and Wear Particles［J］.Wear，2003，255（1／6）：315－326.

［3］郭中領(lǐng)，符素華，張學(xué)會，等.土壤粒徑重量分布分形特征的無標度區(qū)間［J］.土壤通報，2010，41（3）：537－541.Guo Zhongling，F(xiàn)u Suhua，Zhang Xuehui，et al.Scale－free Domain of Fractal Characteristic of the Soil Particle－size Distribution［J］.Chinese Journal of Soil Science，2010，41（3）：537－541.

［4］Tang G J，Du B Q，Wang S L.Scaleless Band Automatic Identification for Fractal Fault Diagnosis of Rotor System［J］.J.Power Eng.，2009，29：440－444.

［5］秦海勤，徐可君，江龍平.分形理論應(yīng)用中無標度區(qū)自動識別方法［J］.機械工程學(xué)報，2006，42（12）：106－109.Qin Haiqin，Xu Kejun，Jiang Longping.Fractal Scaleless Band Automatic Identification for Fractal Theory Application［J］.Chinese Journal of Mechanical Engineering，2006，42（12）：106－109

［6］Fang Zhijun，Zhou Yuanhua，Zou Daowen，et al.Efficient Scheme for Determining Fractal Scaleless Range［J］.Journal of Infrared and Millimeter Waves，2004，23（5）：321－324.

［7］黨建武，施怡，黃建國.分形研究中無標度區(qū)的計算機識別［J］.計算機工程與應(yīng)用，2003，39（12）：25－27.Dang Jianwu，Shi Yi，Huang Jianguo.The Identification of Fractal Scaleless Band in the Study of Fractal with Computers［J］.Computer Engineering And Applications，2003，39（12）：25－27.

［8］巫兆聰.分形分析中的無標度區(qū)確定問題［J］.測繪學(xué)報，2002，31（3）：240－244.Wu Zhaochong.Determination of Fractal Scaleless Range［J］.Acta Geodaetica Et Cartographica Sinica，2002，31（3）：240－244.

［9］龐茂，吳瑞明，謝明祥.關(guān)聯(lián)維數(shù)快速算法及其在機械故障診斷中的應(yīng)用［J］.振動與沖擊，2010，29（12）：106－109.Pang Mao，Wu Ruiming，Xie Jingxiang.Improved Correlation Dimension Algorithm with Its Application to Mechanical Fault Diagnosis［J］.Journal of Vibration and Shock，2010，29（12）：106－109.

［10］欒海軍，汪小欽.多源遙感影像分形特征分析［J］.遙感信息，2010，3：7－12.Luan Haijun，Wang Xiaoqin.Fractal Features Analysis of Multi－source Remote Sensing Images［J］.Remote Sensing Information，2010，3：7－12.

［11］蔡金華，龍毅，毋河海，等.基于反S數(shù)學(xué)模型的地圖目標分形無標度區(qū)自動確定［J］.武漢大學(xué)學(xué)報（信息科學(xué)版），2004，29（3）：249－253.Cai Jinhua，Long Yi，Wu Hehai，et al.Automatic Determination of Fractal Non－scale Interval of Map Objects Based on Inverse‘S’Mathematical Model［J］.Geomatics and Information Science of Wuhan University，2004，29（3）：249－253.

［12］Yang H Y，Ye H，Wang G Z.Identification of Scaling Regime in Chaotic Correlation Dimension Calculation［C］／／3rd IEEE Conferenceon Industrial Electronics and Applications.ICIEA，2008：1383－1387.

［13］Ji C C，Zhu H，Jiang W.A Novel Method to Identify the Scaling Region for Chaotic Time Series Correlation Dimension Calculation［J］.Chinese Science Bull，2011，56：925－932.

［14］王成棟，凌丹，苗強.分形無標度區(qū)的一種自動識別方法［J］.計算機工程與應(yīng)用，2012，48（6）：9－12，27.Wang Chengdong，Ling Dan，Miao Qiang.Automatic Identification Method of Fractal Scaling Region［J］.Computer Engineering and Applications，2012，48（6）：9－12，27.

［15］Ajay K B，Jibitesh M.On Calculation of Fractaldimension of Images［J］.Pattern Recognition Letters，2001，22（6）：631－637.

［16］Kalyani S，Swarup K S.Particle Swarm Optimization Based K－means Clustering Approach for Security Assessment in Power Systems［J］.Expert Systems with Applications，2011，38（9）：10839－10846.

［17］Vasan A，Raju K S.Comparative Analysis of Simulated Annealing，Simulated Quenching and Genetic Algorithms for Optimal Reservoir Operation［J］.Applied Soft Computing，2009，9（1）：274－281.

［18］金元生.鐵譜技術(shù)及在磨損研究中的應(yīng)用［M］.北京：機械工業(yè)出版社，1991.