2006年10月 第41卷　第5期　3山东省东营市中国石油大学(华东)信息与控制工程学院,257061本文于2005年12月21日收到,修改稿于2006年5月12日收到。

本项研究受高等学校博士学科点专项科研基金(No.20020008004)部分资助。

・处理方法・高阶统计量地震子波估计建模戴永寿3①②　郑德玲①　魏　磊②　霍志勇②(①北京科技大学信息工程学院;②中国石油大学(华东)信息与控制工程学院)摘 要戴永寿,郑德玲,魏磊,霍志勇.高阶统计量地震子波估计建模.石油地球物理勘探,2006,41(5):514～518,540 本文在反射系数序列为非高斯、平稳和统计独立的随机过程,地震子波为非因果、混合相位的假设条件下,分别应用滑动平均(MA )和自回归滑动平均(ARMA )模型对地震记录进行建模,并采用运算代价较小的基于高阶累积量的线性化求解方法———累积量矩阵方程法进行了子波提取和模型适应性的研究。

数值模拟结果和实际地震数据处理结果表明:自回归滑动平均(ARMA )模型比滑动平均(MA )模型具有参数节省、模型更为高效的特点;累积量矩阵方程法可以有效地压制加性高斯噪声,但对累积量样本估计的准确性要求较高;如果累积量样本估计的误差和方差适度,结合自回归滑动平均(ARMA )模型描述的累积量矩阵方程法可以高效、准确地估计出地震子波。

关键词　高阶累积量　子波　自回归滑动平均(ARMA )　滑动平均(MA )　建模1　引言作为地震资料反褶积处理、波阻抗反演以及正演模拟的基础工作,准确的地震子波估计对于高分辨率、高信噪比、高保真度的地震勘探数据处理具有极为重要的意义。

统计性子波提取方法的基本原理是首先对反射系数序列的分布做某种假设,然后利用地震记录的统计信息进行子波估计。

在没有任何先验知识的情况下,通常假设反射系数序列为一个非高斯、平稳和统计独立的随机过程,假设子波为一个非因果、非最小相位系统,加性噪声为高斯色噪声。

因此在利用地震记录的统计信息进行子波估计时,其高阶累积量不仅能保留系统的相位信息,而且能较好地压制高斯色噪声,显示出此法的优越性。

近年来,基于高阶累积量的参数化子波估计方法得到了快速发展。

Lazear [1]首先引入滑动平均(MA )模型描述地震记录,然后将子波四阶矩和地震资料的四阶累积量在最小均方误差意义下进行拟合,并用梯度下降法求解目标函数。

随后,Velis 等人[2]及尹成等人[3]试图应用特性更好的全局最优化方法解拟合函数,但求解效率普遍较低。

石殿祥等人[4]基于高阶累积量研究了非最小相位子波提取问题,虽取得了一定的成果,但依然沿用了滑动平均(MA )模型来描述地震记录。

本文分别采用滑动平均(MA )模型和自回归滑动平均(ARMA )模型来描述地震记录,并借助基于高阶累积量的线性化参数估计方法———矩阵方程法求解模型参数,最终精确估计了地震子波。

2　地震记录的滑动平均(MA)模型描述及矩阵方程法子波提取 地震记录y (n )可视为一个零均值的平稳随机过程,且符合如下褶积模型 y (n )=∑qi =0w (i )r (n -i )+v (n )=w (n )3r (n )+v (n )(1)式中:w (n )为地震子波;r (n )为反射系数序列;v (n )为环境噪声。

显然,式(1)符合典型的滑动平均(MA )模型表达式,因此可以把地震记录看作是有限脉冲响应(FIR )系统的含噪输出。

对于上述模型有如下假设:　第41卷　第5期 戴永寿等:高阶统计量地震子波估计建模 (1)反射系数序列r (n )是零均值、独立同分布的非高斯过程,其方差σw 2<∞,且至少存在一k 阶累积量|γkw |<∞;(2)环境噪声v (n )是零均值、因果的加性高斯色噪声,且与r (n )统计独立;(3)地震子波w (n )可以是非因果、非最小相位的,其非因果性表征了检波器或信道引入了失真。

通过上述设定可把褶积模型中子波的求解转化为对一个有限脉冲响应(FIR )系统进行基于高阶累积量的滑动平均(MA )参数的估计。

通过寻找信号累积量C kx (m ,n )和系统响应的线性关系,建立滑动平均(MA )模型参数w (i )与输出信号y (n )高阶累积量的一套矩阵方程,求解该方程组可得到有限脉冲响应(FIR )系统的冲击响应参数,即为子波参数。

根据Zhang 等人的累积量算法[5]中关于累积量切片的表达式∑qi =0h (i )C k-1kx (q ,i -t )=C kx (t ,0)C k-3kx (q ,0)C kx (q ,q )(2)∑qi =0hk-1(i )C kx (q ,i +t )=C kx (t ,0)C kx (q ,q )C kx (q ,0)(3)其中:t =-q ,…,-1,0,1,…,q ,q 为地震子波的长度;k 为选择的累积量的阶数。

并采用总体最小二乘法求解式(2)或式(3)即可得到子波估计。

3　地震记录的自回归滑动平均(AR MA)模型描述及矩阵方程法子波提取 从褶积模型出发,地震记录y (n )还可由如下等式表示y (n )=x (n )+v (n )(4)其中x (n )为理想地震记录,且满足差分方程∑pi =0a (i )x (n -i )=∑qj =0b (j )r (n -j )(5)的非高斯信号,即理想地震记录x (n )可以看作是反射系数序列r (n )通过一个非因果、最小相位的自回归滑动平均(A RMA )系统得到的输出,子波就是该系统的脉冲响应。

式(5)中:a (i )为自回归滑动平均(A RMA )模型的自回归参数;b (j )为自回归滑动平均(ARMA )模型的滑动平均参数。

v (n )是满足∑p vi =0c (i )v (n -i )=∑q vj =0d (j )e (n -j )(6)的有色高斯观测噪声。

其中:c (i )为自回归滑动平均(A RMA )模型的自回归参数;d (i )为此模型的滑动平均参数;e (n )为高斯白噪声。

信号除满足前述假设条件(1)、(2)之外,还应满足如下假设:(1)与式(5)对应的ARMA (p ,q )模型为非因果、非最小相位系统,其中a (0)=1,其传递函数为H (z )=∑∞j =-∞h (i )z -j=B (z )A (z )=∑qj =0b (j )z-j1+∑p i =1a (i )z-i(7)x (n )为稳定信号,所以|z |=1时,A (z )≠0,并且H (z )中无零极点对消;(2)H (z )有p 1个圆内极点,有p 2=p -p 1个圆外极点。

令αj 为H (z )的极点,则有A (z )=A 1(z )A 2(z )g (z )式中A 1(z )=∏p 1j =1(1-αjz-1)=1+∑p 1j =1a1(j )z -1,|αj |<1(8)称之为H (z )的因果自回归(AR )部分,而 A 2(z )=∏p 2j =1(1-βjz )=1+∑p 2j =1a2(j )z jβj =1αj +p 1,|αj +p 1|>1g (z )=z-p 2∏p 2j =1(-αj +p 1)(9)称A 2(z )g (z )=∏p 2j =1(1-αj +p 1z-1)为H (z )的反因果自回归(AR )部分;(3)H (z )有q 1个圆内零点,有q 2=q -q 1个圆外零点。

则与极点情况类似有B (z )=B 1(z )B 2(z )s (z )其中:B 1(z )为最小相位部分;B 2(z )s (z )为最大相位部分。

从上述模型假设可以看出,子波被视为式(5)中ARMA (p ,q )系统的冲激响应,条件((1)～(3))充分考虑到了子波的非因果、非最小相位特性,从而将子波提取转化为非高斯信号激励一个盲非因果、非最小相位自回归滑动平均(ARMA )系统而产生的信号伴有加性高斯色噪声的系统参数辨识问题。

鉴于高阶累积量不仅可以保留自回归滑动平均515 石油地球物理勘探2006年　(A RMA)系统的相位信息,而且可以有效地压制高斯色噪声,所以基于高阶累积量的方法可以用于非因果、非最小相位ARMA系统的建模与参数估计。

Tugnait最早进行了基于高阶累积量的非因果自回归(AR)系统的参数估计方法研究[6]。

在此基础上,G iannakis进一步研究了非因果、非最小相位的ARMA系统,提出了一种基于高阶累积量的线性化参数辨识方法[7]。

通过对上述方法的总结和发展,逐渐形成了一种成熟的三步识别方法[8]。

笔者将上述研究成果应用于子波的估计,求解过程分为三步:(1)求功率谱等价意义下的因果最小相位模型的自回归(A R)参数,由式(4)～式(6)可以得到Y(z)=B(z)A(z)R(z)+D(z)C(z)E(z)(10)其中:Y(z)、R(z)、C(z)、D(z)、E(z)分别是y(n)、r(n)、c(i)、d(j)、e(n)的Z变换,且A(z)=A1(z)××A2(z)g(z)。

由式(10)可知地震记录y(n)的功率谱为P(z)=B(z)B(z-1)A(z)A(z-1)σ2r+D(z)D(z-1)C(z)C(z-1)σ2e=B(z)B(z-1)C(z)C(z-1)σ2r+A(z)A(z-1)D(z)D(z-1)σ2eA(z)A(z-1)C(z)C(z-1)(11) 由于A2(z)的根均在单位圆外,将A2(z)所对应的自回归滑动平均(A RMA)模型的极点取共轭倒数,反演到单位圆内得到 A2(z),则 A2(z)的根均在单位圆内。

因 A2(z)与A2(z)的功率谱等价,现构造 A(z)=A1(z) A2(z),则A1(z)A2(z)与 A(z)功率谱等价。

由功率谱等价原理,并利用Cadzow自相关法估计自回归(AR)参数的奇异值分解—总体最小二乘(SVD2TLS)算法可以估计出A(z)C(z)=A1(z) A2(z)C(z)基于地震记录自相关估计因果信号自回归(A R)参数的算法[9]已经比较成熟,且数值模拟结果证明其具有很高的分辨率和精度。

(2)求H(z)的反因果自回归(AR)部分的参数利用估计出的A1(z) A2(z)C(z)滤波地震记录y(n)得到y1(n),其Z变换为 Y1(z)=A1(z) A2(z)C(z)Y(z)=B(z) A2(z)C(z)A2(z)g(z)R(z)++A1(z) A2(z)D(z)E(z)(12) 可见y1(n)是一反因果信号,上式等号右边第二项为一高斯信号,其高阶累积量为零。

用基于地震记录样本高阶累积量估计反因果系统自回归(A R)参数的奇异值分解—总体最小二乘(SVD2 TL S)算法[8]可估计出反因果自回归(A R)部分A2(z)g(z)的参数 a(1), a(2),…, a(p2)。