一種辨識磨粒群分形無標度區的新算法

張懷亮 鄒佰文 肖 雷

中南大學，長沙，410083

0 引言

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

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

1 磨粒群分形維數計算理論

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

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

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

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

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

圖2 無標度辨識中的局部最優

2 分形無標度區辨識新算法

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

2.1 模擬退火K－means算法

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

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

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

模擬退火K－means算法的實現步驟如下：

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

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

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

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

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

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

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

2.2 無標度區辨識新算法

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

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

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

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

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

圖3 磨粒群分形無標度區辨識過程

2.3 方法檢驗及結果分析

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

表1 測試圖形分形維數計算結果

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

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

3 磨粒群的分形無標度區辨識

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

表2 應用新方法計算磨粒群分形維數結果

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

4 結論

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

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

［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］郭中領，符素華，張學會，等.土壤粒徑重量分布分形特征的無標度區間［J］.土壤通報，2010，41（3）：537－541.Guo Zhongling，Fu 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］秦海勤，徐可君，江龍平.分形理論應用中無標度區自動識別方法［J］.機械工程學報，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］黨建武，施怡，黃建國.分形研究中無標度區的計算機識別［J］.計算機工程與應用，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］巫兆聰.分形分析中的無標度區確定問題［J］.測繪學報，2002，31（3）：240－244.Wu Zhaochong.Determination of Fractal Scaleless Range［J］.Acta Geodaetica Et Cartographica Sinica，2002，31（3）：240－244.

［9］龐茂，吳瑞明，謝明祥.關聯維數快速算法及其在機械故障診斷中的應用［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數學模型的地圖目標分形無標度區自動確定［J］.武漢大學學報（信息科學版），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］王成棟，凌丹，苗強.分形無標度區的一種自動識別方法［J］.計算機工程與應用，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］金元生.鐵譜技術及在磨損研究中的應用［M］.北京：機械工業出版社，1991.