摘 要:基于序列分形自仿射特性,提出一種實現一維信號分形維數估計的方法。按不同尺度將信號序列分段為映射區間和象區間,采用搜索算法確定與各象區間最優匹配的映射區間,并根據迭代函數系統理論估計信號的分形維數。以分形維數已知的Mackey-Glass和Lorenz信號為例,仿真表明提出的方法能準確估計信號的分形維數,對實際應用具有一定的參考價值。
關鍵詞:分形維數; 自仿射; 迭代函數系統; 搜索算法
中圖分類號:TN971+1 文獻標志碼:A 文章編號:1001-3695(2008)08-2521-03
Research and implementation of fractal dimension estimation
HE Taoa, ZHOU Zheng-oub
(a.Institute of Astronautics Aeronautics, b.School of Electronic Engineering, University of Electronic Science Technology of China, Chengdu 610054, China)
Abstract:This paper proposed a method to implement fractal dimension estimation based on the self-affine property of signals. It decomposed the signal series to mapping and mapped segments, and adopted a searching algorithm to find the mapping sections that match best for each mapped sect. It estimated fractal dimensionby the iterative function system. Simulations on Mackey-Glass time series and Lorenz system indicate that the proposed method can estimate the fractal dimension accurately.
Key words:fractal dimension; self-affine; iterative function system; searching algorithm
自從1975年被正式提出,分形理論在數學、物理、化學、生物、經濟、管理、社會等學科[1,2]都得到了廣泛的應用。在信號分析和處理領域,分形研究對象包括各種常見的噪聲與雜波干擾、心電信號、地震波、語音和圖像信號等[3,4]。信號分形分析可以對這些信號進行建模、識別和分割處理。分形維數能反映信號波形的幾何特征信息,對信號的復雜度、不規則性和全局正則性進行定量描述。準確、穩定地估計分形對象的分形維數在信號的分形分析中是十分關鍵的。
根據分形對象的不同性質,研究者提出了各種分形維數估計方法,如計盒法[5]是根據覆蓋分形曲線所需的尺度數與測度之間的指數關系,而功率譜法[6]則是基于分形曲線功率譜具有的指數規律。借助數學形態覆蓋法,Maragos 等人[7]用腐蝕和膨脹算子提出了面積法。Frandrin由分數布朗運動的小波系數方差滿足冪定律的性質出發導出了分形維數的小波分析法[8]。文獻[9]從分數布朗運動特性出發提出了基于離散時間序列的有限長度方差法。這些估計方法在準確度、穩定性、復雜度和對信號長度及信噪比的魯棒性方面具有不同的性能。
分形最基本的特征是自相似性,即局部可以看做是整體的一個小復制品,而分形對象可以看做是通過一系列仿射變換得到的小復制品拼貼而構成。在序列的分形自仿射特性基礎上,本文利用迭代函數系統,根據吸引子定理和拼貼定理建立一維信號序列的自仿射模型,并提出一種由模型參數估計信號分形維數的方法。
1 分形自仿射與迭代函數系統[10]
目前對分形沒有一個嚴格的數學定義,人們通常從分形的特性來對其進行研究。分形是一類復雜性很高、沒有特征長度,但具有一定意義的自相似的圖形和結構的總稱。由于現實世界中的對象通常只是在某些方面表現出分形特性,更多的是考慮研究對象的近似自相似性。一維信號的自仿射分形表現為其結構在新的長度重標下性質不變。將序列某一區間縮小一個因子b-1,接著用一個等于b-H的不同因子重標序列,該區間將重標為其自身,即F(X)≈b-HF(bx)。其中:H>0為某常指數。
迭代函數系統(IFS)能夠僅用少數幾個仿射變換產生復雜的函數和圖形,而形態各異的分形對象可以用IFS進行表示和描述。IFS是建立在壓縮映射基礎上的。
定義1 設(X,d)為一距離空間,一個變換ω:X→X稱為壓縮映射,即存在一個正的常數c<1,使d(ω(x),ω(y))≤c×d(x,y),x,y∈X。其中:c為壓縮因子。
定義2 一個雙曲(hyperbolic)迭代函數系統(IFS)是由一個完備距離空間(X,d)及其上的一組壓縮映射ωi:XX組成,ωi的壓縮因子為ci,且0≤ci<1(i=1,2…,N)。系統表示成{X;ωi,i=1,2,…,N},且系統的壓縮因子為c=max{(ci;i=1,2,…,N)}。
定義3 設(X,d)為一完備距離空間,空間H(X)的元素由X上非空集組成。A,B∈H(X),A、B之間的Hausdorff距離定義為h(A,B)=d(A,B)∨d(B,A)。其中:d(A,B)表示A到B的距離,d(A,B)=max{d(x,B),x∈A},且d(x,B)=min{d(x,y),y∈B}。(H(X),h(d))即Hausdorff度量空間。
吸引子定理給出了IFS吸引子存在的條件:
定理1 設{X;ωi;i=1,2,…,N}為一雙曲IFS,壓縮因子為c,則在(H(X),h(d))上定義的映射ω=∪Ni=1ωi對所有B∈H(X)也是一個壓縮映射,并有壓縮因子c,即對任意的B,C∈H(X),總存在h(ω(B),ω(C))≤c×h(B,C)。且在H(X)上必存在一個惟一吸引子A,B∈H(X),有A=ω(A)=limn→∞ωn(B)。
對于時間序列函數y=f(x)來說,其函數圖A={x,f(x)}是二維的。假定A的一個自仿射劃分為A=∪Mi=1Ai。對每一個子集Ai,均存在一個自仿射變換ωi使Ai=ωi(A)以概率pi(pi>0 且Mi=1 pi=1)成立。ωi的集合ω=∪Mi=1ωi構成一個IFS:{X;ωi,pi,i=1,2,…,N} 。子集Ai相鄰但互不重疊。
信號序列映射區間分別通過IFS的各個自仿射變換映射到對應的象區間,將這些象區間拼貼在一起,所得結果與原序列的誤差描述了IFS產生的吸引子與原序列的相似程度;誤差越小,兩者相似程度越高。拼貼定理給出了誤差的度量:
定理2 設(X,d)為一完備距離空間,給定L∈H(X)和ε≥0,選取一個壓縮因子為c(0≤c<1)的IFS{X;ωi,i=1,2,…,N},使得h(L,∪Ni=1ωi(L))≤ε,則有h(L,A)≤ε/ (1-c)。其中:A為IFS的吸引子;h(·,·)為Hausdorff距離。
上面的敘述為基于序列分形自仿射的分形維數估計提供了理論依據。設IFS由M個自仿射變換構成,每個仿射變換的壓縮因子為di,則序列分形維數D滿足:
D=1+log Mi=1|di| / log M(1)
顯然,對于給定的序列,只要確定相應的仿射變換收縮因子di,即可由式(1)決定的IFS分形維數給出序列分形維數的估計值。
2 算法實現及分析
給定R2中的信號序列為{(un,vn):n=0,1,…,N;un<un+1},{(xi,yi):i=0,1,…,M;M≤N}是序列的一個子集。序列自仿射的IFS由如下形式的M個仿射變換ωi構成:
ωix
y=ai 0
ci dix
y+ei
fi=r
r;i=1,2,…,M(2)
且滿足ωixi,1
yi,1=xi-1
yi-1,ωixi,2
yi,2=xi
yi和(xi,2-xi,1)>(xi-xi-1)。其中:di為壓縮因子,|di|<1。每個仿射變換ωi將映射區間{(xi,1,yi,1),(xi,2,yi,2)}映射到象區間{(xi-1,yi-1),(xi,yi)}。
由于映射沿x軸是收縮的,即沿x軸一個較大區間上的離散數據映射到一個較小的區間,計算映射到象區間上每點的映射值時,采用平均近似算法,即取
r=int(ai×x+ei)(3)
其中:int(·)表示取最靠近的整數。r為變換到點(r,yr)處的所有映射值的平均。
r=1/kk∈[m-1/2,m+1/2]yk(4)
在收縮因子di確定之后,即可由
ai=(xi-xi-1) (/ xi,2-xi,1)(5)
ei=(xi,2xi-1-xi,1xi) /( xi,2-xi,1)(6)
ci=(yi-yi-1) /( xi,2-xi,1) - di ( yi,2-yi,1) /(xi,2-xi,1)(7)
fi=(xi,2 yi-1-xi,1yi) /( xi,2-xi,1) - di(xi,2yi,1-xi,1yi,2) /
( xi,2-xi,1)(8)
確定自仿射ωi。
2.1di的確定
di應使映射區間經ωi映射后得到的結果與象區間盡可能接近,即應將它們之間的距離最小化:
E=xir=xi-1[r-yr]2=i,2j=i,1[ci×j+di×yj+fi-yint[ai×j+ei]]2→min(9)
將式(7)(8)代入(9)并令
λj=(xi,2-j) /( xi,2-xi,1)(10)
Aj=yj-[λjyi,1+(1-λj)yi,2](11)
Bj=yint[ai×j+ei]-[λjyi-1+(1-λj)yi](12)
得E=i,2j=i,1(Ajdi-Bj)2。令E / di=0得
di=∑i,2j=i,1AjBj / i,2j=i,1Aj2(13)
2.2 算法描述
分形自仿射分形模型具有較大的自由度,為了得到實際有效的搜索算法,對映射區間和象區間可能的選擇進行限制。預設映射區間和象區間的寬度分別為η=xi,2-xi,1和δ=xi-xi-1。其中:η和δ為常數且η>δ,各區間僅在區間的端點重合。這實際上是預先把給定的信號序列分成兩個集合:一個集合由寬度為η,共Q段的映射區間組成;另一個集合由寬度為δ共M段的象區間組成(M>Q)。搜索算法就是對每一個象區間尋找對應的與之匹配的映射區間。根據拼貼定理,最優匹配應使映射區間在仿射變換作用下的象與對應象區間的距離最小。
a)選擇η,δ(η>δ);
b)for i=1:M do;
for j=1:Q do
●按式(10)~(13)計算第i個象區間和第j個映射區間伴隨的仿射變換壓縮因子di,j,若|di,j|≥1,下一j;
●按式(5)~(8)計算映射因子,并求出映射區間與象區間之間的Hausdorff距離h(i,j)。存儲di,j和h(i,j),下一j;
●選擇最小的h(i,j),并令對應的di,j=di;
●下一i;
c)按式(1)計算序列分形維數D。
2.3算法分析
信號觀測長度對算法的影響較大,數據越多估計結果越準確,但計算量增大。研究表明,序列長度的選取與使用的方法、信號模型、信號的真實分形維數以及信噪比等基本無關。當觀測長度超過500時,分形維數估計結果基本上保持不變[9]。此外,該算法具有較大的靈活性,步驟a)中η和δ應根據具體情況選擇,通常不宜太大;b)中收縮因子d的判斷保證了序列是自仿射的。
3 仿真結果
為了驗證前面章節所提方法的形象,采用Mackey-Glass時滯序列和Lorenz系統x分量進行仿真。
3.1 Mackey-Glass時滯序列仿真
自從1977年Mackey和Glass發現時滯系統中的混沌現象以來,時滯混沌系統常常用來作為檢驗非線性系統模型性能的標準。該混沌時間序列由時滯微分方程產生:
dx(n) / dn=βx(n)+αx(n-τ) /[ 1+x(n-τ)10]
其中:α=0.2,β=-0.1。
當τ≥17時呈現混沌性,并且τ值越大,混沌程度越高。τ=17時奇怪吸引子及x(n)~n的示意圖如圖1(a)和(b)所示。已知其分形維數約為2.1。
采用基于序列自仿射的分形維數估計算法,取η=50,δ=25。當信號長度為100、250、500、1 000和1 500時,計算得分形維數分別為1.325、1.733、2.089、2.213、2.177??梢姺中尉S數與信號的長度關系比較明顯,當長度超過500時,計算值基本保持恒定,且與已知值(2.1)接近。
3.2 Lorenz系統x分量分形維數估計
Lorenz混沌系統為=σ(x+y)
=-xz+rx-y
=xy-bz,δ=10,r=28,b=8/3,其吸引子結構和x(n)~n序列分別如圖2(a)和(b)所示。已知其分形維數約為2.06。
采用5階龍格—庫塔法采樣得到x分量,并利用x(n)=[x(n)-minn{x(n)}] / [maxn{x(n)}-minn{x(n)}]進行歸一化。
實驗中模型參數取為η=30,δ=15。對不同的數據長度進行計算,當信號長度取500、750、1 000、1 500時,對應的估計分形維數分別為2.437、2.201、2.189和2.134,平均值與真實值的相對誤差為0.07。
仿真還表明,在保證壓縮因子存在的情況下,模型參數η和δ的選取對分形維數計算結果的影響甚小。
4 結束語
本文提出了一種基于序列分形自仿射特性的分形維數估計及其實現方法。從信號序列的自相似特性出發,基于迭代函數系統理論,通過求取序列本身仿射變換因子而獲得分形維數的估計。對已知分形維數的序列進行仿真,結果表明本方法具有較好的估計性能。
參考文獻:
[1]MANDELBROT B B.The fractal geometry on nature[M]. New York:Freeman,1983.
[2]屈世顯,張建華.復雜系統的分形理論與應用[M].西安:陜西人民出版社,1996.
[3]HADJILEONTIADIS L J. Wavelet-based enhancement of lung and bowel sounds using fractal dimension thresholding-part I: methodology [J].IEEE Trans on Biomedical Engineering, 2005(52):1143-1148.
[4]DISTASI R,NAPPI M,RICCIO D.A range/domain approximation error-based approach for fractal image compression [J]. IEEE Trans on Image Processing, 2006(15):89-97.
[5]BISOI A K, MISHRA J. On calculation of fractal dimension of images [J]. Pattern Recognition Letters, 2001,22(6):631-637.
[6]PENTLAND A P. Fractal based description of natural scense[J]. IEEE PAMI, 1984,6(6):661-674.
[7]MARAGOS P, SUN Fang-kuo. Measuring the fractal dimension of signals:morphological covers and iterative optimization[J].IEEE Trans on SP, 1993,41(1):108-121.
[8]ARNEODO A, DECOSTER N, ROUX S G.A wavelet-based method for multifractal image analysis.I. Methodology and test application on isotropic and anisotropic random rough surface[J]. Eur Phy J B,2000,15(13):567-600.
[9]謝文錄,謝維信.時間序列中的分形分析和參數提取[J].信號處理, 1997,13(2): 97-104.
[10]謝和平,薛秀謙.分形應用中的數學基礎[M].北京:科學出版社,1997.
注:本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文