劉 慧, 萬 麗, 曾祥健, 鄧小成

(廣州大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院, 廣東 廣州 510006)

多重分形是由定義在具有自相似性的形態(tài)結(jié)構(gòu)上的無窮多個(gè)標(biāo)度指數(shù)的測(cè)度所組成的度量集合，其譜函數(shù)可以描述分形結(jié)構(gòu)上不同的局域條件或分形結(jié)構(gòu)在演化過程中不同層次所導(dǎo)致的特殊行為與特征[1-2]。常用于計(jì)算多重分形譜及相關(guān)參數(shù)的方法是配分函數(shù)法，但該方法無法準(zhǔn)確反映有趨勢(shì)影響的非平穩(wěn)數(shù)據(jù)的多重特征。近年來，新提出了多重分形降趨波動(dòng)分析法(Multifractal Detrended Fluctuation Analysis, MFDFA)[3]、多重分形降趨移動(dòng)平均分析法(Multifractal Detrending Moving Average Analysis, MFDMAA)[4]、小波領(lǐng)導(dǎo)分析法(Wavelet Leaders Analysis, WLA)[5]和小波模極大值法(Wavelet Transform Modulus Maxima, WTMM)[6]等，其中，MFDFA方法由于運(yùn)算簡(jiǎn)便且易于實(shí)現(xiàn)的優(yōu)點(diǎn)，已證實(shí)是描述非線性序列復(fù)雜性的有效定量化工具之一，被廣泛應(yīng)用在化學(xué)、生物、醫(yī)學(xué)、地質(zhì)和物理等各個(gè)學(xué)科領(lǐng)域，但該方法也存在著產(chǎn)生偽波動(dòng)誤差及多項(xiàng)式擬合階數(shù)的影響較大的缺點(diǎn)[7-10]。Barnsley于1986年首次提出了分形插值技術(shù)[11]，為擬合分形數(shù)據(jù)提供了新思想，分形插值方法具有擬合精度高的優(yōu)點(diǎn)。

為了找出更合適的去趨勢(shì)和去噪聲的方式，本文將分形插值擬合技術(shù)引入到MFDFA方法中，給出了基于分形插值的多重分形降趨波動(dòng)分析法(Fractal Interpolation based Multifractal Detrended Fluctuation Analysis, FI-MFDFA)，并驗(yàn)證其有效性；進(jìn)一步從算法模型異同、數(shù)據(jù)樣本量的變化和多重參數(shù)計(jì)算的統(tǒng)計(jì)精度等方面，對(duì)比分析FI-MFDFA和MFDFA方法的優(yōu)劣，為應(yīng)用這2種方法來分析實(shí)際序列的多重分形特征和長(zhǎng)程相關(guān)性提供理論支持。

1 理論與方法

1.1 分形插值

分形插值方法是建立在迭代函數(shù)系(IFS)理論基礎(chǔ)上，通過給定一組插值點(diǎn){(xj,yj)∈R2,j=0,1,2,…,N}構(gòu)造出滿足插值條件f(xj)=yj(j=0,1,2,…,N)的連續(xù)函數(shù)f:[x0,xN]→R，其中，x0

(1)設(shè)插值區(qū)間為I=[x0,xN]，兩點(diǎn)區(qū)間為Ii=[xi-1,xi],i=1,2,…,N，壓縮變換Li:I→Ii和Fi:K=I×R→R滿足

Li(x0)=xi-1，Li(xN)=xi

(1)

Fi(x0,y0)=yi-1，F(xiàn)i(xN,yN)=yi

(2)

(2)定義仿射變換Wi(xj,yj)=(Li(xj),Fi(xj,yj)),i=1,2,…,N;j=0,1,2,…,N，從而構(gòu)造出迭代函數(shù)系IFS{K;Wi,i=1,2,…,N}，當(dāng)Li(xj,yj)和Fi(xj,yj)是線性函數(shù)時(shí)，可得到如下公式：

i=1,2,…,N;j=0,1,2,…,N

(3)

其中，(x′m,y′m)為進(jìn)行第一次迭代插值之后得到的數(shù)據(jù)長(zhǎng)度為m的序列，在進(jìn)行第二次迭代時(shí)，將(x′m,y′m)作為新的插值點(diǎn)進(jìn)行插值，以此類推。當(dāng)?shù)螖?shù)為k時(shí)，生成的插值數(shù)據(jù)個(gè)數(shù)m和插值點(diǎn)個(gè)數(shù)N之間的關(guān)系為m=(N-1)2k+1。

(3)將式(1)和式(2)代入式(3)中，解得

(4)

(5)

(6)

(7)

式中，di(i=1,2,…,N)為垂直比例因子，是仿射變換Wi的一個(gè)自由變量，滿足|di|<1，代入相應(yīng)的di，系數(shù)ai,ei,ci,fi(i=1,2,…,N)可由插值點(diǎn)數(shù)據(jù)計(jì)算得出，因此，分形插值擬合可以通過改變垂直比例因子di的值，來改變插值函數(shù)方程f(di)的表達(dá)式。

1.2 FI-MFDFA

研究非平穩(wěn)序列的途徑是先通過降趨勢(shì)過程將序列平穩(wěn)化，再對(duì)序列的相關(guān)性特征進(jìn)行分析。降非線性趨勢(shì)的常用方式是多項(xiàng)式擬合，而在實(shí)際應(yīng)用中，無法準(zhǔn)確判斷數(shù)據(jù)具有幾階多項(xiàng)式趨勢(shì)，從而在結(jié)果分析時(shí)容易產(chǎn)生偏差。FI-MFDFA通過分形插值擬合消除序列的分形趨勢(shì)，避免多項(xiàng)式階數(shù)選取的主觀性影響。

對(duì)給定長(zhǎng)度為n的序列{xi}(i=1,2,…,n)，F(xiàn)I-MFDFA的計(jì)算步驟如下：

步驟1：求序列{xi}的累積離差序列y(i)

(8)

步驟3：把序列y(i)從第一個(gè)數(shù)據(jù)開始等長(zhǎng)度地分割成尺度為s的Ns=int(N/s)個(gè)互不相交的數(shù)據(jù)段，由于長(zhǎng)度N經(jīng)常不是s的整數(shù)倍，為了不丟棄尾部剩余部分，從序列最后一個(gè)數(shù)據(jù)重復(fù)這一分割過程，因此,得到2Ns個(gè)區(qū)間。

步驟4：用最小二乘擬合法求均方誤差F2(s,v)。設(shè)fv(i)為第v個(gè)小區(qū)間的分形插值函數(shù)方程，控制每一個(gè)小區(qū)間的函數(shù)方程中的垂直比例因子d不變，均等于第二步計(jì)算出的d。

當(dāng)v=1,2,...,Ns時(shí)，有

(9)

當(dāng)v=Ns+1,Ns+2,…,2Ns時(shí)，有

(10)

步驟5：求序列的q階波動(dòng)函數(shù)(q為整數(shù))

(11)

當(dāng)q=0時(shí)，

(12)

步驟6：確定波動(dòng)函數(shù)的標(biāo)度指數(shù)，根據(jù)Fq(s)與s的關(guān)系

Fq(s)∝sh(q)

(13)

先固定階數(shù)q，作lnFq(s)對(duì)lns的函數(shù)關(guān)系圖，其擬合直線的斜率即為所得的標(biāo)度指數(shù)h(q)。這里h(q)稱為廣義赫斯特(Hurst)指數(shù)，當(dāng)序列是平穩(wěn)時(shí)間序列時(shí)，h(2)稱為Hurst指數(shù)。通常，波動(dòng)函數(shù)值Fq(s)是s的增函數(shù)，廣義Hurst指數(shù)h(q)是隨q變化的單調(diào)減函數(shù)。當(dāng)序列的小波動(dòng)和大波動(dòng)具有不同的標(biāo)度行為時(shí)，h(q)顯著依賴于q，序列表現(xiàn)為多重分形；當(dāng)h(q)獨(dú)立于q為一常數(shù)時(shí)，即廣義指數(shù)函數(shù)h(q)不隨q變化，則序列表現(xiàn)為單一分形。

步驟7：q階廣義Hurst指數(shù)h(q)與質(zhì)量指數(shù)τ(q)的關(guān)系為

τ(q)=qh(q)-1

(14)

根據(jù)legendre變換得到多重分形奇異指數(shù)α和奇異譜函數(shù)f(α)：

α=dτ(q)/dq

(15)

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

(16)

奇異指數(shù)用來描述觀測(cè)序列中不同區(qū)間的奇異程度，奇異譜函數(shù)用于描述不同區(qū)間奇異指數(shù)的分形維數(shù)。當(dāng)f(α)獨(dú)立于α為一常數(shù)時(shí)，序列表現(xiàn)為單一分形特征；當(dāng)f(α)的形狀呈單峰凸分布時(shí)，序列表現(xiàn)為多重分形特征。

2 FI-MFDFA方法檢驗(yàn)

選擇經(jīng)典二項(xiàng)多重分形序列(BMS)模型檢驗(yàn)FI-MFDFA方法的有效性。該模型構(gòu)造如下：

x(i)=pn(i-1)(1-p)n-n(i-1)，i=1,2,…,2n

(17)

式中參數(shù)p的取值范圍是0

2.1 多重分形特征識(shí)別

選取參數(shù)p=0.20、0.25、0.30、0.35和0.40，長(zhǎng)度L=1 024的BMS序列作為分析對(duì)象，運(yùn)用FI-MFDFA方法分析序列在不同參數(shù)下的多重分形特征。無標(biāo)度區(qū)間s的取值范圍從20到L/4，步長(zhǎng)為10；階數(shù)q從-3到3，以0.2為步長(zhǎng)均勻取31個(gè)值。圖1(a)是不同參數(shù)p下的廣義Hurst指數(shù)h(q)隨著q變化的曲線，當(dāng)q從-3增加到3時(shí)，不同參數(shù)p下的h(q)均隨著q的增大而減小，說明FI-MFDFA能識(shí)別出BMS序列的多重分形特征，且Δh=h(-3)-h(3)隨著p增大而減小，說明奇異性越大的序列，多重分形特征越明顯；圖1(b)是不同參數(shù)p下的多重分形譜α-f(α)曲線，從圖形可以看出，隨著p的減小，奇異指數(shù)α增大，即分布奇異性越大的序列多重分形強(qiáng)度越大，說明FI-MFDFA能較好地識(shí)別序列的多重分形性。

圖1 BMS序列不同參數(shù)p的FI-MFDFA計(jì)算結(jié)果Fig.1 FI-MFDFA calculation results of different parameters p in the BMS sequence

2.2 與MFDFA的對(duì)比分析

2.2.1 方法步驟比較

FI-MFDFA與MFDFA主要不同點(diǎn)在第3步。MFDFA方法用多項(xiàng)式擬合法求殘差，當(dāng)擬合的多項(xiàng)式階數(shù)為m，相應(yīng)的方法記作MFDFAm，雖然擬合多項(xiàng)式的階數(shù)可以根據(jù)具體情況靈活設(shè)置，但階數(shù)的確定具有主觀性，選取的階數(shù)過小會(huì)導(dǎo)致不完全去除趨勢(shì)，選取的階數(shù)過大會(huì)引起過擬合；而FI-MFDFA用分形插值擬合法求殘差，可以解決多項(xiàng)式階數(shù)的選取不恰當(dāng)對(duì)分析結(jié)果造成的影響。

2.2.2 參數(shù)的統(tǒng)計(jì)精度比較

采用參數(shù)p=0.3，長(zhǎng)度L=256的BMS序列對(duì)比FI-MFDFA和MFDFA 2種方法的計(jì)算統(tǒng)計(jì)精度。無標(biāo)度區(qū)間s取值范圍為20到L/4，步長(zhǎng)為10，階數(shù)q從-3到3，間隔為0.2均勻取31個(gè)值。圖2(a)和(b)分別是Hurst指數(shù)h(q)和質(zhì)量指數(shù)τ(q)與階數(shù)q的關(guān)系圖；圖2(c)是FI-MFDFA與MFDFA1、MFDFA2及MFDFA3方法分析得到的多重分形譜f(α)和理論值譜線；為了更好地比較4種方法的統(tǒng)計(jì)精度，計(jì)算理論Hurst指數(shù)H(q)與實(shí)際估計(jì)的Hurst指數(shù)h(q)的差值：

Δh(q)=H(q)-h(q)

(18)

圖2(d)是Δh(q)與q的關(guān)系圖。圖2(a)和(b)中h(q)隨著q的變化而變化，且τ(q)與q之間不是直線關(guān)系，故該序列是多重分形序列；FI-MFDFA的q階Hurst指數(shù)更接近理論值，其次是MFDFA1和MFDFA3、MFDFA2的計(jì)算結(jié)果離理論值最遠(yuǎn)。由圖2(c)可以看出，多重分形譜曲線均往左偏離理論值曲線；相比于q<0部分，F(xiàn)I-MFDFA方法計(jì)算出來的分形譜曲線逼近理論值的效果優(yōu)于q>0部分；MFDFA方法隨著多項(xiàng)式擬合階數(shù)的增加，計(jì)算出的多重分形譜與理論值的偏差變化較大，當(dāng)階數(shù)m=1時(shí)，效果優(yōu)于m=2和m=3，擬合階數(shù)過大會(huì)引起過擬合現(xiàn)象，導(dǎo)致多重分形譜的形狀和寬度偏離理論值。圖2(d)中FI-MFDFA的Δh(q)波動(dòng)斜率小于MFDFA，隨著q增大，4種方法下的Δh(q)均有所下降，其中,FI-MFDFA的Δh(q)小于MFDFA，更接近0。

表1是在FI-MFDFA和MFDFA方法下計(jì)算的Hurst指數(shù)h(q)、質(zhì)量指數(shù)τ(q)以及多重分形譜參數(shù)α和f(α)對(duì)于理論值的均方根誤差。4種方法下的參數(shù)h(q)、τ(q)、α和f(α)均方根誤差的大小關(guān)系：FI-MFDFA

表1 多重分形參數(shù)的均方根誤差

2.2.3 樣本量的影響比較

選取p=0.3，樣本量L分別為256、512和1 024的BMS序列，運(yùn)用FI-MFDFA、MFDFA1、MFDFA2和MFDFA34種方法分析序列的多重分形性，控制方法中的無標(biāo)度區(qū)間不變，尺度區(qū)間s取值范圍為20到L/4，步長(zhǎng)為10，階數(shù)q從-3到3，間隔為0.2均勻取31個(gè)值。圖3分別為FI-MFDFA、MFDFA1、MFDFA2和MFDFA34種方法計(jì)算出的多重分形譜，并與理論值對(duì)比。隨著樣本量的增加，4種方法計(jì)算的分形譜右邊部分的波動(dòng)變化均大于左邊部分，MFDFA2和MFDFA3的波動(dòng)明顯大于FI-MFDFA。FI-MFDFA對(duì)于樣本量的變化不是特別敏感，在樣本量為256時(shí)，多重分形譜的形狀在q<0部分略偏離理論值，在樣本量為512和1 024時(shí)，多重分形譜的形狀更接近于理論值；MFDFA方法在樣本量L為256時(shí)，分形譜曲線偏離理論值最遠(yuǎn)，在樣本量為1 024時(shí)，分形譜曲線最接近理論值，其中MFDFA2和MFDFA3受小樣本量的影響較大。

圖3 不同樣本量序列的多重分形譜Fig.3 Multifractal spectrum of different sample series

表2是不同樣本量序列在4種方法下的Hurst指數(shù)的均方根誤差。FI-MFDFA方法在計(jì)算同一樣本量序列的實(shí)際Hurst指數(shù)與理論值的偏差均小于MFDFA。4種方法計(jì)算的h(q)均方根誤差均隨著樣本量L的增大而逐漸減小，其中FI-MFDFA方法對(duì)于樣本量的增加，Hurst指數(shù)偏離理論值的變化波動(dòng)最小，當(dāng)L=512和L=1 024時(shí)的h(q)與理論值的均方根誤差值分別為0.074 1和0.066 9，當(dāng)L=256時(shí)，均方根誤差為0.101 5，均屬于可以接受的誤差范圍；MFDFA方法的多項(xiàng)式擬合階數(shù)取1時(shí)，計(jì)算3個(gè)樣本量序列的均方根誤差均小于MFDFA3，而MFDFA2最大，只有MFDFA1在L=1 024時(shí)均方根誤差值小于0.1。綜合分析，使用MFDFA方法對(duì)序列進(jìn)行分析最少需要1 024個(gè)數(shù)據(jù)點(diǎn)，而FI-MFDFA方法取512個(gè)數(shù)據(jù)點(diǎn)就可滿足計(jì)算精度，對(duì)于小樣本量也可達(dá)到滿意的精度。

表2 不同樣本量序列的Hurst指數(shù)均方根誤差

3 結(jié) 論

將分形插值技術(shù)與降趨勢(shì)波動(dòng)方法相結(jié)合給出基于分形插值的降趨勢(shì)波多重分形方法(FI-MFDFA)，并利用BMS模型對(duì)該方法進(jìn)行了檢驗(yàn)，從方法的算法步驟、參數(shù)統(tǒng)計(jì)精度和樣本容量的敏感性3個(gè)方面，對(duì)比分析了FI-MFDFA方法與MFDFA方法的優(yōu)劣性。分析結(jié)果顯示：FI-MFDFA方法能有效識(shí)別多重分形強(qiáng)度，其多重參數(shù)h(q)、τ(q)、α和f(α)計(jì)算結(jié)果的均方根誤差均小于MFDFA方法的對(duì)應(yīng)值，且受數(shù)據(jù)量大小的影響也較小，表明了FI-MFDFA方法要明顯優(yōu)于MFDFA方法，還能避免多項(xiàng)式擬合階數(shù)的變化對(duì)多重參數(shù)計(jì)算結(jié)果的影響，為進(jìn)一步應(yīng)用該方法研究實(shí)際數(shù)據(jù)的多重分形特征提供了理論支持。該方法對(duì)二維和高維數(shù)據(jù)的應(yīng)用有待進(jìn)一步研究。