基于分形插值的多重分形降趋波动分析法及其有效性检验

刘 慧, 万 丽, 曾祥健, 邓小成

(广州大学 数学与信息科学学院, 广东 广州 510006)

多重分形是由定义在具有自相似性的形态结构上的无穷多个标度指数的测度所组成的度量集合，其谱函数可以描述分形结构上不同的局域条件或分形结构在演化过程中不同层次所导致的特殊行为与特征[1-2]。常用于计算多重分形谱及相关参数的方法是配分函数法，但该方法无法准确反映有趋势影响的非平稳数据的多重特征。近年来，新提出了多重分形降趋波动分析法(Multifractal Detrended Fluctuation Analysis, MFDFA)[3]、多重分形降趋移动平均分析法(Multifractal Detrending Moving Average Analysis, MFDMAA)[4]、小波领导分析法(Wavelet Leaders Analysis, WLA)[5]和小波模极大值法(Wavelet Transform Modulus Maxima, WTMM)[6]等，其中，MFDFA方法由于运算简便且易于实现的优点，已证实是描述非线性序列复杂性的有效定量化工具之一，被广泛应用在化学、生物、医学、地质和物理等各个学科领域，但该方法也存在着产生伪波动误差及多项式拟合阶数的影响较大的缺点[7-10]。Barnsley于1986年首次提出了分形插值技术[11]，为拟合分形数据提供了新思想，分形插值方法具有拟合精度高的优点。

为了找出更合适的去趋势和去噪声的方式，本文将分形插值拟合技术引入到MFDFA方法中，给出了基于分形插值的多重分形降趋波动分析法(Fractal Interpolation based Multifractal Detrended Fluctuation Analysis, FI-MFDFA)，并验证其有效性；进一步从算法模型异同、数据样本量的变化和多重参数计算的统计精度等方面，对比分析FI-MFDFA和MFDFA方法的优劣，为应用这2种方法来分析实际序列的多重分形特征和长程相关性提供理论支持。

1 理论与方法

1.1 分形插值

分形插值方法是建立在迭代函数系(IFS)理论基础上，通过给定一组插值点{(xj,yj)∈R2,j=0,1,2,…,N}构造出满足插值条件f(xj)=yj(j=0,1,2,…,N)的连续函数f:[x0,xN]→R，其中，x0

(1)设插值区间为I=[x0,xN]，两点区间为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，Fi(xN,yN)=yi

(2)

(2)定义仿射变换Wi(xj,yj)=(Li(xj),Fi(xj,yj)),i=1,2,…,N;j=0,1,2,…,N，从而构造出迭代函数系IFS{K;Wi,i=1,2,…,N}，当Li(xj,yj)和Fi(xj,yj)是线性函数时，可得到如下公式：

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

(3)

其中，(x′m,y′m)为进行第一次迭代插值之后得到的数据长度为m的序列，在进行第二次迭代时，将(x′m,y′m)作为新的插值点进行插值，以此类推。当迭代次数为k时，生成的插值数据个数m和插值点个数N之间的关系为m=(N-1)2k+1。

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

(4)

(5)

(6)

(7)

式中，di(i=1,2,…,N)为垂直比例因子，是仿射变换Wi的一个自由变量，满足|di|<1，代入相应的di，系数ai,ei,ci,fi(i=1,2,…,N)可由插值点数据计算得出，因此，分形插值拟合可以通过改变垂直比例因子di的值，来改变插值函数方程f(di)的表达式。

1.2 FI-MFDFA

研究非平稳序列的途径是先通过降趋势过程将序列平稳化，再对序列的相关性特征进行分析。降非线性趋势的常用方式是多项式拟合，而在实际应用中，无法准确判断数据具有几阶多项式趋势，从而在结果分析时容易产生偏差。FI-MFDFA通过分形插值拟合消除序列的分形趋势，避免多项式阶数选取的主观性影响。

对给定长度为n的序列{xi}(i=1,2,…,n)，FI-MFDFA的计算步骤如下：

步骤1：求序列{xi}的累积离差序列y(i)

(8)

步骤3：把序列y(i)从第一个数据开始等长度地分割成尺度为s的Ns=int(N/s)个互不相交的数据段，由于长度N经常不是s的整数倍，为了不丢弃尾部剩余部分，从序列最后一个数据重复这一分割过程，因此,得到2Ns个区间。

步骤4：用最小二乘拟合法求均方误差F2(s,v)。设fv(i)为第v个小区间的分形插值函数方程，控制每一个小区间的函数方程中的垂直比例因子d不变，均等于第二步计算出的d。

当v=1,2,...,Ns时，有

(9)

当v=Ns+1,Ns+2,…,2Ns时，有

(10)

步骤5：求序列的q阶波动函数(q为整数)

(11)

当q=0时，

(12)

步骤6：确定波动函数的标度指数，根据Fq(s)与s的关系

Fq(s)∝sh(q)

(13)

先固定阶数q，作lnFq(s)对lns的函数关系图，其拟合直线的斜率即为所得的标度指数h(q)。这里h(q)称为广义赫斯特(Hurst)指数，当序列是平稳时间序列时，h(2)称为Hurst指数。通常，波动函数值Fq(s)是s的增函数，广义Hurst指数h(q)是随q变化的单调减函数。当序列的小波动和大波动具有不同的标度行为时，h(q)显著依赖于q，序列表现为多重分形；当h(q)独立于q为一常数时，即广义指数函数h(q)不随q变化，则序列表现为单一分形。

步骤7：q阶广义Hurst指数h(q)与质量指数τ(q)的关系为

τ(q)=qh(q)-1

(14)

根据legendre变换得到多重分形奇异指数α和奇异谱函数f(α)：

α=dτ(q)/dq

(15)

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

(16)

奇异指数用来描述观测序列中不同区间的奇异程度，奇异谱函数用于描述不同区间奇异指数的分形维数。当f(α)独立于α为一常数时，序列表现为单一分形特征；当f(α)的形状呈单峰凸分布时，序列表现为多重分形特征。

2 FI-MFDFA方法检验

选择经典二项多重分形序列(BMS)模型检验FI-MFDFA方法的有效性。该模型构造如下：

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

(17)

式中参数p的取值范围是0

2.1 多重分形特征识别

选取参数p=0.20、0.25、0.30、0.35和0.40，长度L=1 024的BMS序列作为分析对象，运用FI-MFDFA方法分析序列在不同参数下的多重分形特征。无标度区间s的取值范围从20到L/4，步长为10；阶数q从-3到3，以0.2为步长均匀取31个值。图1(a)是不同参数p下的广义Hurst指数h(q)随着q变化的曲线，当q从-3增加到3时，不同参数p下的h(q)均随着q的增大而减小，说明FI-MFDFA能识别出BMS序列的多重分形特征，且Δh=h(-3)-h(3)随着p增大而减小，说明奇异性越大的序列，多重分形特征越明显；图1(b)是不同参数p下的多重分形谱α-f(α)曲线，从图形可以看出，随着p的减小，奇异指数α增大，即分布奇异性越大的序列多重分形强度越大，说明FI-MFDFA能较好地识别序列的多重分形性。

图1 BMS序列不同参数p的FI-MFDFA计算结果Fig.1 FI-MFDFA calculation results of different parameters p in the BMS sequence

2.2 与MFDFA的对比分析

2.2.1 方法步骤比较

FI-MFDFA与MFDFA主要不同点在第3步。MFDFA方法用多项式拟合法求残差，当拟合的多项式阶数为m，相应的方法记作MFDFAm，虽然拟合多项式的阶数可以根据具体情况灵活设置，但阶数的确定具有主观性，选取的阶数过小会导致不完全去除趋势，选取的阶数过大会引起过拟合；而FI-MFDFA用分形插值拟合法求残差，可以解决多项式阶数的选取不恰当对分析结果造成的影响。

2.2.2 参数的统计精度比较

采用参数p=0.3，长度L=256的BMS序列对比FI-MFDFA和MFDFA 2种方法的计算统计精度。无标度区间s取值范围为20到L/4，步长为10，阶数q从-3到3，间隔为0.2均匀取31个值。图2(a)和(b)分别是Hurst指数h(q)和质量指数τ(q)与阶数q的关系图；图2(c)是FI-MFDFA与MFDFA1、MFDFA2及MFDFA3方法分析得到的多重分形谱f(α)和理论值谱线；为了更好地比较4种方法的统计精度，计算理论Hurst指数H(q)与实际估计的Hurst指数h(q)的差值：

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

(18)

图2(d)是Δh(q)与q的关系图。图2(a)和(b)中h(q)随着q的变化而变化，且τ(q)与q之间不是直线关系，故该序列是多重分形序列；FI-MFDFA的q阶Hurst指数更接近理论值，其次是MFDFA1和MFDFA3、MFDFA2的计算结果离理论值最远。由图2(c)可以看出，多重分形谱曲线均往左偏离理论值曲线；相比于q<0部分，FI-MFDFA方法计算出来的分形谱曲线逼近理论值的效果优于q>0部分；MFDFA方法随着多项式拟合阶数的增加，计算出的多重分形谱与理论值的偏差变化较大，当阶数m=1时，效果优于m=2和m=3，拟合阶数过大会引起过拟合现象，导致多重分形谱的形状和宽度偏离理论值。图2(d)中FI-MFDFA的Δh(q)波动斜率小于MFDFA，随着q增大，4种方法下的Δh(q)均有所下降，其中,FI-MFDFA的Δh(q)小于MFDFA，更接近0。

表1是在FI-MFDFA和MFDFA方法下计算的Hurst指数h(q)、质量指数τ(q)以及多重分形谱参数α和f(α)对于理论值的均方根误差。4种方法下的参数h(q)、τ(q)、α和f(α)均方根误差的大小关系：FI-MFDFA

表1 多重分形参数的均方根误差

2.2.3 样本量的影响比较

选取p=0.3，样本量L分别为256、512和1 024的BMS序列，运用FI-MFDFA、MFDFA1、MFDFA2和MFDFA34种方法分析序列的多重分形性，控制方法中的无标度区间不变，尺度区间s取值范围为20到L/4，步长为10，阶数q从-3到3，间隔为0.2均匀取31个值。图3分别为FI-MFDFA、MFDFA1、MFDFA2和MFDFA34种方法计算出的多重分形谱，并与理论值对比。随着样本量的增加，4种方法计算的分形谱右边部分的波动变化均大于左边部分，MFDFA2和MFDFA3的波动明显大于FI-MFDFA。FI-MFDFA对于样本量的变化不是特别敏感，在样本量为256时，多重分形谱的形状在q<0部分略偏离理论值，在样本量为512和1 024时，多重分形谱的形状更接近于理论值；MFDFA方法在样本量L为256时，分形谱曲线偏离理论值最远，在样本量为1 024时，分形谱曲线最接近理论值，其中MFDFA2和MFDFA3受小样本量的影响较大。

图3 不同样本量序列的多重分形谱Fig.3 Multifractal spectrum of different sample series

表2是不同样本量序列在4种方法下的Hurst指数的均方根误差。FI-MFDFA方法在计算同一样本量序列的实际Hurst指数与理论值的偏差均小于MFDFA。4种方法计算的h(q)均方根误差均随着样本量L的增大而逐渐减小，其中FI-MFDFA方法对于样本量的增加，Hurst指数偏离理论值的变化波动最小，当L=512和L=1 024时的h(q)与理论值的均方根误差值分别为0.074 1和0.066 9，当L=256时，均方根误差为0.101 5，均属于可以接受的误差范围；MFDFA方法的多项式拟合阶数取1时，计算3个样本量序列的均方根误差均小于MFDFA3，而MFDFA2最大，只有MFDFA1在L=1 024时均方根误差值小于0.1。综合分析，使用MFDFA方法对序列进行分析最少需要1 024个数据点，而FI-MFDFA方法取512个数据点就可满足计算精度，对于小样本量也可达到满意的精度。

表2 不同样本量序列的Hurst指数均方根误差

3 结 论

将分形插值技术与降趋势波动方法相结合给出基于分形插值的降趋势波多重分形方法(FI-MFDFA)，并利用BMS模型对该方法进行了检验，从方法的算法步骤、参数统计精度和样本容量的敏感性3个方面，对比分析了FI-MFDFA方法与MFDFA方法的优劣性。分析结果显示：FI-MFDFA方法能有效识别多重分形强度，其多重参数h(q)、τ(q)、α和f(α)计算结果的均方根误差均小于MFDFA方法的对应值，且受数据量大小的影响也较小，表明了FI-MFDFA方法要明显优于MFDFA方法，还能避免多项式拟合阶数的变化对多重参数计算结果的影响，为进一步应用该方法研究实际数据的多重分形特征提供了理论支持。该方法对二维和高维数据的应用有待进一步研究。