基于L0范数稀疏约束的地震数据反褶积

梁东辉，陈生昌

(浙江大学地球科学系，浙江杭州310027)

随着油气勘探进程的不断推进，需要寻找地质构造复杂、规模较小及埋藏较深的油气藏，这就要求地震勘探资料要有更高的分辨率。反褶积是地震数据处理中的一个基本环节，它的作用是压缩地震记录中的地震子波，同时可以压制鸣震和多次波，提高地震记录的纵向分辨率[1]。实际地震记录中不可避免的包含有未知噪声，所以反褶积是一个病态的反问题，需要采取正则化方法来求解。求解此类问题时常用的阻尼最小二乘法(预白化法)在子波自相关矩阵的对角线上加一个小的扰动量，等价于零阶二次正则化，虽然这样能够很好地提高算法的稳定性，但这是一种较宽泛的二次约束，求得的解是连续的、光滑的，分辨率不够高。为了提高地震勘探资料的分辨率，可运用稀疏优化手段，使反褶积结果稀疏化，这就是稀疏约束反褶积。传统的反褶积方法需要假设子波为最小相位和反射系数是白噪声，而实际地震资料中子波通常不是最小相位的，反射系数也是非白噪的[2-3]，这就导致传统的反褶积效果不理想。稀疏约束反褶积则避免了这两种假设，它利用反射系数由稀疏脉冲序列组成这一先验信息，从带有噪声的地震记录中提取出反射系数的时间位置和振幅值，同时去除噪声,可得到比常规反褶积具有更高分辨率的反演结果。

稀疏约束反褶积方法有很多，前人已做了很多研究，如L1范数稀疏约束反褶积[4-5]，Lp范数反褶积[6]，L1-L2范数联合约束稀疏脉冲反演的应用[7]，基于柯西准则的预条件共轭梯度法[8-11]，地震反褶积中的重加权策略[12]等。柯西准则约束稀疏反褶积对于较大的反射系数能准确提取，但对于较小的反射系数会有压制作用[13]。L1范数约束放松了对解的稀疏性的约束，所求的并不是最稀疏的解，并且对于较长的信号其计算复杂度较高[14]。

稀疏约束反褶积由两部分组成。一是对噪声的约束:实际地震记录中的噪声一般认为是高斯分布的，子波与所求得的反射系数褶积得到正演地震记录，对噪声的约束即是使实际地震记录与正演地震记录的残差平方和最小。二是对反射系数的约束:这部分用来限制反射系数的先验概率分布情况。对于一个向量来说，L0范数就是其中非零元素的个数，因此是向量稀疏性的最佳度量，L0范数稀疏优化已在地震数据压缩重构中得到很好的应用[15-16]。层状地层模型的地下地层反射系数是一系列稀疏脉冲，从理论上来讲，对反射系数采取L0范数约束能求得最为稀疏的解，可以有效地提高地震资料的分辨率，有利于后期的处理和解释。

我们将L0范数稀疏约束引入到地震资料反褶积处理中，并通过两组模型试验，对比柯西准则，L1范数和L0范数稀疏约束反褶积的效果。

1 方法原理

由震源激发的脉冲信号经过大地滤波器的作用变成一个长度为几十毫秒的波形，称为地震子波。根据褶积模型，地震记录相当于子波和反射系数序列的褶积，又因为野外实际地震数据采集时不可避免地会含有噪声，故地震记录可表示为

(1)

式中:d表示地震记录;ω0表示地震子波;r表示反射系数序列;n表示噪声。可写成如下矩阵形式:

(2)

这里ω是由地震子波ω0组成的褶积矩阵。

由于地震子波的影响，来自相邻界面的反射波重叠在一起，难以区分开来，因此需要压缩子波来提高地震勘探资料的纵向分辨率。传统的反褶积方法是设计一个与子波性质相反的滤波因子，然后与地震记录褶积，得到期望的窄脉冲，这个过程称为反褶积。

1.1 约束优化

稀疏约束反褶积运用正则化方法直接求解方程(2)，假设子波已知或已准确提取出，优化目标函数的一般形式可表示为

(3)

其中，λ>0，是反射系数稀疏性与准确性之间的权衡因子。该目标函数前一部分(J0)对残差进行约束，保证所求的结果与地震记录数据符合，也就是对噪声的约束,这部分一般采用最小平方约束;后一部分(λJr)称为正则化项或罚函数，是为了对反射系数进行稀疏约束，保证求解结果稀疏，不同的稀疏约束方法选取不同的罚函数。

也可以从贝叶斯定理出发，给定噪声的分布概率和反射系数的先验分布概率，然后推导出相应的稀疏约束形式。贝叶斯定理可以表示为

(4)

其中，P(r|d)表示反射系数的后验概率分布;P(d|r)表示给定反射系数时观测数据的符合程度，称为似然函数，等于噪声的分布概率;P(r)表示反射系数的先验概率分布。

假设噪声满足高斯分布，且均值为零，标准差为σn，单独一个噪声分量的概率密度函数为

(5)

噪声整体的分布概率为

(6)

又因为

(7)

所以似然函数

(8)

(9)

反射系数总体的分布概率为

(10)

根据(4)式，反射系数的后验概率分布为

(11)

两边同时取对数，得

(12)

使后验概率分布最大，等价于求解优化函数:

(13)

(13)式就是柯西准则稀疏约束反褶积的优化目标函数。

同理，如果我们假设噪声满足高斯分布，均值为0，反射系数满足均值为0的拉普拉斯分布，可根据贝叶斯定理推导出L1范数稀疏约束反褶积的优化目标函数:

(14)

1.2 L0范数约束反褶积

稀疏约束反褶积希望求得稀疏的反射系数，即较少的非零元素值。反射系数序列的L0范数就是其中非零元素的个数，因此采用L0范数稀疏约束来求解反褶积问题是一种理想的策略，理论上能求得最为稀疏的解。L0范数要解决的优化问题是

(15)

其中，δ是与噪声有关的参数，表示对噪声水平的估计量。(15)式表示在残差满足一定条件时，要求反射系数序列中非零元素数量最少。通过拉格朗日乘子法可将(15)式转化为无约束问题，构建L0范数稀疏约束反褶积优化目标函数:

(16)

如果按照贝叶斯定理，假设噪声满足高斯分布，反射系数为稀疏分布，同样可推导出(16)式。解此优化问题的方法很多，如迭代重加权最小二乘(IRLS)算法[17]、基追踪(BP)算法[18]、贪婪算法[19-20]、迭代硬阈值法等[21]。IRLS算法构建加权的二范数来近似零范数，需要进行矩阵求逆运算，不适合大规模数据的处理;BP算法将L0范数最小优化转化为L1范数最小优化问题，放松了对稀疏性的限制并且计算量较大;贪婪算法，如匹配追踪、正交匹配追踪等，求解精度较低，且不适合处理欠稀疏数据;迭代硬阈值计算形式简单，计算量小，有利于处理大数据量的地震勘探资料。本文选取迭代硬阈值求解L0范数稀疏约束优化问题，其迭代格式为

医学和药学是实践性非常强的学科，课堂学生内需要接受的知识繁多枯燥，以教师、课堂、教材为中心的教学方式无法满足学生可持续发展的需要。如果没有课外科技活动作为学习内容和形式的补充，学生的学习积极性会受到一定的影响。因此，在医药院校中，学生课外科技活动在人才培养中占有重要地位，它能在校内提供学生理论联系实际的机会，培养学生的创新思维，提高人才培养质量，是教育行业“供给侧结构性改革”的途径之一。

(17)

其中，μ为步长。为了保证迭代收敛[22]，步长的选取范围要满足

(18)

Ha为硬阈值函数:

(19)

a>0表示阈值，计算时先给定一个较大的阈值:

(20)

每次迭代后降低阈值，再进行下一次迭代。

因为实际地震记录中的噪声可认为是均值为0，标准差为σn的高斯分布噪声。

(21)

迭代终止的条件可设定为

(22)

其中，Nd表示观测值d的个数。

L0范数稀疏约束反褶积的计算步骤如下:

1) 给定反射系数初值r0为0，估计噪声的标准差σn，迭代次数n=0;

2) 选取合适的步长μ及阈值初值;

3) 按照(17)式由rn计算rn+1;

4) 阈值下降，迭代次数n增加1;

5) 重复步骤3)和步骤4),直至满足迭代终止条件，输出反褶积结果。

2 模型试验

选取延迟30ms，主频为25Hz的雷克子波(图1)作为下面两个模型试验中的地震子波，采样间隔为2ms。

图1 试验中所使用的子波

2.1 模型试验一

构建一个简单的稀疏分布反射系数序列(图2a)，与图1中的子波褶积并加入标准差为0.01的高斯噪声，得到相应的合成地震记录(图2b)。

分别采用柯西准则约束、L1范数约束和L0范数约束进行反褶积计算，结果如图3所示。对比图3a,图3b,图3c和图3d可见，上面两个椭圆中有几个较弱的反射系数，柯西准则约束反褶积结果未能显示，L1范数约束能显示但振幅值偏弱，L0范数约束全部显示且较清晰;下面两个椭圆中各有1个较弱的反射系数，柯西准则和L1范数约束反褶积结果未能显示，L0范数约束结果中有显示。说明L0范数对弱反射系数保护能力最强。

图2 简单的稀疏反射系数模型(a)及其合成地震记录(b)

图3 稀疏反射系数序列(a)和柯西准则(b)，L1范数(c)，L0范数(d)稀疏约束反褶积结果

2.2 模型试验二

构建一个复杂的反射系数模型(图4a)，与图1所示子波褶积并加入标准差为0.01的高斯噪声，得到相应的地震记录(图4b)。

分别采用柯西准则约束、L1范数约束和L0范数约束进行反褶积计算，结果如图5a，图6a和图7a 所示。红色表示正的反射系数，蓝色表示负的反射系数。3张图中黑色方框对应的放大图示分别为图5b，图6b和图7b,图中箭头所指处有3条较弱的反射系数同相轴，在图5b中基本没有显示，在图6b 和图7b中有显示，说明柯西准则约束对弱反射系数有压制作用，L1范数和L0范数约束能够较好地保护弱反射系数。

图5b，图6b和图7b中黑色椭圆内有两条蓝色的反射系数同相轴，在图5b中显示不全，图6b中显示的振幅较弱，在图7b中清晰显示。图8为3种稀疏约束反褶积结果的波形剖面。图8a,图8b 和图8c中的红色椭圆内有两条负的(向左)反射系数同相轴，可以明显看出，柯西准则稀疏约束反褶积(图8a)只显示了左边一部分;L1范数稀疏约束反褶积(图8b)显示较全，但幅值较弱;L0范数稀疏约束反褶积(图8c) 清晰显示了这两条负的反射系数同相轴。

表1为模型试验二中3种反褶积方法的迭代次数与计算时间。可以看到，L0范数约束具有迭代次数少、计算时间短的优点，有利于处理大规模地震勘探数据。柯西准则约束只针对一维的地震记录，因此需要逐道进行反褶积运算，而利用L1范数和L0范数约束进行反褶积运算时，可以直接处理一个二维的地震记录，各道之间互不干扰。

表1 模型试验二3种方法的迭代次数与计算时间

图4 复杂的稀疏反射系数模型(a)及其合成地震记录(b)

图5 柯西准则稀疏约束反褶积a 反褶积结果; b 图5a黑色方框部分的放大显示

图6 L1范数稀疏约束反褶积a 反褶积结果; b 图6a黑色方框部分的放大显示

图7 L0 范数稀疏约束反褶积a 反褶积结果; b 图7a黑色方框部分的放大显示

图8 复杂的反射系数模型3种稀疏约束反褶积结果的波形剖面a 柯西准则约束; b L1范数约束; c L0范数约束

3 结束语

我们将L0范数稀疏约束引入地震勘探数据的反褶积处理，并与柯西准则约束和L1范数约束反褶积结果进行对比。柯西准则约束反褶积对于较弱的反射系数有压制作用；L1范数约束反褶积对弱反射系数的保护比柯西准则强，但是求得的反射系数振幅可能会偏弱，且迭代次数较多。L0范数约束反褶积对弱反射系数的保护最好，有利于提高地震勘探资料的精度，并且具有迭代次数少、计算时间短的优点。地震勘探需要处理的数据量很大，因此在进行反褶积运算时以采用L0范数约束可以提高处理效果和计算效率。

本文的模型试验是对合成地震记录进行的，子波是已知的。在实际数据的反褶积处理过程中，一个至关重要的问题是如何从含噪地震道中提取准确的高精度地震子波，这将是下一步的研究方向。

参 考 文 献

[1] 牟永光.地震数据处理方法[M].北京:石油工业出版社,2007:57-79

Mou Y G.Processing methods of seismic data[M].Beijing:Petroleum Industry Press,2007:57-79

[2] 吴常玉,杨瑞娟,鲍峥,等.基于负熵的地震盲反褶积方法及其应用[J].石油物探,2009,48(3):232-238

Wu C Y,Yang R J,Bao Z,et al.Seismic blind deconvolution method based on negative entropy and its application[J].Geophysical Prospecting for Petroleum,2009,48(3):232-238

[3] 邬世英,孙赞东,朱兴卉.动态反褶积中的反射系数序列时频特征研究[J].石油物探,2011,50(4):324-330

Wu S Y,Sun Z D,Zhu X H.Research on the time-frequency spectrum of reflectivity sequence based on dynamic deconvolution[J].Geophysical Prospecting for Petroleum,2011,50(4):324-330

[4] Taylor H L,Banks S C,Mccoy J F.Deconvolution with the1norm[J].Geophysics,1979,44(1):39-52

[5] Dossal C,Mallat S.Sparse spike deconvolution with minimum scale[C]∥Signal Processing with Adaptive Sparse Structured Representations.Rennes,France:SPARS workshop,2005:1-4

[6] Debeye H,Riel V P.Lp‐norm deconvolution1[J].Geophysical Prospecting,1990,38(4):381-403

[7] 王宇,韩立国,周家雄,等.L1-L2范数联合约束稀疏脉冲反演的应用[J].地球科学:中国地质大学学报,2009,34(5):835-840

Wang Y,Han L G,Zhou J X,et al.Application of combined norm constrained sparseness spike inverse[J].Earth Science—Journal of China University of Geosciences,2009，34(5):835-840

[8] 刘喜武,刘洪.实现稀疏反褶积的预条件双共轭梯度法[J].物探化探计算技术,2003,25(3):215-219

Liu X W,Liu H.The preconditional dual conjugate gradient algorithm for sparse deconvolution[J].Computing Techniques for Geophysical and Geochemical Exploration,2003,25(3):215-219

[9] 朱振宇,刘洪.稀疏反褶积方法及其应用[J].石油大学学报:自然科学版,2005,29(6):20-22

Zhu Z Y,Liu H.Sparse deconvolution method and its application[J].Journal of the University of Petroleum,2005,29(6):20-22

[10] 孟小红,吴何珍,刘国峰.盲源反褶积方法与应用研究[J].石油地球物理勘探,2005,40(6):642-645

Meng X H,Wu H Z,Liu G F.Study of blind deconvolution and application of method[J].Oil Geophysical Prospecting,2005,40(6):642-645

[11] 刘喜武,宁俊瑞,张改兰.Cauchy稀疏约束Bayesian估计地震盲反褶积框架与算法研究[J].石油物探,2009,48(5):459-464

Liu X W,Ning J R,Zang G L.Cauchy sparse constrained bayesian estimation based seismic blind deconvolution frame and algorithm[J].Geophysical Prospecting for Petroleum,2009,48(5):459-464

[12] Sacchi M D.Reweighting strategies in seismic deconvolution[J].Geophysical Journal International,1997,129(3):651-656

[13] 张繁昌,刘杰,印兴耀,等.修正柯西约束地震盲反褶积方法[J].石油地球物理勘探,2008,43(4):391-396

Zang F C,Liu J,Yin X Y,et al.Modified cauchy-constrained seismic blind deconvolution[J].Oil Geophysical Prospecting,2008,43(4):391-396

[14] 焦李成,杨淑媛,刘芳,等.压缩感知回顾与展望[J].电子学报,2011,39(7):1651-1662

Jiao L C,Yang S Y,Liu F,et al.Development and prospect of compressive sensing[J].Chinese Journal of Electronics,2011,39(7):1651-1662

[15] 曹静杰,王彦飞,杨长春.地震数据压缩重构的正则化与零范数稀疏最优化方法[J].地球物理学报,2012,55(2):596-607

Cao J J,Wang Y F,Yang C C.Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization[J].Chinese Journal of Geophysics,2012,55(2):596-607

[16] 陈国新,陈生昌,王汉闯,等.基于 L0 范数最小化的地球物理数据稀疏重构[J].应用地球物理，2013,10(2):181-190

Chen G X,Chen S C,Wang H C,et al.Geophysical data sparse reconstruction based on L0-norm minimization[J].Applied Geophysics,2013,10(2):181-190

[17] Wohlberg B,Rodr guez P.An iteratively reweighted norm algorithm for minimization of total variation functionals[J].Signal Processing Letters,IEEE,2007,14(12):948-951

[18] Chen S S,Donoho D L,Saunders M A.Atomic decomposition by basis pursuit[J].SIAM Journal on Scientific Computing,1998,20(1):33-61

[19] Mallat S G,Zhang Z.Matching pursuits with time-frequency dictionaries[J].IEEE Transactions on Signal Processing,1993,41(12):3397-3415

[20] Tropp J A,Gilbert A C.Signal recovery from random measurements via orthogonal mathing pursuit[J].IEEE Transactions on Information Theory,2007,53(12):4655-4666

[21] Blumensath T,Davies M E.Iterative thresholding for sparse approximations[J].Journal of Fourier Analysis and Applications,2008,14(5):629-654

[22] Beck A,Teboulle M.A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J].SIAM Journal on Imaging Sciences,2009,2(1):183-202