基于Bregman迭代的不...值模型地震数据重构方法分析_赵虎.pdf

2023 年第 38 卷 第1期2023，38(1):0409-0418地球物理学进展Progress in Geophysicshttp:/wwwprogeophyscnISSN 1004-2903CN 11-2982/P赵虎，赵子涵，陈伟，等 2023 基于 Bregman 迭代的不同阈值模型地震数据重构方法分析 地球物理学进展，38(1):0409-0418，doi:10 6038/pg2023FF0310ZHAO Hu，ZHAO ZiHan，CHEN Wei，et al 2023 Analysis of seismic data reconstruction methods for different threshold models based on Bregmaniteration Progress in Geophysics(in Chinese)，38(1):0409-0418，doi:106038/pg2023FF0310基于 Bregman 迭代的不同阈值模型地震数据重构方法分析Analysis of seismic data reconstruction methods for different threshold modelsbased on Bregman iteration赵虎1，2，赵子涵2，陈伟3，邸志欣4，韩嵩5，胥良君6，杨暾2ZHAO Hu1，2，ZHAO ZiHan2，CHEN Wei3，DI ZhiXin4，HAN Song5，XU LiangJun6，YANG Tun2收稿日期2021-10-19;修回日期2022-08-31投稿网址http:/www progeophys cn基金项目中国石油-西南石油大学创新联合体科技合作项目(2020CX010201)和国家自然基金青年科学基金(41704134)联合资助第一作者简介赵虎，男，1983 年生，博士，教授，现在西南石油大学从事地震采集方法与解释教学和科研工作E-mail:cumtzhaohu163 com1 西南石油大学天然气地质四川省重点实验室，成都6105002 西南石油大学地球科学与技术学院，成都6105003 中石油西南油气田公司勘探事业部，成都6100414 中石化石油工程地球物理有限公司科技研发中心，南京2100005 中石油西南油气田公司勘探开发研究院，成都6100416 中石油西南油气田分公司重庆气矿，重庆4007071 Sichuan Province University Key Laboratory of Natural Gas Geology，SWPU，Chengdu 610500，China2 School of Geoscisence and Technology，Southwest Petroleum University，Chengdu 610500，China3 Petrochina Southwest Oil and Gas field Company Exploration Division，Chengdu 610041，China4 Sinopec Geophysical Corporation，esearch Center，Nanjing 210000，China5 Exploration and Development esearch Institute，PetroChina Southwest Oil and Gasfield Company，Chengdu 610041，China6 PetroChina Southwest Oil and Gasfield Company Chongqing Gas District，Chongqing 400707，China摘要传统的地震数据采样必须严格遵循 Nyquist采样定理，而野外实际数据的采集可能由于施工条件或者地表障碍物的限制，不一定能记录到完整的地震波场，所以地震资料处理中的数据重建是非常重要的问题 压缩感知理论最先来自信号处理领域，它所包括的问题类型有信号的稀疏表征和数学组合优化，它给地震数据重建这类问题指明了思考方向 而其中如何选择最优的迭代算法是数据重建中的关键问题 本文将地震数据插值问题归纳到约束最优化问题，选择能有效稀疏表征地震波场的傅里叶变换，对于压缩感知理论框架下的混合范数反问题，再用 Bregman 迭代方法去求解，在地震数据的重建过程中，传统的阈值参数收敛慢，为了降低迭代次数并且提高地震数据恢复的AbstractTraditional seismic data sampling must strictlyfollow the Nyquist sampling theorem，and actual field dataacquisition may not be able to record the complete seismicwave field due to construction conditions or obstacles onthe ground Therefore，data reconstruction in seismic dataprocessing is an important factor problem Compressedsensingtheoryoriginatesfromthefieldofsignalprocessing It includes sparse representation of signals andmathematicalcombinationoptimizationproblemsItprovides a new solution for the reconstruction of seismicdata And how to choose the optimal iterative algorithm isthekeyissueindatareconstructionThispapersummarizes the seismic data interpolation problem into aconstrainedoptimizationproblem，choosesFouriertransform，which can effectively represent the seismic wavefield sparsely，and then applies the Bregman iterativemethod to solve the mixed norm inverse problem under the地球物理学进展www progeophys cn2023，38(1)精度，总结出改进型指数衰减规律的阈值参数，选择用硬阈值算子来重建恢复地震数据 通过对理论模型和实际地震资料的处理结果表明该方法可以快速、有效的恢复地震波场的缺失数据关键词地震数据插值;Bregman 迭代;阈值参数;压缩感知中图分类号P631文献标识码Adoi:10 6038/pg2023FF0310frameworkofcompressedsensingtheoryInthereconstructionprocess，forthedisadvantageofslowconvergenceoftraditionalthresholdparameters，thethreshold parameters of the improved exponential decay laware summarizedIn order to reduce the number ofiterations andimprovetheaccuracyofseismicdatareconstruction，the hard threshold operator is selected toperform the two-dimensional seismic dataebuild andrestore The processing results of theoretical model andactual seismic data show that this method can quickly andeffectively restore the missing data of seismic wave fieldKeywordsSeismic data interpolation;Bregman iteration;Threshold parameter;Compressed sensing0引言随着我国油气勘探程度的提高，如何提高复杂勘探区的地震勘探效果成为下一步工作的关键 在野外地震数据采集的进程中，因为受限于采集环境和经济条件，导致某些位置地震数据缺失，使得整个地震剖面不完整 从而会对室内资料处理和解释造成不良影响，因此如何进行地震数据规则化重建显得尤为关键地震数据重建问题，最早由 Larner 等(1981)提出 他最开始对不完整地震道的恢复进行深入的研究，之后陆续出现了不同的地震数据重建方法;onen(1987)提出把缺失道视为零道再代入波动方程部分偏移的叠前地震数据重构方法;Spitz(1991)提出了线性同相轴在 f-k 域是能够被预测的地震数据重建方法;随后 Porsani(1999)对这一方法提出了改进，得出 f-k 域半步长预测滤波的地震道插值方法;Naghizadeh 和 Sacchi(2007)提出了基于多步自回归预测滤波的不规则缺失道插值重建;Kaplan 等(2010)采用最小平方反演方法进行数据重建 不难发现以上方法基本都是针对均匀采样数据进行重建而近年来提出的压缩感知理论对于解决此类问题 指 明 了 新 的 方 向，压 缩 感 知 由 EmmanuelCandes、omberg 和 Donoho 等人在 2004 年正式提出 基于压缩感知的地震勘探技术，目前国内外已经做了一些研究，Herrmann 等(2008)、Herrmann(2010)提出了一种基于 Curvelet 的多尺度非线性地震数据处理方法 Candes 等(2006)基于鲁棒不确定性原理对高度不完整频率信息的精确信号重建进行了研究 Hennenfent 和 Herrmann(2008)提出了基于随机抖动的欠采样方案，对相邻观测点之间的最大距离进行控制，便于实际施工 Wu 等(2009)用 dreamlet 方 法 在 压 缩 域 进 行 成 像Naghizadeh 和 Sacchi(2010)研究了采样函数和傅里叶重建方法 Mosher 等(2012b)首次提出并且完善了一种根据约束条件选择最佳炮点和检波点位置的非均匀采样方法，用于观测系统设计，取得了显著的成效 Mosher 等(2012a)对压缩感知成像这一方向有一个更深层次的探究 Brown 等(2017)在2015 年将非规则排列观测系统应用在阿拉斯加地区的陆上冰雪覆盖区域，通过利用数据重建的方法，得到一个相比传统采集方法更好的高分辨率地震资料数据重构中最为关键的问题是求解策略，常用方法是基于某种变换，Xu 等(2005)提出了基于反泄露傅里叶变换的地震数据重建方法 Jin(2010)提出了基于阻尼最小范数 Fourier 反演，实现了五维地震数 据 的 重 建 Herrmann 在 2008 年 提 出 基 于Curvelet 数据重建的稀疏促进反演方法 Liu 和Fomel(2010)提出稀疏多尺度变换域中的 OC-seislet变换重建缺失的地震数据并消除随机噪声 Yang 和Fomel(2015)阐述了用快速 Fourier 反演策略代替非均匀快速 Fourier 变换的观点，采用共轭梯度法来优化最优化问题，用来提高计算的效率 还有复值曲波变换(徐卫等，2016)和 Shearlet 稀疏变换基(王常波，2018)对于地震数据重建的迭代算法，主要有凸集投影算法，最小加权范数法，不准确 Uzawa 方法，两步法以及迭代阈值法(IST)等 但是，已有的迭代算法往往难以同时兼顾计算精度和效率 为此，本文拟引用 Bregman 迭代算法，分析不同迭代算法的重建特性，通过提出改进型指数衰减规律的阈值参数来改善在地震数据恢复过程中，传统的阈值参数收敛得较为缓慢这一不佳的表现，运用理论模型和实际地震资料经过方法验证，有效解决地震数据缺失数据重建问题0142023，38(1)赵虎，等:基于 Bregman 迭代的不同阈值模型地震数据重构方法分析(www progeophys cn)1方法原理1.1理论基础压缩感知理论和传统的 Nyquist 采样不同，它是同时在进行着采样和压缩这两个步骤，不被Nyquist 采样定理所约束，而是受控于稀疏性和非相干性这两个基本准则 采用稀疏变换去构建数据空间和模型空间之间的联系 对地震数据进行重构，就是把残缺的地震数据构建成完整的数据，模型如下:d=m，(1)其中，d 为采集到的地震数据，m 为插值重建后的地震数据，为稀疏变化算子 利用信号具有稀疏性的先验条件，通过一些数学变换(Curvelet 变换、Seislet 变换和 Fourier 变换等)来构建数据空间和模型空间的联系 由于公式(1)为数学 g 问题，要通过增加另外的正则化条件来求得合理的 m1.2Bregman 迭代重建算法根据上述分析可得出，压缩感知方法是利用最小化策略这一关键的技术来对前面的欠定问题进行求解的 与 L2准则不同，压缩感知理论是通过稀疏反变换算子来构建一个稀疏的数据空间 m(其中m=Tm)，用来解决下面的约束最优化问题:minmm1s tKTm=d(2)该模型是利用求解带约束的 L1范数最小化问题，来求取线性系统 KTm=d 的稀疏解 通常时候，考虑用到最小二乘法来求解此问题:m=argminmKTm d22(3)将公式(2)的带约束最小化问题转化为公式(3)的无约束问题后，需要在公式(3)的右边添加正则项 m1:minm m1+12KTm d22，(4)其中 为惩罚参数 但事实上，公式(4)与公式(2)并不等价，为了满足 KTm=d 的条件，需要取为较小值Bregman 迭代方法是 Osher 等(2005)于图像处理领域提出来的方法，该方法的迭代框架为:dk+1=dk+(d KTmk)mk+1=argminm m1+12KTm dk+122(5)求解公式(5)里面的非约束问题能够采用迭代硬阈值方法，再结合稀疏正变换算子，用于解决数据插值重建问题的 Bregman 迭代方法为:dk+1=dk+(d KTmk)mk+1=S T1dk+1，(6)其中，m0=0，d0=0，S 为硬阈值算子 这样，最终插值结果 m 即可通过 m=Tmk+1来求得 应用Bregman 迭代方法时，关键问题在于公式(6)中阈值参数的选取，接下来，我们讨论该参数的选取方法图 1阈值模型迭代次数变化图Fig 1Threshold model iteration number change graph图 2四种阈值参数的评价函数 J 曲线Fig 2Evaluation function J curve of four threshold parameters图 3重构数据的评价函数 J 值1 为线性、2 为指数型、3 为反比例型、4 为改进型指数Fig 3Evaluation function J value of reconstructed data1 is linear，2 is exponential，3 is inverse proportional，4 is improved exponent114地球物理学进展www progeophys cn2023，38(1)图 4(a)理论模型;(b)理论数据 FK 谱Fig 4(a)Theoretical model;(b)The FK spectrum of theoretical data图 5(a)50%随机欠采样数据;(b)欠采样数据 FK 谱Fig 5(a)50%random undersampled data;(b)The FK spectrum of undersampled data图 6(a)线性阈值参数重建结果;(b)重建结果 FK 谱Fig 6(a)Linear threshold parameter reconstruction results;(b)The FK spectrum of reconstruction result2142023，38(1)赵虎，等:基于 Bregman 迭代的不同阈值模型地震数据重构方法分析(www progeophys cn)图 7(a)指数型阈值参数重建结果;(b)重建结果 FK 谱Fig 7(a)Exponential threshold parameter reconstruction results;(b)The FK spectrum of reconstruction result图 8(a)反比例型阈值参数重建结果;(b)重建结果 FK 谱Fig 8(a)econstruction results of inverse proportional threshold parameters;(b)The FK spectrum of reconstruction result图 9(a)改进型指数阈值参数重建结果;(b)重建结果 FK 谱Fig 9(a)Improved index threshold parameter reconstruction results;(b)The FK spectrum of reconstruction result314地球物理学进展www progeophys cn2023，38(1)图 10不同阈值参数重建结果误差图(a)线性阈值参数重建误差;(b)指数型阈值参数重建误差;(c)反比例型阈值参数重建误差;(d)改进型指数阈值参数重建误差Fig 10Error graph of reconstruction results with different threshold parameters(a)Linear threshold parameter reconstruction error;(b)Exponential threshold parameter reconstruction error;(c)Inverse proportionalthreshold parameter reconstruction error;(d)Improved exponential threshold parameter reconstruction error1.3阈值模型分析合理的阈值参数，可以提高计算效率，以下分别对传统阈值参数和改进型指数阈值参数进行讨论(1)线性阈值模型，表达式为:Sk=Smax(k 1)S，k 1，2，N，(7)其中 Smax和 Smin分别为不规则地震数据在频率域中设定的最大和最小傅里叶系数，S1=Smax，SN=Smin保证 Smin0，步长 S=(Smax Smin)/(N 1)N 为最大迭代次数(2)指数型阈值模型，表达式为:Sk=Smaxeb(k 1)，k=1，2，N，(8)其中 b=In(Smax/Smin)/(N 1)(3)反比例型阈值模型，表达式为:Sk=aSmax/k+b，k=1，2，N，(9)其中 a=N(Smax Smin)/Smax(N 1)，b=(NSminSmax)/(N 1)(4)改进型指数阈值参数，表达式为:Sk=Smaxebi 1N1，k=1，2，N，(10)其中 b=In(Smax/Smin)图 1 分别表示四种阈值参数的迭代次数变化图，可以看出四种阈值参数随着迭代次数的增加收敛速度不同为了充分论证各阈值模型的收敛效果和精度，这里来构造一个评价函数 J 进行分析:Jk=dk dk12dk12，k=1，2，N，(11)4142023，38(1)赵虎，等:基于 Bregman 迭代的不同阈值模型地震数据重构方法分析(www progeophys cn)图 11(a)含噪地震数据;(b)含噪地震数据 FK 谱Fig 11(a)Noisy seismic data;(b)FK spectrum of noisy seismic data图 12(a)60%随机欠采样数据;(b)欠采样数据 FK 谱;(c)重建数据;(d)重构数据的 FK 谱Fig 12(a)60%random undersampled data;(b)FK spectrum of undersampled data;(c)econstructed data;(d)FK spectrum of reconstructed data514地球物理学进展www progeophys cn2023，38(1)图 13非均匀采样数据重构(a)复杂原始数据;(b)非均匀采样数据;(c)重建后数据Fig 13Non-uniformly sampled data reconstruction(a)Complex raw data;(b)Non-uniformly sampled data;(c)econstructed data式中:dk 1为第 k 1 次迭代的重构数据，N 为迭代总次数，dk为第 k 次迭代的重构数据利用该函数评价阈值模型的收敛速度和收敛精度如图 2 和图 3 所示 通过分析可知，本文所提出的改进型指数阈值参数重构误差较少，精度较高，而且计算效率也更快，在处理海量数据时更有优势各阈值模型对算法的影响总结如表 1 所示表 1各阈值参数的效果对比Table 1Comparison of the effect of eachthreshold parameter阈值类型 线性指数型反比例型改进型指数收敛速度 最慢一般，优于线性最快低于反比例型收敛精度 一般差于改进型指数 一般，优于线性最佳2模型数据分析可以把不规则地震数据分成两种类型:(1)均匀道缺失地震数据，它具有空间上规则采样，但存在地震道随机缺失现象的特点;(2)非均匀缺失地震数据，它在空间上是不规则采样，且存在道缺失现象 首先讨论均匀道缺失地震数据的重建 图 4a、b是用频率为 30 Hz 的雷克子波合成的 300 道二维地震数据及其 FK 谱，每道有 5000 个采样点，其中采样间隔是 0.2 ms，数据合成时每一层的反射能量都存在一些差异 图 5a、b 是对理论数据进行 50%随机欠采样后的不完整地震记录及其 FK 谱，可以看出欠采样数据 FK 谱中产生了很多噪声对于以上理论二维地震数据，拟通过采用硬阈6142023，38(1)赵虎，等:基于 Bregman 迭代的不同阈值模型地震数据重构方法分析(www progeophys cn)图 14实际采集地震数据重构对比(a)原始数据;(b)60%随机欠采样数据;(c)重建后数据;(d)重建误差Fig 14Comparison of reconstruction of actual acquired seismic data(a)aw data;(b)60%random undersampled data;(c)econstructed data;(d)econstruction error值算子进行重构 迭代次数为50 次 图6 图9 分别为采用线性、指数型、反比例型、改进型指数阈值重建的结果以及其 FK 谱 图 10 为各阈值参数的数据重构结果和原始理论数据的误差图 四种阈值模型所求结果的信噪比如表 2 所示表 2不同阈值模型重建的信噪比Table 2Signal-to-noise ratio reconstructedby different threshold models阈值模型信噪比/dB线性15.3指数型20.2反比例型17.1改进型指数22.5从上述的结果能够得出，改进型指数阈值参数重建误差相对较少，计算效率也较快，这与前文阈值参数的讨论结果一致，改进型指数阈值参数更具有优势因此以下重建过程都采用该阈值模型进行计算 为了检验该阈值模型的抗噪能力，在理论模型(图 4a所示)里面添加了随机噪声，如图11 所显示，接着再加入抗噪模拟实验 图 12a、b 为 60%随机欠采样的实验数据及其 FK 谱;图 12c、d 为重构结果及其 FK谱，重构后数据信噪比为 12.6 dB 通过对比两幅图可以看出，含噪不规则缺失地震数据的所有信息都有了一个很好的重构，前后的有效信息损伤较小，消除了 FK 谱图中存在的假频噪声，这表明数据重构的效果较好，同时具有良好的抗噪能力以上分析了均匀采样的不规则数据重构，而对714地球物理学进展www progeophys cn2023，38(1)于非均匀采样的不规则地震数据本文方法是否适应，建立了如图 13a 所示的理论模型，共 400 道，每道共 1024 个采样点，道间距是 10 m，采样率是2 ms 图 13b 为非均匀采样数据，其最大道间距为40 m，最小道间距为 10 m 利用本文方法进行重建重构结果如图13c 所示，信噪比8.6 dB 看图能够得出经过重建后的效果变好，同相轴变得更加连续，有效波能量损失少，重构后的地震数据与原始数据的差异小，验证了方法的正确性3实际数据分析图 14a 为新疆某地区实际采集地震数据，该野外地震数据的道间距为 25 m，采样间隔为 1 ms，660道接收，时间采样长度为 3 s 图 14b 为 60%随机欠采样条件下的地震数据 对缺失地震道数据采用本文提出的方法来进行重建处理 重建结果如图 14c所示，其信噪比为8.3 dB 图14d 为其重建误差剖面图 从图中我们能够得出，在 60%随机欠采样的条件下，本文采用的重建方法具有良好的效果，缺失地震道数据恢复得较好，重建后的数据同相轴基本与原始数据一致，可以达到后续处理和解释的要求4结论与认识本文主要在压缩感知理论的基础上，研究了基于Bregman 迭代的地震数据重建方法，形成如下认识:(1)满足压缩感知稀疏采样的不规则地震数据，能够突破所受到的传统 Nyquist 采样定理的局限，利用基于 Bregman 迭代的重建方法，能够较理想的重建出缺失的地震数据，并且该方法也具有一定的抗噪声能力(2)阈值模型是影响地震数据重建结果的重要因素 阈值模型不相同会获得不一样的重建效果，通过比较，本文提出的改进型指数阈值模型重建误差较少，精度较高，且计算效率也最快，在处理海量地震数据时具有优势致谢感谢审稿专家提出的修改意见和编辑部的大力支持!eferencesBrown L，Mosher C C，Li C B，et al 2017 Application of compressiveseismic imaging at Lookout Field，Alaska The Leading Edge，36(8):670-676，doi:101190/tle36080670 1Candes E J，omberg J，Tao T 2006 obust uncertainty principles:Exactsignalreconstructionfromhighlyincompletefrequencyinformation IEEE Transactions on Information Theory，52(2):489-509，doi:101109/TIT2005862083Hennenfent G，Herrmann F J2008Simply denoise:wavefieldreconstruction via jittered undersamplingGeophysics，73(3):V19-V28，doi:101190/12841038Herrmann F J 2010 andomized sampling and sparsity:Getting moreinformation from fewer samplesGeophysics，75(6):WB173-WB187，doi:101190/13506147Herrmann F J，Wang D L，Hennenfent G，et al 2008 Curvelet-basedseismic data processing:A multiscale and nonlinear approachGeophysics，73(1):A 1-A5，doi:101190/12799517Jin S D 2010 5D seismic data regularization by a damped least-normFourier inversion Geophysics，75(6):WB103-WB111，doi:101190/13505002Kaplan S T，Naghizadeh M，Sacchi M D 2010 Data reconstruction withshot-profile least-squares migration Geophysics，75(6):WB121-WB136，doi:101190/13478375Larner K，Gibson B，othman D 1981 Trace interpolation and thedesign of seismic surveys Geophysics，46(4):407-415Liu Y，Fomel S 2010 OC-seislet:seislet transform construction withdifferential offset continuationGeophysics，75(6):WB235-WB245，doi:101190/13479554Mosher C C，Kaplan S T，Janiszewski F D 2012b Non-uniform optimalsampling for seismic survey design/74th EAGE Conference andExhibition Incorporating EUOPEC 2012 European Association ofGeoscientists Engineers，doi:103997/2214-460920148781Mosher C C，Keskula E，Kaplan S T，et al 2012a Compressive seismicimaging/SEGTechnicalProgramExpandedAbstracts 2012Society of Exploration Geophysicists，1-5NaghizadehM，SacchiMD2007Multistepautoregressivereconstruction of seismic records Geophysics，72(6):V111-V118，doi:10 1190/12771685Naghizadeh M，Sacchi M D 2010Beyond alias hierarchical scalecurvelet interpolation of regularly and irregularly sampled seismicdata Geophysics，75(6):WB189-WB202，doi:10 1190/13509468Osher S，Burger M，Goldfarb D，et al 2005 An iterative regularizationmethodfortotalvariation-basedimagerestorationMultiscaleModeling and Simulation，4(2):460-489，doi:10 1137/040605412Porsani M J 1999 Seismic trace interpolation using half-step predictionfilters Geophysics，64(5):1461-1467，doi:101190/11444650onen J 1987 Wave-equation trace interpolation Geophysics，52(7):973-984，doi:101190/1 1442366Spitz S1991SeismictraceinterpolationintheF-XdomainGeophysics，56(6):785-794，doi:101190/11443096Wang C B 2018 Compressed sensing seismic data reconstruction withShearlet transformation Progress in Geophysics(in Chinese)，33(6):2441-2449 doi:106038/pg2018BB0499Wu S，Wu B Y，Geng Y 2009 Imaging in compressed domain usingdreamlets/Beijing 2009 International Geophysical Conference andExposition Beijing，China:Society of Exploration Geophysicists，212Xu S，Zhang Y，Pham D，et al 2005 Antileakage Fourier transform forseismic data regularization Geophysics，70(4):V87-V95，doi:101190/11993713Yang P L，Fomel S 2015 Seislet-based morphological component analysisusing scale-dependentexponentialshrinkageJournalofAppliedGeophysics，118:66-74，doi:101016/jjappgeo201504003附中文参考文献王常波 2018 基于 Shearlet 稀疏变换基的压缩感知重建技术 地球物理学进展，33(6):2441-2449 doi:10 6038/pg2018BB0499徐卫，张华，张落毅 2016 基于复值曲波变换的地震数据重建方法 物探与化探，40(4):750-756814