基于纯相位滤波器的子波相位校正方法研究

张亚南，戴永寿，丁进杰，张漫漫，王蓉蓉

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

传统的地震反褶积处理常常将地震子波假设为因果最小相位的，但在实际应用过程中，地震波在地层中传播受各种复杂地质因素的影响，子波的最小相位假设很难满足;且由于接收检波器等仪器性能的限制将非因果成分引入到地震子波中，地震子波的因果性也不能严格保证。故传统反褶积结果中往往残留子波相位成分，从而导致分辨率无显著提高或者地层假象的出现。

为解决最小相位反褶积后相位残余问题，许多学者假设子波振幅已经准确提出，对反褶积后的剖面应用相位校正方法。周兴元[1]在最小相位反褶积处理后的剖面中对子波剩余相位进行常相位校正，一定程度上提高了剖面的信噪比和分辨率。郭向宇等[2]提出了一种估算剩余混合相位子波相位的算法，对地震记录中的子波进行相位校正，可使子波接近或达到零相位，从而消除剩余相位的影响，提高了地震记录的分辨率。王君等[3]提出了一种在最小相位反褶积后估算子波相位补偿因子的方法，并对地震记录进行相位校正，可使剩余相位接近或达到零相位，从而提高资料的分辨率。但以上方法由于假设条件的苛刻性、模型描述的局限性和选取的补偿因子的特殊性，其效果并不十分理想，且不能涵盖所有因果和相位性的地震子波相位残余。

我们采用自回归滑动平均(ARMA)模型和单位化振幅谱的方式来构造纯相位滤波器，以描述任意因果和任意特性的地震子波相位残余，并在最大方差模准则约束下，采用改进的粒子群算法对子波相位残余进行非线性寻优，最终实现子波相位的精确校正。通过正演模拟数据和实际地震资料处理结果验证了所提方法的有效性和实用性。

1 纯相位滤波器的构造

在子波提取的过程中往往假设子波为最小相位，故反褶积结果中子波残余相位有可能是非因果的。为准确消除子波残余相位影响，我们采用单位化ARMA模型振幅谱构造纯相位滤波器，以描述任意特性的子波残余相位。相对于其它参数化模型，ARMA模型具有参数吝啬的特点，能够采用较少的参数描述一个精确的滤波器。ARMA模型可以表示为

(1)

式中：x(n)为滤波器输出;e(n)为滤波器输入信号;自回归(AR)部分的参数ai和滑动平均(MA)部分的参数bk均为实数;p为AR阶数，q为MA阶数，不失一般性，令a0=b0=1。若将地震子波视为一个滤波器系统，Robinson褶积模型[4]也可由(1)式表示。

假设ARMA描述下滤波器系统z域变换的系统函数W(z)无零极点对消，则

(2)

式中：w(m)为滤波器响应的时间域序列；W(z)为其z域表示。由于ai和bk均为实系数，因此系统零点ck与极点di均为实根或共轭复根。

(2)式所示的ARMA滤波器模型中，当零点全部在单位圆内时为零相位滤波器，内外兼有为混合相位滤波器，全部在单位圆上及单位圆外时为最大相位滤波器;当极点全部在单位圆内时为因果滤波器，内外兼有为混合因果滤波器，全部在单位圆上及单位圆外时为反因果滤波器。故ARMA模型可以描述任意因果特性、任意相位特性的滤波器系统。

令z=ejω，(2)式中滤波器的z域表示转换为频率响应，即

(3)

纯相位滤波器的时域为纯相位序列，舒立华[5]在其研究纯相位序列的能量传递性质的过程中，给出了纯相位序列的验证方法：一个序列g(n)为纯相位序列的充要条件是其z变换G(z)z=ejω能表示成ejφ(ω)的形式。则根据(3)式可构造纯相位滤波器：

(4)

由(4)式可以看出，任意滤波器均可以通过频率谱单位化来构造纯相位滤波器。为避免计算错误，当(4)式分母中某些频率成分为零值时，相除结果规定为单位值1。

相对于其它模型，由ARMA模型单位化振幅谱构造出的纯相位滤波器可以描述任意特性的子波相位残余，包括最小相位、混合相位、最大相位，以及因果、混合因果、反因果等。

2 子波相位校正

反褶积结果中的相位残余可表示为

(5)

式中：D(ejω)和R(ejω)分别为反褶积结果和反射系数序列的频域表示;φR(ω)为子波残余的相位谱。采用(4)式中构造的纯相位滤波器对(5)式进行相位校正处理，

(6)

式中：y(n)和Y(ejω)分别为相位校正结果的时间域和频率域表示;F-1为逆傅里叶变换算子。当构造的纯相位滤波器和子波相位残余具有相同相位谱时，(6)式可表示为

(7)

此时地震子波被校正为零相位，反褶积结果被校正为y0(n)，即反射系数序列r(n)。

若反射系数序列为稀疏序列，则相位校正结果具有最大的纵向分辨率，本文采用输出结果的最大方差模准则[6]来检验相位校正效果。最大方差模准则表示为

(8)

其中，L为相位校正输出的数据长度。采用最大方差模准则Ψ(y)对相位校正结果y(n)进行评价，当Ψ(y)取最大值时，可以认为构造出的纯相位滤波器的相位谱最接近于子波残余的相位谱。

对目标函数Ψ(y)求取最大值的过程是一个多参数、多极值的非线性优化问题，对该目标函数应用的优化算法需具有参数向量整体的全局随机搜索能力和对单个参数的深度搜索能力。杨文采[7]、师学明等[8-9]针对地球物理学中的非线性寻优问题引入了多种非线性优化算法，并对此类算法进行了分析和探讨[10-11]。在此我们采用改进的粒子群算法对参数进行寻优，改进的粒子群算法继承了基本粒子群算法的优点，并具有全局寻优能力强、收敛速度快、可避免陷入局部最优值和自适应学习的优点[12]。

基于纯相位滤波器的子波相位校正方法的实现步骤为：

1) 采用ARMA模型描述子波残余，并且通过相关函数法确定子波残余的AR阶数p，采用信息量准则和高阶累积量相结合的方法确定子波残余的MA阶数q[13];

2) 初始化AR参数和MA参数;

3) 通过(4)式构造纯相位滤波器，通过(6)式进行相位校正;

4) 通过(8)式的最大方差模准则对相位校正结果进行评价，将评价值与历史最优进行比较，若优于历史最优，则将当前参数和评价值设定为历史最优;

5) 判断是否达到最大迭代次数，若达到则转向6);若未达到，则粒子飞翔，生成新的参数，转向3);

6) 读取历史最优参数，算法结束。

3 正演模拟数据测试

为验证本文所提方法的有效性，首先对合成地震记录进行子波相位校正测试。合成地震记录中反射系数序列为独立同分布(IID)的随机过程，且服从Bernoulli-Gaussian分布，序列采样间隔为1ms，长度为1000ms，如图1a所示。为了体现ARMA模型描述下子波的一般性，选用混合因果混合相位的地震子波合成地震记录。合成地震记录所采用的原始地震子波在ARMA模型描述下差分方程形式为

x(t)-4.02x(t-1)+8.43x(t-2)-8.15x(t-3)+2.86x(t-4)=

r(t)-0.8r(t-1)+0.2r(t-2)-0.82r(t-3)

(11)

其z域系统函数形式为

(12)

子波的时间域波形、合成地震记录和子波相位谱如图1b,图1c和图1d所示，其中相位谱的取值范围为(-π,π]。

针对合成地震记录进行地震子波提取，在模拟过程中构造与原始合成地震子波具有相同振幅谱的因果最小相位地震子波进行反褶积处理。反褶积处理所用最小相位地震子波时间域波形和相位谱如图2a和图2b所示。应用构造的子波进行反褶积处理，反褶积结果如图2c所示，反褶积结果中所含有的子波相位残余如图2d所示。

采用本文所提方法构造纯相位滤波器，并采用最大方差模准则进行约束，通过改进的粒子群算法针对图2c中的反褶积结果进行相位校正，相位校正结果如图3所示。

图1 合成地震记录a 反射系数序列; b 原始地震子波; c 合成的地震记录; d 原始地震子波相位谱

图2 子波估计与反褶积结果a 构造的因果最小相位地震子波时域波形; b 构造的因果最小相位地震子波相位谱; c 反褶积结果; d 子波相位残余

图3 子波相位校正结果a 提取的子波相位残余; b 子波相位校正结果

对比图2d与图3a可以看出，提取的子波相位残余和理论子波相位残余略有差异，但整体相似度较高，可以认为本文所提方法能够较为准确地提取出子波相位残余。对比图3b与图1a可以看出，由于提取的子波相位残余与理论子波相位残余的差异，子波相位校正结果与原始反射系数序列不能完全吻合，但二者仅在细节方面存在差异，整体吻合度较高，子波相位校正结果能够准确反映出原始反射系数序列的各个反射界面。模拟数据测试结果表明本文所提方法能够实现子波相位的校正，可以较为准确地恢复出原始的反射系数序列，证明该方法是有效可行的。

4 实际地震资料处理

为进一步验证本文所提方法对子波相位校正和提高地震剖面分辨率的作用，应用该方法对实际地震资料进行相位校正处理。图4为胜利油田某区块经反褶积处理后的地震剖面，图5为应用本文方法进行相位校正后的反褶积剖面。

由图4和图5对比可以看出，经相位校正后的剖面质量得到明显改善：同相轴连续性变好，一些弱反射变得清晰。可以认为，本文所提基于纯相位滤波器的子波相位校正方法可以实现提高反褶积后剖面分辨率的作用，具有一定的实际应用价值。

图4 胜利油田某区块经反褶积处理后的剖面

图5 胜利油田某区块经本文方法相位校正后的反褶积剖面

5 结束语

我们提出了一种基于纯相位滤波器的地震子波相位校正方法，通过ARMA模型和振幅谱单位化的方式构造纯相位滤波器以描述地震子波相位残余，在最大方差模准则约束下采用改进的粒子群算法针对残余的子波相位进行准确寻优，以实现反褶积结果中子波残余相位的校正。正演模拟数据和实际地震资料的应用结果表明该方法是有效的，并具有一定的实际应用价值。

将反褶积中子波残余相位校正问题转化为信号的智能寻优问题是提高地震记录分辨率的一种新途径，下一步研究重点是进一步提高相位谱的寻优精度。

参 考 文 献

[1] 周兴元.常相位校正[J].石油地球物理勘探,1989,24(2):119-129

Zhou X Y.Constant phase correction [J].Oil Geophysical Prospecting,1989,24(2):119-129

[2] 郭向宇,周兴元,董敏煜.混合相位子波的相位估算及校正[J].石油地球物理勘探,1998,33(2):214-221

Guo X Y,Zhou X Y,Dong M Y.Estimation and correction of mixed-phase wavelet phase [J].Oil Geophysical Prospecting,1998,33(2):214-221

[3] 王君,周兴元.子波相位校正及效果分析[J].石油地球物理勘探,2008,43(增刊2):142-145

Wang J.Zhou X Y.Wavelet phase correction and analysis of effect [J].Oil Geophysical Prospecting,2008,43 (S2):142-145

[4] Robinson E A.Predictive decomposition of time series with application to seismic exploration [J].Geophysics,1967,32(3):418-484

[5] 舒立华.纯相位序列的能量传递性质[J].数学学报,1974,17(1):20-27

Shu L H.The energy transfer properties of phase-only sequence[J].Acta Mathematica Sinica,1974,17(1):20-27

[6] 王有新,周兴元.最小熵反褶积和集约化准则[J].石油地球物理勘探,1993,28(6):678-684

Wang Y X,Zhou X Y.Minimum entropy deconvolution and parsimonious criterion [J].Oil Geophysical Prospecting,1993,28(6):678-684

[7] 杨文采.地球物理反演的遗传算法[J].石油物探,1995,34(1):116-122

Yang W C.Genetic algorithm for geophysical inversion [J].Geophysical Prospecting for Petroleum,1995,34(1):116-122

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

Shi X M,Wang J Y,Yi Y Y,et a1.A study on the simulated atomic transition algorithm for geophysical inversion [J].Chinese Journal of Geophysics,2007,50(1):305-312

[9] 师学明,王家映,张胜业,等.多尺度逐次逼近遗传算法反演大地电磁资料[J].地球物理学报,2000,43(1):122-130

Shi X M,Wang J Y,Zhang S Y,et a1.Multi-scale Genetic Algorithm and its application in magnetotelluric sounding data inversion [J].Chinese Journal of Geophysics,2000,43(1):122-130

[10] 杨文采.非线性地球物理反演方法:回顾与展望[J].地球物理学进展,2002,17(2):255-261

Yang W C.Non linear geophysical inversion methods:review and perspective [J].Process in Geophysics,2002,17(2):255-261

[11] 王家映.地球物理资料非线性反演方法讲座(一):地球物理反演问题概述[J].工程地球物理学报,2007,4(1):1-3

Wang J Y.Lecture on non-linear inverse methods in geophysics(1):introduction to geophysical inverse problems [J].Chinese Journal of Engineering Geophysics,2007,4(1):1-3

[12] 戴永寿,牛慧,彭星,等.基于自回归滑动平均模型和粒子群算法的地震子波提取[J].中国石油大学学报(自然科学版),2011,35(3):47-50

Dai Y S,Niu H,Peng X,et al.Seismic wavelet extraction based on auto-regressive and moving average model and particle swarm optimization [J].Journal Of China University of Petroleum (Edition Of Natural Science),2011,35(3):47-50

[13] 张亚南,戴永寿,王少水,等.高效ARMA模型高分辨率地震子波提取方法[J].石油地球物理勘探,2011,46(5):686-694

Zhang Y N,Dai Y S,Wang S S,et al.High resolution wavelet estimation by ARMA modeling [J].Oil Geophysical Prospecting,2011,46(5):686-694