郭鵬程，孫龍剛，李輝，袁江霞

(西安理工大學(xué) 水利水電學(xué)院，陜西 西安 710048)

分形一詞來源于拉丁文Fractus，是用來描述非常不規(guī)則以至不能視為經(jīng)典幾何研究對象的物體。

對分形集一般基于它的屬性來定義，具有精細(xì)、不規(guī)則、自相似、遞歸的性質(zhì)[1]。一般情況下，具有分形特征的圖形，其分形維數(shù)高于其拓?fù)渚S數(shù)。

分形理論在包括信號處理在內(nèi)的許多領(lǐng)域都得到了廣泛的應(yīng)用。研究過程發(fā)現(xiàn)，旋轉(zhuǎn)機(jī)械的振動(dòng)信號是一種典型的分形信號[2]，因此可將分形理論應(yīng)用到旋轉(zhuǎn)機(jī)械故障診斷中來。

單一維數(shù)由于對盒子尺度的大小依賴性較強(qiáng)，對信號的局部特征描述性不夠，而多重分形可以有效的彌補(bǔ)這一問題。

在用分形維數(shù)對分形體的分類研究中，出現(xiàn)這樣的現(xiàn)象：具有不同紋理表面的兩個(gè)分形體可能有相近或相同的分形維數(shù)。

基于這一事實(shí)，在如何提高這一領(lǐng)域的分形幾何的算法方面人們做了大量工作。Evertsz和Mandelbrot于1992年提出，在描述一個(gè)分形集的幾何特性方面，多重分形維數(shù)可能是一個(gè)比分形維數(shù)更合適的參數(shù)[3]。對于多重分形的研究表明，多重分形是一個(gè)比分形維數(shù)(盒維數(shù))更適用的參數(shù)。文獻(xiàn)[4] 基于數(shù)學(xué)形態(tài)學(xué)的分維數(shù)計(jì)算方法，采用形態(tài)學(xué)操作計(jì)算分形維數(shù)，提出了一種快速有效判斷滾動(dòng)軸承故障狀態(tài)的新方法，對比分析結(jié)果表明該方法計(jì)算速度快，計(jì)算結(jié)果準(zhǔn)確穩(wěn)定。文獻(xiàn)[5] 將多重分形方法應(yīng)用于故障信號診斷中，并進(jìn)行了實(shí)例測試和分析，證明了該方法在實(shí)際應(yīng)用中的可行性；文獻(xiàn)[6] 提出了廣義維數(shù)最小二乘法的計(jì)算公式，并對轉(zhuǎn)子系統(tǒng)實(shí)測時(shí)域故障信號進(jìn)行了廣義維數(shù)計(jì)算，結(jié)果顯示各分形維數(shù)能較好地診斷和識別故障狀態(tài)。文獻(xiàn)[7] 采用基于多分辨率分析的小波系數(shù)去噪方法，通過在頻域上對各層的小波分解系數(shù)的畸變量進(jìn)行抑制以實(shí)現(xiàn)信號去噪，并通過試驗(yàn)驗(yàn)證了算法的有效性。文獻(xiàn)[8] 提出了一種改進(jìn)的多重分形維數(shù)算法，改變了傳統(tǒng)多重分形維數(shù)對q維特征進(jìn)行累加的計(jì)算方法，在保證算法計(jì)算復(fù)雜度基本不變的情況下增加了信號特征的規(guī)律性和類內(nèi)聚集度。

本文介紹了基本的多重分形維數(shù)的計(jì)算方法，研究了影響其計(jì)算準(zhǔn)確性的因素，發(fā)現(xiàn)在實(shí)際的計(jì)算中，除了盒子的尺度，盒子不同的運(yùn)動(dòng)方向也會(huì)導(dǎo)致多重分形維數(shù)的不穩(wěn)定性。針對此問題，本文展開了有關(guān)的研究，提出一種改進(jìn)的多重分形算法并進(jìn)行了仿真計(jì)算。

1 盒子覆蓋法

標(biāo)準(zhǔn)的盒子計(jì)算方法一般用來分析點(diǎn)集，每一組點(diǎn)集用廣義維數(shù)或多重分形維數(shù)來描述。

廣義維數(shù)通過改變一系列的值來定義函數(shù)進(jìn)行計(jì)算，然后通過變換計(jì)算。

盒子計(jì)算法在1989年被提出，由Vicsek在1990年完善了其應(yīng)用。其計(jì)算方法如下。

數(shù)學(xué)上，多重分形是根據(jù)其具有的分區(qū)功能來定義的：把所研究的對象分為N個(gè)小區(qū)域，分形體生長界面在第i個(gè)區(qū)域的生長概率為Pi(ε)，以盒維數(shù)計(jì)算為基礎(chǔ)，當(dāng)盒子尺寸為ε時(shí)，定義Si(ε)為第i個(gè)小盒子內(nèi)所有像素點(diǎn)灰度的數(shù)值之和，則第i個(gè)小盒子的平均灰度值可表示為：

(1)

不同的小區(qū)域生長的概率不同，用不同的標(biāo)度指數(shù)αi表示。

(2)

在多重分形中，αi用來表示分形體小區(qū)域的分維，因?yàn)檫@些小區(qū)域數(shù)目趨于無窮，因此可以用f(αi)表示這個(gè)由不同αi組成的無窮序列構(gòu)成的譜。f(α)和α是描述多重分形的一組參數(shù)。f(α)又被稱為奇異譜。由于直接計(jì)算f(αi)和αi比較困難，因此一般選用另一組參數(shù)D(q)和q來描述。

首先定義一個(gè)配分函數(shù)Xq(ε)，對概率Pi(ε)用q次方進(jìn)行加權(quán)求和:

(3)

式中q為權(quán)重因子，τ(q)為質(zhì)量指數(shù)，其表達(dá)式為：

(4)

定義q次信息維D(q)為:

(5)

兩套參數(shù)α,f(α)和q,D(q)之間可以用Legendre法則進(jìn)行變換：

(6)

或：

f(α)=qα-τ(q)

(7)

其中:

τ(q)=(q-1)D(q)

(8)

由式(6)可知，若已知α與其譜f(α)，即可求出D(q)。反之，若測得Pi可得出D(q)。

式(8)兩邊同時(shí)對q求導(dǎo)，得：

(9)

根據(jù)式(7)可求出f(α)。

這種方法的優(yōu)點(diǎn)是盒子一般在結(jié)構(gòu)的中心，因此，即使有的盒子包含的信息很少，也不會(huì)溢出邊界，當(dāng)q<0時(shí)，包含信息少的盒子對結(jié)果也有很大的貢獻(xiàn)。

盒子覆蓋方法能夠解決邊界問題，當(dāng)q<0時(shí)，允許重構(gòu)多重分形譜，在計(jì)算時(shí)，需要二進(jìn)制圖像。運(yùn)用小波模極大值計(jì)算有一定的限制，雖然其在選擇尺度方面有一定的自由度，但計(jì)算比盒子覆蓋法難以實(shí)現(xiàn)。

因此，本文在編制程序過程中，選用盒子覆蓋法進(jìn)行計(jì)算。

2 多重分形維數(shù)值的影響因素及其改進(jìn)算法

2.1 尺度ε的取值

如上所述，選用盒子覆蓋法計(jì)算多重分形譜，在計(jì)算過程中主要考慮計(jì)算概率測度時(shí)的盒子大小ε和配分函數(shù)中的q。

由式(1)可知，概率測度的大小取決于盒子大小ε，即ε的變化是影響Pi(ε)的主要因素。

概率測度基礎(chǔ)是盒維數(shù)，其計(jì)算的思想與盒維數(shù)類似：用直徑為ε的正方形盒子去覆蓋一個(gè)分形體(信號)，選擇不同的ε會(huì)得到一系列的Si(ε)，兩者之間用式(10)表示的線性關(guān)系就是盒維數(shù)。

(10)

式中，D表示盒維數(shù)。

由上式可知，在盒維數(shù)的計(jì)算過程中，影響其結(jié)果的因素主要是盒子大小ε和Si(ε)，而ε的大小決定Si(ε)的大小。這與多重分形中計(jì)算概率測度的性質(zhì)一樣。

在一般的盒維數(shù)計(jì)算中，一般只根據(jù)ε大小與Si(ε)的關(guān)系，并對一系列ε和Si(ε)線性關(guān)系進(jìn)行判別，在實(shí)際的應(yīng)用中筆者發(fā)現(xiàn)，即使對同一個(gè)分形體，在ε大小相同的情況下，得到的Si(ε)也可能是不同的。

如圖1所示，盒子初始的位置左右相反，由相同的ε得出的Si(ε)分別是12和14。即在盒子的直徑和分形體相同的情況下，由于盒子移動(dòng)的初始位置不同，覆蓋分形體所需要的盒子總數(shù)Si(ε)也是不同的。

因此，盒維數(shù)的計(jì)算結(jié)果也取決于盒子的相對位置，這個(gè)是多重分形譜計(jì)算過程中同樣需要考慮的。

圖1 相同的ε對應(yīng)不同的Si(ε)

針對此問題，本文采用如圖2所示算法流程，首先計(jì)算盒維數(shù)，在盒子ε變化的情況下，同時(shí)考慮盒子位置的變化，定義1～12個(gè)初始位置，分別計(jì)算其盒維數(shù)。

對得到的數(shù)值進(jìn)行平均，最終以得出的均值維數(shù)作為實(shí)際的結(jié)果。

圖2 算法流程

2.2 權(quán)重因子q的選擇

基于統(tǒng)計(jì)物理方法計(jì)算的多重分形譜的配分函數(shù)Xq(ε)可以表示為式(3)形式。由式(3)可知，q值的不同表征了不同子集在配分函數(shù)中的重要作用。若q→+∞，則最大概率起決定作用，此時(shí)Xq(ε)反映了概率最大的子集的性質(zhì)；相反地，若q→-∞，則表明最小概率起決定作用，即Xq(ε)反映了概率最小的子集的有關(guān)性質(zhì)。多重分形就是以不同的q值作為邊界對象，從而將分形體分成具有不同層次的區(qū)域來研究。

在理論上，q取值越大越好(-∞

圖3(a)和(b)分別畫出了生成元為0.4/0/0.6和0.498/0/0.502的Cantor集的多重分型譜，圖中的實(shí)線和點(diǎn)分別表示解析法計(jì)算的f(α)曲線和配分函數(shù)得出的f(α)曲線，標(biāo)注為相應(yīng)的q值。由圖可以看出，隨著|q|的增加，α和f(α)的值逐漸接近理論極限值。

圖3(c)和(d)分別畫出了0.4/0/0.6的Cantor集和0.2/0.6/0.2的Cantor集的α隨q的變化，可見隨著|q|增大到一定程度時(shí)，α基本上不再隨|q|的增加而變化。

圖3 Cantor集的多重分形譜

2.3 無標(biāo)度區(qū)間的確定

分形理論基于自相似原理，在任何標(biāo)度下系統(tǒng)局部均相似于整體，這里的標(biāo)度表示廣義上的尺度。然而對于實(shí)際的情況，這一相似特性不可能在任何區(qū)間都成立，而只能在某一特定的區(qū)間內(nèi)成立，即所謂的無標(biāo)度區(qū)間。

2.3.1 相空間重構(gòu)

實(shí)現(xiàn)自動(dòng)確定的第一步就是由觀測序列生成m維向量集合，該過程稱為相空間重構(gòu)。此時(shí)m表示嵌入維數(shù)。

由嵌入定理可知，若嵌入維數(shù)m≥2D+1，那么在重構(gòu)空間Rm中系統(tǒng)結(jié)構(gòu)同構(gòu)于原空間。具體實(shí)現(xiàn)過程如下：假設(shè)已獲得某一原始的觀測序列{xi,i=1,2,…,N}，選取合適的嵌入維數(shù)m和時(shí)間延遲τ進(jìn)行相空間重構(gòu)，將會(huì)得到空間中的一個(gè)點(diǎn)集Jm，該點(diǎn)集中的元素表示為：

Xi=[xi,xi+k,xi+2k, …,xi+(m-1)k]

i=1,2,…,Nm

(11)

式中，Nm=N-(m-1)τ，為點(diǎn)集Jm中向量的個(gè)數(shù)；k是時(shí)間延遲數(shù)，若Δt為采樣間隔，則τ=k·Δt，對離散時(shí)間序列而言，k也就是尺度r。

2.3.2 無標(biāo)度區(qū)間的自動(dòng)選擇

從點(diǎn)集Jm中任選一個(gè)任意元素Xi，然后計(jì)算Xi到其余的Nm-1個(gè)元素的歐式距離，表示為：

(12)

對所有的Xi重復(fù)這一過程，其中，當(dāng)j=i時(shí)，rij=0。然后在相空間中計(jì)算q次關(guān)聯(lián)積分，由于式(13)無法直接計(jì)算，重新定義q次關(guān)聯(lián)積分如下：

(13)

式中，H(r-rij)是Heaviside階躍函數(shù)，并且有如下關(guān)系：

(14)

當(dāng)嵌入維數(shù)m足夠大時(shí)，可以用q次關(guān)聯(lián)積分得到廣義分形維數(shù)Dq。

(15)

在各種尺度ε下按式(13)計(jì)算可得到Xq(ε)，在對數(shù)坐標(biāo)中做一元線性回歸分析，擬合lnXq(ε)～lnε，所得直線的斜率即為廣義分形維數(shù)Dq。

在實(shí)際計(jì)算中，標(biāo)度ε滿足的最大范圍應(yīng)該在兩點(diǎn)間距的最小值與最大值之間，即:

min{‖Xi-Xj‖}≤ε≤max{‖Xi-Xj‖}

i,j∈Nm,且i≠j

(16)

當(dāng)標(biāo)度ε較小時(shí)，關(guān)聯(lián)積分Xq(ε)趨于零并且與后面的關(guān)聯(lián)積分值有較低的線性相關(guān)程度，這樣的ε值應(yīng)該舍去；反之，如果標(biāo)度ε很大，關(guān)聯(lián)積分Xq(ε)將會(huì)達(dá)到飽和，會(huì)超出無標(biāo)度區(qū)間的最大范圍，這樣的ε值也應(yīng)該舍去。

使關(guān)聯(lián)積分Xq(ε)的值飽和或者接近于零的ε值稱為野值。

記Un=PilnPi，舍掉使NUM小于2的所有ε值。NUM表達(dá)式如下：

(17)

具體算法流程如圖4所示。

圖4 無標(biāo)度區(qū)間算法流程圖

從圖4可看出，確定無標(biāo)度區(qū)間的基本過程如下：首先進(jìn)行q次關(guān)聯(lián)積分的計(jì)算，以q是否為1為界，若q=1，則直接進(jìn)行關(guān)聯(lián)積分的計(jì)算并確定尺度ε的取值范圍，而當(dāng)q≠1時(shí)，按照給定的尺度范圍計(jì)算出中間變量，隨后兩種情況均對野值進(jìn)行剔除操作，尋找相應(yīng)區(qū)間的上下限或計(jì)算出相關(guān)系數(shù)，最終確定無標(biāo)度區(qū)間的精確范圍。

3 仿真結(jié)果與分析

為了驗(yàn)證改進(jìn)多重分形算法的有效性，本文選用Weierstrass函數(shù)進(jìn)行驗(yàn)證仿真分析。

Weierstrass函數(shù)是數(shù)學(xué)上一個(gè)非常重要的函數(shù)，該函數(shù)具有處處連續(xù)但處處不可微的性質(zhì)，是德國科學(xué)家Karl Weierstrass在1872年正式創(chuàng)造發(fā)表的。該函數(shù)的數(shù)學(xué)表達(dá)式為：

(18)

式中，0<α<1,λ>1。

Weierstrass函數(shù)具有分形性質(zhì)，任何局部的任意放大都與整體相似。其在數(shù)學(xué)上有嚴(yán)格意義上的分形盒維數(shù)，數(shù)學(xué)上對其證明的結(jié)果為：

DB=2-α

(19)

設(shè)a=λ-α=0.5，b=λ=3時(shí)，則函數(shù)可表達(dá)為：

(20)

這里的盒維數(shù)理論值為:

DB=2-α=2-0.631=1.369

在多重分形計(jì)算中，廣義維數(shù)譜中D0即為盒維數(shù)。

圖5為上式Weierstrass函數(shù)的曲線，選定參數(shù)如下：盒子尺度ε∈[1,45]，權(quán)重因子q∈[-20,20]，運(yùn)用初始不考慮其盒子的相對位置算法得出的維數(shù)值D0=1.209 7。對其用改進(jìn)的多重分形方法進(jìn)行計(jì)算，初始參數(shù)不變，設(shè)定盒子的初始位置為4個(gè)。圖6給出了維數(shù)計(jì)算的結(jié)果。表1為改進(jìn)前后兩種算法計(jì)算結(jié)果對比。

圖5 Weierstrass 函數(shù)曲線

圖6 廣義維數(shù)譜

表1 結(jié)果對比

從表1可看出，用改進(jìn)后的方法計(jì)算得到的結(jié)果其準(zhǔn)確性更高。但是估計(jì)值始終與理論值存在差異，原因是盒維數(shù)實(shí)際上是點(diǎn)集(lnε, lnN(ε))的斜率,而擬合直線始終會(huì)存在著斜距。

4 結(jié) 論

本文對多重分形維數(shù)的計(jì)算方法作了基本介紹，并對影響其計(jì)算的因素-盒子尺度 和權(quán)重因子 進(jìn)行了研究和分析。研究發(fā)現(xiàn)，除了盒子尺度 ，盒子的運(yùn)動(dòng)方向也對結(jié)果影響明顯。因而，在編制計(jì)算程序時(shí)，加入了此影響因素并對算法進(jìn)行了改進(jìn)。實(shí)例結(jié)果表明：考慮盒子的相對位置得到的分形維數(shù)值更為準(zhǔn)確和穩(wěn)定。

參考文獻(xiàn):

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

[2] 徐玉秀, 鐘建軍, 劉薇, 等. 基于廣義分形的旋轉(zhuǎn)機(jī)械

故障診斷識別與分類[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào), 2005, 24(4):591-594.

Xu Yuxiu, Zhong Jianjun, Liu Wei, et al. Application of multifactal theory to couping fault diagnosis[J]. Journal of Liaoning Technixal University, 2005, 24(4):591-594.

[3] Jelinek H, Karperien A, Cornforth D, et al. Micromod - an l-systems approach to neural modelling[C]// Sixth Australia-Japan joint workshop on intelligent and evolutionary systems. Canberra, Australia, 2002.

[4] 李兵, 張培林, 任國全, 等, 基于數(shù)學(xué)形態(tài)學(xué)的分形維數(shù)計(jì)算及在軸承故障診斷中的應(yīng)用 [J]. 振動(dòng)與沖擊, 2010, 29(5): 191-194.

Li bing, Zhang Peilin, Ren Guoquan, et al. Calculation of the fractal dimension based on mathematical morphology and application in bearing fault diagnosis[J]. Journal of Vibrarion and Shock, 2010, 29(5): 191-194.

[5] Chunsheng E, Jiang Dongxiang, Ni Weidou. Fractal geometry and its applications in vibration failure symptom abstraction of the steam turbine [J]. Turbine Technology,1999,41(6):321-324.

[6] 徐玉秀, 鐘建軍, 聞邦椿. 旋轉(zhuǎn)機(jī)械動(dòng)態(tài)特性的分形特征及故障診斷[J]. 機(jī)械工程學(xué)報(bào), 2005, 41(12): 186-189.

Xu Yuxiu, Zhong Jianjun, Wen Bangchun. Fractal fault diagnosis and classification to modal characteristic of rotor system[J]. Chinese Journal of Mechanical Engineering, 2005, 41(12): 186-189.

[7] 王振朝, 岳瑩昭, 師潔, 等. 基于多分辨率分析的小波系數(shù)壓擴(kuò)去噪算法[J].中國電機(jī)工程學(xué)報(bào), 2008, 28(10): 76-81.

Wang Zhenchao, Yue Yingzhao, Shi Jie, et al. Wavelet cofficients companding ce-noising algorithm based on multiresolution analysis[J]. Proceeding of the CSEE, 2008, 28(10): 76-81.

[8] 寧宇. 一種簡化的多重分形維數(shù)算法[J]. 價(jià)值工程, 2014, (9):181-182.

Ning Yu. A simplified multi-fractal dimension algorithm[J]. Valve Engineering, 2014,(9):181-182.