地球物理学进展  2017, Vol. 32 Issue (1): 210-216   PDF    
基于偏斜度准则的地震子波相位估计方法特征研究
刘俊州1, 马铭2, 时磊1, 夏红敏1     
1. 中国石油化工股份有限公司石油勘探开发研究院, 北京 100083
2. 中国石油大学 (北京) 油气资源与探测国家重点实验室, 北京 102249
摘要:本文研究了基于偏斜度准则的常相位旋转地震子波估计方法的优势,并分析了偏斜度准则与峰度准则的差异性及其适用条件.主要通过数值模拟试验充分说明偏斜度准则相对于峰度准则具有更高的动态范围及稳定性.偏斜度准则克服了分析信号带宽必须大于中心频率的条件,降低了反射系数对于子波相位估计的影响.相对于峰度准则,偏斜度准则在处理低信噪比的地震资料过程中表现出了更好的稳定性和准确性.通过蒙特卡洛试验进一步说明基于偏斜度准则的常相位子波估计能够准确地估计出相位为常数的地震子波,并且子波相位估计结果对地震资料噪声水平的强弱不敏感.
关键词地震子波相位估计    常相位    偏斜度    峰度    
Characteristics study of seismic wavelet phase estimation based on the skewness criterion
LIU Jun-zhou1 , MA Ming2 , SHI Lei1 , XIA Hong-min1     
1. SINOPEC Exploration & Production Research Institute, Beijing 100083, China
2. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China
Abstract: We make comparison analysis of skewness and kurtosis criteria for constant-phase wavelet estimation, which plays a significant part in high-resolution seismic data processing and interpretation. Synthetic examples illustrate that skewness criterion has higher dynamic range and better stability than kurtosis criterion. As a high-order statistics, skewness criterion breaks the limitation that the bandwidth of input signal must be greater than center frequency and reduces the influence of noise on wavelet phase estimation. Like kurtosis criterion, the estimated phase corresponding to maximum skewness is also affected by the distribution of reflectivity. Through the Monte Carlo simulation it can be illustrated that the result of wavelet phase estimation based on the skewness criterion is not sensitive to the noise level.
Key words: seismic wavelet phase estimation     constant phase     skewness     kurtosis    
0 引言

叠后地震资料可以认为是地震子波和地下反射系数的褶积运算结果 (Robinson and Treitel, 1980; Robinson, 1984; 奈德尔, 1986).地震勘探的一个主要目的是得到地下有效的反射信息.地震资料处理中消除地震子波影响,使得叠后地震剖面直接反映地下地层反射系数,在地震剖面上显示出更多的细节,从而提高纵向分辨率,这个过程即为反褶积.如何精确地估计地震子波并且应用于高分辨率处理和精细解释领域是当前地震勘探中的重要研究内容.地震子波振幅谱影响的消除可以拓宽地震资料的频带,压缩地震子波;而相位谱的消除在改善同相轴的横向连续性方面起到重要的作用 (Yuan and Wang, 2011; 于永才, 2012; Yu et al., 2012; Yuan and Wang, 2013; 张亚南等, 2013).地震子波的估计方法可以分为确定性方法和统计性方法 (杨培杰和印兴耀, 2008; 高少武等, 2009).在井资料可用的区域,可以使用井资料获得反射系数序列,然后结合井旁道地震记录采用确定性方法进行子波估计.而缺少井资料或井资料无法使用的探区,统计性子波估计方法则只需要从地震资料本身出发来计算子波 (Ulrych, 1971; Lazear, 1993; 王嘉松和曹桂荣, 1998; 戴永寿等, 2008).传统的维纳反褶积被广泛地用于消除最小相位子波,在一定程度上提高了地震道的分辨率 (Leinbach, 1995).但是对于反射系数是白噪以及地震子波为最小相位的假设条件在实际资料中并不完全满足.针对地震子波振幅谱估计方面,由于大地滤波作用,导致地震记录的振幅谱不光滑,因此可以通过平滑地震记录振幅谱的方法得到地震子波的振幅谱.脉冲反褶积可以准确提供一个子波的振幅谱,并且消除了与真实混合相位子波对等的一个最小相位子波成分,而剩余相位会残存在反褶积结果中.这些剩余相位特别是通频带内的剩余相位经常会导致脉冲反褶积结果的波形发生严重畸变,进一步可能在解释过程中错误地判定地下地质异常体的几何形状和位置 (刘玉金等, 2014).因此,从地震剖面上消除子波的正确相位在提高纵向分辨率及横向连续性方面有着重要的作用.

针对估计地震子波相位的方法,常相位旋转和移根算法都能取得良好的效果.常相位旋转方法的假设条件是地震子波的相位为固定的常数.通过给定范围内子波相位的扫描得到使得准则函数达到极值的相位, 该相位即为地震子波的相位.根变换估计方法是基于子波Z变换的根进行对称变换时子波相位发生改变的现象,通过在相同振幅谱不同相位谱的子波库中选择使得准则函数达到极值的子波.无论是常相位旋转还是移根算法,估计出的子波相位的准确性与地震资料的性质和质量密切相关,并且受到准则函数的影响.一般在实际资料中采用井震匹配作为准则函数来确定最佳的子波相位.但是测井资料并不总是可用的.在地震子波的相位谱不依赖于频率变化的条件下,Levy和Oldenburg (1987)Longbottom等 (1988)White (1988)、以及van der Baan (2008)探索了基于峰度准则的地震子波相位谱估计,通过常相位旋转方式和峰度最大值准则确定子波的相位,峰度准则属于稀疏准则.随后,很多学者研究峰度准则的适用条件,并且进行了许多改进,取得了良好的效果.Sacchi等 (1994)Lu (2005)Lu和Liu (2007)分别对峰度准则进行了局部改进,减小了传统峰度对相对较大的反射系数的敏感性.Parsimony准则是Claerbout (1977)提出一种反褶积约束准则,该准则对强弱反射系数都具有较高的灵敏度,并且能很好的保存弱反射信息.Ooe和Ulrych (1979)针对峰度准则对噪声的不敏感性,提出了指数变换准则, 其改善了噪声压制与恢复小反射系数序列之间的平衡.Sacchi (2002)研究了鲁棒统计学中的Sech准则和Cauchy准则,并分别将其作为正则化因子引入到稀疏脉冲反褶积中.两种准则在抗噪性能和寻找强反射体上具有较大的优势,但对弱反射信息有压制作用.Fomel和van der Baan (2014)使用局部偏斜度对地震资料局部相位进行估计,并且将该准则拓展到空间域,定义了偏斜度的局部属性,对非稳态资料的处理得到了很好的效果.偏斜度相对于峰度具有较大的动态范围,并且在估计地震子波相位的应用中具有独特地性质.

本文分析基于偏斜度和峰度准则的常相位旋转方法估计地震子波相位方面的特点和差异,并进行了一系列的地震数据测试,通过蒙特卡洛试验验证基于偏斜度准则的常相位地震子波估计方法的稳定性和准确性.得出偏斜度准则函数用于指示地震子波相位的特征总结.

1 峰度和偏斜度准则函数

通过地震记录直接估计地震子波的相位是统计性地震子波估计的重要步骤.通过常相位旋转方法将地震记录的相位校正为零,则旋转相位的大小即为地震子波的相位.而一般具有零相位的信号会使得准则函数值达到极大或极小.常见的准则函数为峰度准则.它测量了信号偏离高斯性的程度 (Wiggins, 1978),其表达式为

(1)

其中N表示采样点数.对于高斯信号的峰度值为3,峰度值偏离3越大,表示信号偏离高斯性的程度越大.当信号具有零相位时,其峰度值会达到最大.但只有在信号的带宽大于中心频率时, 上述条件才能够满足.Fomel和van der Baan (2014)提出了不同于峰度准则的测量方法,即偏斜度准则.对于信号的偏斜度定义为

(2)

在统计学中,偏斜度是用来测量概率分布的不对称性 (Bulmer, 1979).通过简单的公式推导可以得到偏斜度函数值可以使用信号自相关的平方表示.即:

(3)

具有零相位的信号由于和其自身的平方具有较高的相关性,因此其偏斜度也高.不同于测量信号高斯程度的峰度函数,偏斜度与对称性密切联系.两种测量函数均能在信号为零相位时达到最大值,但是通过实际的测试可以看出:偏斜度准则函数相对于峰度准则具有较大的动态范围.并且当信号的带宽小于中心频率时,其仍然能保持在零相位处取得极大值.

2 带限信号测试

本文通过带限相移滤波器测试偏斜度准则函数值随频率变化的影响,并与峰度准则函数进行对比,验证在不同相位下函数最大值及最小值随信号中心频率和带宽的变化.考虑带限相移滤波器(White, 1988),公式为

(4)

其中,B=(fH-fl) 表示带宽,f0=(fH+fL)/2表示中心频率,fHfL分别表示高截频率和低截频率.θ表示相移因子,即滤波器的相位谱.该信号的谱可以看作地震道经过叠后反褶积、带通滤波、有效频带内的谱白化之后得到的谱.以该滤波器为基础,研究偏斜度准则和峰度准则的适用性条件.

图 1给出了低截频率和高截频率分别为20 Hz和120 Hz、相移因子为0°的带限相移滤波器的脉冲响应及其振幅谱,振幅谱出现跳动的原因是Gibbs效应引起的.通过改变频宽与中心频率的关系,可以观察到偏斜度和峰度值的变化情况,图 2给出了三个频带的带限相移滤波器的相移因子从-π到π对应的偏斜度值以及峰度值.带限相移滤波器的频率组合分别为15~25 Hz、15~65 Hz、15~125 Hz.三个频带对应的频宽和中心频率分别为 (10, 20) Hz、(50, 40) Hz、(110, 70) Hz.其中,第一个滤波器的带宽小于中心频率,第二个和第三个滤波器则相反.从图 2a图 2b可以看出,对于15~65 Hz和15~125 Hz的相移滤波器,偏斜度准则函数值和峰度准则函数值均在零相位处取得极大值.但15~25 Hz的滤波器对应的峰度准则并未有强烈的数值变化,并且零相位处达到极小值,因此无法从极大值中识别出零相位.而偏斜度准则函数值依然在零相位处取得极大值,整体值的变化趋势并未发生改变.通过对比偏斜度和峰度我们可以看出:偏斜度准则相对于峰度准则在中心频率大于频宽时仍然能保持在零相位时取得极大值.但准则函数值的变化范围已经非常小,为10-5.

图 1 fL=20 Hz, fH=120 Hz, θ=0°的带限相移滤波器的脉冲响应 (a) 及振幅谱 (b) Figure 1 Impulse response (a) and amplitude spectrum (b) of a band-limited phase-shifter with fL=20 Hz, fH=120 Hz, and θ=0°

图 2 不同带限相移滤波器的相移因子从-π到π对应的偏斜度值 (a) 以及峰度值 (b) Figure 2 The change of skewness (a) and kurtosis (b) in three band-limited phase-shifters with different fL and fH

固定带限相移滤波器的低截频率为15 Hz,高截频率以1 Hz的间隔从16 Hz增加到200 Hz.该过程类似于反褶积方法中拓宽频带的过程,通过该实验研究偏斜度和峰度与高频拓宽之间的关系.图 3分别为高截频率对应的最大最小偏斜度值和峰度值变化,其中的圆圈代表中心频率B与频宽f0相等的点.从图中可以看出,在该点的右端B < f0,随着高截频率的逐渐增大 (带宽越来越大于中心频率),偏斜度和峰度的变化范围逐渐增大.而该点处,B=f0时,偏斜度的极大值与极小值仍可以在图中区分,而峰度的极大值与极小值变化已经难以分辨.因此可以得出认识:偏斜度准则函数值相对于峰度准则函数值在一定程度上适应了信号中心频率大于带宽条件下的指示结果,并且其在不同频率下的变化范围明显大于峰度值的变化范围.

图 3 固定带限相移滤波器的低截频率为15 Hz,增大高截频率时偏斜度值范围 (a) 和峰度值范围 (b) 的变化情况 Figure 3 The change of skewness (a) and kurtosis (b) with increased high cut-off frequency, with low cut-off frequency fixed at 15 Hz
3 理论合成数据测试:反射系数及噪声影响

对不同类型的反射系数与零相位雷克子波褶积合成的地震记录经过常相位旋转形成相位旋转地震道.通过偏斜度准则函数和峰度准则函数估计地震子波的相位.图 4分别给出了稀疏反射系数模型和主频为30Hz的零相位雷克子波合成地震记录进行常相位旋转得到的旋转地震道集以及其对应的偏斜度值和峰度值,其中第26道是原始地震道 (旋转相位为0°).通过偏斜度函数值和峰度函数值可以看出,当反射系数为稀疏时,通过偏斜度极大值确定的地震子波相位约为5°,与真实相位之间存在微弱偏差.而峰度极大值确定地震子波为零相位子波.同样的步骤用于处理基于非稀疏反射系数模型的合成资料,即反射系数序列不再稀疏时,峰度准则通过最大值确定的相位不再为零相位,出现了较大的偏差,偏斜度极大值也没有出现在零相位处,如图 5所示.但相对于峰度准则,其对应的相位要更接近于零相位.当反射系数非稀疏时,峰度准则函数极大值对应的相位与真实子波相位与真实的零相位之间大约存在28°的偏差,而偏斜度极大值对应的相位与真实相位保持在5°左右.通过对比可以发现,使用常相位旋转方法估计地震子波相位时,合成地震记录与原始地震子波的相位并不相同.即褶积运算会改变地震子波所对应地震记录的相位.而准则函数需克服反射系数的影响,从地震记录中得到正确的相位.偏斜度准则函数相对于峰度准则函数在消除反射系数对子波相位估计的影响方面更加优越,对于非稀疏反射系数,其效果更加显著.

图 4 稀疏反射系数模型测试 (a) 稀疏反射系数模型;(b) 合成地震记录的相位旋转地震道;(c) 不同旋转相位下的偏斜度值;(d) 不同旋转相位下的峰度值. Figure 4 Test of the wavelet estimation based on skewness and kurtosis criteria applied to a synthetic record with sparse reflectivity (a) Sparse reflectivity; (b) Synthetic record rotated through various phases; (c) Skewness as function of the phase rotation angle; (d) Kurtosis as function of the phase rotation angle.

图 5 非稀疏反射系数模型测试 (a) 稀疏反射系数模型;(b) 合成地震记录的相位旋转地震道;(c) 不同旋转相位下的偏斜度值;(d) 不同旋转相位下的峰度值. Figure 5 Test of the wavelet estimation based on skewness and kurtosis criteria applied to a synthetic record with non-sparse reflectivity (a) Non-sparse reflectivity; (b) Synthetic record rotated through various phases; (c) Skewness as function of the phase rotation angle; (d) Kurtosis as function of the phase rotation angle.

为了考虑噪声对准则函数识别子波相位的影响,通过零相位的雷克子波与非稀疏反射系数合成地震记录.向地震记录中加入不同程度的高斯随机噪声.得到合成地震记录的信噪比 (SNR) 分别为30、20、15、10.计算其在不同相位旋转下的偏斜度值和峰度值,通过偏度值或峰度值最大来估计地震子波的相位.为了更加准确地反映两种相位指示准则函数的特征,本文进行了1000次蒙特卡洛模拟实验.最终统计结果如图 6图 7所示.通过测试可以看到:随着含噪程度的增加,峰度准则的极大值对应的相位更加偏离正确相位.而对于偏斜度准则函数值,其整体趋势并未受到噪声的影响.极大值对应的相位更加接近于零相位.因此可以得出认识:偏斜度准则相对于峰度准则在噪声方面表现出较强的低敏感性.不同的信噪比对偏斜度函数值的影响有限.表现出较强的鲁棒性.

图 6 不同信噪比条件下基于偏斜度准则子波相位估计结果的统计特征 (a) SNR=30; (b) SNR=20; (c) SNR=15; (d) SNR=10. Figure 6 The statistical distributions of estimated phases based on skewness with SNR=30 (a), 20 (b), 15 (c) and 10 (d) applied on the seismic records

图 7 不同信噪比条件下基于峰度准则子波相位估计结果的统计特征 (a) SNR=30; (b) SNR=20; (c) SNR=15; (d) SNR=10. Figure 7 The statistical distributions of estimated phases based on kurtosis with SNR=30 (a), 20 (b), 15 (c) and 10 (d) applied on the seismic records
4 结论

偏斜度和峰度计算量相对较小,处理实际资料中通常采用井震匹配估计零相位子波,而通过合适的准则函数评判可以得到最佳的常相位子波或者混合相位子波.通过带限相移滤波器和合成地震记录对偏斜度和峰度函数的测试,可以得出如下结论:

(1) 偏斜度准则函数相对于峰度准则函数具有较大的动态范围和稳定性.

(2) 当信号的中心频率大于频宽时 (一定范围内), 偏斜度函数值仍然保持在零相位处取得极大值.

(3) 反射系数对地震子波相位估计存在一定的影响,会导致偏斜度准则函数和峰度准则函数极大值发生偏转,在非稀疏反射系数假设条件下,偏斜度准则函数表现要优于峰度准则函数.

(4) 常规高斯噪声存在情况下,偏斜度准则函数在抗噪效果方面表现要优于峰度准则函数.但当噪声为一定频率范围内有色噪声时,偏斜度准则函数和峰度准则函数值在相位表现方面都出现较大的偏差.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持!
参考文献
[] Bulmer M. 1979. Principles of Statistics[M]. New York: Dover Publications.
[] Claerbout J F. 1977. Parsimonious deconvolution[J]. Stanford Exploration Project, 13: 1–9.
[] Dai Y S, Wang J L, Wang W W, et al. 2008. Seismic wavelet extraction via cumulant-based ARMA model approach with linear and nonlinear combination[J]. Chinese Journal of Geophysics (in Chinese), 51(6): 1851–1859. DOI:10.3321/j.issn:0001-5733.2008.06.027
[] Fomel S, van der Baan M. 2014. Local skewness attribute as a seismic phase detector[J]. Interpretation, 2(1): SA49–SA56. DOI:10.1190/INT-2013-0080.1
[] Gao S W, Zhao B, He Z H, et al. 2009. Research progress of seismic wavelet extraction[J]. Progress in Geophysics (in Chinese), 24(4): 1384–1391. DOI:10.3969/j.issn.1004-2903.2009.04.028
[] Lazear G D. 1993. Mixed-phase wavelet estimation using fourth-order cumulants[J]. Geophysics, 58(7): 1042–1051. DOI:10.1190/1.1443480
[] Leinbach J. 1995. Wiener spiking deconvolution and minimum-phase wavelets:A tutorial[J]. The Leading Edge, 14(3): 189–192. DOI:10.1190/1.1437110
[] Levy S, Oldenburg D W. 1987. Automatic phase correction of common-midpoint stacked data[J]. Geophysics, 52(1): 51–59. DOI:10.1190/1.1442240
[] Liu Y J, Li Z C, Wu D, et al. 2014. Well constrained non-stationary phase correction[J]. Chinese Journal of Geophysics (in Chinese), 57(1): 310–319. DOI:10.6038/cjg20140127
[] Longbottom J, Walden A T, White R E. 1988. Principles and application of maximum kurtosis phase estimation[J]. Geophysical Prospecting, 36(2): 115–138. DOI:10.1111/gpr.1988.36.issue-2
[] Lu W K. 2005. Non-minimum-phase wavelet estimation using second-and third-order moments[J]. Geophysical Prospecting, 53(1): 149–158. DOI:10.1111/gpr.2005.53.issue-1
[] Lu W K, Liu D Q. 2007. Frequency recovery of band-limited seismic data based on sparse spike train deconvolution[C].//77th Annual International Meeting, SEG. Expanded Abstracts, 1977-1981.
[] Neidel. 1986. The Convolutional Model (in Chinese)[M]. Niu Y Q Trans. Beijing:Petroleum Industry Press.
[] Ooe M, Ulrych T J. 1979. Minimum entropy deconvolution with an exponential transformation[J]. Geophysical Prospecting, 27(2): 458–473. DOI:10.1111/gpr.1979.27.issue-2
[] Robinson E A. 1984. Seismic Inversion & Deconvolution:Classical Methods[M]. London-Amsterdam: Geophysical Press.
[] Robinson E A, Treitel S. 1980. Geophysical Signal Analysis[M]. Englewood Cliffs, NJ: Prentice-Hall Inc.
[] Sacchi M D. 2002. Statistical and transform methods in geophysical signal processing[D]. Edmonton, Canada:University of Alberta.
[] Sacchi M D, Velis D R, Cominguez A H. 1994. Minimum entropy deconvolution with frequency-domain constraints[J]. Geophysics, 59(6): 938–945. DOI:10.1190/1.1443653
[] Ulrych T J. 1971. Application of homomorphic deconvolution to seismology[J]. Geophysics, 36(4): 650–660. DOI:10.1190/1.1440202
[] van der Baan M. 2008. Time-varying wavelet estimation and deconvolution by kurtosis maximization[J]. Geophysics, 73(2): V11–V18. DOI:10.1190/1.2831936
[] Wang J S, Cao G R. 1998. The method of fractal impulse deconvolution[J]. Acta Geophysica Sinica (in Chinese), 41(1): 99–108.
[] White R E. 1988. Maximum kurtosis phase correction[J]. Geophysical Journal International, 95(2): 371–389. DOI:10.1111/j.1365-246X.1988.tb00475.x
[] Wiggins R A. 1978. Minimum entropy deconvolution[J]. Geoexploration, 16(1-2): 21–35. DOI:10.1016/0016-7142(78)90005-4
[] Yang P J, Yin X Y. 2008. Summary of seismic wavelet pick-up[J]. Oil Geophysical Prospecting (in Chinese), 43(1): 123–128.
[] Yu Y C. 2012. Method Study of seismic wavelet estimation and seismic wave attenuation parameter estimation (in Chinese)[Ph. D. thesis]. Beijing:China University of Petroleum.
[] Yu Y C, Wang S X, Yuan S Y, et al. 2012. Phase spectrum estimation of the seismic wavelet based on a criterion function[J]. Petroleum Science, 9(2): 170–181. DOI:10.1007/s12182-012-0197-6
[] Yuan S Y, Wang S X. 2011. Influence of inaccurate wavelet phase estimation on seismic inversion[J]. Applied Geophysics, 8(1): 48–59. DOI:10.1007/s11770-011-0273-5
[] Yuan S Y, Wang S X. 2013. Spectral sparse Bayesian learning reflectivity inversion[J]. Geophysical Prospecting, 61(4): 735–746. DOI:10.1111/gpr.2013.61.issue-4
[] Zhang Y N, Dai Y S, Chen J, et al. 2013. The research on the influence of wavelet phase on the inversion results of reflection coefficient sequences by using the ARMA model of symmetrical mapping Pole-Zeros[J]. Chinese Journal of Geophysics (in Chinese), 56(6): 2043–2054. DOI:10.6038/cjg20130625
[] 戴永寿, 王俊岭, 王伟伟, 等. 2008. 基于高阶累积量ARMA模型线性非线性结合的地震子波提取方法研究[J]. 地球物理学报, 51(6): 1851–1859. DOI:10.3321/j.issn:0001-5733.2008.06.027
[] 高少武, 赵波, 贺振华, 等. 2009. 地震子波提取方法研究进展[J]. 地球物理学进展, 24(4): 1384–1391. DOI:10.3969/j.issn.1004-2903.2009.04.028
[] 刘玉金, 李振春, 吴丹, 等. 2014. 井约束非稳态相位校正方法[J]. 地球物理学报, 57(1): 310–319. DOI:10.6038/cjg20140127
[] 奈德尔. 1986.褶积模型[M].牛毓荃译.北京:石油工业出版社.
[] 王嘉松, 曹桂荣. 1998. 分形脉冲反褶积方法[J]. 地球物理学报, 41(1): 99–108.
[] 杨培杰, 印兴耀. 2008. 地震子波提取方法综述[J]. 石油地球物理勘探, 43(1): 123–128.
[] 于永才. 2012.地震子波估计与地震波吸收衰减参数估计方法研究[博士论文].北京:中国石油大学.
[] 张亚南, 戴永寿, 陈健, 等. 2013. 用对称映射ARMA模型的零极点研究子波相位对反射系数序列反演的影响[J]. 地球物理学报, 56(6): 2043–2054. DOI:10.6038/cjg20130625