基于聚類分析理論的磨粒群無標度區自動辨識

0無標度區的辨識在機器的操作中產生的磨粒群包含豐富的摩擦磨損信息，可以提取磨粒群的特征，可用于設備的狀態監測和故障排除。研究表明,磨粒群在一定尺度上表現出分形特性,傳統基于歐氏幾何的參數無法準確描述磨粒群特征［１－２］。分形維數是表征磨粒群分形特性的一個重要參數,因而提高磨粒群分形維數的計算精度,對于磨粒群特征的提取和設備磨損狀態監測有重要意義。由于產生機制復雜,磨粒群的分形自相似特性是近似或統計意義上的,即在一定尺度范圍內的自相似,這個尺度范圍即稱為無標度區,在無標度區外討論磨粒群的分形特性沒有意義,因此,無標度區的準確辨識是磨粒群分形維數計算中必須解決的關鍵問題之一。辨識無標度區的常用方法有以下幾類:1基于r~N(r)(尺度~測度)雙對數曲線的人工識別法［３］;2基于r~N(r)點對相關性檢驗的排序逼近法［４－８］;3基于三折線、曲線-直線-曲線、反S等模型的模型擬合法［９－１１］。人工識別法主觀性強,排序逼近法在r~N(r)點對較多時計算量很大且較易出現局部最優,模型擬合法在模型的選取方面缺乏理論依據,應用時有一定局限性。因而,Yang等［１２］利用無標度區一階導數近似不變的特點提出應用聚類分析理論辨識無標度區的算法;Ji等［１３］在此基礎上結合K-means和點斜率誤差法提高了算法的精度;王成棟等［１４］根據無標度區二階導數在0附近波動的特點提出了一種自動辨識方法,這類利用r~N(r)雙對數曲線導數特點的辨識算法幾何意義明確,計算量較小,但大多對初始值敏感,容易出現局部最優,可能導致辨識出來的無標度區雖然相關性檢驗參數符合要求,但并不是完整的無標度區,甚至不在真實無標度區內,即無標度區的局部最優。本文結合模擬退火和K-means算法對雙對數曲線的一階、二階導數進行全局聚類,可有效避免陷入局部最優,從而提高磨粒群分形維數的計算精度。1圖像分形特征表征磨粒群分散于設備的潤滑油中,通過鐵譜儀可將潤滑油樣中的磨粒群分離出來,制成鐵譜片,經鐵譜顯微鏡拍照成像,得到磨粒群的圖像(圖1a),原始圖像背景復雜,通過背景去噪、閾值分割等圖像預處理過程,得到磨粒群的二值圖像(圖1b),磨粒群分形維數的計算過程基于其二值圖像。實際應用中,為保證狀態監測數據的可比性,獲取磨粒群圖片時應遵守同樣的取樣倍數及視場選取規范。分形維數的定義包括相似維、Hausdorff維、信息維、計盒維數等多種,其中,計盒維數(box-countingdimension)由于其物理意義明確,便于程序化計算等優點,廣泛應用于圖像分形特征的表征。計盒維數的計算原理如下:將邊長為r的方形網格覆蓋所需測量的圖形,測算與圖形重疊的網格數目為測度N(r),不斷改變網格邊長,重復上述過程,得到一組不同尺度r下的測度N(r),那么計盒維數即為根據計盒維數的定義,磨粒群分形維數計算的一般過程如下:假設磨粒群圖像的分辨率為M×M,以邊長為r個像素的網格不斷覆蓋磨粒群圖像并統計重疊網格數目N(r),其中r的取值范圍為{r|1≤r<M/2,r=1,2,4,…,2n}［１５］,根據r~N(r)點對序列繪出雙對數曲線,將曲線中的直線部分作為無標度區間,在無標度區間內對r~N(r)點對進行線性回歸,取分形維數D=-α,其中α為回歸直線的斜率。在磨粒群分形維數的計算過程中,無標度區的辨識結果直接影響磨粒群分形維數的計算精度,辨識出來的無標度區若出現偏離、過短或過長等誤差,會導致無法正確評估磨粒群的分形維數。而應用已有的無標度區辨識方法計算磨粒群的分形維數時,往往會出現局部最優的現象,如圖2所示。圖2中的數據點為磨粒群圖像(圖1)應用計盒維數算法得到的r~N(r)點對,A、B、C、D4個區為不同的無標度區辨識結果,其中,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,差異很大,所以,無標度區的準確辨識是磨粒群分形維數計算的關鍵環節。2k-melus聚類算法的辨識為避免無標度區辨識算法陷入局部最優,將K-means聚類算法與全局搜索算法模擬退火算法相結合,從r~N(r)雙對數曲線中準確辨識磨粒群無標度區。2.1模擬退火算法K-means算法是一種基于距離的聚類算法,該算法以距離為相似性指標,將觀察對象分為k類［１６］。基本過程如下:首先從研究對象中隨機選擇k個對象作為聚類中心,其余的對象根據它們到聚類中心的距離劃分到最近的聚類。然后,以各聚類中所有對象的平均值作為新的聚類中心,重復上述劃分過程,直到準則函數收斂。應用K-means算法對雙對數曲線的導數進行聚類,可分離出無標度區,但是從K-means的聚類過程可以看出,算法對初始聚類中心的選擇依賴很大,初始聚類中心選擇不當,可能出現局部最優,將它應用到磨粒群無標度區的識別可能會出現無標度區過短等現象。模擬退火算法是由Kirkpatrick提出的一種啟發式、隨機優化算法［１７］。該算法與初始值無關,算法求得的解與初始解狀態無關,具有漸近收斂性,已在理論上被證明是一種以概率1收斂于全局最優解的全局優化算法。因而,利用模擬退火算法能夠彌補K-means算法依賴初始聚類中心的缺點,跳出局部最優。模擬退火K-means算法以K-means聚類的結果作為模擬退火算法的初始解,通過隨機改變聚類中某個元素的類別產生新解,目標函數為當前聚類劃分的總類間離散度:其中,X為研究對象,w為聚類劃分,為第i個聚類的中心,為對象到對應聚類中心距離,聚類準則函數Jｗ為各樣本到對應聚類中心距離的總和。模擬退火K-means算法的實現步驟如下:(1)將數據樣本的K-means聚類結果作為初始解w,根據式(2)計算目標函數值Jｗ。(2)初始化溫度t０,t０=Jｗ,初始化退火速度a和最大退火次數。(3)在某一溫度t,迭代進行步驟(4)~步驟(6),直到最大迭代次數跳到步驟(7)。(4)通過隨機改變一個聚類樣本的當前所屬類別產生新解w′,計算新的目標函數值J′ｗ。(5)判斷J′ｗ是否為最優目標函數值,若是,則保存聚類劃分w′為最優聚類劃分、J′ｗ為最優目標函數值;否則跳到下一步。(6)計算新的目標函數值與當前目標函數值的差ΔJ。判斷ΔJ是否小于0:若ΔJ<0,則接受新解,即將新解作為當前解;若ΔJ≥0,則根據Metropolis準則,以概率p(p=eΔＪ／（ｋｔ）)接受新解。(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=-α,作為磨粒群的分形維數。2.3圖形無標度區的識別效果為檢驗新方法對無標度區識別的準確性及精度,選取Koch曲線、Sierpinski三角形、Sierpinski方毯三個已知理論分形維數的圖形作為測試對象,三個圖形的大小一致,采用計盒法獲取三個圖形的r~N(r)雙對數曲線,利用新方法識別相應的無標度區,對無標度區內的點進行最小二乘擬合,確定圖形的分形維數。同時,為了對比研究本文無標度區識別方法的精度,同時用人工識別法和文獻提出的二階導數法對三個圖形無標度區進行了識別,計算結果如表1所示,圖4示出了三種方法計算結果與理論值的誤差。應用新方法辨識三個圖形無標度區的過程中均出現局部最優現象,圖4表明,人工識別法在三個圖形的無標度區識別中都取得了較好的效果,平均誤差在3%~5%之間,但是人工識別法需與計算機進行交互操作,無法進行自動識別,且不同觀察者一般得到不同的結論,主觀性強;二階導數法對Koch曲線的識別效果很好,但對其他兩個圖形的無標度區識別效果較差;本文提出的識別方法雖對Koch曲線的識別效果與人工識別法效果相當,但對Sierpinski三角形和Sierpinski方毯的識別誤差均在1%以下,由于Sierpinski三角形和Sierpinski方毯的特點與磨粒群的特點相似,都是一種多孔結構,因而應用本文提出的方法對磨粒群分形無標度區的識別是合理的,分形維數的計算結果是可信的。3磨粒群分形維數為分析新方法對磨粒群無標度區的辨識能力,應用新方法計算了正常磨損、切削磨損、滾動疲勞磨損、滾動和滑動復合磨損及嚴重滑動磨損的磨粒群分形維數,磨粒群圖片來自文獻,圖片大小為400像素×400像素,盒子尺寸增長方式為{r|1≤r<M/2,r=1,2,4,…,2n},計算程序基于MATLAB開發,計算機CPU主頻為3.00GHz,計算結果如表2所示,其中,無標度區中的數值范圍表示r~N(r)雙對數曲線中識別的無標度區取點范圍。計算結果表明,在圖片大小為400像素×400像素,盒子尺寸增長方式為1,2,4,…,2n時,磨粒群的分形無標度區大約為1~20個點之間,取此區間內的點對進行擬合,可以得到較精確的磨粒群分形維數,這段區間內,點對的線性相關系數大于0.99,同時區間長度占整個雙對數曲線的50%左右,未出現無標度區辨識的局部最優。另外,試驗中不同磨損狀態的磨粒群分形維數差異較大,從1.44到1.71,這表明磨粒群的分形維數可作為摩擦副磨損狀態識別的一個判據。新算法對給定條件下磨粒群無標度區進行識別所耗用