基于自回归滑动平均模型和粒子群算法的地震子波提取

戴永寿，牛 慧，彭 星，王少水

(中国石油大学信息与控制工程学院，山东东营 257061)

基于自回归滑动平均模型和粒子群算法的地震子波提取

戴永寿，牛 慧，彭 星，王少水

(中国石油大学信息与控制工程学院，山东东营 257061)

基于自回归滑动平均(ARMA)模型理论，对地震子波进行参数化建模，采用累积量拟合法精确估计参数，使地震子波提取问题最终归结为一个多参数、多极值的非线性函数优化问题。对基本粒子群算法进行改进，通过自适应参数调整和边界约束，克服基本粒子群算法易陷入局部极值的缺陷，同时提高算法寻优精度和计算效率。仿真数据试验结果验证了改进的粒子群算法在地震子波提取方法中的有效性和稳定性。

地震数据处理;自回归滑动平均模型;地震子波;系统辨识;累积量拟合;粒子群算法

地震子波提取是地震资料反褶积处理、波阻抗反演和正演模型的基础［1］。Lazear［2］对基于四阶累积量的混合相位子波提取方法进行研究，提出将子波的四阶矩和地震数据的四阶累积量进行拟合优化以提取子波，但由于参数模型的未知性，确定模型初始参数范围较为困难，严重影响算法的寻优效率。粒子群算法是基于群体智能的多峰值寻优算法，概念简单且收敛速度快，但算法寻优精度与参数选择有着密切的关系。笔者基于自回归滑动平均(ARMA)模型理论，将累积量拟合法和参数化辨识方法相结合，将子波提取问题最终归结为一个多参数、多极值的非线性函数优化问题［3］。

1 地震子波ARMA模型和累积量拟合函数

通常地震激励假设为一零均值平稳随机过程y(n)，用如下褶积模型［3］描述一道地震记录：

式中，h(n)为地震子波;lc和lnc分别为子波因果、非因果部分的长度;r(n)为发射系数序列;v(n)为加性噪声。模型(1)一般有如下假设：①反射系数序列r(n)为零均值、独立同分布的非高斯过程［4］，其方差＜∞，且至少存在一k阶累积量满足＜∞;②E(θ)≤4为环境噪声、仪器噪声及激发产生的多次反射噪声等的合成噪声信号，为一与r(n)统计独立的随机过程，一般为加性色噪声［5-6］，且其高斯成分远大于非高斯部分;③h(n)为非因果、混合相位的地震子波，其非因果性表征了检波器或信道引入的失真。

根据文献［4］，将地震子波的褶积模型转换为

其中，x(n)为平稳随机过程，为一ARMA模型在输入为反射系数序列r(n)时的响应，该ARMA模型满足以下差分方程：

式中，p和q分别为ARMA模型的AR(自回归)和MA(滑动平均)部分的阶数，模型对应冲激响应为式(1)中的地震子波h(n)。

由于奇数阶累积量恒为零，偶数阶累积量不为零，故应采用四阶累积量进行地震子波估计。对地震记录求四阶累积量，应用Bartlett-Brillinger-Rosenblatt公式，由假设条件结合累积量定义和性质，得

式中，C4y(t1，t2，t3) 为地震数据的四阶累积量;m4b(t1，t2，t3)为估计的地震子波四阶矩;r4r为反射系数序列峰度。r4r为一标量，在不影响所求取子波形态的情况下，可作为系数被m4b(t1，t2，t3)吸收，即地震记录的四阶累积量与子波的四阶矩函数之间仅差一个标量因子，可通过计算地震记录的四阶累积量来估算子波。由于高阶累积量关于其变元是对称的，因此在最小平方误差意义下可建立目标函数［5］

式中，L为地震子波的长度。目标函数φ为地震子波的函数，当φ达到最小值时，所对应的一组变量即为子波的估计值。此目标函数的求解是一个非线性多参数多极值的优化问题，需要应用非线性优化算法进行求解。通过对优化算法进行研究，提出了改进的粒子群算法，并将其运用到模型求解当中，以改善累积量拟合优化的计算效率和寻优精度。

2 改进的粒子群算法

粒子群算法［6］概念简明，实现方便，参数设置少和收敛速度快，但算法精度与参数选择有着密切的关系，不恰当的参数选择，易导致算法结果发散，寻优精度降低［7-8］。

基本粒子群算法的数学表达式如下：

式中，ω为惯性权重，控制前一代对下一代速度的影响程度;c1和c2为学习因子;r1和r2为在［0，1］范围内服从均匀分布的随机变量;t为迭代次数。

基本粒子群算法存在如下缺陷［9］：①算法精度同参数选择有着密切的关系，若参数选择不当，算法精度大大降低;②在整个实现过程中，参数均为定值，而不考虑具体的优化模型和迭代过程，影响了算法精度。为此，本文中提出一种改进的粒子群算法：令惯性权重ω、学习因子c1和c2随着迭代次数的增加线性递减或线性递加，从而保证算法初期个体能搜索整个空间而不陷入局部最优值，后期能朝着全局最优值收敛而找到全局最优值;为避免粒子位置超出给定位置，对粒子边界条件进行设置，保证寻优解在有效的范围内，从而应用到地震子波提取中。

2.1 惯性权重

Shi和 Eberhart［10］提出随着迭代次数增加，惯性权重线性递减的策略，表达如下：

式中，ωmax为最大惯性权重;ωmin为最小惯性权重;ts为算法迭代总次数;t*为算法当前迭代次数。

2.2 学习因子

c1和c2分别表征搜索局部极值和全局极值的能力。算法初期，较大的c1和较小的c2使粒子在整个可行空间进行搜索，保证地震子波初始范围的多样性;算法后期，较小的c1和较大的c2使粒子快速收敛于全局最优值，找到地震子波的准确解。根据分析，数学表达式如下：

式中，c1start和c2start为学习因子的初始值;c1end和c2end表示两个学习因子的终值。

2.3 粒子边界条件

当粒子的位置超出了给定位置时，重新设置粒子的当前位置能够使逃逸出去的粒子又回到种群中，而且其值也不再单一地确定为边界值，确保了如果粒子过多发生逃逸时新种群的多样性，保证地震子波提取的有效性。

3 地震子波提取的具体实现步骤

(1)数据生成。用待估计ARMA子波模型与满足反射系数序列假设的随机序列合成地震数据。

(2)参数预估计。运用奇异值分解(SVD)法、Zhang算法、基于奇异值分解的总体最小二乘法(SVD-TLS)等模型辨识方法确定待估计模型的阶数p、q和自回归参数、滑动平均参数。

(3)确定参数搜索空间。在初步估计模型参数的基础上，确定模型参数向量θ的搜索范围。

(4)拟合优化。在给定参数向量的搜索范围内用改进的粒子群算法对累积量拟合目标函数E(θ)进行参数精确估计，寻找最优解x。

(5)评价函数分析。计算最优解x的评价函数值，当其明显增加或降低时继续，否则转至步骤(7)。

(6)阶数扰动。根据拟合误差，调整上述方法的阈值，对模型的阶数进行扰动，生成一组新的p、q值，然后转至步骤(3)。

(7)结束。将当次扰动前所得阶数p、q及精确参数估计值x视为模型最优解。

4 合成地震记录处理

通过合成地震记录进行仿真试验来验证方法的有效性和实用性。合成地震道所需的地层反射系数序列为独立同分布(IID)的随机过程，且服从伯努力分布;地震子波为因果混合相位的，其ARMA模型的差分方程为

地震子波ARMA模型的零极点分布如图1所示。图中两个极点全部在单位圆内，根据子波分类标准可知，子波是因果的。同理，图中两个零点，一个在单位圆内，一个在单位圆外，子波为混合相位的，即子波为因果混合相位的(子波波形见图2)。

地震记录中的噪声一般假设为高斯有色噪声，若地震记录表示为y(t)=s(t)+v(t)，其中s(t)为有效波，v(t) 为噪声，定义噪信比 RNS［11］为

在频域定义噪信比为

式中，ω为地震记录在频域中的频率;Avv(w)和Ass(w)分别为噪声和有效波的振幅谱;Rvv(w)和Rss(w)分别为噪声和有效波的功率谱。本文中采用式(13)形式定义噪信比。

采用基于相关函数的SVD-TLS算法估计AR参数，同时采用累积量法对MA参数进行估计，建立初始子波寻优范围，并用改进的粒子群算法进行优化，在不同数据长度、噪声的情况下，对该记录运用本文方法进行子波提取，验证方法的有效性和实用性。

4.1 数据长度对子波ARMA模型参数估计的影响

每次试验均在噪信能量比为50%的情况下，随机生成长度分别为10、8、5和3 s的合成地震记录。

由表1可知，AR和MA参数在长数据情况下均能得到较好的估计结果，E(θ)也满足E(θ)≤4的标准，而在短数据(5 s和3 s)时估计结果与实际值有一定偏差，不过在不同数据长度情况下提取的子波波形具有较好的连续性，如图3(a)所示。

4.2 环境噪声对子波ARMA模型参数估计的影响

为考察环境噪声对模型参数估计方法的影响，每次试验均生成长度为5 s的合成地震记录，对该记录分别加入噪信能量比为20%、50%、100%、200%的高斯色噪声后，利用本文方法对子波模型进行参数估计，所得结果见表2，提取的子波波形见图3(b)。由表2和图3(b)可知，随着噪声强度的增加，模型参数估计误差略有增加，即使噪声与信号的能量比达到1∶1，模型参数估计仍得到很好的结果，验证了模型参数估计和粒子群算法既有很好的抗噪容噪能力，又能在适宜的噪声环境中应用，也说明了高阶累积量的适用性。

表1 不同数据长度下地震子波参数估计结果Table 1 Seismic wavelet parameter estimation results with different length

图3 数据长度和高斯色噪声对提取的子波波形的影响Fig.3 Effects of trace length and Gaussian noise on waveforms of extracted wavelets

表2 不同强度高斯色噪声环境下地震子波参数估计结果Table 2 Seismic wavelet parameter estimation results with different noise-signal-ratio

5 结论与建议

(1)用改进的粒子群算法进行寻优，假设地震子波为因果混合相位时的仿真结果表明，所提方法能够有效地确定模型参数，在不同数据长度和适宜噪声环境下仍能得到很好的结果。

(2)改进的粒子群算法提高了算法收敛的稳定性和收敛质量，但是迭代次数较高，一般迭代到80～100次时才会满足目标函数E(θ)≤4的条件。因此，如何在保证地震子波提取精度的同时，提高粒子群算法寻优的效率，将是以后研究的重点。

(3)实际地震记录是非因果、混合相位的，因此下一步的工作目标是采用改进的粒子群算法对非因果混合相位子波和实际地震记录进行优化，并与前期的优化算法，如遗传算法等优化结果进行比较，以得到较好的优化结果。

［1］ROBINSON E A，TREITEL Sven.Geophysical signal analysis［M］.Englewood Cliffs：Prentice-Hall，1980.

［2］LAZEAR G D.Mixed-phase wavelet estimation using fourthorder cumulant［J］.Geophysics，1993，58：1042-1051.

［3］师学明，王家映，易远元，等.一种新的地球物理反演方法——模拟原子跃迁反演法［J］.地球物理学报，2007，50(1)：305-312.

SHI Xue-ming，WANG Jia-ying，YI Yuan-yuan，et al.A study on the simulated atomic transition algorithm for geophysical inversion［J］.Chinese Journal Geophysics，2007，50(1)：305-312.

［4］ 戴永寿，王俊岭，王伟伟，等.基于高阶累积量ARMA模型线性非线性结合的地震子波提取方法研究［J］.地球物理学报，2008，51(6)：1851-1859.

DAI Yong-Shou，WANG Jun-ling，WANG Wei-wei，et al.Seismic wavelet extraction via cumulant-based ARMA model approach with linear and nonlinear combination［J］.Chinese Journal Geophysics，2008，51(6)：1851-1859.

［5］ADNAN Al-Smadi.Cumulant-based order selection of non-Gaussian autoregressive moving average models：the corner method［J］.Signal Processing，2005，85(3)：449-456.

［6］KENNEDY J，EBERHART C.Particle swarm optimization：proceeding of IEEE International Conference on Neural Networks［C］.Piscataway，NJ：IEEE，c1995：1942-1948.

［7］KENNEDY J，EBERHART R C.A new optimizer using particle swarm theory：proceedings of the Sixth International Symposium on Micro Machine and Human Science［C］.Nagoya：IEEE，c1995：1045-1049.

［8］NIU Hui，DAI Yongshou，PENG Xing.Particle swarm optimization with adaptive parameters and boundary constraints：proceedings of the 3rd International Conference on Networks Security，Wireless Communications and Trusted Computing［C］.Wuhan：IEEE，c2011：1082-1085.

［9］CLERC M，KENNEDY J.The particle swarm-explosion，stability，and convergence in a multidimensional complex space［J］.IEEE Transactions on Evolutionary Computation，2002，6(1)：58-73

［10］KANNAN S，SLOCHANAL S M R，SUBBARAJ P，et al.Application of particle swarm optimization technique and its variants to generation expansion planning problem［J］.Electric Power Systems Research，2004，70(3)：203-210.

［11］YIN P Y.A discrete particle swarm algorithm for optimal polygonal approximation of digital curves［J］.Journal of Visual Communication and Image Representation，2004，15(2)：241-260.

Seismic wavelet extraction based on auto-regressive and moving average model and particle swarm optimization

DAI Yong-shou，NIU Hui，PENG Xing，WANG Shao-shui
(College of Information and Control Engineering in China University of Petroleum，Dongying 257061，China)

A seismic wavelet parametric model was developed based on auto-regressive and moving average(ARMA)model theory.The model parameters were accurately determined based on cumulant fitting method.So the seismic wavelet can be a multi-parameters，multi-extremes nonlinear functional optimization problem.An improved particle swarm optimization with adaptive parameters and boundary constraints was proposed for the local extreme value defects of elementary particle swarm optimization.The optimization accuracy and computation efficiency are also improved.Simulation results show that the method has good applicability and stability in seismic wavelet extraction.

seismic data processing;auto-regressive and moving average(ARMA)model;seismic wavelet;system identification;cumulant fitting;particle swarm optimization

P 631

A

10.3969/j.issn.1673-5005.2011.03.009

1673-5005(2011)03-0047-04

2011－04－02

国家自然科学基金项目(40974072);山东省自然科学基金项目(ZR2010DM14)

戴永寿(1963－)，男(汉族)，安徽巢湖人，教授，博士，博士生导师，从事信号与信息处理和计算机测控技术等教学研究工作。

(编辑 修荣荣)