地球物理学进展  2014, Vol. 29 Issue (4): 1780-1784   PDF    
引用本文
王万里, 杨午阳, 魏新建, 桂金咏. 随机稀疏脉冲非线性反褶积[J]. 地球物理学进展, 2014, 29(4): 1780-1784, doi: 10.6038/pg20140439  
WANG Wan-li, YANG Wu-yang, WEI Xin-jian, GUI Jin-yong. Stochastic sparse spike nonlinear deconvolution. Progress in Geophysics, 2014, 29(4): 1780-1784, doi: 10.6038/pg20140439   
随机稀疏脉冲非线性反褶积
王万里, 杨午阳, 魏新建, 桂金咏    
中国石油勘探开发研究院西北分院, 兰州 730020
收稿日期: 2013-11-11 修回日期: 2014-5-5 .
基金项目：大型油气田及煤层气开发国家重大科技专项“基于GeoEast的多波地震、测井和地质信息的裂缝综合储层综合分析软件集成”子课题（2011ZX05019-008-005）资助.
作者简介：王万里，男，1986年生，硕士，工程师，中石油勘探勘探开发研究院西北分院，主要从事复杂构造地震资料处理和储层预测研究工作.（E-mail：wang_wl@petrochina.com.cn）
摘要：常规的反褶积方法通过线性褶积压缩子波提高地震记录的分辨率，其能力受到有效信号频带的限制.随机稀疏脉冲非线性反褶积方法将传统的以子波压缩为核心理念的反褶积方法转移到反射系数位置和大小的检测上来，它直接从地震记录中通过非线性反演方法得到反射系数的位置和大小，突破了地震资料有效频带的限制，能够较大幅度提高地震记录的分辨率.同时通过对反射系数统计特征的有效约束，减小了反褶积结果的多解性.模型实验表明，随机稀疏脉冲反褶积对噪声和子波的敏感性较小，能够较好的保护弱反射信号.在模型实验的基础上，利用随机稀疏脉冲反褶积对实际地震资料进行了实验处理，有效的改善了地震资料的分辨率.
关键词： 随机稀疏     正则化     反褶积     分辨率    
Stochastic sparse spike nonlinear deconvolution
WANG Wan-li, YANG Wu-yang, WEI Xin-jian, GUI Jin-yong    
Research Institute of Petroleum Exploration & Development-northwest, PetroChina, Lanzhou 730020, China
Abstract: Deconvolution is an effective and essential method to improve the resolution of seismic data. Conventional deconvolution methods assume that the wavelet is minimum phase and the reflection coefficients are white noise distribution, which improve the resolution of seismic records by wavelet compression and linear convolution, but their abilities are limited by the frequency band of effective signal. So the detailed study of stochastic sparse spike nonlinear deconvolution is provided. Instead of compressing the wavelet, stochastic sparse spike deconvolution improves the resolution of seismic records by detecting the location and size of reflection coefficients, and which directly gets the location and size of reflection coefficients by nonlinear inversion from the seismic records, thus breaking through the frequency band limit of seismic data and improving the resolution of seismic records greatly. While combining the effective constraints of statistical characteristics of reflection coefficients, the multi solutions of deconvolution are reduced. The experimental results demonstrate that stochastic sparse spike nonlinear deconvolution is less sensitive to noise and wavelet and protects the weak reflection better than other deconvolution methods. Based on the model experiments, stochastic sparse spike nonlinear deconvolution is applied to real seismic data and improves the resolution of seismic data effectively.
Key words: stochastic sparse spike     regularization     deconvolution     resolution    
0 引 言

随着油气勘探开发工作的不断深入，地震资料分辨率与储层预测精度的矛盾日益突出.迫切需要提高地震资料的分辨率和成像精度，依托高分辨率地震资料对储层的结构特征和沉积层序进行精细解释和深入描述(李振春等，2012；刘志伟等，2013).

反褶积技术是目前提高地震资料分辨率的主要方法，目前工业界应用较多的反褶积方法有脉冲反褶积、预测反褶积、谱白化、谱模拟和同态反褶积等(戴永寿等，2011；蔡连芳等， 2012).尽管反褶积的方法较多，但大多都需要某些数学假设才能成立，高斯白噪反射系数和最小相位子波是两个最为常用的反褶积假设.随着储层预测精度对地震资料分辨率要求的不断提高，反褶积的数学假设与实际地层统计特征的矛盾已经成为制约进一步提高地震资料分辨率的主要技术障碍(张繁昌等，2008；安西峰等，2010).

董敏煜(1991)、董恩清(2001)等通过求解卡曼滤波方程进行自适应反褶积，克服了提高分辨率而信噪比降低的问题.王嘉松(1998)通过研究证明了反射系数并不都是白噪，而是某种分形噪声，并提出了分形反褶积，较大程度上提高了地震资料的分辨率.李海山(2012)通过对反射系数的概率密度统计特征分析，得出反射系数更趋近于超高斯分布特征，并实现了地震子波和反射系数同时估计的地震盲反褶积.高静怀(2012)指出地震资料中子波非最小相位的，同时给出了地震记录变子波模型的一种近似表达式.通常，已有的反褶积方法没有考虑或者弱化了子波带限、子波的相位特性、反射系数的特征、分辨率和信噪比的关系等问题.鉴于此，本文试图寻求一些相对于常规的反褶积方法而言，能够识别微小结构，特别是薄层，并具有更合理的理论和更可靠的结果的反褶积方法.

随机稀疏脉冲反褶积摒弃了这两种假设，它直接从含噪信号(地震记录)中估算脉冲(反射系数)的振幅及其位置.其主要目标是寻找最少数量的脉冲，它们与估计的有限带宽地震子波褶积时，褶积结果与地震资料间的误差满足给定的要求.脉冲时间(位置)的确定是个高度非线性的最优化问题，可用快速模拟退火(SA)来求解.每次SA迭代时可用线性最小二乘方法方便地计算出脉冲振幅.反褶积方法中可引入阻抗约束，而二次型正则化的使用保证了解的稳定性，也减少非唯一性，同时使得求解更具有物理意义(Walden， 1985；Sacchi，1994；Okaya，1995). 1 基本原理

地震记录的褶积模型可以表示为

式中s(t)表示地震记录，w(t)是地震子波，r(t)表示反射系数，n(t)表示噪声. 稀疏反射系数序列可以表示为如下的形式：

式中M是非零脉冲的个数，M< j是脉冲的时延，αj是脉冲的振幅，δj是单位脉冲函数. 将(2)式代入(1)式中得到：

式中，st为地震道，wt为地震子波，nt为随机噪音. 将(3)式写成矩阵形式：

式中A 是一个N*M矩阵，其中Aij=ωi－τj+1.它是一个N维非线性方程组，包含2M个未知数，可以通过实际地震记录与合成地震记录的最小平方差来求解.因此，可定义如下目标函数

使用快速SA(模拟退火)求取时间延迟的同时将Ja最小化.振幅通过在每次SA迭代中用线性最小二乘法(已知时间延迟)估计出：

得到

式中Aij=ωi－τj+1. 为使反演结果更稳定，反演时需要增加正则化准则，尤其是当M较大的时候.考虑对待求反射系数的幅值α进行约束，从而保证解的稀疏性，实现最少的非零脉冲(反射系数)同子波的合成数据与实际地震数据的匹配差最小.稀疏解的特征，可以选择采用一范数约束，二次型正则化约束以及柯西范数约束，本文采用了二次型约束. 如果存在阻抗资料时，还可以加入阻抗约束.假设已知阻抗为ξ，新的目标函数可以定义为：

其中Jα=是解的L2范数形式，μ为阻抗权系数，C 为下三角全为1的方阵.在实际中，令β=β0f0，f0=maxi Fii，β0是一些正的小数.然而，对于一个给定的M，只要反演出来的结果尊重地震数据和约束条件，希望β0尽可能的大. 当已知时间位置后，新的解可以写成以下形式

式中β是衰减因子(即预白化)，I是单位矩阵，F=A T A+μC T C .

随机稀疏脉冲反褶积的求解是一个迭代过程，其中的每次SA迭代都涉及到估算一系列M延迟(位置)，矩阵和C、F，同时通过计算误差Ja来检验其收敛性.当达到期望误差或者用户定义的最高迭代数itmax时，SA迭代停止.主要过程如下：

(1)固定M，设置阻抗约束项权系数μ和二次型正则化项的阻尼系数β为

(2)随机给定反射系数所在的初始时间位置，j=1…M；

(3)对于每次SA(模拟退火)迭代中的第l次扰动时，按下述过程操作：

(a)更新反射系数的时间位置 ，j=1…M；

(b)计算新的矩阵 A、C，和振幅α ；

(c)计算出J；

(d)根据Metropolis准则(SA算法中的概率接受准则)，判断是否接受反射系数的时间位置.

(e)同时计算，如果>σ或者l<itmax，返回；

(4)达到收敛(或者最大迭代数)，收敛或 SA(模拟退火)降温次数达最大，停止SA迭代.

第一步到第三步可以使用不同的脉冲数或者平衡参数重复运行.同样，一个或者多个脉冲位置上的参数也可以很容易的运行.这就允许引入一些从测井资料或者数据处理中得到的先验知识，尤其是对于一些很重要的反射层.因为每一个未知的参数(地震道的样本指数)都有一个相关的搜索范围，在SA迭代过程中很容易的固定一个或者多个脉冲的位置(或者缩小它们的搜索范围).整个过程会减少解的不唯一性和反射系数位置的不确定性. 2 模型实验

如图 1a为一稀疏的反射系数序列(包含12脉冲)，图 1b为一混合相位的地震子波，主频为50 Hz，图 1c为合成地震记录.图 2a为M=15时随机稀疏脉冲反褶积的结果，图 2b为M=20时随机稀疏脉冲反褶积的结果，反演时假设子波是已知的，且和原始子波是相同的.图 2a是经过了100次SA运算，图 2b是经过了120次SA运算，从反演的结果可以看出，真正的反射系数恢复的较好，并且存在较低的不确定性(小的脉冲).M=20时和M=15时相比，他们之间的不同之处不是很明显，这表明M的选择对于随机稀疏脉冲反褶积并没有决定性作用.

图 1 混合相位子波和稀疏反射系数合成地震记录 (a)稀疏反射系数序列；(b)地震子波；(c)合成记录. Fig. 1 The seismograms synthesized with mix phase wavelet and sparse reflection coefficients (a)The sparse reflection coefficients；(b)Seismic wavelet； (c)Seismic record.

图 2 随机稀疏反演的结果 (a)反演得到的反射系数(M=15)； (b)反演得到的反射系数(M=20). Fig. 2 The result of stochastic sparse spike inversion (a)The reflection coefficients got from inversion(M=15)； (b)The reflection coefficients got from inversion(M=20).

对图 1c中的地震记录加上高斯噪声，来考查噪声对随机稀疏脉冲反褶积的结果的影响.如图 3a为加上S/N=10的高斯噪声的地震记录，图 3b为M=20时随机稀疏脉冲反褶积的结果.从反褶积结果上可以看出，噪声对反褶积的结果有一定的影响，反演结果出现小的振铃现象，但是影响结果不是很大.

图 3 含噪地震记录及其反演结果 (a)含噪地震记录；(b)反演得到的反射系数. Fig. 3 The seismic record with noise and its inversion result (a)The seismic record with noise； (b)The reflection coefficients got from inversion.

下面考查地震子波对反演结果的影响.把真实子波进行30度相位旋转后将其作为输入地震子波进行反演.图 4为M=20时随机稀疏脉冲反褶积的结果.从反演的结果上看，虽然原始模型的反射系数大体上恢复出来，但反演结果出现了一些原始反射系数中没有的分量，可见子波对随机稀疏脉冲反褶积有一定的影响，但是影响较小.

图 4 地震子波的相位存在30度误差时的反演结果 Fig. 4 The inversion results when the seismic wavelet has 30 degree error
3 实际资料处理

利用随机稀疏脉冲反褶积对叠后数据进行了随机稀疏脉冲反褶积处理，如图 5所示.可以看出，反褶积之后的分辨率得到明显改善，沉积内幕和接触关系清晰自然.尤为难得的是，由于消除了反褶积之后剩余相位对弱反射信号的干涉影响，1700~1900 ms之间的波阻抗差异较小的弱反射信号得到了有效恢复，为精细地质解释和储层描述提供更加可靠的基础数据.

图 5 随机稀疏脉冲反褶积前(上)后(下)比较 Fig. 5 The contrast of seismic profile before(up) and after(down)stochastic sparse spike deconvolution
4 结 论

随机稀疏脉冲反褶积假设地质模型是稀疏块状的，它通过建立了一个动态的褶积模型，通过采用模拟退火算法和最小二乘法来分别反演反射系数的位置和大小.实质上，该方法的最终目标是估计出一个最简结构的反射系数脉冲.

随机稀疏脉冲反褶积可以看作是一个探测或者估计问题.可以使用VFSA来探测脉冲位置这个高度非线性问题，当然这其中涉及到正则化问题.正则化不仅使得数值解稳定，而且保证反演出来的反射系数序列更具有物理意义.同时，该方法对反演时子波相位不准确时，其反演结果的误差是在允许的范围内的.当子波不准时，可以通过常相位扫描的方法来确定一个准确的子波，来改善反褶积的效果.

该反褶积方法的最大优点是通过使用不同的种子多次运行SA去减少反演的不确定性，并且通过反演初值的临道借用来保证横向连续性和准确性.该方法在反演过程中需要对反演反射系数脉冲有先验认识，这些信息可以通过测井统计信息提供.

致 谢 感谢审稿专家和编辑部老师的指导和帮助.
参考文献
[1] An XF, Liu P, Zhang X.2010.Application of surface consistent deconvolution to seismic prospecting and its effect.Progress in Geography 25(6):2125-2129, doi:10.3969/j.issn.1004-2903.2010.06.032.
[2] Antoine, G, and William, WS.2003.Robust inversion of seismic data using the Huber norm.Geophysics, 68(4):1310-1319.
[3] Cai LF, Tian XM.2012.doi:10.6038/j.issn.1004-2903.2012.03.036.
[4] Dai YS, Peng X, Niu H.2011.Research on high precision seismic estimation based on information feedback.Progress in Geography, 26(3):1004-1009,doi:10.3969/j.issn.1004-2903.2011.03.027.
[5] Gao JH, Zhang M, Zhang B.2012.An improved blind deconvolution based on mutual information rate.Progress in Geography, 26(2):0540-0548, doi:10.3969/j.issn.1004-2903.2011.02.019.
[6] Gao JH, Wang LL, Zhao W.2009.Enhancing resolution of seismic traces based on the changing wavelet model of the seismogram.Chinese J.Geophys, 52(5):1289-1300.A blind seismic deconvolution method based on Particle Swarm Optimization.Progress in Geography, 27(3):1116-1122,
[7] Liu ZW, Wang YC, Zhao HX, et al.2013.High-resolution processing methods of thin interbeds imaging.doi:10.6038/cjg20130429.,
[8] Zuo BX, Hu XY, Han B.2012.The geophysical model enhancement based on the convolution model. Chinese J.Geophys, 55(12):4058-4068, doi:10.6038/j.issn.0001-5733.2012.12.018.,
[9] Dong EQ, Liu GZ, Zhang ZP.2001.Fast implementation technique of adaptive kalman filtering deconvolution via dyadic wavelet transform. Chinese J.Geophys, 44(2):255-261.Chinese J.Geophys, 56(4):1350-1359,
[10] Wang JS, Cao GR.1998.The method of fractal impulse deconvolution.Chinese J.Geophys, 41(1):99-108.
[11] Dong MY, Wang ZH.1991.Adaptive Time-variable deconvolution for seismic data using kalman predictor model.Chinese J.Geophys, 34(2):248-258.
[12] Kaaresen, KF, and Taxt, TM.1998.Multichannel blind deconvolution of seismic signals.Geophysics, 63(6):2093-2017.
[13] Li GF, Qin DH, Peng GX et al.2013.Experimental analysis and application of sparsity constrained deconvolution.Applied Geophysics, 2(10):191-200.
[14] Li GF, Peng GX, Wang WL, et al.2012.Signal-purity-spectrum-based colored deconvolution.Applied Geophysics, 3(9):333-340.
[15] Li ZC, Li D, Wang DY, et al.2013. doi:10.6038/pg20130133.
[16] Li HS, Wu GC, Yin XY.2012.Seismic blind deconvolution method based on generalized Gaussian distribution.Progress in Geography, 27(3):939-944, doi:10.6038/j.issn.1004-2903.2012.03.014.Signal to noise spectrum-constrained adaptive modeling deconvolution.Progress in Geography, 28(1):301-309,
[17] Okaya DA.1995.Spectral properties of the earth's contribution to seismic resolution.Geophysics, 60(1):241-251.
[18] Oldengurg D, Scheuer T, Levy S.1983.Recovery of the acoustic impedance from reflection seismograms.Geophysics, 48(2):1318-1337.
[19] Peng GX, Wang WL, Chen M, et al.2013.Analysis on statistical characteristics of reflection coefficient and its application in high resolution seismic data processing.Journal of oil and gas Technology, 35(4):55-58.
[20] Sacchi MD, Velis DR, Cominguez AH.1994.Minimum entropy deconvolution with frequency-domain constraints.Geophysics, 59(6):938-945.
[21] Sacchi, MD, Ulrych, TJ, and Walker, C.1998.Interpolation and extrapolation using a high resolution discrete Fourier transform.IEEE Trans on Signal Processing, 46(1):31-38.
[22] Sacchi, MD., and Ulych, TJ.2000.Nonminimum-phase wavelet estimation using higher order statistics.The Leading Edge, 19(1):80-83.
[23] Tarantola, A.1984.Inversion of seismic reflection data in acoustic approximation.Geophysics, 49(6):1259-1266.
[24] Walden AT, Hosken JW.1985.A investigation of the spectral properties of primary reflection coefficients.Geophysical Prospecting, 33(2):400-435.
[25] Wang W L, Li GF, Gui JY.2013.Colored deconvolution of mixed phase wavelet.Lithologic reservoirs, 25(3):82-86.
[26] Zhang FC, Liu J, Yin X Y, et al.2008.Modified Cauchy-constrained seismic blind deconvolution.Oil Geophysical Prospecting, 43(4):391-396.
[27] 蔡连芳,田学民.2012.基于PSO的地震盲反褶积.地球物理学进展, 27(3):1116-1122, doi:10.6038/j.issn.1004-2903.2012.03.036.
[28] 戴永寿,彭星,牛慧.2011.一种基于信息反馈的高精度地震子波提取方法.地球物理学进展,26(3):1004-1009,doi:10.3969/j.issn.1004-2903.2011.03.027.
[29] 高静怀,张明,张兵.2012. 一种改进的基于互信息率的盲反褶积方法.地球物理学进展,26(2):540-548, doi:10.3969/j.issn.1004-2903.2011.02.019.
[30] 高静怀,汪玲玲,赵伟.2009.基于反射地震记录变子波模型提高地震记录分辨率.地球物理学报,52(5):1289-1300.
[31] 刘志伟,王彦春,赵会欣,等.2013.薄互层地震成像中高分辨率处理方法.地球物理学报,56(4):1350-1359, doi:10.6038/cjg20130429.
[32] 左博新,胡祥云,韩波.2012.基于褶积模型的地球物理反演模型增强.地球物理学报,55(12):4058-4068, doi:10.6038/j.issn.0001-5733.2012.12.018.
[33] 董恩清,刘贵忠,张宗平.2001.自适应Kalman滤波反褶积的快速实现方法.地球物理学报,44(2):255-261.
[34] 王嘉松,曹桂荣.1998.分形脉冲反褶积方法.地球物理学报,41(1):99-108.
[35] 董敏煜,王振华.1991.地震资料自适应变卡尔曼反褶积.地球物理学报,34(2):248-258.
[36] 李振春,李栋,王德营,等.2013.信噪比谱约束的自适应谱模拟反褶积方法研究.地球物理学进展,28(1):301-309, doi:10.6038/pg20130133.
[37] 李海山,吴国忱,印兴耀.2012.基于广义高斯分布的地震盲反褶积方法研究.地球物理学进展,27(3):939-944, doi:10.6038/j.issn.1004-2903.2012.03.014.
[38] 彭更新,王万里,陈猛,等.2013.反射系数统计特征分析及其在高分辨率地震资料处理中的应用.石油天然气学报,35(4):55-58.
[39] 王万里,李国发,桂金咏.2013.混合相位子波有色反褶积.岩性油气藏,25(3):82-86.
[40] 张繁昌,刘杰,印兴耀,等.2008.修正柯西约束地震盲反褶积方法.石油地球物理勘探,43(4):391-396
[41] 安西峰,刘平,张旭.2010.地表一致性反褶积在地震勘探中的应用及效果.地球物理学进展,25(6):2125-2129, doi:10.3969/j.issn.1004-2903.2010.06.032.