测井数据的多重分形特征

测井数据的多重分形特征

0ctal的时间序列分析在应用差分分析技术时，对测量数据的应用取得了很大进展，主要体现在应用差分分析技术可以帮助开展测量数据重建。提高测量数据的分辨率。对测量数据进行校正。沉积相分析、裂缝研究、储层研究、岩石力学领域的井测量数据分析等。重标度极差分析(R/S分析)和功率谱分析是单一分形(monofractal)时间序列分析中最常用的2个方法。R/S分析的结果是Hurst指数(Hu),其数值范围为0<Hu<1;功率谱分析求得全域标度指数(H),而H的数值范围在理论上为-∞<H<+∞。笔者在实际应用中发现,Hu和H往往被误认为是可以互相替代或可对比的,这显然是不恰当的,不但给初学者带来混淆,还为正确认识时间序列的复杂性带来麻烦和困难,这是因为Hu和H是2个物理意义基本不同的参量,尽管在特定的情况下它们具有相似的功能,即代表了时间序列的复杂性和相关性。另外一个非常重要的问题是如何看待R/S分析方法的适用性以及如何在R/S分析中采取恰当的数据处理步骤。不恰当的分析结果常常会导致Hu趋近于1,而且似乎对几乎所有的地球物理时间序列都是如此,这就使得R/S分析失去其应用意义。本文对R/S分析和功率谱分析方法进行了仔细的研究和对比,提出了工作中的问题及可能的合理解决方案,为进一步的工作提供依据。文章中还将简单讨论具1/fβ现象的分形特征在测井数据中普遍存在的成因问题。1时间序列预处理功率谱分析方法几乎对所有复杂的时间序列都是适用的,因为对任何时间序列人们都可以通过傅立叶变换来获得其功率谱,然后分析其全域标度指数的大小,进而可以估计时间序列的分数维。然而R/S分析方法的直接应用只对静态(stationary)数据(比如以分形高斯噪声fGn为模型的时间序列)适用。R/S分析的数学公式为R(τ)S(τ)=max[1≤t≤τ∑i=1S(ξi-〈ξ〉τ)]-min[1≤t≤τ∑i=1S(ξi-〈ξ〉τ)][1ττ∑i=1(ξi-〈ξ〉τ)2]1/2=(aτ)Ηu(1)R(τ)S(τ)=max[∑i=11≤t≤τS(ξi−⟨ξ⟩τ)]−min[∑i=11≤t≤τS(ξi−⟨ξ⟩τ)][1τ∑i=1τ(ξi−⟨ξ⟩τ)2]1/2=(aτ)Hu(1)式中,1≤t≤τ∑i=1S∑i=11≤t≤τS表示从1到τ时间段内的逐项累积相加;〈ξ〉τ=1ττ∑i=1ξi⟨ξ⟩τ=1τ∑i=1τξi,代表了时间序列ξi在τ时间段内的平均值。依据式(1)列出R/S分析的6个主要步骤。①取初始τ值,求得时间序列ξi在τ时间段内的平均值〈ξ〉τ,然后观测时间序列相对于平均值的摆动情况ξi-〈ξ〉τ;②ξi-〈ξ〉τ逐项相加累积,得到时间序列ξi相对于其平均值〈ξ〉τ的累积偏离(accumulateddeparture)程度1≤t≤τ∑i=1S(ξi-〈ξ〉τ)∑i=11≤t≤τS(ξi−⟨ξ⟩τ)。显然,如果ξi是一随机时间序列(如高斯白噪随机序列),其累积偏离中不可能会有明显的趋向;相反,如果ξi是相关的分形高斯噪声,且具有持续或反持续的空间分布,其累积偏离的程度就会随着累积点数或区间τ的增加而益趋明显,所以累积偏离是对时间序列的持续或反持续程度的一种度量;③在区间τ范围内确定累积偏离的极大值与极小值之间的极差R;④为了对具不同幅度的时间序列分析结果作对比,应用区间τ范围内的标准差S将极差R归一化,得到比值R/S;⑤选取不同的τ值,重复步骤①到④;⑥最后在(R/S)-τ的双对数坐标下做线性回归分析来估计Hu。步骤①要求先去掉时间序列平均值,而对于非静态(或被称为伪平稳)时间序列(如分形布朗运动fBm,其点的数值代表某一时刻粒子运动所在的位置,其算术平均值没有清楚的物理意义)用式(1)来求得的Hu很难代表任何物理意义。实际研究表明,如果不对非静态时间序列做任何预处理,直接应用R/S方法最终算得的Hu数值常常偏大,并且非常接近甚至大于1,这就是为什么在文献中会看到Hu大于1的情况。从许多学者的研究成果分析表明,在全域一阶近似上,原始测井数据是类似于分形布朗运动的非静态时间序列。这样,应用R/S方法对原始测井数据直接进行分析往往会带来很大误差。许多学者在研究类似的时间序列中也注意到了这些问题并对R/S分析技术做了不同程度的改进[19,20,21,22,23,24]。为了消除这些问题,在R/S方法的实际操作中,人们尝试将原始的测井数据进行一些预处理,主要目的是为了消除时间序列中的短周期(或波长)变化,使原始的非高斯分布、非静态序列变为接近具高斯分布的静态序列,以此来希望提高R/S分析的准确性。这些预处理手段包括趋向消除和归一化。然而在这些方法的具体操作中,人为因素常常会由于参数选择的不同而使最终分析结果差别较大。实际上,这些预处理手段所隐含的目的或最终结果是将非静态的具fBm特征的时间序列转化为相应静态的具fGn特征的时间序列,然后在新的时间序列上采用R/S技术。那么,这些预处理手段与直接采用某一时间序列相对应的静态的具fGn特征的时间序列所得的结果之间是否可类比?如果可以,哪种方法更准确可靠一些?分析结果表明,采用某一时间序列相对应的静态的具fGn特征的时间序列(逐项增量序列successiveincrement)所得结果从(R/S)-τ的双对数图上看,具有更长的可用于回归分析的τ值范围和更高的回归相关系数(见图1),从而使得Hu的分析结果更准确。图1中的垂直虚线表示了不同预处理方法过渡带(Transientzone)的大小。另外,在数据处理上直接求得fBm序列中相邻2点的差值(即逐项增量序列)也比趋向消除和归一化等方法更客观,没有人为的受参数选择的影响而影响分析结果。更重要的是,采用某一时间序列相对应的静态的具fGn特征的逐项增量时间序列所求得的Hu值与应用功率谱分析对原始fBm时间序列的直接分析而得到的H结果具有较高的可类比性,从而使不同分析方法对同一数据分析所得到的结果可以互相对比并用以分析物理意义。在具体应用中,有些学者也是这样处理的,但是由于事先并不知道所要分析的时间序列是否具fGn特征还是fBm特征,如此的分析往往具有一些盲目性,且没有明确的标准来遵循。另外,这种R/S分析中的做法除了为了提高分析精度和可应用性外,有许多经验性成份。考虑到R/S分析中诸多误差影响和缺乏明确的理论依据,建议在实际应用中尽量避开R/S分析方法,而采取功率谱分析来估计全域标度指数H。2时间序列的物理含义在功率谱分析方法中,首先需要明确这种分析是对时间序列本身还是对时间序列的逐项累积(successiveaccumulation)而产生的序列。我们已经看到在R/S分析方法中,为了观测静态时间序列相对于其平均值的累积偏离,对时间序列逐项依次累积相加,得到与原来时间序列相同长度的新的时间序列。如果这种做法对于R/S分析方法还可以理解,若对功率谱分析则不应该采取同样的处理,因为这样的处理不但缺乏理论基础,而且还会影响人们对时间序列本身物理含义的正确认识。对于任何时间序列F(t),若有F(σt)=σΗF(t)(2)那么F(t)的功率谱S(f)与其频率f成幂指数关系S(f)∝f-β(3)可以证明β=2Η+1(4)需要注意的是式(4)对符合式(2)的时间序列本身,或是其逐项累积或逐项相减而产生的新的序列同样是成立的。由于H是通过傅里叶变换后在频率域中求得的,而功率谱反映的是时间序列的整体信息,H显然是对时间序列整体特征复杂性的度量,所以笔者称H为全域标度指数。与全域标度指数相对应用来描述信号局部标度(scaling)变化特征的量是Hölder指数α,又被称为Lipschitz指数。因为α刻划了某一点上或在某一点周围很小区域内的奇异性,α又被称为奇异指数,其数值的估计一般是在小波域中实现的,因为与傅立叶变换相比,小波变换具有优异的局部分析功能。对于单一分形时间序列,由于Hölder指数a在所有点上变化很小,可以认为全域标度指数H是所有点的Hölder指数的算术平均值Η≈1ΝΝ∑i=1(5)认识到H和α之间的密切联系,不仅它们在定义上接近,在算法上相似,而且这种联系揭示了H的物理意义,即H表征了时间序列整体上的复杂性。同时还说明H的估计应该是在原始时间序列上直接进行的,而不应该也没有必要在时间序列逐项累积序列上进行。与Hu不同,H的范围并不局限在,而主要是在[-3,3]的数值范围内,这样H比Hu更精确地反映了时间序列的复杂性。3静态fag类型的时间序列如果对fBm时间序列直接应用R/S方法或者用时间序列的逐项累积序列做功率谱分析,所求得的Hu或H往往会接近于1,这会使人们错误地认为几乎所有的地球物理数据都有Hu或H→1。那么对于任一地球物理测井数据,正确地应用R/S方法分析它的分形特征,其合理的步骤如下。①首先确定待分析的时间序列是静态fGn类型的时间序列,还是非静态fBm类型的数据。②应用功率谱分析技术估计时间序列的全域标度指数H,若H>0则表明该时间序列为fBm类型的数据;相反,若H<0则为fGn类型的时间序列(见图2)。③如果待分析的时间序列是非静态fBm类型的数据,应该先通过逐项相减而获得其相应的静态fGn类型的时间序列,然后对新的时间序列作R/S分析。如果待分析的时间序列已经是静态fGn类型的数据,可以直接对其做R/S分析。④比较2种估计方法所得到的Hu和H的大小。注意只有当待分析的时间序列是非静态fBm类型的数据时,按上述步骤所得到的Hu和H才具有可对比性。如果待分析的时间序列是静态的fGn类型的数据,则Hu>0,而H<0。由于所分析的测井数据几乎都是非静态fBm类型的数据,按上述步骤所得到的Hu和H应该接近,具有可对比性。同时认识到,只用其中的一种方法做分形分析有很多潜在的误差因素,所以应该采用至少2种不同的相对独立的分析方法,通过互相验证来确定所分析的结果。4多重分形谱的特征前面提到,原始测井资料在一阶近似上是fBm类型的数据,即通过傅里叶变换分析所得到的全域标度指数H>0。这样在一阶近似上从测井数据的整体上看,测井数据具有单一分形的特征。而在小波域中的分析结果表明,表征局部标度特征的Hölder指数在空间上变化很大。从图3可以看出,Hölder指数在不同点上变化明显。这说明测井数据并不只是简单的单一分形,而且还具有多重分形(multifractal)的特征。对数据的多重分形分析主要是通过对其多重分形谱的数值和形状的研究分析来实现。分析多重分形谱时需要注意分析的几个量是,(1)αmax和αmin;(2)αmax-αmin;(3)α0,f(α),在α0点达到最大值;(4)多重分形谱在某一点的斜率。多重分形谱在平面图上表现为向上突起的拱形,纵坐标f(α)代表了具不同奇异值数据点相对含量或分形维,其横坐标的数值分布范围或其宽度(αmax-αmin)则指示了数据中奇异值的分布范围。宽度越大表明了测井数据中的奇异值大小差别较大,宽度小表明了测井数据中的奇异值分布比较窄。αmax和αmin是数据中奇异性最小和最大点的奇异值,而时间序列中似乎存在一条在每一点上奇异值都为α0的虚拟序列,所以f(α0)=1。多重分形谱在某一点切线的斜率反映了具不同奇异值点的相对含量是如何变化的,在某奇异值点上对应的斜率越大,表明具有此奇异值的点的相对含量向相邻奇异值的点的相对含量有较快速的变化,反之则表明变化较慢。与具对称多重分形谱的二项乘积过程不同,实测的测井数据多重分形谱很少是对称的(见图4)。一般是奇异值小于α0的分支比奇异值大于α0的分支斜率要大,显示了较奇异的点的维数(或相对含量)随着奇异值的变化而变化较快。51沉积过程对总体成分的影响许多人对不同的地球物理数据的观测和分析结果表明,地球物理数据具有普遍的1/fβ现象,即在频率域内功率谱与频率的β次方成反比。这一现象表明不同尺度或量级的事件或物理参数(譬如波速度、密度、电导率,反射系数等)以不同的频率重复出现或集中,而且不同尺度或量级的事件或物理参数在频率-功率谱幅度的双对数坐标下呈现线性改变关系。一般原始地球物理测井数据具有更多的低频成份,与之相反,反射系数频率分布则具有更多的高频成份。这样当地震波在介质中传播时,地震波的高频成份与反射系数成更多更活跃的相互作用,这或许可以部分解释高频地震波衰减较快的原因。对于测井数据中普遍存在的1/fβ现象,不同学者对其成因有不同的观点。主要有两大派:地层成因和裂隙成因。笔者认为无论是地层成因,还是裂隙成因还是多因素的复合,都可能会导致1/fβ现象,只倾向于某一成因而摈弃其他成因的态度是不科学的。当然对于某一特定的研究目标,某些因素可能较其他因素占有更重要的地位。笔者在研究江西北部上寒武统西阳山组碳酸盐岩岩石学与地层分布的特征中发现,整套地层由一系列不同级别的旋回层组成。小尺度的旋回层镶嵌于大尺度的旋回层内,而大尺度的旋回层又镶嵌于更大尺度的旋回层内。剖面中泥、炭质含量的变化有其独特的规律性,无论是整个剖面,还是大旋回层或小旋回层,泥、岩质含量自底到顶都呈增加趋向,直到最后碳酸盐岩消失而被页岩所代替。另外,就整个剖面而言,自底至顶,每个大尺度旋回层的厚度变小,旋回交替加快;就大尺度旋回层而言,自底到顶,其中的小旋回层的厚度也变小,旋回交替也加快。这说明这一地层的分布既不是随机的,又不是周期性的重复,而是具有内在规律的分形分布。这一碳酸盐岩沉积层序的奇特的分布特征很难用传统的马尔可夫链-蒙托卡罗方法(Markovchain-MonteCarlo)来模拟,但是可以很好地用二项乘积过程(binomialmultiplicativeprocess)来模拟,而二项乘积过程是产生多重分形的最基本的方法。这一简单的例子表明,单纯的沉积过程中各因素间复杂的相互作用而形成的沉积层本身就会具有很好的自相似和多重分形特征。当然其他因素的叠加,譬如裂隙的发育,有可能会进一步导致测井数据的1/fβ特征。6数据处理及分析地下介质物性参数的分布不是随机的,而是具有复杂的分形和多重分形的特征。这就需要精确刻划和描述地下介质,从而帮助物性参数(如孔隙度、含油性等)的预测和统计。下面以具体的测井数据分析结果说明分形分析方法正确合理运用的重要性。我们选择美国堪萨斯州西南部某油田的γ射线测井数据,该测井长度为1296.162m,测井数据间隔为0.125m,主要岩性为灰岩与页岩互层。对该组数据分别做R/S分析和功率谱分析,发现当对数据直接进行R/S分析时,Hu=0.91。这一数值明显偏高。然后对测井数据进行一些基本的预处理(见图1),包括归一化、趋势消除和计算逐项增量序列等。预处理后的数据Hu值在0.38～0.54之间。而功率谱分析的结果是H=0.27,显然与R/S分析的结果相差很大。类似的由于不同的分析技术对同一数据的分析结果的明显差异也在前人的研究中出