基于扩展伪多道匹配的保幅型多次波压制方法_李振春

２０２０年６月第５５卷　第３期　＊山东省青岛市市南区宁夏路３０８号青岛大学电子信息学院，２６６０７１。

Ｅｍａｉｌ：ｔｈｕｌｚｘ＠１６３．ｃｏｍ本文于２０１９年１０月７日收到，最终修改稿于２０２０年３月２３日收到。

本项研究受国家自然科学基金项目“基于卷积神经网络的多次波自适应相减方法”（４１８０４１１０）、国家重点研发计划项目“超深层弱信号增强、速度建模与保幅偏移技术研究”（２０１６ＹＦＣ０６０１１０５０１）和中石油重大科技项目“塔里木盆地深层复杂高陡构造与碳酸盐岩储层地震速度建模及成像关键技术研究”（ＺＤ２０１９－１８３－００３）联合资助。

·处理技术·文章编号：１０００－７２１０（２０２０）０３－０５３０－１１基于３Ｄ匹配滤波器和伪地震数据算法的多次波自适应相减方法李钟晓＊①　高好天①　陈鑫泽①　李永强②　李振春②（①青岛大学电子信息学院，山东青岛２６６０７１；②中国石油大学（华东）地球科学与技术学院，山东青岛２６６５８０）李钟晓，高好天，陈鑫泽，李永强，李振春．基于３Ｄ匹配滤波器和伪地震数据算法的多次波自适应相减方法．石油地球物理勘探，２０２０，５５（３）：５３０－５４０．摘要　多次波自适应相减是预测减去法压制多次波的关键步骤。

为进一步去除残余多次波，基于常规２Ｄ匹配滤波方法，文中引入３Ｄ匹配滤波器，同时利用多个预测多次波道集以匹配原始数据。

针对３Ｄ匹配滤波器可能造成的一次波损伤现象，利用相同的３Ｄ匹配滤波器同时拟合多个原始数据道集；同时，引入伪地震数据算法求解对一次波施加Ｈｕｂｅｒ范数最小化约束的优化问题，不需满足一次波与多次波正交的假设，能有效分离一次波与多次波。

另外，在整个迭代过程中，伪地震数据算法只需利用Ｃｈｏｌｅｓｋｙ分解算法进行一次矩阵分解，计算效率较高。

模型和实际数据的处理结果表明，与基于一次波能量最小化的３Ｄ匹配滤波器方法和基于伪地震数据算法的２Ｄ匹配滤波器方法相比，所提方法能更好地均衡一次波保护与多次波分离。

各向异性高分辨率Radon变换压制多次波范景文;李振春;宋翔宇;张凯;周隶华【期刊名称】《石油地球物理勘探》【年(卷),期】2016(051)004【摘要】针对VTI介质和海洋大炮检距地震数据反射波同相轴不能严格满足双曲时距关系,一次波和多次波能量在Radon域中聚焦性较差,多次波的压制效果不理想的问题,通过在双曲Radon变换积分路径中引入各向异性参数η,描述VTI介质和海洋大炮检距地震数据的时距曲线关系,在Radon域中有效聚焦一次波和多次波能量.同时,引入基于共轭导向梯度(CGG)算法的各向异性高分辨率Radon变换,基于L1范数对模型空间进行稀疏约束,提高Radon变换的分辨率.模型数据和海上大炮检距实际资料测试表明,该方法能够很好地压制多次波.【总页数】5页(P665-669)【作者】范景文;李振春;宋翔宇;张凯;周隶华【作者单位】中国石油大学(华东)地球科学与技术学院,山东青岛266580;中国石油大学(华东)地球科学与技术学院,山东青岛266580;中海油研究总院,北京100028;中国石油大学(华东)地球科学与技术学院,山东青岛266580;新疆油田分公司采油一厂地质研究所,新疆克拉玛依834000【正文语种】中文【中图分类】P631【相关文献】1.混合域高分辨率双曲Radon变换及其在多次波压制中的应用 [J], 巩向博;韩立国;王升超2.各向异性Radon变换及其在多次波压制中的应用 [J], 巩向博;韩立国;李洪建3.各向异性Radon变换在深水地震多次波压制中的应用 [J], 梁全有;韩立国;巩向博4.基于低频约束的高分辨率Radon变换多次波压制方法研究 [J], 刘仕友;马继涛;孙万元;应明雄5.基于SLA范数的高阶高分辨率λ-f域Radon变换多次波压制方法（英文） [J], 孙文之;李振春;曲英铭;李志娜因版权原因，仅展示原文概要，查看原文内容请购买。

东北石油大学学报第42卷第2期 2018年4月JOURNAL OF NORTHEAST PETROLEUM UNIVERSITY Vol.42 No. 2 Apr.2018DOI 10. 3969/j.issn. 2095 —4107. 2018. 02. 010基于频率拓展的表面多次波压制方法李星缘，王维红，王海娇，张文武(东北石油大学地球科学学院，黑龙江大庆163318 )摘要：为消除多次波预测中多余的子波效应，拓宽地震数据的频带，使预测多次波的走时、振幅和相位与实际的多次波相吻合，根据原始数据，衍生频率拓展道等多道地震数据；根据最小二乘原则，得到自适应滤波器和压制多次波后的有效波数据n结果表明：基于频率拓展的多次波压制方法，能够有效压制多次波，提高地震资料的信噪比；对水平层状模型测试结果频谱分析，自适应相减后得到的有效波同预测的多次波相比，频带拓宽，高频能量增强；通过对SM AART模型测试，基于频率拓展的多次波压制方法在噪声较强的情况下，也能得到很好的多次波压制效果n对于地质结构较复杂的情况，与拋物Radon变换滤波法相比，该方法多次波压制效果更好，并且不损害有效波能量。

关键词：表面多次波；SRME;预测；自适应相减；频率拓展；压制中图分类号：TE132. 1;P631 文献标识码：A文章编号=2095 - 4107(2018)02 - 0088 - 070引百简单地区油气勘探基本完毕，随着油气需求量的增大，复杂介质的地震资料处理方法被广泛关注。

尤 其是存在大量多次波的地震剖面，会产生构造假象、解释困难等问题，多次波的压制一直是地震数据处理，尤其是表面多次波丰富的海洋勘探亟待解决的问题之一。

常规的多次波压制方法主要包括滤波法和预测相减法。

滤波法利用多次波和一次波之间的运动学差 异压制多次波，预测相减法是以波动方程为基础，利用一次波和多次波之间的动力学差异压制多次 波[1^3]。

在简单的油气勘探区域，地质条件良好，一次波和多次波有较大时差，滤波方法可有效压制多次 波，且具有计算成本低的特点；与滤波法相比，预测相减法的优点是无需地下介质的先验信息，可更好地适 应构造复杂的地区[4]。

基于共散射点道集的多次波压制方法研究谢玉洪;李列;刘玉金;李振春;李林;袁全社【摘要】多次波压制是地震数据处理中的重要环节.一种常用的多次波压制方法是根据一次波和多次波在成像空间的曲率差异并借助Radon变换进行分离和压制.成像空间既可以是叠前时间域也可以是叠前深度域.考虑到对复杂介质的适应性以及叠前处理对计算效率的要求,提出基于共散射点(CSP)道集的高分辨率Radon变换多次波压制方法.该方法只需要一个简单的初始速度场,就可将常规共中心点道集(CMP)映射到覆盖次数更高、炮检距覆盖范围更广的共散射点道集,然后在该道集上应用高分辨率双曲Radon变换,可较好地分离一次波和多次波.相对于以水平层状介质为假设条件的CMP道集,CSP道集更适应复杂地质构造,且时距关系满足双曲规律.模型和实际资料测试结果表明,该方法可以较好地实现速度谱能量团的聚焦,有利于较复杂地质条件下的多次波压制.【期刊名称】《石油地球物理勘探》【年(卷),期】2015(050)002【总页数】6页(P232-237)【关键词】共散射点道集;多次波压制;高分辨率Radon变换【作者】谢玉洪;李列;刘玉金;李振春;李林;袁全社【作者单位】中海油湛江分公司,广东湛江524000;中海油湛江分公司,广东湛江524000;中国石油大学(华东),山东青岛266580;中国石油大学(华东),山东青岛266580;中海油湛江分公司,广东湛江524000;中海油湛江分公司,广东湛江524000【正文语种】中文【中图分类】P6311 引言多次波压制一直是勘探地球物理领域的热点问题，多次波属于相干噪声，其直接影响一次波成像的真实性以及地震资料解释的准确性。

由于一次波和多次波往往耦合在一起，因此多次波一直是地震资料处理中难以克服的问题。

目前常用的多次波压制方法根据压制空间主要可以分为两类。

一是基于数据空间多次波压制方法，主要包括预测反褶积法［1］、F-K 滤波法［2，3］、聚束滤波法［4］、表层多次波压制（SRME）［5］方法等。

基于反馈迭代模型的多次波压制方法综述
贺紫林;李振春;李志娜;徐夷鹏
【期刊名称】《物探与化探》
【年(卷),期】2022(46)2
【摘要】在油气勘探中,多次波通常被视为影响地震数据处理和解释结果准确性的相干噪声,因此,多次波压制方法研究一直是地球物理勘探领域中的重要研究方向。

鉴于自由表面多次波是地震勘探尤其是海上地震勘探资料中最为发育的多次波类型,因此,对自由表面多次波压制方法的研究是十分必要的。

目前国内外精度较高且应用广泛的自由表面多次波压制方法即为基于波动方程的反馈迭代法,该方法以地震原始数据为预测因子,不需要其他先验条件,可以在数据驱动下实现高精度的自由表面多次波压制。

在该方法理论基础上,近年来发展了一系列更为先进的方法技术。

本文首先简要介绍了反馈迭代法的理论基础——反馈迭代模型;然后阐述了基于反馈迭代模型发展的SRME、反数据域法、EPSI、CL-SRME四种方法压制自由表面多次波的机理;最后对比了四种方法各自的优缺点,展望了基于反馈迭代模型的多次波压制方法的研究前景。

【总页数】10页(P275-284)
【作者】贺紫林;李振春;李志娜;徐夷鹏
【作者单位】中国石油大学(华东)地球科学与技术学院
【正文语种】中文
【中图分类】P631
【相关文献】
1.反馈迭代法在自由表面多次波压制中的应用
2.基于相关迭代的非因果匹配滤波器多次波压制方法
3.基于反馈迭代和独立分量分析的多次波压制技术研究
4.基于迭代反演的层间多次波压制方法
5.反馈迭代法压制表面多次波效果分析
因版权原因，仅展示原文概要，查看原文内容请购买。

一种数据驱动的VSP多次波压制方法X王　鹏(同济大学海洋与地球科学学院,上海　200092) 摘　要:VSP资料用于提供地下结构的高分辨率信息,但其成像质量通常受到多次波的干扰,影响了地下真实反射界面的识别。

本文介绍一种数据驱动的VSP多次波压制方法,该方法基于反馈迭代原理,综合利用地面和井中地震数据,资料中的多次波进行预测并基于最小能量准则实现自适应相减;同时,利用VSP资料中直达波,可以预测出下行波场,进而提取上行波记录。

简单模型结果表明,该方法可以有效实现VSP资料多次波压制及上下行波场分离,且无需地下模型先验信息。

关键词:VSP数据;多次波;数据驱动 中图分类号:P631.8+1 文献标识码:A 文章编号:1006—7981(2012)05—0027—02 VSP资料中的多次波作为一种规则干扰,影响其成像质量,造成虚假同相轴继而降低了解释的可靠度。

因此,VSP资料的多次波压制是VSP资料数据处理中极为关键的一步。

传统的VSP资料多次波压制方法是利用上下行波之间的反褶积运算来实现的,该方法仅适用于零偏VSP数据。

Ikelle给出了基于逆散射级数的VSP多次波压制方法,该方法综合利用地面及井中地震测量数据,并可用于实现波场分离[1]。

Jiang提出了一种模型驱动的,基于双程波动方程的VSP资料多次波压制方法,该方法需要地下模型的先验信息,且只能压制特定界面的多次波[2]。

此后,Jiang又介绍了一种数据驱动的VSP多次波压制方法,通过地面地震记录和VSP数据进行褶积运算进行多次波预测并匹配相减,同时,利用VSP资料中的直达波预测下行波分量,但没有给出计算实例[3]。

Ma借助CFP概念,把SRME技术应用于VSP 资料的多次波预测中,该方法无需地下介质的先验信息,结合地面及井中地震记录,预测所有与自由表面相关多次波;此外,基于VSP资料中的直达波信息,预测下行波场进而实现波场分离[4]。

本文在详细分析SRME原理的基础上,通过地面地震和VSP资料接收方式的对比,认为将SRME 技术中多次波预测算子(即不含或部分衰减多次波的地震记录)的共接收点道集替换为VSP数据的共接收点道集,可以实现VSP资料多次波的预测。

一种基于高效迭代解法的频率域全波形反演解飞;黄建平;李振春【摘要】巨大的计算量是制约全波形反演(FWI)生产实用化的难题之一.为此,本文提出了一种高效的波场迭代解法,将其应用于频率域常密度声波方程FWI,并给出了详细的反演流程.通过建立用于波场迭代的目标函数,推导相应梯度、步长公式,新方法将反演中波场正传和残差波场反传过程转化为无约束优化问题,从理论上分析了新方法的计算效率显著高于常规FWI.在数值试验中,本文方法通过几次迭代便能获得高精度的正传、残差反传波场,收敛速度明显高于未经预处理的GMRES方法.进一步引入高效编码策略,新方法的计算时间约为常规编码FWI的1/8,与理论分析结果吻合(波场迭代次数为8,模型未知量个数约为7万),且波场迭代次数为6时,反演效果已与常规编码FWI相近.【期刊名称】《地球物理学报》【年(卷),期】2018(061)008【总页数】10页(P3346-3355)【关键词】全波形反演;计算量;迭代解法;编码【作者】解飞;黄建平;李振春【作者单位】中国石油大学(华东)地球科学与技术学院地球物理系,青岛 266580;海洋国家实验室海洋矿产资源评价与探测技术功能实验室,青岛 266071;中国石油大学(华东)地球科学与技术学院地球物理系,青岛 266580;海洋国家实验室海洋矿产资源评价与探测技术功能实验室,青岛 266071;中国石油大学(华东)地球科学与技术学院地球物理系,青岛 266580;海洋国家实验室海洋矿产资源评价与探测技术功能实验室,青岛 266071【正文语种】中文【中图分类】P6310 引言全波形反演(FWI)作为一种非线性反演方法，通过优化迭代的方式使模拟数据和观测数据之差最小化来估计地下物性参数(Lailly, 1983; Tarantola, 1984; Pratt et al., 1998).由于该方法基于双程波动方程，可以模拟所有波形(包括层间多次波)，因此相比于常规基于射线和单程波动方程的成像方法具有更高的成像精度.然而，目前FWI并未得到广泛应用，特别是在高密度地震采集中，巨大的计算量成为一个瓶颈.目前，编码技术是解决多炮高效成像和建模的有效方法之一.Romero等(2000)在叠前深度偏移中使用了相位编码方法，并测试了多种编码方式的效果.Krebs等(2009)将随机相位编码策略应用于时间域FWI，在每次迭代中随机生成新的编码矩阵以压制非相干正、反传残差波场产生的串扰噪音.Ben-Hadj-Ali等(2011)在频率域FWI中测试了随机相位编码的反演效果，编码后对随机噪音更为敏感并且敏感程度取决于组合炮的数量.Anagaw和Sacchi(2014)研究发现在每次迭代中改变编码函数和震源位置可以高效地压制串扰噪音.此外，使用不同的射线参数将炮集记录转化为平面波记录也可以显著降低计算量，适用于海洋观测系统且对随机噪音不敏感(Liu et al., 2006; Tao and Sen, 2013).Vigh和Starr(2008)将时间域平面波FWI应用于合成数据和实际3D海洋数据反演，而Tao和Sen(2013)表示在频率域可以通过相移合成平面波记录，并在实现2D频率域平面波FWI的基础上分析了射线参数的选取对反演结果的影响.与上面不同的是，Godwin和Sava(2010)提出了一种确定性编码方法称为截断奇异值编码(TSV)并用矩阵的形式表示了常规和联合震源成像的过程，对加窗后的近似单位矩阵进行奇异值分解以获得编码矩阵. 除编码策略外，一些学者从其他方面对减少FWI计算量进行了有益的研究.Wang 等(2013)将L-BFGS方法和不精确牛顿方法相结合更高效地估计Hessian算子的逆，从而加快收敛以减少迭代次数.曲英铭等(2016)在时间域将双变网格技术应用于FWI，在不同区域采用不同网格间距，并采取局部变时间采样来减少计算量.Van Leeuwen和Herrmann(2013)提出了一种分批次增加震源的反演方法，反演初期参与的震源较少，具有较高的迭代效率，随着反演的进行，与常规优化方法越来越相近，因而收敛速度加快.2D频率域FWI中，直接解法(如LU分解)由于具有求解多炮的优势因而得到广泛应用，其计算量主要集中在大型稀疏矩阵的分解上.Abubakar等(2009)将有限差分对比源反演(FDCSI)方法应用于频率域FWI用以减少矩阵分解的次数.该方法只需在每个频率开始时进行一次阻抗矩阵分解，其后的迭代中保持背景模型不变，矩阵分解的结果可以重复使用，因而具有很高的迭代效率.He等(2015)利用两步反演策略将FDCSI方法扩展到纵波速度、密度双参数反演.随后通过对多分量弹性波场进行线性变换，形成拟守恒弹性波动方程，将FDCSI方法应用于弹性波方程反演(He et al., 2016).由于FDCSI方法的第一步是计算每炮的对比源，当炮数过多时，其计算效率将受到挑战；其次，该反演方法形式比较特殊，与常规FWI方法相差较大，所以可移植性差，从而制约了其发展.借鉴FDCSI方法的思想，考虑到方法的计算量与可移植性问题，在常规频率域FWI理论的基础上，本文发展了一种基于波场迭代解法的高效FWI方法.该方法同样只需要在每个频率开始时进行一次阻抗矩阵的分解运算，其后迭代中波场的正传和残差波场反传都是基于共轭梯度法的矩阵回带运算.配合高效的编码策略，理论分析与实际测试都证明了新方法的计算效率要显著高于常规编码FWI.1 方法原理1.1 非线性最小二乘反演在频率域，二维各向同性常密度声波方程可以表示为(Marfurt,1984)(1)其中，r=(x,z)是笛卡尔坐标系下的地下空间位置；rs为震源空间位置；ω是角频率；表示拉普拉斯算子；f(ω)代表频率域的震源项；δ(r-rs)为狄拉克函数；u(r,rs,ω)表示频率域波场；k(r,ω)=ω/v(r)代表波数，v(r)是地下介质的纵波速度. 作为一个最小二乘反演问题，全波形反演正是通过优化的方法使模拟数据和观测数据的误差最小来对模型进行迭代更新(Lailly, 1983; Tarantola, 1984).其基于L2范数的目标函数定义如下：(2)其中，δd(rg,rs,ω)表示检波器位置处观测数据uobs(rg,rs,ω)和模拟数据usyn(rg,rs,ω)之差组成的列向量，rg指检波器所在位置；p(r)是模型参数(如纵波速度)；上标t和*分别表示转置和共轭.目标函数的梯度表达式为{Jtδd(rg,rs,ω)*}，(3)式中，Re表示取实部；矩阵J是Fréchet导数矩阵.当模型参数个数m较大时，直接计算Fréchet导数矩阵的计算量巨大.将Tarantola的残差波场反传理论引入到频率域波形反演，可以高效地得到目标函数对第i个模型参数的梯度(Pratt, 1999; 曹书红和陈景波, 2014)：}，(4)其中，A是与模型参数、频率、离散格式以及边界条件有关的阻抗矩阵；是由δd(rg,rs,ω)增广所得的残差列向量(非检波器位置为零).在获取梯度以后，可以使用最速下降法(Choi et al., 2005)或进一步利用共轭梯度法(Kamei et al., 2014)及截断牛顿法(王义和董良国, 2015)等构造优化方向q(r).模型参数的更新迭代公式如下：pn+1(r)=pn(r)+αnqn(r)，(5)上式中，n表示迭代次数；步长α可以通过线搜索的方式获得(Gauthier et al., 1986; Pica et al., 1990).1.2 迭代法生成波场本文将总波场u分解为背景波场ubac和散射场usct之和的形式(u=ubac+usct)，并且背景波场满足方程：(6)其中，kb表示背景介质的波数.从方程(1)中减去方程(6)可以得到散射场满足以下方程：(7)上式中，η(r,ω)=(r,ω)-k2(r,ω).为了简化表示，引入线性算子Lb，则散射场可以表示为usct(r,rs,ω)=Lb[η(r,ω)u(r,rs,ω)].由此可以定义用于波场迭代的L2型目标函数(为了方便，下文中不再刻意强调坐标位置r、角频率ω等)：C(u)=‖u-ubac-Lb[η u]‖2，(9)其中，线性算子Lb由背景模型的阻抗矩阵分解结果确定并且在同一个频率内保持不变.选定需要进行波场模拟的模型后，通过对(9)式最小化可以求得对应波场u，从而避免了对模型阻抗矩阵的重复分解计算(直接解法).本文采用共轭梯度法对(9)式进行迭代更新，其梯度公式如(10)式所示：}，(10)r(u)=u-ubac-Lb[η u]，(11)式中，表示Lb的伴随算子.通过梯度信息构造Polak-Ribière共轭梯度方向w，并令∂ C(u+β w)/∂ β=0得到步长β的计算公式如下：(12)其中，〈〉表示内积运算.由于步长β可以利用公式(12)显式表达，省去了对步长的估算，从而进一步提高了方法的效率.上面的推导过程中解决了正传波场的求解问题，而全波形反演中另一个重要部分就是残差波场的反传.用残差波场的反传结果δu替换u，得到下面的目标函数：D(δu)=‖δu-δubac-Lb[η δu]‖2，(13)上式中，的求解过程与u基本一致，在此不再赘述.值得注意的是，此时相当于震源项并且在模型更新中不断变化，所以每次更新需要重新计算δubac.图1给出的是基于波场迭代方法的FWI反演流程，其中抛物插值法计算步长α的正演流程与红色方框中一致.图1 单次迭代反演流程Fig.1 Flow chart of a single iterative inversion1.3 计算量分析上文中，完成了波场正传以及残差波场反传的计算，线性算子Lb和只与背景模型有关(背景模型保持不变)，因此，采用直接解法(LU分解)是一种理想的选择.在这过程中只需对背景模型进行一次阻抗矩阵的分解，其余计算主要是基于共轭梯度法的矩阵乘法与加法运算.假设离散的未知量个数为N，震源数量为NS，模型迭代次数为Niter，用于生成正传波场和残差反传波场的迭代次数均为niter，LU分解计算量的数量阶为O(2N1.5)，向前向后回带的计算量为O(N logN)(Davis and Duff, 1997).本文采用Polak-Ribière共轭梯度法对模型参数进行更新，并利用抛物插值法估算步长α(韩雨桐等, 2016).由于每次模型更新中至少额外进行两次正演模拟，所以新方法与传统CG方法的计算量分别为Tnew≈O(2N1.5)+NiterniterO(4NSNlogN)，(14)TCG≈Niter{3×O(2N1.5)+O(4NSNlogN)}，(15)由于直接法求解波动方程的回带阶段所需要的时间不到LU分解阶段CPU计算时间的5%(Pratt and Worthington, 1990)，所以本文将计算量主要集中于回带阶段可以获取较高的计算收益.另外，新方法只需几次迭代就可以合成出高精度的波场(niter较小)，再结合高效的编码策略和硬件并行加速方法进一步压缩NS，可以获得更高的反演效率.图2中绘制了理论计算量随模型参数个数N变化的关系曲线(Niter、niter和NS分别取20、8和1)，本文方法的计算量显著低于常规FWI，并且在N较大时，优势更为明显.图2 理论计算量对比Fig.2 Comparison of the theoretical computation of the new method and the conventional FWI2 数值试验为了验证本文提出的波场迭代解法的有效性，采用了简单洼陷模型进行波场正传和残差波场反传测试，并与GMRES方法进行对比以突出高效性.然后将新方法与截断奇异值编码策略(TSV)(Godwin and Sava, 2010)进行了结合，利用Marmousi模型测试了常规编码FWI和niter取不同值时新方法的反演效果.2.1 迭代法生成波场真实洼陷模型及其初始模型(利用Madagascar的sfsmooth函数平滑所得，平滑窗口为20)如图3所示，参数如下：水平网格数nx=150，垂直网格数nz=100，水平和垂直网格间距均为10 m.检波器设置为地表全接收，震源位于x=750 m、z=20 m处，震源子波采用主频为10 Hz的雷克子波.利用图3a中的真实模型进行LU分解计算得到图4a中的正传波场(8 Hz实部数据，下同)作为标准波场，同样，图4b是由图3b中平滑初始模型进行LU分解计算所得.采用本文的波场迭代方法，以图4b作为背景波场进行迭代计算，10次迭代后的结果如图4c所示.图4d是真实洼陷模型使用GMRES方法求解的波场结果(以图4b作为初始值)，其迭代时间与图4c一致(约0.7 s).图4中的第3行(e，f，g)是第2行(b，c，d)分别与图4a作差的结果，从中可以看出GMRES方法迭代以后的波场残差(图4g)相比于图4e有所减小，而本文方法迭代后的残差(图4f)趋近于零，与标准波场基本一致.图5a是检波器接收到的残差波场在真实洼陷模型中的反传结果(8 Hz实部数据，下同)，而图5b是在初始模型中的反传结果.此外，图5中图形的排列方式与图4完全一致，可以看出经过迭代以后，GMRES方法的波场残差(图5g)明显大于本文方法(图5f).由此说明，在提供一个背景波场的前提下，借助背景模型的阻抗矩阵分解结果，新方法可以高效地获得已知速度场的正传和残差反传波场，从而为下一步计算梯度打下基础.图3 洼陷模型(a)以及初始平滑模型(b)Fig.3 The sag model (a) and initial smoothing model (b)由于本文提出的波场迭代方法需要配合震源编码技术来压缩NS，才能充分发挥新方法的计算效率.而当NS很小时，常用的迭代解法(如：GMRES)就成为了一种可以考虑的选择.为了比较新方法和GMRES方法的迭代效果，在图4和图5中对比了两种方法在相同时间下的迭代结果.图6 (a—b)分别是正传和反传两种情况下，两种迭代方法的归一化波场误差随迭代时间变化曲线，本文方法的收敛速度明显高于GMRES方法(新方法不包含背景模型矩阵分解时间，因为一个频率反演中只需分解一次).主要原因是新方法是一种基于CG方法的无约束优化方法，具有明确的优化方向和步长值，而GMRES方法需要进行合适的预处理，而这通常是难以满足的(Virieux and Operto, 2009).图6 (c—d)分别是正、反传时新方法的波场误差随迭代次数变化曲线，大约进行6次迭代后残差基本收敛.为了测试新方法对初始模型的适应性，分别采用图7a中平滑程度更高的初始模型(平滑窗口为30)和图8a中线性梯度初始模型进行波场合成.图7(b—d)分别是初始模型的模拟波场以及2次和10次迭代后的波场结果，而图7 (e—g)则分别是图7 (b—d)与图4a中标准波场的差值.图8的排列方式与图7一致，并且可以看出，新方法对初始模型具有较好的适应性，随着迭代的进行，残差稳定减小.值得注意的是，本次测试中是以初始模型的波场作为迭代的初始值，与真实波场差异较大，实际反演中相邻迭代间模型差异较小，并以上一次迭代的波场作为初始值，因此需要的迭代次数通常更少.图4 8 Hz正传波场实部(a) 真实模型正演波场； (b) 初始模型正演波场； (c) 新方法迭代生成的正演波场； (d) GMRES方法迭代生成的正演波场，与(c)中计算时间一致； (e—g)分别是(b—d)与(a)之差.Fig.4 The real part of forward wavefields of 8 Hz data(a) The forward wavefield of true model； (b) The forward wavefield of initial model； (c) The forward wavefield obtained by the new method； (d) The forward wavefield obtained by the GMRES method with same computing time as (c)； (e—g) The difference between (b—d) and (a), respectively.图5 8 Hz残差数据反传波场实部(a) 真实模型反传波场； (b) 初始模型反传波场；(c) 新方法迭代生成的反传波场； (d) GMRES方法迭代生成的反传波场，与(c)中计算时间一致； (e—g)分别是(b—d)与(a)之差.Fig.5 The real part of back-propagated wavefields of 8Hz residual data(a) The back-propagated wavefield of true model； (b) The back-propagated wavefield of initial model； (c) The back-propagated wavefield obtained by the new method；(d) The back-propagated wavefield obtained by the GMRES method with same computing time as (c)； (e—g)The difference between (b—d) and (a), respectively.图6 归一化波场误差收敛曲线(a) 正传波场误差随时间变化曲线； (b) 反传波场误差随时间变化曲线； (c) 正传波场误差随迭代次数变化曲线； (d) 反传波场误差随迭代次数变化曲线.Fig.6 The history of normalized wavefields errorThe variation of forward wavefields error with (a) time and (c) iteration； The variation of back-propagated wavefields error with (b) time and (d)iteration.图7 平滑初始模型及波场模拟(a) 平滑初始模型； (b) 初始模型正演波场； (c) 2次迭代后的波场； (d) 10次迭代后的波场； (e—g) 分别是(b—d)与图4a之差.Fig.7 The smooth initial model and wavefields modeling(a) The smooth initial model； (b) The forward wavefield of initial model； (c) The wavefield after 2 iterations； (d) The wavefield after 10 iterations； (e—g) The difference between (b—d) and Fig.4a, respectively.图8 梯度初始模型及波场模拟(a) 梯度初始模型； (b) 初始模型正演波场； (c) 2次迭代后的波场； (d) 10次迭代后的波场； (e—g) 分别是(b—d)与图4a之差.Fig.8 The gradient initial model and wavefields modeling(a) The gradient initial model； (b) The forward wavefield of initial model； (c) The wavefield after 2 iterations； (d) The wavefield after 10 iterations； (e—g) The difference between (b—d) and figure 4a, respectively.2.2 合成数据反演为了验证本文方法的反演效果，采用Marmousi模型(Anagaw and Sacchi, 2014)进行合成数据测试.所使用的真实和初始模型如图9所示(初始模型由真实模型平滑所得)，网格大小为126×384，横、纵向空间采样间隔均为15 m.震源和检波器位于地表，第一炮位于15 m 处，炮间距为30 m，共191炮，检波器设置为地表全接收，震源子波为主频15 Hz的雷克子波.图9 Marmousi真实模型(a)以及平滑初始模型(b)Fig.9 The Marmousi model (a) and initial smoothing model (b)图10 串扰噪音矩阵Fig.10 Cross-talk matrix obtained by multiplying encoding matrix with its transposition反演中所使用的频率为4、6、8、10、12、15、18、21、24、28、32和36 Hz共12个频率，每个频率迭代20次，低频的反演结果作为相邻高频的输入.本文使用截断奇异值编码(TSV)方法，将191炮编码为15炮，并采取并行反演策略，图10是相应的串扰噪音矩阵.常规编码FWI反演结果如图11a所示，niter分别取10、8和6时，新方法的反演结果对应于图11(b—d).由于串扰噪音的影响，反演结果的浅部低速层受影响比较明显，但总体上说，图11反演结果较好并且比较接近.图12给出了图11中4个反演结果的归一化收敛曲线，常规编码FWI和新方法niter等于10时的收敛效果相差不大，并且常规编码FWI的模型误差最小，而新方法随着niter增大收敛性有所改善.图13记录了反演时间随迭代次数变化曲线，常规编码FWI的计算用时要显著高于本文方法，并且niter越小新方法计算时间越少.本次测试中模型四周设有20个网格的PML边界，实际计算用时与图2中的理论用时比较吻合，常规编码FWI计算时间约为niter等于8时新方法的8倍(N≈7×104).3 结论本文提出了一种高效的迭代解法生成频率域波场方法，将其应用于声波常密度方程FWI，并给出了详细的反演流程.在实现算法的基础上，通过洼陷模型和Marmousi模型测试，得到如下几点认识：(1) 新方法将频率域波场计算转化为一种无约束优化问题，通过几次迭代便能生成高精度的正传和残差反传波场，其误差收敛速度显著高于未经预处理的GMRES 方法.(2) 基于新的波场迭代法的FWI，每个频率只需进行一次阻抗矩阵的分解，其计算量与波场迭代次数niter成正比.配合高效的震源编码技术，新方法的计算时间相比常规编码FWI有明显降低(niter=8、N≈7×104时，计算时间约为常规编码FWI 的1/8).图11 最终反演结果(a) 常规FWI； (b) 新方法niter=10； (c) 新方法niter=8；(d) 新方法niter=6.Fig.11 The final inversion results(a) Conventional FWI；New method with (b) niter=10，(c) niter=8， (d) niter=6.图12 归一化模型误差收敛曲线Fig.12 The normalized model fit图13 计算时间与迭代次数关系曲线Fig.13 The relationship between the inversion time and the number of iterations(3) niter越大，新方法的反演效果与常规编码FWI越接近，但是当生成的波场精度达到一定条件后，niter的增大对反演结果的改善并不明显.因此，需要在计算量与反演效果之间做好权衡.虽然本文提出的基于波场迭代法FWI可以获得与常规FWI基本一致的反演效果，并且在计算效率上有显著优势，但是该方法还不是十分完善，也存在一些需要改进的地方，比如：(1)当NS较大时，新方法的计算效率将受到严重影响，因此必须配合高效的震源编码技术，这样不可避免地会引入一些串扰噪音；(2)新的波场迭代法本质上是一种基于直接解法的迭代方法，当模型维度过大时，同样存在存储量的问题.ReferencesAbubakar A, Hu W Y, Habashy T M, et al. 2009. Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data. Geophysics, 74(6): WCC47-WCC58.Anagaw A Y, Sacchi M D. 2014. Comparison of multifrequency selection strategies for simultaneous-source full-waveform inversion. Geophysics, 79(5): R165-R181.Ben-Hadj-Ali H, Operto S, Virieux J. 2011. An efficient frequency domain full-waveform inversion method using simultaneous encoded sources.Geophysics, 76(4): R109-R124.Cao S H, Chen J B. 2014. Studies on complex frequencies in frequency domain full waveform inversion. Chinese Journal of Geophysics (in Chinese), 57(7): 2302-2313, doi: 10.6038/cjg20140724.Choi Y, Shin C, Min D J, et al. 2005. Efficient calculation of the steepest descent direction for source-independent seismic waveform inversion: An amplitude approach. Journal of Computational Physics, 208(2): 455-468. Davis T A, Duff I S. 1997. An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM Journal on Matrix Analysis and Applications, 18(1): 140-158.Gauthier O, Virieux J, Tarantola A. 1986. Two-dimensional nonlinear inversion of seismic waveforms: Numerical results. Geophysics, 51(7): 1387-1403.Godwin J, Sava P. 2010. Simultaneous source imaging by amplitude encoding, Technical Report CWP-645, Center for Wave Phenomena, Colorado School of Mines.Han Y T, Huang J P, Qu Y M, et al. 2016. Encoded full waveform inversion method by an optimum inversion step length for aliasing data. Progress in Geophysics (in Chinese), 31(1): 198-204, doi: 10.6038/pg20160123.He Q L, Han B, Chen Y, et al. 2015. Application of the finite-difference contrast source inversion method to multiparameter reconstruction using seismic full-waveform data. Journal of Applied Geophysics, 124: 4-16.He Q L, Chen Y, Han B, et al. 2016. Elastic frequency-domain finite-difference contrast source inversion method. Inverse Problems, 32(3):035009.Kamei R, Pratt R G, Tsuji T. 2014. Misfit functionals in laplace-fourier domain waveform inversion, with application to wide-angle ocean bottom seismograph data. Geophysical Prospecting, 62(5): 1054-1074.Krebs J R, Anderson J E, Hinkley D, et al. 2009. Fast full-wavefield seismic inversion using encoded sources. Geophysics, 74(6): WCC177-WCC188. Lailly P. 1983. The seismic inverse problem as a sequence of before stack migration. ∥SIAM Conference on Inverse Scattering: Theory and Applications. Philadelphia: Society for Industrial and Applied Mathematics, 206-220.Liu F Q, Hanson D W, Whitmore N D, et al. 2006. Toward a unified analysis for source plane-wave migration. Geophysics, 71(4): S129-S139.Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equation. Geophysics, 49(5): 533-549.Pica A, Diet J P, Tarantola A. 1990. Nonlinear inversion of seismic reflection data in a laterally invariant medium. Geophysics, 55(3): 284-292.Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography. Geophysical Prospecting, 38(3): 287-310.Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362.Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model. Geophysics, 64(3): 888-901.Qu Y M, Li Z C, Huang J P, et al. 2016. Full waveform inversion based on multi-scale dual-variable grid in time domain. Geophysical Prospecting for Petroleum (in Chinese), 55(2): 241-250.Romero L A, Ghiglia D C, et al. 2000. Phase encoding of shot records in prestack migration. Geophysics, 65(2): 426-436.Tao Y, Sen M K. 2013. Frequency-domain full waveform inversion with plane-wave data. Geophysics, 78(1): R13-R23.Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266.VanLeeuwen T, Herrmann F J. 2013. Fast waveform inversion without source-encoding. Geophysical Prospecting, 61(S1): 10-19.Vigh D, Starr E. 2008. 3D prestack plane-wave, full-waveform inversion. Geophysics, 73(5): VE135-VE144.Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26.Wang Y, Dong L G, Liu Y Z. 2013. Improved hybrid iterative optimization method for seismic full waveform inversion. Applied Geophysics, 10(3): 265-277.Wang Y, Dong L G. 2015. Multi-parameter full waveform inversion for acoustic VTI media using the truncated Newton method. Chinese Journal Geophysics (in Chinese), 58(8): 2873-2885, doi: 10.6038/cjg20150821.附中文参考文献曹书红, 陈景波. 2014. 频率域全波形反演中关于复频率的研究. 地球物理学报,57(7): 2302-2313, doi: 10.6038/cjg20140724.韩雨桐, 黄建平, 曲英铭等. 2016. 基于优化步长的混叠数据编码全波形反演方法. 地球物理学进展, 31(1): 198-204, doi: 10.6038/pg20160123.曲英铭, 李振春, 黄建平等. 2016. 基于多尺度双变网格的时间域全波形反演. 石油物探, 55(2): 241-250.王义, 董良国. 2015. 基于截断牛顿法的VTI介质声波多参数全波形反演. 地球物理学报, 58(8): 2873-2885,doi: 10.6038/cjg20150821.【相关文献】Abubakar A, Hu W Y, Habashy T M, et al. 2009. Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data. Geophysics, 74(6): WCC47-WCC58.Anagaw A Y, Sacchi M D. 2014. Comparison of multifrequency selection strategies for simultaneous-source full-waveform inversion. Geophysics, 79(5): R165-R181.Ben-Hadj-Ali H, Operto S, Virieux J. 2011. An efficient frequency domain full-waveform inversion method using simultaneous encoded sources. Geophysics, 76(4): R109-R124. Cao S H, Chen J B. 2014. Studies on complex frequencies in frequency domain full waveform inversion. Chinese Journal of Geophysics (in Chinese), 57(7): 2302-2313, doi: 10.6038/cjg20140724.Choi Y, Shin C, Min D J, et al. 2005. Efficient calculation of the steepest descent direction for source-independent seismic waveform inversion: An amplitude approach. Journal of Computational Physics, 208(2): 455-468.Davis T A, Duff I S. 1997. An unsymmetric-pattern multifrontal method for sparse LU factorization. SIAM Journal on Matrix Analysis and Applications, 18(1): 140-158. Gauthier O, Virieux J, Tarantola A. 1986. Two-dimensional nonlinear inversion of seismic waveforms: Numerical results. Geophysics, 51(7): 1387-1403.Godwin J, Sava P. 2010. Simultaneous source imaging by amplitude encoding, Technical Report CWP-645, Center for Wave Phenomena, Colorado School of Mines.Han Y T, Huang J P, Qu Y M, et al. 2016. Encoded full waveform inversion method by an optimum inversion step length for aliasing data. Progress in Geophysics (in Chinese), 31(1): 198-204, doi: 10.6038/pg20160123.He Q L, Han B, Chen Y, et al. 2015. Application of the finite-difference contrast source inversion method to multiparameter reconstruction using seismic full-waveform data. Journal of Applied Geophysics, 124: 4-16.He Q L, Chen Y, Han B, et al. 2016. Elastic frequency-domain finite-difference contrast source inversion method. Inverse Problems, 32(3): 035009.Kamei R, Pratt R G, Tsuji T. 2014. Misfit functionals in laplace-fourier domain waveform inversion, with application to wide-angle ocean bottom seismograph data. Geophysical Prospecting, 62(5): 1054-1074.Krebs J R, Anderson J E, Hinkley D, et al. 2009. Fast full-wavefield seismic inversion using encoded sources. Geophysics, 74(6): WCC177-WCC188.Lailly P. 1983. The seismic inverse problem as a sequence of before stack migration.∥SIAM Conference on Inverse Scattering: Theory and Applications. Philadelphia: Society for Industrial and Applied Mathematics, 206-220.Liu F Q, Hanson D W, Whitmore N D, et al. 2006. Toward a unified analysis for source plane-wave migration. Geophysics, 71(4): S129-S139.Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equation. Geophysics, 49(5): 533-549.Pica A, Diet J P, Tarantola A. 1990. Nonlinear inversion of seismic reflection data in a laterally invariant medium. Geophysics, 55(3): 284-292.Pratt R G, Worthington M H. 1990. Inverse theory applied to multi-source cross-hole tomography. Geophysical Prospecting, 38(3): 287-310.Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newtonmethods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362.Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model. Geophysics, 64(3): 888-901.Qu Y M, Li Z C, Huang J P, et al. 2016. Full waveform inversion based on multi-scale dual-variable grid in time domain. Geophysical Prospecting for Petroleum (in Chinese), 55(2): 241-250.Romero L A, Ghiglia D C, et al. 2000. Phase encoding of shot records in prestack migration. Geophysics, 65(2): 426-436.Tao Y, Sen M K. 2013. Frequency-domain full waveform inversion with plane-wave data. Geophysics, 78(1): R13-R23.Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266.VanLeeuwen T, Herrmann F J. 2013. Fast waveform inversion without source-encoding. Geophysical Prospecting, 61(S1): 10-19.Vigh D, Starr E. 2008. 3D prestack plane-wave, full-waveform inversion. Geophysics, 73(5): VE135-VE144.Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26.Wang Y, Dong L G, Liu Y Z. 2013. Improved hybrid iterative optimization method for seismic full waveform inversion. Applied Geophysics, 10(3): 265-277.Wang Y, Dong L G. 2015. Multi-parameter full waveform inversion for acoustic VTI media using the truncated Newton method. Chinese Journal Geophysics (in Chinese), 58(8): 2873-2885, doi: 10.6038/cjg20150821.曹书红, 陈景波. 2014. 频率域全波形反演中关于复频率的研究. 地球物理学报, 57(7): 2302-2313, doi: 10.6038/cjg20140724.韩雨桐, 黄建平, 曲英铭等. 2016. 基于优化步长的混叠数据编码全波形反演方法. 地球物理学进展, 31(1): 198-204, doi: 10.6038/pg20160123.曲英铭, 李振春, 黄建平等. 2016. 基于多尺度双变网格的时间域全波形反演. 石油物探, 55(2): 241-250.王义, 董良国. 2015. 基于截断牛顿法的VTI介质声波多参数全波形反演. 地球物理学报, 58(8): 2873-2885,doi: 10.6038/cjg20150821.。

扩展多道匹配法压制海上多次波的开题报告一、研究背景海洋勘探一直是石油勘探的重要领域之一，其中含有大量的复杂地质构造，如海底隆起、漫扇、坳陷等，容易产生多次波，影响勘探的准确性和精度。

因此，如何进行多次波的压制，成为了海洋地震勘探研究的重要领域。

目前海上多次波压制的常用方法是多道匹配法。

多道匹配法是一种基于记录道之间相互关系的波形压制方法，它可以利用记录道之间的相关性对多次波进行抑制，从而提高地震数据的信噪比和分辨率。

但是，传统的多道匹配法往往不能很好地处理复杂的海底地形和结构，所以需要进一步研究和改进。

二、研究目的本研究旨在利用多道匹配法来压制海上多次波，特别是针对复杂地质构造和海底地形情况下的多次波进行研究。

具体研究目的如下：1. 探究复杂海底地形和结构下多次波形成的原因和规律。

2. 研究常规多道匹配法在海上多次波压制中存在的问题。

3. 提出针对海上多次波的改进多道匹配法，并对改进方法进行实验验证和性能评估。

三、研究方法本研究将采用以下方法和步骤：1. 对已有的多次波压制方法进行调研和比较，分析其优缺点和适用范围。

2. 分析复杂海底地形和结构下多次波形成的机理和规律。

3. 提出改进的多道匹配法并进行实验验证和性能评估。

4. 对结果进行分析和总结，并提出下一步的研究方向。

四、研究意义本研究的意义在于：1. 提高海洋地震勘探数据的质量和可靠性，对于石油勘探具有重要的指导意义。

2. 探究多次波形成的规律和机理，可以对海洋地质构造的研究提供一定的帮助。

3. 对多道匹配法进行改进，可以为其他领域的波形压制提供一些启示和借鉴。