2013, Vol. 56
地震子波提取是地震资料处理和地震解释的重要组成部分，但子波提取过程中子波相位往往估计不准确，严重影响地震反演处理和地震解释.针对地震子波相位提取的问题，相关领域的学者做了大量研究.唐斌等[1]将子波相位分成最小相位和最大相位两个部分，通过高阶累积量方法将二者恢复，从而提取子波相位.伊振林等[2]利用延迟因子构建混合相位滤波器进行混合相位反褶积，通过最小熵准则和Lp范数约束反褶积结果，进而获取最佳滤波器和子波相位.崔庆辉等[3]通过地震道自相关提取子波振幅，利用双谱法提取子波相位，最终重构出混合相位的地震子波.但在实际应用过程中，以上方法很难得到精确结果，以至于造成反演结果的失真.为改善反演结果，国内外科研人员提出了多种相位校正方法.Levy等[4-5]假设提取的子波与真实子波的差异为不依赖于频率的常数，采用常相位校正方法对反演结果进行校正，以消除子波相位估计不准的影响. Wang[6-7]提出了一种基于波场下延理论的常Q反滤波方法，通过两项指数对大地Q滤波的相位和振幅进行校正，并给出其稳定性的证明.到目前为止，人们提出的各种相位校正方法并不能完全有效地消除反演结果中子波相位估计不准的影响，因此子波相位对反演影响的研究将成为今后的攻关课题和研究热点.Yuan等[8]针对滑动平均（Movingaverage，MA）模型描述下地震子波相位估计不准对反演结果影响的问题进行了讨论，并给出了多种方法评价反演结果.但由于其模型所含参数较多，造成多种相位可能性；且所用子波的长度过短，与实际子波相比具有一定差异[9].因此，本文在ARMA模型描述子波的基础上，忽略噪声等因素的影响，重点探讨子波相位对反射系数序列反演结果的影响，并通过理论分析和仿真实验研究其影响规律.
本文首先采用在z域对称映射ARMA模型零极点的方式构造出一系列相同振幅谱、不同相位谱的子波，进而对人工合成地震记录进行反射系数序列反演并对结果进行分析，通过仿真实验和实际地震资料处理来验证理论分析结果.本文采用两种方法对反演结果进行评价，以确定出真实或准确的反射系数序列.
2 构造不同相位谱的地震子波戴永寿等[10]将ARMA模型引入到地震子波的估计当中，ARMA模型相对于MA模型具有参数吝啬的性质，用较少的参数即可描述一个精确的子波.Robinson褶积模型[11]可用ARMA模型表示为
|
(1) |
其中，x（n）为地震记录；r（n）为反射系数序列；子波AR部分的参数ai和MA部分的参数bk均为实数；p为AR阶数，q为MA阶数.
假设ARMA描述下地震子波z域变换的系统函数W（z）无零极点对消，则
|
(2) |
其中w（m）为地震子波的时间域序列，W（z）为其z域表示.由于ai和bk均为实系数，因此系统零点ck与极点di均为实根或者共轭复根[12]，A为增益常数.
当z=ejω时，式（2）中子波的z域表示转换为频率响应，即
|
(3) |
其中
|
(4) |
对Ck（ejω）求取振幅谱Ck（ω）和相位谱φCk（ω），可得
|
(5) |
|
(6) |
同理，令di=si+jli，由式（5）和式（6）可得子波ARMA模型的振幅谱和相位谱
|
(7) |
|
(8) |
子波系统函数零极点分布决定着子波的因果特性和相位特性，采用对称映射其零极点方法可以改变其相位特性和因果特性，从而构造出一系列具有相同振幅谱的、涵盖各种不同因果性和相位特性的子波.将式（2）中的Ck（z）=1-ckz-1由C′k（z）=
|
(9) |
其振幅谱和相位谱为：
|
(10) |
|
(11) |
由式（5）、式（6）、式（10）和式（11）可知，当Ck（z）被对称映射到C′k（z）时，地震子波的振幅谱不变，当|ck|=1（即xk2
+yk
2=1）时，Ck（z）与C′k（z）完全相同.同理，当Di（z）被对称映射到D′i（z）=
在已经构造出一系列相同振幅谱、不同相位谱子波的基础上研究子波相位不准对反射系数反演的影响，首先需要找寻一个合适的反演方法.现有的反演方法已经十分成熟，如最小二乘反褶积、最小熵反褶积[13]、稀疏脉冲反褶积[14]、多分辨率地震信号反褶积[15]、基追踪反演[16]等.由于ARMA模型描述下的子波在时域表现为无限长脉冲响应序列，时域反褶积处理不可避免地要对子波进行截断.为减少截断误差对反演过程影响，我们选择谱除法在频域进行反褶积处理，处理结束后将频域结果映射到时域进行比较.
将Robinson褶积模型经傅里叶变换可得
|
(12) |
其中X（ejω）、W（ejω）和R（ejω）分别为地震记录、地震子波和反射系数序列的频域表示，则反射系数序列的估计可表示为
|
(13) |
其中
由于所估计的地震子波和真实地震子波具有相同的振幅谱，褶积过程在频域内应为频域相乘的形式；反褶积过程去除子波的影响，为频域相除的形式.W（ejω）与
|
(14) |
其中
由式（14）可知，子波相位估计不准的反褶积结果中会残留一个纯相位滤波器，即反射系数序列反演的结果为真实反射系数序列与一个纯相位滤波器的褶积，且此纯相位滤波器的相位谱为真实子波和估计子波相位谱之差.估计的反射系数序列应当和真实反射系数序列具有相同的振幅谱和能量，仅相位谱存在差异.由于纯相位滤波器的影响，估计的反射系数序列出现了相位偏移，不符合真实反射系数序列的分布规律，有可能造成假象的出现和分辨率的降低[8].
4 数据仿真验证为了验证ARMA模型描述下子波相位估计不准对反射系数反演结果的影响，本文使用具有四个反射系数脉冲的稀疏模型进行实验，其长度为1000ms，采样频率为1kHz，如图 1所示.其中稀疏反射系数脉冲的位置和大小分别为（200ms，1.0）、（300 ms，-0.8）、（600ms，0.8）和（800ms，-1.0）.
|
图 1 反射系数序列模型的（a）时域波形图和（b）振幅谱 Fig. 1 The time domain waveform (a) and amplitude spectrum (b) of reflection coefficient sequence model |
本实验中采用因果混合相位模型的真实子波，其ARMA描述下的差分形式为
|
(15) |
其z域系统函数形式为
|
(16) |
分别求取式（16）中AR和MA部分的根可得真实子波的极点：α1=0.8643+0.1666j，α2=0.8643-0.1666j，α3=0.3107 + 0.4177j，α4=0.3107 -0.4177j；真实子波的零点：β1=1.2015，β2=-0.2008+0.8013j，β3=-0.2008-0.8013j.子波的时域和频域表示如图 2所示，其中时域波形进行了截断处理.文中相位谱中令相位角φ ∈（-π，π].
|
图 2 真实子波的（a）时域波形图，（b）振幅谱和（c）相位谱 Fig. 2 The time domain waveform (a), amplitude spectrum (b) and phase spectrum (c) of real wavelet |
将真实子波与反射系数序列作褶积处理合成地震记录，合成的地震记录如图 3所示.由于所采用的子波为ARMA模型描述，在时间域具有很强的延续性，所以图 3a中会出现较长的拖尾现象.图 3b与图 1相比，合成地震记录的能量大部分集中在了低频区域.
|
图 3 合成的地震记录（a）时域波形图和（b）振幅谱 Fig. 3 The time domain waveform (a) and amplitude spectrum (b) of synthetic seismogram |
按照第2节中所介绍的方法构造出一系列相同振幅谱、不同相位谱的子波，其z域中零极点的分布如图 4所示.其中“○”代表零点，“×”代表极点.第一行（4a-4d）、第二行和第三行（4e-4l）、第四行（4m-4p）分别为因果、混合因果、反因果子波的零极点分布，第一列、第二列和第三列、第四列分别为最小相位、混合相位、最大相位子波的零极点分布.图 4b为真实的因果混合相位子波在z域的零极点分布.
|
图 4 构造出的不同相位谱子波z域零极点分布图 Fig. 4 Pole-Zeros of constructed seismic wavelets with different phase spectrums |
图 5所示为相同振幅谱不同相位谱的地震子波的时域序列和相位谱.由图 5可知，时间域各道子波具有相同振幅谱，由于相位谱的不同，在时域中表现出了不同的性质.图 5b中的每道波形为图 5a中相应地震子波的相位谱，相位角φ∈（-π，π].由图 5a可以看出，1-4道子波为因果子波，其时间域响应均在零点之后；5-12道子波为混合因果子波，在零点前后均有时间域响应；13-16道子波为反因果子波，其时间域响应均在零点之前，其中第2道为真实子波.在实际的地震勘探过程中，子波应当是因果混合相位的，但由于检波器等原因会引入少量的反因果成分，故在实际处理过程中不能单纯地将子波看作仅含因果成分.由图 4和图 5a可以看出，在因果稳定系统的情况下，零点均在单位圆内为最小相位，内外均有为混合相位，均在单位圆外为最大相位.当系统为反因果时则与因果系统完全相反，混合因果时需要分情况讨论，在此不再赘述.
|
图 5 构造出的相同振幅谱不同相位谱子波的（a）时域波形图和（b）相位谱 Fig. 5 The time domain waveform (a) and phase spectrum (b) of constructed seismic wavelets with the same amplitude spectrum and different phase spectrums |
应用构造出的相同振幅谱不同相位谱子波对合成的地震记录进行反射系数序列反演，即反褶积处理.16个不同相位谱子波反褶积处理后的结果如图 6所示，其中第2道为真实子波反演结果.
由图 6可以看出，仅真实子波能够完全恢复真实的反射系数序列，其余相位不准的子波反演结果均在时域呈现出反射系数拓展的形式，不满足真实反射系数序列的稀疏性质.如果本实验结果满足第3节中的结论---在子波相位估计不准的情况下反射系数序列反演的结果为反射系数序列与一个纯相位滤波器的褶积，则对所有估计子波进行反射系数序列反演的结果应当具有相同的振幅谱，即相同的能量.根据能量在时域和频域等价原则，计算图 6中所有反射系数序列的时域和频域能量，其结果如表 1所示，第2道为真实子波反演的能量值.
|
图 6 不同相位谱子波反射系数反演结果 Fig. 6 The reflection coefficient sequences inversion results of wavelets with different phase spectrums |
由表 1可知，不同相位子波的反演结果所含能量略有差异，其根本原因在于采用的ARMA模型为无限长脉冲响应系统，由于实验用反射系数序列的长度为1000 ms，在实际处理过程中进行了截断处理，从而产生截断误差；数值计算存在精度误差.但估计子波与真实子波反演结果的能量值最大差异不超过0.031%，属可接受范围，可认为估计子波的反演结果具有相同的能量，可初步验证文中第3节中的结论.
|
表 1 不同子波反演结果的能量值 Table 1 The energy inversion results with different wavelets |
为进一步验证第3节中的结论，我们将求取出纯相位滤波器.由式（14）可知，纯相位滤波器P(ejω)=
|
(17) |
若上式中分母中所估计的子波频率响应的某些频率为零值时，将其用极小实数值代替以避免计算错误.
纯相位滤波器应当具有相同的单位振幅谱，并且在时域和频域应当具有单位能量.图 7所示为由式（17）所计算出的纯相位滤波器，每道中各个频率的幅值均为1，图 7a为其时间域波形，图 7b为其相位谱.对图 7中所示的16个纯相位滤波器进行时域和频域能量的求取，结果如表 2所示.
|
图 7 纯相位滤波器的（a）时域波形图和（b）相位谱 Fig. 7 The time domain waveform (a) and phase spectrum (b) of pure-phase filter |
|
|
表 2 不同纯相位滤波器的能量值 Table 2 The energy of different pure-phase filters |
由表 2可以看出，数值计算误差和截断误差在0.03%以内，属可接受范围，可认为所求得的结果在时域和频域均具有相同的单位能量.由图 7和表 2可以验证由式（17）计算的结果为纯相位滤波器.
令真实反射系数序列激励各纯相位滤波器，所得结果如图 8所示.由图 8所得结果与图 6反射系数序列反演结果比较，二者十分相似，为进一步验证本文结果，将图 6与图 8数据根据相同道号求相似度，其结果如表 3所示.
|
图 8 反射系数序列激励不同纯相位滤波器的输出结果 Fig. 8 The output of reflection coefficient sequences incentive different pure-phase filters |
|
|
表 3 反射系数序列反演结果与纯相位滤波结果的相似度 Table 3 Similarity between inversion results of reflection coefficient sequences and pure-phase filter results |
由表 3可以看出，反射系数序列反演的结果与真实反射系数序列激励纯相位滤波器输出结果的相似度非常高，其中的误差可以认为是数据截断误差和数值运算精度误差.由此可以认为在子波相位估计不准的情况下反射系数序列反演的结果为反射系数序列与一个纯相位滤波器的褶积.
为验证纯相位滤波器的相位谱为真实子波与估计子波的相位谱之差，将图 5中真实子波与各道子波的相位谱做相减运算，其结果如图 9所示.
|
图 9 真实子波与估计子波相位谱之差 Fig. 9 The phase spectrum difference between real wavelet and estimated wavelets |
由图 9和图 7b对比可看出二者十分相似，为进一步验证二者的相似性，将图 9和图 7b中的数据按相同道求取相似度，其结果如表 4所示.由表 4可以看出纯相位滤波器的相位谱与真实子波-估计子波的相位谱之差相一致，可以认为，纯相位滤波器的相位谱为原始子波与估计子波的相位谱之差.
|
|
表 4 纯相位滤波器的相位谱与原始-估计子波相位谱之差的相似度 Table 4 Phase spectrums similarity between pure-phase filters and the difference of real-estimated wavelets |
以上实验结果表明，子波相位估计不准对反射系数序列反演具有影响，反演结果为真实反射系数序列与一个纯相位滤波器的褶积，且此纯相位滤波器的相位谱为真实子波和估计子波的相位谱之差.
5 反射系数序列反演结果的评价在实际地震信号处理过程中，反射系数序列和真实子波均为未知，真实子波的振幅谱较为容易求得，但其相位谱往往估计不准，造成反射系数序列反演结果的相位偏移.假设真实反射系数序列是稀疏脉冲序列，时域能量集中在数个脉冲内；子波相位不准的反演结果虽然具有相同的能量，但由于改变了其相位谱，时域能量则会分散.因此应用如下评价方法对反演结果进行评价，以确定出真实或准确的反射系数序列.
|
(18) |
|
(19) |
其中式（18）和式（19）分别称为变分[17]和丰度[18]方法，应用二者对本文的反射系数序列反演结果进行评价，得到如图 10所示结果.
|
图 10 反射系数序列的（a）变分评价结果和（b）丰度评价结果 Fig. 10 Evaluation results of reflection coefficient sequences with Variation (a) and Kurtosis (b) |
由图 10可以看出，变分方法在第2道取最小值，丰度方法在第2道取最大值，表明二者均可以确定出准确的反射系数序列.由于真实反射系数序列包含在反演结果中，则确定出的准确反射系数序列即为真实反射系数序列.
为验证以上方法应用的普遍性，假设真实子波为因果最小相位子波（第1道）、混合因果混合相位子波（第11道）、反因果最大相位子波（第16道），其反演结果与评价结果如图 11所示.
|
图 11 不同因果性和相位性子波的反演及评价结果 Fig. 11 Inversion and evaluation results of wavelets with different causal and phase characteristic |
图 11表明，在真实子波为任意相位和任意因果性的条件下，均可以通过变分和丰度方法对反演结果进行评价，从而确定出真实或准确的反射系数序列.同样，采用验证时域能量集中性的其它评价方法也可得到相应的结果，在此不作一一表述.
6 实际地震资料处理为验证子波相位不准对反射系数序列反演的影响规律，我们将对实际地震资料和测井所得反射系数序列进行处理.图 12b为胜利某区块的地震记录，地震记录的采样速率为2 ms，井位于第289道，图 12a为根据该井的测井数据所转换的时域反射系数序列.
|
图 12 实际地震资料 （a）第289道反射系数序列；（b）胜利某区块的地震记录. Fig. 12 Actual seismic data (a) Reflection coefficient sequences of the 289th seismic trace, and (b) seismogram of Shengli Oil Field. |
对图 12所示的实际地震资料进行地震子波提取.采用张广智等[19]提出的井旁道地震子波精细提取方法结合测井资料提取井旁道子波，采用原始的统计性子波提取方法提取最小相位的地震子波，二者具有相同的振幅谱.子波提取结果如图 13所示.
|
图 13 实际地震资料子波提取的（a）时域波形图和（b）相位谱 Fig. 13 Time domain waveforms (a) and phase spectrums of wavelet estimation in actual seismic data |
结合井旁道测井资料所提取的地震子波为混合因果混合相位的，而传统方法所提取的地震子波为因果最小相位的，如图 13所示.二者具有相同的振幅谱，仅相位谱存在差异，为验证本文所提理论的有效性，对二者进行相位谱差值化处理，以构造纯相位滤波器.构造出的纯相位滤波器如图 14所示.
|
图 14 构造出纯相位滤波器的（a）时域波形图和（b）相位谱 Fig. 14 The time domain waveform (a) and phase spectrum of constructed pure-phase filter |
应用提取出的最小相位地震子波对实际地震记录作反射系数反演处理，将反演结果与实际反射系数序列激励纯相位滤波器的输出进行比较，若二者相一致则可证明本文所提理论的有效性.测井所得实际反射系数序列与反演处理结果如图 15所示.
|
图 15 反射系数序列与反演结果 （a）测井所得反射系数序列；（b）最小相位子波反演结果；（c）反射系数序列纯相位滤波结果. Fig. 15 Reflection coefficient sequences and inversion results (a) Reflection coefficient sequences of logging; (b) Minimum phase wavelet inversion results; and (c) Pure-phase filter results of reflection coefficient sequences |
对比图 15a和图 15b可知，当子波相位估计不准时，得到的反演结果与真实的反射系数序列相差甚远，不能够正确反映地层的实际信息；由于提取的子波振幅谱是准确的，故地震记录中的低频成分经反演后压制的较为理想.对比图 15b和图 15c可以看出，反射系数序列激励图 14中的纯相位滤波器的输出与最小相位子波反演的结果较为相似.对最小相位子波反演结果和反射系数序列纯相位滤波结果计算相似度，其相似度为73.74%，二者具有较好的相似性.考虑到测井资料和地震道中的噪声、深时转换的误差等原因，给子波反演结果和反射系数序列求取结果带来一定误差，二者相似度在可接受范围内，从而表明实际数据处理结果可以验证文中结论-子波相位估计不准时反射系数序列反演结果中会残留一个纯相位滤波器，此纯相位滤波器的相位谱为真实子波和估计子波的相位谱之差.
通过丰度和变分方法可以在图 15所示的三个结果中辨识出真实的反射系数序列，验证了评价方法的有效性.丰度和变分方法评价结果如表 5所示，*表示为最佳结果.
|
|
表 5 实际资料处理的变分和丰度评价结果 Table 5 Variation and Kurtosis evaluation results of actual seismic data |
本文用ARMA模型描述地震子波，提出采用对称映射z域中ARMA模型的零极点的方式构造出相同振幅谱、不同相位谱的子波，并对合成的地震记录进行反射系数序列反演.反演结果表明，当估计子波和真实子波具有相同振幅谱、不同相位谱时，反射系数序列反演的结果为真实反射系数序列与一个纯相位滤波器的褶积，此滤波器的相位谱为真实子波和估计子波相位谱之差.由于纯相位滤波器的作用，反褶积结果破坏了真实反射系数序列的稀疏性（时域能量集中性）.通过两种方法对反演结果进行了评价，评价结果表明二者均可确定出真实或准确的反射系数序列.本文针对子波相位估计不准对反射系数序列反演的影响进行了理论分析、实验仿真和实际地震资料处理，仿真和实际处理结果可以验证理论结果，为后续进一步提高反射系数序列反演结果精度指明了研究方向.
致谢感谢中国石油大学（华东）地球科学与技术学院印兴耀教授和张广智教授的大力支持，感谢中国石油大学（北京）袁三一老师的指导与建议.
| [1] | 唐斌, 尹成. 基于高阶统计的非最小相位地震子波恢复. 地球物理学报 , 2001, 44(3): 404–410. Tang B, Yin C. Non-minimum phase seismic wavelet reconstruction based on higher order statistics. Chinese J. Geophys. (in Chinese) , 2001, 44(3): 404-410. |
| [2] | 伊振林, 王润秋. 一种新的混合相位反褶积方法. 石油地球物理勘探 , 2006, 41(3): 266–270. Yi Z L, Wang R Q. A new mixed phase deconvolution method. Oil Geophysical Prospecting (in Chinese) , 2006, 41(3): 266-270. |
| [3] | 崔庆辉, 芮拥军, 尚新民, 等. 混合相位地震子波提取及应用. 石油物探 , 2011, 50(5): 481–486. Cui Q H, Rui Y J, Shang X M, et al. Mixed-phase seismic wavelet extraction and its application. Geophysical Prospecting for Petroleum (in Chinese) , 2011, 50(5): 481-486. |
| [4] | Levy S, Oldenburg D W. The deconvolution of phase-shifted wavelets. Geophysics , 1982, 47(9): 1285-1294. DOI:10.1190/1.1441388 |
| [5] | Levy S, Oldenburg D W. Automatic phase correction of common-midpoint stacked data. Geophysics , 1987, 52(1): 51-59. DOI:10.1190/1.1442240 |
| [6] | Wang Y H. A stable and efficient approach of inverse Q filtering. Geophysics , 2002, 67(2): 657-663. DOI:10.1190/1.1468627 |
| [7] | Wang Y H. Inverse Q-filter for seismic resolution enhancement. Geophysics , 2006, 71(3): V51-V60. DOI:10.1190/1.2192912 |
| [8] | Yuan S Y, Wang S X. Influence of inaccurate wavelet phase estimation on seismic inversion. Applied Geophysics , 2011, 8(1): 48-59. DOI:10.1007/s11770-011-0273-5 |
| [9] | 梁光河. 地震子波提取方法研究. 石油物探 , 1998, 37(1): 3l–39. Liang G H. On the methods of seismic wavelet extraction. Geophysical Prospecting for Petroleum (in Chinese) , 1998, 37(1): 3l-39. |
| [10] | 戴永寿, 王俊岭, 王伟伟, 等. 基于高阶累积量ARMA模型线性非线性结合的地震子波提取方法研究. 地球物理学报 , 2008, 51(6): 1851–1859. Dai Y S, Wang J L, Wang W W, et al. Seismic wavelet extraction via cumulant-based ARMA model approach with linear and nonlinear combination. Chinese J. Geophys. (in Chinese) , 2008, 51(6): 1851-1859. |
| [11] | Robinson E A. Predictive decomposition of time series with application to seismic exploration. Geophysics , 1967, 32(3): 418-484. DOI:10.1190/1.1439873 |
| [12] | 吴镇扬. 数字信号处理. 北京: 高等教育出版社, 2004 . Wu Z Y. Digital Signal Processing (in Chinese). Beijing: Higher Education Press, 2004 . |
| [13] | Wiggins R A. Minimum entropy deconvolution. Geoexploration , 1978, 16(1-2): 21-35. DOI:10.1016/0016-7142(78)90005-4 |
| [14] | Sacchi M D. Re-weighting strategies in seismic deconvolution. Geophysical Journal International , 1997, 129: 651-656. DOI:10.1111/gji.1997.129.issue-3 |
| [15] | 章珂, 李衍达, 刘贵忠, 等. 多分辨率地震信号反褶积. 地球物理学报 , 1999, 42(4): 529–535. Zhang K, Li Y D, Liu G Z, et al. Multiresolution seismic signal deconvolution. Chinese J. Geophys. (in Chinese) , 1999, 42(4): 529-535. |
| [16] | Zhang R, Castagna J. Seismic sparse-layer reflectivity inversion using basis pursuit decomposition. Geophysics , 2011, 76(6): R147-R158. DOI:10.1190/geo2011-0103.1 |
| [17] | Yuan S Y, Wang S X. Noise attenuation without spatial assumptions about seismic coherent events.//80th Ann. Internat. Mtg., Soc. Explor. Geophys., Expanded Abstracts, 2010, 3524-3528. http://www.oalib.com/references/18987224 |
| [18] | Wiggins R. Entropy guided deconvolution. Geophysics , 1985, 50(12): 2720-2726. DOI:10.1190/1.1441892 |
| [19] | 张广智, 刘洪, 印兴耀. 井旁道地震子波精细提取方法. 石油地球物理勘探 , 2005, 40(2): 158–162. Zhang G Z, Liu H, Yin X Y. Method for fine picking up seismic wavelet at up-hole trace. Oil Geophysical Prospecting (in Chinese) , 2005, 40(2): 158-162. |