高阶交错网格有限差分法纵横波波场分离数值模拟

高阶旋转交错网格有限差分方法模拟TTI介质中横波分裂李敏;刘洋【期刊名称】《物探与化探》【年(卷),期】2012(36)6【摘要】笔者给出了一种能够模拟弹性波在任意各向异性介质中传播的二维三分量高阶有限差分算法.相对于常规交错网格有限差分方法,旋转交错网格有限差分方法在介质具有强差异性时能更精确地模拟地震波的传播,避免常规交错网格中因对弹性系数进行插值而带来的误差.采用高阶旋转交错网格有限差分方法模拟并分析了零偏移距横波分裂现象随裂缝介质方位角和倾角变化的响应特征.结果表明:结合完全匹配层(PML)吸收边界条件的高阶旋转交错网格有限差分方法能获得高精度的地震波场模拟数据,并且在边界具有良好的吸收效果；横波分裂现象主要受裂缝走向与波的极化方向之间的夹角影响,受裂缝倾角影响较小,且快慢横波的能量也跟裂缝走向与波极化方向间的夹角有关.具有倾斜对称轴的横向各向同性(TTI)介质倾角的变化可能会导致记录中波到达时的变化,影响快慢横波的时差.利用横波分裂的能量分布和方位各向异性特征,可以帮助检测裂缝的方位角和倾角.横波在多层TTI介质中传播时会发生多次分裂的现象.%A high-order rotated staggered grid scheme ( RSG) has been implemented to simulate the shear-wave splitting in tilted transversely isotropic (TTI) media. The high-order RSG can simulate wave propagation in media that contain high-contrast discontinuities like cracks more precisely than the standard staggered grid scheme (SSG) by avoiding the unstableness of the staggered grid scheme (SSG). The authors conducted a study of zero-offset S-wave splitting withthe high-order RSG. The S-wave splitting study was mainly focused on fractured media which, on the scale of seismic wavelength, could be regarded as transversely isotropic (TI) media. The results of numerical modeling show that the high-order RSG scheme can be used to simulate waves' propagation in general anisotropic media. The perfect matched layer (PML) absorbing boundary condition combined with the high order RSG scheme can well attenuate reflections from the artificial boundary. The S-wave splitting is mainly affected by the angle between polarization direction of incoming wave and strike of the TTI media, and the energy of fast and slow shear waves is also associated with this angle. The dipping angle of ITI media may affect time lag between the fast and slow waves, which may result in variation of arrival time of waves from the same interface. Thus, the analysis of energy distribution of the fast and slow waves and the variation of arrival time may help detect the strike and dipping angle of the fracture. Besides, when propagating in the media that contain more than one layer of TTI media, the S-wave splitting will occur more than once.【总页数】7页(P934-940)【作者】李敏;刘洋【作者单位】中国石油大学油气资源与探测国家重点实验室,北京102249;中国石油大学CNPC物探重点实验室,北京102249;中国石油大学油气资源与探测国家重点实验室,北京102249;中国石油大学CNPC物探重点实验室,北京102249【正文语种】中文【中图分类】P631.4【相关文献】1.黏弹TTI介质旋转交错网格微地震波场模拟 [J], 姚振岸;孙成禹;谢俊法;唐杰2.黏弹TTI介质中旋转交错网格高阶有限差分数值模拟 [J], 严红勇;刘洋3.旋转交错网格在横波分裂和再分裂模拟与分析中的应用 [J], 张建利;刘志斌;周超;田小波;李维新;王赟4.基于卷积完全匹配层的旋转交错网格高阶差分法模拟弹性波传播 [J], 冯德山;王向宇5.基于CPML-RML组合边界条件粘弹TTI介质旋转交错网格有限差分正演模拟[J], 张奎涛;顾汉明;刘少勇;刘春成;陈宝书;张立;肖逸飞因版权原因，仅展示原文概要，查看原文内容请购买。

基于交错网格有限差分弹性波正演模拟及波场特征分析【摘要】为研究和认识多种储层中弹性波的波场特征，以利于多波地震资料解释，高精度数值模拟是有效的方法之一。

本文在弹性波方程基础上，采用高阶交错网格有限差分技术模拟地震波在各向同性介质和各向异性介质中的传播，可得到不同类型介质的弹性波场。

同时，文中也分析了各向异性系数对多波波场特征的影响。

通过对高精度数值模拟得到的波场快照对比研究表明，该方法可有效获得高精度弹性波正演结果，为研究各种复杂介质中弹性波的波场特征和传播规律奠定了基础。

【关键词】多波多分量波场特征各向异性弹性波正演1 引言随着油气田勘探技术的不断发展[1][2]，人们对地震资料的认识也不断加深，纵波地震资料在含油气的显示上存在一些不确定性，单一纵波资料解释的多解性问题尤为突出。

在地震勘探领域中，过去一直把各向同性弹性体理论作为研究地下介质的前提，但是在实际地层中普遍存在各向异性，地下介质的各向异性（如周期薄互层引起的各向异性、以及裂隙引起的各向异性）产生的弹性波场与各向同性介质产生的弹性波场存在着不可忽略的差异。

由此，多波地震勘探作为油储地球物理的主要方法之一应运而生。

在多波资料解释过程中，要求搞清楚储层的岩性与多波的波场特征之间的关系，因此，多波波场数值模拟技术显得非常重要。

高精度数值模拟技术是联系地震、地质、测井以及油藏工程的纽带，其作用主要体现在提高人们对各种复杂介质中地震波传播规律的认知，并可为新技术、新方法提供试验数据，以满足方法技术研究的需要，同时也可以检验解释结果的正确性。

弹性波波动方程高精度数值模拟可以得到全波场信息，包含了地震波的动力学和运动学特点，为准确描述地震波场特征和波的传播规律奠定基础，本文在弹性波方程基础上，采用高阶交错网格有限差分技术模拟地震波在各向同性介质和各向异性介质中的传播，比较地震波在各向同性介质和各向异性介质中的波场响应异同，并分析了各向异性系数对多波波场特征的影响，这对研究各种复杂介质中弹性波的波场特征和传播规律有着重要的意义。

地震波场的高阶交错网格有限差分模拟霍凤斌;李振鹏;徐发;张涛【摘要】This paper analyzes the stability and convergence of the seismic wave ifeld by using the high-order staggered-grid limited differential method of joining the absorbing boundary condition and attenuating zone to simulate the elastic wave equation. The results of the simulation of both isotropic-and anisotropic-medium models show that the grid frequency dispersion of the high-order differential wave equation simulation is smaller and more accurate. Therefore, this method should improve the efifciency of seismic prospecting and of the associated data interpretation.%应用高阶交错网格有限差分算法，并加入吸收边界条件和衰减带，对弹性波方程进行模拟，分析了其稳定性和收敛性。

通过对各向同性和各向异性介质模型的模拟表明，高阶差分波动方程模拟的网格频散较小、精度较高、效果较好，可为地震勘探及其资料解译提供技术手段。

【期刊名称】《上海国土资源》【年(卷),期】2014(000)001【总页数】4页(P97-100)【关键词】地震波场;波动方程;有限差分;边界条件;交错网格【作者】霍凤斌;李振鹏;徐发;张涛【作者单位】中海石油中国有限公司上海分公司，上海200030;中海石油中国有限公司上海分公司，上海200030;中海石油中国有限公司上海分公司，上海200030;中海石油中国有限公司上海分公司，上海200030【正文语种】中文【中图分类】P315.01随着地震波动理论在天然地震和油气地震中的应用，以及计算机技术的飞速发展，在现代地震数值模拟领域逐渐形成了有限差分法、有限元法、虚谱法和积分方程法等求解波动方程的方法。

地震波动方程方向行波波场分离正演数值模拟与逆时成像陈可洋;陈树民;李来林;吴清岭;范兴才;刘振宽【摘要】为了进一步提高对地震波传播规律的认识,将波印廷矢量引入到二维地震波动方程方向行波波场分离正演数值模拟中.根据地震波波印廷矢量的波场数值特征,实现了对全地震波场中左行波、右行波、上行波和下行波波场的自动识别与分离.以均匀介质模型、倾斜界面模型以及Marmousi模型为例,开展了相应的数值模拟实验与逆时成像处理.计算结果表明,该方法准确有效,能够实现任意时刻波场快照中方向行波的波场分离,并合成分别由左行波、右行波、上行波和下行波形成的波场快照与数值模拟记录.该方法简单易行,计算量较小,对实际地震资料中方向行波波场的识别、分离、成像及验证具有一定的应用价值.【期刊名称】《岩性油气藏》【年(卷),期】2014(026)004【总页数】7页(P130-136)【关键词】地震波动方程;正演模拟;单程波;波场分离;波印廷矢量;逆时成像【作者】陈可洋;陈树民;李来林;吴清岭;范兴才;刘振宽【作者单位】中国石油大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712;中国石油大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712;中国石油大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712;中国石油大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712;中国石油大庆油田有限责任公司勘探开发研究院,黑龙江大庆163712;中国石油大庆油田有限责任公司勘探事业部,黑龙江大庆163453【正文语种】中文【中图分类】P631.40 引言地震波正演数值模拟技术一直是国际地球物理学界的热点内容之一。

随着地震波动理论和计算机技术的不断发展，自20世纪60年代以来该项技术便得到了飞速发展，其中采用波动方程的地震波数值模拟技术能更好地保持地震波的几何学、运动学和动力学等特征，因此可达到精确模拟地震波传播规律的目的［1］。

二阶弹性波动方程高精度交错网格波场分离数值模拟陈可洋;杨微;刘洪林;吴清岭【摘要】给出了一种等价的二阶弹性波动方程,以解决弹性波场中完全弹性波动方程不能完全分离耦合的纵、横波波场问题.应用高阶交错网格有限差分法求解该波动方程,并使用通量校正技术(FCT)进一步压制频散,采用均匀介质模型和层状介质模型进行波场分离数值试验,精确得到了混合波场、完全分离的纯纵波及纯横波波场.数值结果分析表明,本文方法在均匀介质情况下准确可靠,在分离后的纯纵、横波波场中可观察到较为丰富的能量转换信息,这对认识复杂弹性波的传播规律及弹性波理论具有重要意义.【期刊名称】《物探与化探》【年(卷),期】2009(033)006【总页数】4页(P700-703)【关键词】地震波场分离;高阶交错网格;等价二阶弹性波动方程;数值模拟;通量校正技术【作者】陈可洋;杨微;刘洪林;吴清岭【作者单位】大庆石油学院,地球科学学院,黑龙江,大庆,163318;大庆石油学院,地球科学学院,黑龙江,大庆,163318;大庆石油学院,地球科学学院,黑龙江,大庆,163318;大庆油田有限责任公司,勘探开发研究院,黑龙江,大庆,163712【正文语种】中文【中图分类】P631.4地震波场数值模拟技术一直是勘探地球物理领域内最为活跃的研究内容之一。

利用精确的波动方程数值解模拟复杂地下波场，为研究地震波传播机理、地震资料的特殊处理方法以及复杂地层的解释等许多方面提供更为科学的数学物理依据［1－3］。

在实际勘探中仍然存在许多复杂问题值得深入探讨，这对研究波动理论，指导地震资料的特殊处理与复杂构造解释具有重要意义。

例如，多波多分量地震记录的每一个分量均包含不同的波型，简单地把垂直分量看成P波，水平分量看成S波是很不合理的。

常规地震资料处理总是希望处理单一分量的标量波场，因此从混合波场中解析出纯P波和纯S波场是进行速度分析、偏移成像等地震资料数据处理的前提。

马德堂等［4］提出满足P波为无旋场、S波为无散场的等价方程思路，采用虚谱法来实现波场分离数值模拟［5－6］，得到较好的数值模拟结果，但是虚谱法很难处理吸收边界条件和自由表面边界条件，同时计算时间和存储量的代价很大，计算网格节点数必须满足2N，且可能存在着Gibbs效应。

三维似French模型高精度弹性波正演数值模拟及其波场特征刘阿男;刘思彤;陈可洋;杨微【摘要】基于三维一阶双曲型各向异性介质速度-应力弹性波动方程,给出各向异性介质向各向同性介质转化的条件,采用高阶交错网格有限差分法和PML吸收边界条件,实现高精度三维弹性波正演数值模拟.以均匀介质为例,研究弹性波在三维空间的波场特征;以三维似French模型为例,合成不同时刻的三分量波场快照和共炮点记录,并对比分析三维各向同性介质和各向异性介质在该模型中的波场响应特征.数值结果验证了方法的正确性和有效性,同时具有精度高和边界吸收效果明显的优点.【期刊名称】《东北石油大学学报》【年(卷),期】2010(034)004【总页数】6页(P45-50)【关键词】三维弹性介质;高精度正演数值模拟;PML吸收边界条件;似French模型;高阶交错网格有限差分法【作者】刘阿男;刘思彤;陈可洋;杨微【作者单位】东北石油大学,地球科学学院,黑龙江,大庆,163318;东北石油大学,地球科学学院,黑龙江,大庆,163318;大庆油田有限责任公司,勘探开发研究院,黑龙江,大庆,163712;大庆油田有限责任公司,第六采油厂,黑龙江,大庆,163712【正文语种】中文【中图分类】TE19弹性波正演数值模拟在多分量地震资料处理、观测系统采集设计、处理参数试验、解释方案验证等方面发挥越来越重要作用,是地震学和地震勘探学研究领域中的前沿课题之一,也是难度较大的课题之一[1].近些年,随计算机计算能力的增强该技术得到飞跃发展.弹性波正演数值模拟技术大多是基于二维波场观测研究地震勘探问题,原因在于二维弹性波场的计算速度快、算法易实现、剖面易显示,而实际的地球介质是三维空间,二维波场观测忽略三维构造体的侧向效应等很多重要信息,无法准确、真实地刻画三维空间的弹性波场响应特征.只有通过三维弹性波正演数值模拟,才能全面地认识弹性介质与弹性波场之间的耦合关系及弹性波场的空间变化特征,后者对于研究和认识地球介质的各向异性特征尤为重要,然而三维弹性波场正演问题研究困难,其原因在于:三维空间的计算变量较多;空间网格节点配置关系复杂;边界吸收问题处理较难;剖面立体显示困难;计算量和存储量巨大,三维声波正演数值模拟的算法要比弹性波正演数值模拟的算法更为简单、易实现.目前,对三维声波正演问题已得到丰富的研究成果[2-7],如频散问题[8-15](精度问题)、边界吸收问题[16]、稳定性问题[17]等,而对三维弹性波正演问题的研究相对较少,且大多是基于均匀单层介质研究弹性波场的传播特征(波场响应),如IgelH等[18]讨论三维各向异性介质中弹性波数值模拟的有限差分实现方法;张文生等[19]对三维正交各向异性介质弹性波场进行有限差分模拟;裴正林[20-21]采用高阶交错网格有限差分法对三维各向异性介质中弹性波场进行高精度数值模拟;董良国等[17]提出三维任意各向异性介质弹性波动方程的稳定性条件;夏凡等[22-23]提出三维均匀介质弹性波动方程的吸收边界条件及其高精度数值计算方法.笔者借助Linux集群计算机的高速运算能力和大容量存储能力,采用高阶交错网格有限差分法求解三维各向异性介质弹性波动方程,通过数值实例分析弹性波在三维空间经典似French速度模型中的三分量波场响应特征,给出三维各向异性介质向各向同性介质转化的条件,同时使用任意方位波场切片技术(包括波场快照切片、模拟记录切片等),实现三分量弹性波场全方位可视化显示,研究三分量弹性波场的传播规律,以此验证方法的准确有效性.通常三维二阶各向异性介质(VTI介质)弹性波动方程可表示为式中:U、V、W分别为x、y、z方向质点的振动位移;C为三维弹性张量(矩阵组合形式C也可称为物性矩阵,一般情况下沿主对角线方向是对称正定的);ρ为弹性介质密度;t为时间.参数均为空间坐标的函数.当各向异性介质对称轴与观测系统不一致时,需要进行观测坐标旋转或者物性矩阵旋转,以达到合理观测目的.此时需要用到Bond矩阵旋转变换,对物性矩阵旋转运算方程为C′=MCMT,其中M为正交正定矩阵,C′为变换后的物性矩阵.根据Thomsen参数的定义,准纵波相速度Vp0=C33/ρ;准横波相速度Vs0=C44/ρ;度量准纵波各向异性强度的参数ε=(C11-C33)/(2C33);度量准横波各向异性强度的参数γ=(C66-C44)/(2C44);影响垂直于VTI介质对称轴方向附加准纵波速度的参数δ=[(C13+C44)2-(C33-C44)2]/[2C33(C33-C44)].为实现三维各向异性介质向各向同性介质转化,借用各向同性介质中的泊松比ν概念,建立三维VTI介质中准纵波和准横波相速度之间的关系式,即同时可推得当三维各向异性介质中的3个Thomsen参数ε=γ=δ=0时,可转化为三维各向同性介质情况,此时有C33=C11=λ+2μ,C44=C66=μ,C13=λ,其中λ和μ为拉梅系数.于是有各向同性介质纵波速度Vp=为便于应用高阶交错网格有限差分法,将三维二阶各向异性弹性波动方程(式(1))变换为一阶双曲型速度—应力弹性波动方程的微分形式,即式中:u、v、w分别为x、y、z方向质点的振动速度正应力;τyz、τxz、τxy分别为沿x、y、z方向的切应力.采用高阶交错网格有限差分法(时间2阶、空间10阶差分精度)对式(2)进行高精度差分离散(形成差分行式),对离散方程中的各变量进行交错网格节点配置,可得到三维高精度弹性波正演递推算子.在人工截断边界处,采用统一格式的最佳匹配层吸收边界条件(吸收层厚度取为20个节点),消除或削弱截断边界处的反射波.要保证边界计算过程稳定,还需满足差分迭代所需的稳定性条件[17].三维均匀各向同性介质模型尺寸为500m×500m×500m,3个方向空间网格长度为5m,采用纯纵波震源(在3个方向正应力场中加入等量外力扰动)并置震源于模型中央,其最大频率为40Hz(波形为雷克子波),介质的纵波速度为2000m/s;泊松比为0.25;密度为1g/cm3;时间步长取0.5ms.满足稳定性条件,波场快照的记录时间为0.125s.波场任意方向切片显示方式的定义:(x=A,yoz)代表沿x方向深度A处且在yoz平面内的波场切片(其中y为图片的横坐标;z为图片的纵坐标);(t=B,yoz)代表在时间B节点处且在yoz平面内的波场切片,其他情况依此类推.三维任一分量弹性波场在3个相互垂直方向上的波场快照切片见图1.由图1可知:在纯纵波震源激励下,x方向质点速度场中纵波波形与yoz平面对称;y方向质点速度场中纵波波形与xoz平面对称;z方向质点速度场中纵波波形与yox平面对称,且对称面两侧极性相反.这验证在均匀各向同性介质中,三维纯纵波波场具有球对称的特点.纯横波三分量弹性波场在任意方向的质点振动速度和与其相垂直的平面存在对称关系,且极性相一致(见图2,在x分量质点的振动速度处加入外力扰动).数值计算结果表明计算精度较高、频散现象很弱.该三维三分量弹性波场特征分析方法同样适合于表征三维各向异性介质准纵波和准横波波场的特点,受各向异性参数的影响,准纵波和准横波波形不再是球对称关系,而是椭圆形甚至是具有不规则特征的对称关系.采用似French模型的速度模型,包含1个水平层界面、1个倾斜断层、1个隆起构造、1个凹陷构造(见图3).模型尺寸为500m×500m×500m,3个方向空间网格长度为5m,震源位于模型地表中央位置,其最大频率为50Hz,每层密度均为1g/cm3,时间步长为0.5ms,满足稳定性条件,接收时间为1s.对于三维各向同性介质模型,上层介质纵波速度为2000m/s,下层介质纵波速度为2500m/s,泊松比为0.25;对于三维各向异性介质模型,上层介质准纵波速度为2000m/s,下层介质准纵波速度为2500m/s,泊松比为0.25,ε为0.1,δ为0.2,γ为0.在0.20,0.25s时,三维各向同性介质和各向异性介质波场快照中的x、y、z方向的中间位置,且沿平行于yoz、xoz、xoy平面的波场快照切片(只给出三维垂直分量的波场快照)见图4和图5.考虑三维似French速度模型与三维各向同性介质弹性波场快照,经过对比容易识别直达P波、倾斜界面的反射PP波和反射PS波、水平层界面引起的反射PP波和反射PS波、隆起构造的反射PP波和反射PS波、凹陷构造引起的回转反射PP波和反射PS波;弹性波传播至下层介质时形成的透射PP 波和透射PS波,以及弹性波在隆起构造和凹陷构造之间相互作用形成的复杂波(包括转换波能量).与三维各向异性介质弹性波场快照对比,弹性波场更为复杂,波场快照中包含由各向异性特性引起的直达S波,模型中构造形成的反射SP波、反射SS 波、透射SP波、透射SS波,以及这些波在构造之间相互作用形成的复杂波形,直达S波的波场能量比其他反射波能量要强.由图4可知,在不同的切片中,凹陷构造和隆起构造的侧向传播效应较为明显.三维各向同性介质和各向异性介质数值模拟记录的在不同位置处的波场切片(见图5),两者最明显的特征在于各向异性介质中存在直达S波及其转换波,其波场情况明显比各向同性介质更为复杂.由此可知,三维弹性波在三维空间中进行传播,仅考虑二维空间的弹性波传播不够准确和全面.另外,地震波到达模型边界时无边界反射波形成,三维弹性波场的传播过程中无明显的频散现象,说明数值计算方法具有较高精度.(1)采用高阶交错网格有限差分法求解三维一阶双曲型各向异性介质速度—应力弹性波动方程及统一格式的PML吸收边界条件方程构建三维高精度弹性波正演递推算子,给出三维VTI介质与三维各向同性介质之间的转换关系,实现了三维均匀介质和似French模型的弹性波高精度正演数值模拟.(2)使用任意方位弹性波场切片技术实现三维三分量弹性波场全方位可视化,研究三分量弹性波场在三维空间中的波场特征及其传播规律,纵波或准纵波具有平面对称、极性相反的特点;横波或准横波具有平面对称、极性相同的特点.(3)三维弹性波场正演数值模拟技术可以更加全面地表征弹性波在全三维空间中的传播过程,计算方法准确有效.【相关文献】[1] 裴正林,牟永光.地震波传播的数值模拟[J].地球物理学进展,2004,19(4):933-941.[2] 刘金平,张向君,董建平.波动方程显式短算子三维正演模拟[J].吉林大学学报:地球科学版,2004,34(1):138-141.[3] 王润秋,张明,费建博,等.精细积分法三维地震波正演模拟[J].勘探地球物理进展,2006,29(6):394-397.[4] 何兵寿,魏修成,刘洋.三维波动方程的数值频散关系及其叠前和叠后数值模拟[J].石油大学学报:自然科学版,2001,25(1):67-71.[5] 熊晓军,贺振华,黄德济.三维波动方程正演及模型应用研究[J].石油物探,2005,44(6):554-556.[6] 熊高君,张学工,贺振华,等.三维定位原理与三维反射波场模拟[J].矿物岩石,2002,22(3):93-97.[7] 熊高君,贺振华,张毅祥.三维混合延拓一步法波动方程正演模拟[J].成都理工学院学报,1999,26(1):48-51.[8] Alford R M,Kelly K R,Boore D M.Accuracy of finite-difference modeling of the acoustic wave equation[J].Geophysics,1974,39 (6):834-842.[9] Marfurt KJ.Accuracy of finite difference and finite element modeling of the scalar and elastic equations[J].Geophysics,1984,49(5): 533-549.[10] Dablain M A.The application of high-differencing to the scalar waveequation[J].Geophysics,1986,51(1):54-66.[11] Fei T,Larner K.Elimination of numerical dispersion in finite-difference modeling and migration by flux-corrected transport[J].Geophysics,1995,60(6):1830-1842.[12] 董良国,李培明.地震波传播数值模拟中的频散问题[J].天然气工业,2004,24(6):53-56.[13] 吴国忱,王华忠.波场模拟中的数值频散分析与校正策略[J].地球物理学进展,2005,20(1):58-65.[14] 宁刚,熊章强,陈持.波动方程有限差分正演模拟误差来源分析[J].物探与化探,2008,32(2):203-206.[15] 董良国,马在田,曹景忠,等.一阶弹性波方程交错网格高阶差分解法[J].地球物理学报,2000,43(3):411-419.[16] 陈可洋,杨微.优化的三维地震波旁轴近似吸收边界条件[J].勘探地球物理进展,2009,32(3):179-181,206.[17] 董良国,马在田,曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究[J].地球物理学报,2000,43(6):856-864.[18] Igel H,Morap,Riolletb.Anisotropic wave propagation through finite-differencegrids[J].Geophysics,1995,60(4):1203-1216.[19] 张文生,宋海斌.三维正交各向异性介质三分量高精度有限差分正演模拟[J].石油地球物理勘探,2001,36(4):422-432.[20] 裴正林.三维各向异性介质中弹性波方程交错网格高阶有限差分法数值模拟[J].石油大学学报:自然科学版,2004,28(5):23-29.[21] 裴正林.三维各向同性介质弹性波方程交错网格高阶有限差分法模拟[J].石油物探,2005,44(4):308-315.[22] Xia F,Dong L G,Ma Z T.Absorbing boundary conditions in 3D elastic-wave numerical modeling[J].Chinese J.Geophys,2004,47 (1):132-136.[23] Xia F,Dong L G,Ma Z T.The numerical modeling of 3-D elastic wave equation using a high-order,staggered-grid,finite difference scheme[J].Applied Geophysics,2004,1(1):38-41.。

第39卷第2期2021年4月海洋科学进展A D V A N C E S I N MA R I N E S C I E N C EV o l .39 N o .2A pr i l ,2021天然气水合物赋存地层的地震波场模拟及对比范 金1,2,刘洋廷1,2,张晓波2,3,刘 凯1,2,郑彦鹏1,2*(1.自然资源部第一海洋研究所,山东青岛266061;2.青岛海洋科学与技术试点国家实验室海洋地质过程与环境功能实验室,山东青岛266237;3.国家深海基地管理中心,山东青岛266237)收稿日期:2020-01-06资助项目:国家重点研发计划项目 南海多类型天然气水合物成藏原理与开采基础研究(2017Y F C 0307302);中央级公益性科研院所基本科研业务费专项资金资助项目 黏弹性海底沉积物的A V O 响应及物性反演方法研究(2018Q 02);泰山学者工程专项经费项目(T S P D 20161007);国家自然科学基金委员会-山东省人民政府联合基金项目 海洋地质过程与环境(U 1606401);国家自然科学基金青年基金项目 基于柔性边界模型的南海北部陆坡峡谷区海底滑动运动定量研究(41606084)作者简介:范 金(1995 ),男,硕士研究生,主要从事海洋地球物理方面研究.E -m a i l :898109884@q q.c o m *通信作者:郑彦鹏(1972 ),男,研究员,博士,主要从事海底构造与海洋地球物理方面研究.E -m a i l :z h e n g y p @f i o .o r g.c n (陈 靖 编辑)摘 要:通过分析天然气水合物在海洋中的6种主要赋存状态类型,总结了每种赋存状态之间的相互转化关系及其物性参数计算方法,并应用到地震波场的正演模拟中㊂对比研究了声波模型㊁弹性波模型和双相介质模型对各种水合物赋存地层的响应特征,结果表明:1)当地层中存在孔隙充填型水合物且下伏地层不含游离气时,双相介质模拟的含水合物层底界表现负极性特征;当充填结核型水合物时,弹性介质和双相介质模拟的水合物底界反射呈负极性;2)当地层充填颗粒包裹型水合物且下伏地层含游离气时,无论是低频(25H z )条件还是提高子波主频(40H z ),3种介质模拟水合物的地震响应特征都很明显,但水合物层底界反射振幅随偏移距变化的关系存在差异;3)当沉积薄层中充填颗粒间胶结型水合物且下伏地层含游离气时,弹性介质和双相介质模拟水合物薄层底界的反射振幅随偏移距的增大而减小;将水合物类型改为颗粒支撑型并提高子波主频,声波介质和弹性介质模拟水合物层底界的反射振幅随偏移距的增大而减小㊂关键词:水合物;赋存状态;正演模型;波场模拟;对比分析中图分类号:P 738.3 文献标志码:A 文章编号:1671-6647(2021)02-0290-14d o i :10.3969/j.i s s n .1671-6647.2021.02.012引用格式:F A NJ ,L I U Y T ,Z HA N GXB ,e t a l .S i m u l a t i o n a n d c o m p a r i s o no f s e i s m i cw a v e f i e l d s i nn a t u r a l g a s h yd r a tef o r m a t i o n s [J ].A d v a n c e s i nM a r i n e S c i e n c e ,2021,39(2):290-303.范金,刘洋廷,张晓波,等.天然气水合物赋存地层的地震波场模拟及对比[J ].海洋科学进展,2021,39(2):290-303.在中低纬度地区,天然气水合物是一种由充足的气源㊁水源在高压(一般大于10M P a)和低温(一般小于10ħ)条件下络合而成的笼形结晶化合物[1-2]㊂因水合物的气源成分以C H 4为主,另有少量C O 2和H 2S 等气态物质,它又被称为 甲烷水合物 ㊂天然气水合物主要分布在陆地的永久冻土带和大陆边缘水深超过300m 的陆坡带[3-4],在我国典型的矿藏区分别是祁连山冻土区和南海北部陆坡区㊂国际地质调查部门将其认定为地球上尚未大规模开发的可替代性清洁能源,对全球气温变化和碳循环过程有至关重要的作用[5]㊂由于地质和地球物理条件的不同,天然气水合物以不同的状态赋存在地层中㊂在研究不同区域的水合物时,根据水合物与骨架颗粒的接触关系㊁沉积物的粒度等要素,赋存状态有多种划分方法[6-9]㊂胡高伟等[10]发现沉积孔隙中的水合物呈混合分布,在不同的形成阶段会以某种分布状态为主导㊂近年来,学者发现标志天然气水合物存在似海底反射(B o t t o m S i m u l a t i n g Re f l e c t i o n s ,B S R )㊁振幅空白㊁速度异常和极性反2期范 金,等:天然气水合物赋存地层的地震波场模拟及对比291转等地震异常信息[11]㊂高精度地震勘探技术已成为圈定天然气水合物成矿区的主导技术,它的基本任务是通过人工地震激发的各种地震波数据反映含水合物沉积层的地质信息,实现对水合物矿藏的精准定位㊂地震波场的正演模拟是地震勘探的基础㊂通过对不同正演模拟方法的改进,D a b l a i n [12]㊁何兵寿等[13]㊁L i u 和S e n [14]的研究成果不仅保证了计算精度,而且降低了计算成本;基于B i o t -G a s s m a n n 理论,L e e [15]提出一种可预测未固结的含水合物沉积层物性特征的方法;C o r d o n 等[16]利用正演模拟方法研究水合物储层中的衰减特征,发现衰减导致的振幅减小可能会影响水合物体积预测的程度㊂随着探测方法和技术的不断提高,天然气水合物在沉积层中的赋存状态不断被重新认识;然而,前人并没有开展不同正演模拟方法识别不同赋存状态水合物差异性的研究㊂通过研究天然气水合物在地层中的赋存状态,本文总结了水合物6种不同的赋存状态:颗粒间胶结型㊁颗粒包裹型㊁颗粒支撑型㊁孔隙充填型㊁基质与填充物共存型和结核/裂缝型,并讨论了其相互转化的关系,概述了相应物性参数的计算方法;介绍模拟含水合物沉积层地震波场的声波模型㊁弹性波模型和双相介质模型,根据水合物6种赋存状态设定适合的孔隙度㊁饱和度等,基于物性参数计算方法求出地层的弹性参数并分别应用到3种正演模型中;改变沉积层的厚度㊁子波主频信息,借助波场快照与地震记录直观地反映含水合物沉积层的地震特征,对不同赋存状态水合物地震响应的差异性进行对比和分析,为识别水合物的勘探研究以及指导地震波反演㊁成像等提供依据㊂图1 天然气水合物在海底的6种赋存状态[7,10,21-23]F i g .1 T h e s i xm a i no c c u r r e n c e s o f n a t u r a l ga s h yd r a t ebe l o wt h e s e af l o o r [7,10,21-23]1 水合物的赋存状态概述天然气水合物的微观分布主要指水合物和沉积物颗粒在地层中的接触关系[17]㊂微观分布又称为赋存状态㊂天然气水合物在地层中的6种赋存状态及其关系转化见图1㊂状态1和状态2属于胶结状态,将岩石颗粒视为随机填充的球体,水合物发生在接触点(状态1)或在颗粒周围生长(状态2)㊂当仅存少量的水合物时,预测的弹性性能大幅增加,随水合物浓度的逐步增加而保持相对接近[7]㊂在状态2中,水合物和沉积物胶结良好,这在降低孔隙度的同时,极大改变了沉积物的刚度[18-20]㊂水合物在状态3中起着颗粒支撑作用,孔隙一般被水和/或游离气等非黏滞性流体充填,沉积物的剪切性质不受影响;水合物在固体颗粒边缘沉积并向孔隙空间内部生长,影响岩石的刚度和热传导性质㊂值得注意的是,随着饱和度的逐步增加(过程一),水合物可能会由状态1转变成状态2,甚至是状态3㊂水合物在状态4中起填充孔隙的作用,充填物包括水合物㊁水或游离气,每种组分的弹性性质各有不同,所以诱发的孔隙压力也不尽相同;当水合物的饱和度较高时,水合物颗粒与相邻的固体矿物颗粒相桥连,成为固体基质的一部分,增强固体骨架的力学稳定性,状态4就可能变为状态3,导致地层刚度变弱[21]㊂状态5是分别将骨架颗粒和水合物作为基质和填充物处理,作为夹杂物嵌入均匀介质中㊂状态1~状态5都被认为天然气水合物在沉积物中均匀分布,而水合物取芯的证据表明,水合物通常作为浅层泥质沉积物中的结核和裂缝填充物存在,该几何形态在状态6中示出[7]㊂产出于裂缝内的水合物迫使地层岩石张开形成裂缝并充填其中,一般呈网状㊁结核状或脉状,肉眼可观测到呈乳白色晶体状㊂含水合物层不仅有较大的裂缝倾角,更受区域构造主应力控制[22]㊂水合物处于状态3且在局部的饱和度逐渐升高时(过程二),就可能转变成状态6㊂如图1中蓝色箭头所示,在实际的形成过程中,天然气水合物可能首先是胶结沉积物颗粒(状态1),中期以接触沉积物颗粒或悬浮在孔隙中(状态3或状态4)为主,后期又以胶292海洋科学进展39卷结形式(状态2)为主,即多种赋存状态共存㊂学者认为游离气的存在和水合物体量增大可能是造成这种变化产生的原因[10,23]㊂不同类型水合物的存在对沉积层的声学特性有直接的影响,会产生不同的地震波场特征㊂含不同赋存状态水合物沉积层的物性参数计算流程见图2㊂当沉积层中不含水合物时,按照图中主干部分(黑色边框)的计算流程,代入相关的参数值,求出地层整体的弹性模量㊁密度㊁速度;当沉积层中储存水合物时,根据不同的赋存状态,选择图中分支部分(彩色边框)相应的计算流程,求出含不同赋存状态水合物地层的弹性模量㊁密度㊁速度㊂然后,根据不同的正演模型(图2中黑色实心箭头),确定相应的参数㊂因篇幅所限,本文对计算公式不做具体介绍,详见文献[6,22,24-25]㊂图2含水合物沉积层的物性参数计算流程F i g.2 F l o wc h a r t f o r c a l c u l a t i n gp h y s i c a l p a r a m e t e r s o f h y d r a t e-b e a r i n g s e d i m e n t a r y l a y e r2正演模拟的理论模型地震波传播理论是各类地震学科的基础,核心内容之一就是地震波传播和介质模型[26]㊂弹性波动方程是地震波传播的动力学核心,解释了地震波在复杂地层中的传播机理㊂不同的模型中介质的速度㊁相态㊁密度等参数的计算方法各有不同,主要通过孔隙度和饱和度来界定,进而影响地震波的传播㊂本文主要使用二2期范 金,等:天然气水合物赋存地层的地震波场模拟及对比293维各向同性的声波模型㊁弹性波模型和双相介质(即孔隙流体介质)模型,其波动方程详见文献[27-29]㊂在地震波场的模拟方法中,有限差分法因实现简单㊁计算效率高等优势被应用于多项研究㊂董良国等[30]将高阶差分法和交错网格相结合,这一方法目前在地震波场数值模拟研究中被广泛使用㊂除此之外,差分格式的相容性㊁收敛性㊁稳定性和边界条件等问题也已被前人研究[27-28,30-33]㊂图3 双层模型F i g .3 T w o -l a ye rm o d e l 参考E c k e r 对海底沉积物的参数设定[20],本文模拟的沉积层中各矿相的物性参数值见表1,其中沉积骨架由5%的石英㊁35%的方解石和60%的黏土共同构成[34]㊂为对比3种正演模型的差异,设计了一种双层模型(图3),物性参数见表2;其中,上部是孔隙度为35%的饱水沉积层,下部是孔隙度为20%的饱水沉积层㊂计算模拟的波场空间为1500mˑ1500m ,网格大小为5mˑ5m ,时间步长为0.5m s㊂采用交错网格有限差分法,精度为时间2阶㊁空间10阶㊂震源采用主频为25H z 的雷克子波㊂表1 沉积物中各矿相的物性参数T a b l e 1 P h y s i c a l p r o p e r t yp a r a m e t e r s o f v a r i o u sm i n e r a l ph a s e s i n s e d i m e n t s 成 分纵波速度V p /(m ㊃s-1)横波速度V s /(m ㊃s -1)密度ρ/(k g ㊃m -3)体积模量K /G P a 剪切模量G /G P a沉积骨架4370260026303513.8海水(卤水)1500010002.250自由气84903000.2160纯水合物330016809005.62.4表2 双层模型的物性参数表T a b l e 2 P h y s i c a l p r o p e r t yp a r a m e t e r s o f t w o -l a ye rm o d e l 层 位纵波速度V p固相参数ANρ11横波速度V s流相参数Rρ22密度ρ耦合参数Qρ12ρ21上层18964.710.9518437201.0474818431.22-283-283下层24295.742.42250410250.3760023041.33-400-400注:V p 和V s 单位为m ㊃s -1;A ,N ,R 和Q 单位为(ˑ109k g ㊃m s -2);ρ11,ρ22,ρ12,ρ21和ρ单位为k g ㊃m -3图4是0.25s 时刻地震波在3种模型中传播的网格化波场快照,展现了清晰的差异性:激发纯纵波在声波介质中波阻抗界面只存在透射和反射;当地震波到达弹性介质的波阻抗界面时,除了产生透射和反射,还会产生转换波;区别于弹性介质,震源在双相介质中会产生与快纵波相互伴生的慢纵波㊂慢纵波在实际地层传播过程中具有极大的衰减性,在波阻抗界面也会产生更多类型的地震波㊂观察波场快照还可知,基于等效弹性模量模型和B i o t 理论计算的双相介质参数,模拟的纵波走时与声波介质㊁弹性介质中纵波走时存在明显的差异,主要是由不同的拉梅系数计算方法导致的㊂294 海 洋 科 学 进 展39卷注:1为直达快纵波;2为反射快纵波;3为快纵波产生的转换反射横波;4为快纵波产生的转换透射横波;5为透射快纵波;6为快纵波产生的转换反射慢纵波;7为快纵波产生的转换透射慢纵波;8为直达慢纵波;9为慢纵波产生的转换反射快纵波;10为慢纵波产生的转换反射横波;11为慢纵波产生的转换透射横波;12为慢纵波产生的转换透射快纵波图4 0.25s 时刻地震波在3种模型中的网格化波场快照F i g .4 M e s h e d s n a ps h o t s o f s e i s m i cw a v e s i n t h e t h r e em o d e l s a t t i m e 0.25s 3波场模拟及对比分析图5 正演模拟与对比所需的层状模型F i g .5 F o r w a r d s i m u l a t i o na n d l a ye r e dm o d e lf o r t h e p u r p o s e o f c o m p a r i s o n 根据4个不同的层状模型(图5),设计了合理的地层孔隙度㊁不同赋存状态下水合物的饱和度以及游离气的饱和度,并依据图2的流程来计算相应层位的物性参数(表3)㊂采用与前文相同的波场空间㊁采样间隔和震源类型㊂在海水中激发纯纵波,炮点位于近海面处,检波器排列于炮点下方,61道接收,道间距为25m ㊂采用P M L 吸收边界,只考虑地震波在介质中的传播㊂对含不同赋存状态水合物的沉积层进行了大量的正演模拟分析,并挑选出几种含水合物沉积层的波场快照和地震记录,展开对比与分析㊂层状模型1从上往下分别是海水层㊁饱水沉积层㊁含孔隙充填型水合物沉积层和饱水沉积层(图5a )㊂图6~图8给出了低频(子波主频采用25H z)条件下,在声波介质㊁弹性介质和双相介质分别模拟不含游离气的孔隙充填型水合物沉积层的典型波场快照㊁去除直达波的单炮地震记录和含水合物层底界振幅随偏移距变化的曲线㊂1)对比波场快照(图6)可知,地震波在声波介质中的反射波能量要明显强于弹性介质和双相介质,主要是因为纵波在弹性介质中向上部地层反射时,遇到了向下部地层透射的横波,产生能量耗散;在双相介质中,地震波不仅发生能量耗散,且孔隙流体的存在导致地震波将能量按2期范 金,等:天然气水合物赋存地层的地震波场模拟及对比295孔隙度大小分配给了固㊁流相,反射振幅进一步减小㊂2)对比地震记录(图7)可知,当含水合物层下伏的沉积层中不含游离气时,声波介质和弹性介质模拟含水合物层底界只表现出与海底反射极性相同㊁与水合物层顶界反射能量相近的反射;双相介质模拟含水合物层顶界的反射微弱,但底界的负极性反射异常明显㊂3)从含水合物层底界振幅随偏移距变化的曲线(图8)可知,由于所提取曲线受到了其他波的干扰,造成个别位置产生局部跳动,但不影响整体的变化㊂随着偏移距的增大,声波介质模拟的反射振幅明显增大,弹性介质模拟的反射振幅整体上减小,而双相介质模拟的底界振幅是先增大后减小㊂表3 基于等效弹性模量模型和B i o t 理论计算的层状模型物性参数T a b l e 3 P h y s i c a l p r o p e r t yp a r a m e t e r s o f l a ye r e dm o d e l c a l c u l a t e db a s e do n e q u i v a l e n t e l a s t i cm o d u l u sm o d e l a n dB i o t t h e o r y层 位纵波速度V p固相参数ANρ11横波速度V s流相参数Rρ22密度ρ耦合参数Qρ12ρ21海水层15002.25010000100010000饱水沉积层18964.710.9518437201.0474818431.22-283-283含水合物层状态124833.583.08218412000.6763721351.38-343-343状态228533.963.08217712080.7662421081.56-336-336状态326328.321.1521437361.1526021230.17-140-140状态420243.551.1621847360.6663721351.47-343-343状态540278.885.71146520090.756501414-0.46-350-350状态620213.201.1621917350.5965021411.32-350-350含游离气层14890.990.7918576390.4465119390.65-279-279饱水沉积层24295.742.42250410250.3760023041.33-400-400注:V p 和V s 单位为m ㊃s -1;A ,N ,R ,Q 单位为(ˑ109k g ㊃m s -2);ρ11,ρ22,ρ12,ρ21和ρ单位为k g ㊃m -3图6 模拟孔隙充填型水合物的网格化波场快照F i g .6 S n a p s h o t s o f g r i d d e dw a v e f i e l d s t h a t s i m u l a t e p o r e -f i l l i n g h yd r a t es 图7 模拟孔隙充填型水合物的去除了直达波的单炮地震记录F i g .7 As i n g l e s h o t s e i s m i c r e c o r d t h a t s i m u l a t e s p o r e -f i l l i n g h yd r a te sw i t hd i r e c tw a v e s r e m o v e d296海洋科学进展39卷图8模拟含孔隙充填型水合物层底界振幅随偏移距变化的曲线F i g.8 T h e v a r i a t i o n c u r v e o f t h eb a s e a m p l i t u d e o f t h e s i m u l a t e d p o r e-f i l l i n g h y d r a t e l a y e rw i t ho f f s e t层状模型2从上往下分别是海水层㊁饱水沉积层㊁含颗粒包裹型水合物沉积层㊁含游离气层和饱水沉积层(图5b)㊂低频(子波主频采用25H z)条件下,在3种介质分别模拟含游离气的颗粒包裹型水合物沉积层的正演结果㊂图9~图11分别展示了正演模拟得到的波场快照㊁单炮地震记录和含水合物层底界振幅随偏移距变化的关系㊂1)从波场快照(图9)可知,地震波在3种介质中反射的能量大小关系与图6相似,且含水合物高速层与含游离气低速层的界面反射与海底反射极性相反㊂2)对比地震记录(图10)可以看出,3种介质模拟水合物底界都有清晰的B S R响应[35],与波场快照的现象吻合;弹性介质和双相介质模拟含水合物层底界的振幅在红色边框处突然变小,而声波介质模拟不存在这一现象㊂究其原因,猜测是由于临界角的变化或者是极性相反的波形相叠加造成的㊂3)从振幅随偏移距变化的曲线(图11)可知,随着偏移距的增大,声波介质模拟水合物层底界的反射振幅整体上先保持稳定,后快速增大;弹性介质模拟水合物底界的振幅先缓慢减小,后产生震荡式变化;双相介质模拟水合物底界振幅整体上呈震荡式减小㊂与图6~图8相比,当含水合物层下伏沉积层中含游离气时,无论使用哪一种正演模型,水合物层底界的负极性反射都尤为明显㊂区别于图9~图11,将子波主频提高为40H z后,地震波形变窄㊂声波介质模拟水合物底界的反射振图9模拟颗粒包裹型水合物的网格化波场快照F i g.9 S n a p s h o t s o f g r i d d e dw a v e f i e l d s t h a t s i m u l a t e p a r t i c l e-e n c a p s u l a t e dh y d r a t e s图10模拟颗粒包裹型水合物的去除了直达波的单炮地震记录F i g.10 As i n g l e s h o t s e i s m i c r e c o r d t h a t s i m u l a t e s p a r t i c l e-e n c a p s u l a t e dh y d r a t e sw i t hd i r e c tw a v e s r e m o v e d2期范金,等:天然气水合物赋存地层的地震波场模拟及对比297幅连续性较好,弹性介质和双相介质模拟水合物层底界的同相轴出现振幅突然减小的特征㊂相比于图9~图11,提高子波主频对3种正演模型识别含水合物沉积层的影响不大㊂图11模拟含颗粒包裹型水合物层底界振幅随偏移距变化的曲线F i g.11 T h e v a r i a t i o n c u r v e o f t h eb a s e a m p l i t u d e o f t h e s i m u l a t e d p a r t i c l e-e n c a p s u l a t e dh y d r a t e l a y e rw i t ho f f s e t层状模型3从上往下分别是海水层㊁饱水沉积层㊁含颗粒间胶结型水合物薄层㊁含游离气层和饱水沉积层(图5c)㊂低频(子波主频采用25H z)条件下,在3种介质分别模拟含游离气的颗粒间胶结型水合物薄层的正演结果㊂图12~图14分别展示了正演模拟得到的波场快照㊁单炮地震记录和含水合物层底界振幅随偏移距变化的关系㊂1)从波场快照(图12)可以看出,透射纵波在声波介质中传播时的能量最强,弹性介质中次之,双相介质中最弱㊂2)由地震记录(图13)可知,当含水合物沉积层的厚度较薄时,声波介质模拟含水合物层顶界的反射较为明显;在弹性介质中可以清晰地看到薄层顶底面的反射,且纵波走时相差很小,两种极性相反的同相轴相互干涉;而双相介质中能在含水合物薄层与含游离气层之间看到明显的负极性反射,却几乎看不到薄层顶界的反射㊂3)从振幅随偏移距变化的曲线(图14)可以看出,随着偏移距的增大,声波介质模拟薄层底界的振幅逐渐增大,而弹性介质和双相介质模拟的振幅都逐渐减小㊂图12模拟颗粒间胶结型水合物的网格化波场快照F i g.12 S n a p s h o t s o f g r i d d e dw a v e f i e l d s t h a t s i m u l a t e i n t e r p a r t i c l e c e m e n t a t i o nh y d r a t e s图13模拟颗粒间胶结型水合物的去除了直达波的单炮地震记录F i g.13 As i n g l e s h o t s e i s m i c r e c o r d t h a t s i m u l a t e s i n t e r p a r t i c l e c e m e n t a t i o nh y d r a t e sw i t hd i r e c tw a v e s r e m o v e d298海洋科学进展39卷图14模拟含颗粒间胶结型水合物层底界振幅随偏移距变化的曲线F i g.14 T h e v a r i a t i o n c u r v e o f t h eb a s e a m p l i t u d e o f t h e s i m u l a t e d i n t e r p a r t i c l e c e m e n t a t i o nh y d r a t e l a y e rw i t ho f f s e t将层状模型3(图5c)中含水合物薄层中水合物的类型改为颗粒支撑型后进行分析,得到高频(子波主频采用40H z)条件下,在3种介质中分别模拟含游离气的颗粒支撑型水合物薄层的正演结果(图15~图17)㊂1)从波场快照(图15)可知,在提高子波主频后,地震波形变窄㊂声波介质和弹性介质模拟水合物薄层顶底面的反射可以被区分㊂与之相比,双相介质模拟水合物薄层顶底界的反射能量要弱很多,不易被识别㊂2)从地震记录(图16)可知,声波介质和弹性介质模拟的水合物薄层出现了相互叠加的同相轴,但不影响薄层顶底面的识别㊂3)从振幅随偏移距变化的曲线(图17)可知,随着偏移距的增大,声波介质和弹性介质模拟薄层底界的振幅整体上逐渐减小,而双相介质模拟的振幅整体是增加的㊂与图12~图14相比,图15~图17中声波模型对水合物薄层顶底面的识别能力在高频条件下得以增强,弹性波模型识别水合物薄层的能力差别不大,而双相介质模型对其识别能力在高频条件下得以削弱㊂图15模拟颗粒支撑型水合物的网格化波场快照F i g.15 S n a p s h o t s o f g r i d d e dw a v e f i e l d s t h a t s i m u l a t e p a r t i c l e-s u p p o r t e dh y d r a t e s图16模拟颗粒支撑型水合物的去除了直达波的单炮地震记录F i g.16 As i n g l e s h o t s e i s m i c r e c o r d t h a t s i m u l a t e s p a r t i c l e-s u p p o r t e dh y d r a t e sw i t hd i r e c tw a v e s r e m o v e d2期范金,等:天然气水合物赋存地层的地震波场模拟及对比299图17模拟含颗粒支撑型水合物层底界振幅随偏移距变化的曲线F i g.17 T h e v a r i a t i o n c u r v e o f t h eb a s e a m p l i t u d e o f t h e s i m u l a t e d p a r t i c l e-s u p p o r t e dh y d r a t e l a y e rw i t ho f f s e t层状模型4从上往下分别是海水层㊁饱水沉积层㊁含结核型水合物层和饱水沉积层(图5d)㊂低频(子波主频采用25H z)条件下,在3种介质分别模拟含结核型水合物沉积层的正演结果(图18~图20)㊂多个聚集成块状体的纯水合物(以下称 结核 )分布在沉积层中,形成多个波阻抗界面㊂1)从波场快照(图18)可知,3种介质模拟波场的能量大小关系和前文模拟的相似,在水合物结核处出现明显的振幅空白特征,主要是因为水合物结核内部阻抗差异较小,反射较弱;各种地震波之间相互干涉,波场异常复杂㊂2)从地震记录(图19)可知,由于水合物下部不存在游离气,声波介质模拟水合物底界的反射呈正极性,与图7a模拟的结果一致;弹性介质模拟水合物底界的反射呈负极性,这与图7b模拟的结果相反;双相介质模拟水合物底部的反射与弹性介质模拟的规律一致,与图7c模拟的结果相似㊂3)从振幅随偏移距变化的曲线(图20)可知,3种介质模拟水合物底界的振幅变化有着相似的规律,且除了双相介质模拟的在大偏移距处能量较小之外,其他位置基本一样㊂图18模拟结核型水合物的网格化波场快照F i g.18 S n a p s h o t s o f g r i d d e dw a v e f i e l d s t h a t s i m u l a t en o d u l eh y d r a t e s图19模拟结核型水合物的去除了直达波的单炮地震记录F i g.19 As i n g l e s h o t s e i s m i c r e c o r d t h a t s i m u l a t e sn o d u l eh y d r a t e sw i t hd i r e c tw a v e s r e m o v e d300 海 洋 科 学 进 展39卷图20 模拟含结核型水合物层底界振幅随偏移距变化的曲线F i g .20 T h e v a r i a t i o n c u r v e o f t h eb a s e a m p l i t u d e o f t h e s i m u l a t e dn o d u l eh y d r a t e l a y e rw i t ho f f s e t 分析上述试验,对比结果见表4㊂其中,3种介质模拟水合物层底界的反射振幅大小整体高于10-3数量级的标记为强振幅,否则,为弱振幅㊂在利用声波介质与弹性介质模拟上层含20%水合物且下层不含游离气的波场时,孙小芳[36]发现含水合物层底界基本不产生负极性的B S R 现象且振幅较小;相同的地质模型下,吴志强等[37]也认为很难看到明显的似海底反射现象㊂本文图6~图8中设定的地质模型和水合物饱和度与之相似,得到的试验结果也基本一致㊂当下伏介质含少量游离气且上层介质含适量的水合物时,含水合物层底界就会产生强振幅的负极性反射[36],与本文图9~图11中得到的试验结果基本一致㊂吴志强等[37]在研究水合物饱和度的变化对沉积层顶底界振幅的影响时,发现在其他条件相同的情况下,水合物饱和度达到40%以上将使含水合物层的顶界出现明显的反射㊂同时,该研究还认为提高子波主频对分辨薄层水合物十分有效㊂这两点在本文图12~图14和图15~图17的对比中都得到了验证㊂表4 三种介质模拟不同类型水合物地震波场的对比T a b l e 4 C o m p a r i s o no f t h r e em e d i a s i m u l a t i n g d i f f e r e n t t y p e s o f h yd r a te s e i s m i cw a v ef i e l d s 位 置正演模型基于不同地质条件得到的地震记录特征图6~图8图9~图11图12~图14图15~图17图18~图20含水合物层底界声波介质正极性㊁弱振幅负极性㊁强振幅极性记录不明显㊁强振幅负极性㊁强振幅正极性㊁弱振幅弹性介质正极性㊁弱振幅负极性㊁强振幅负极性㊁强振幅负极性㊁强振幅负极性㊁弱振幅双相介质负极性㊁弱振幅负极性㊁弱振幅负极性㊁弱振幅负极性㊁弱振幅负极性㊁弱振幅 注:极性分为正极性㊁负极性两种;振幅分为强振幅和弱振幅两种4 结 论本文总结介绍了天然气水合物在海洋中的6种主要赋存状态类型,概述了每种赋存状态之间的相互转化关系及其物性参数计算方法,并基于这些方法,求取了层状模型的各项参数㊂通过改变沉积层的厚度㊁子波主频信息,利用3种正演模拟方法对含不同赋存状态的水合物沉积层进行波场模拟,对比分析波场快照㊁地震记录的差异性,得到以下结论:1)当地层中存在孔隙充填型水合物且下伏地层不含游离气时,声波介质和弹性介质模拟含水合物层底界表现出与海底反射极性相同的弱反射,而双相介质模拟含水合物层底界表现出负极性特征㊂当地层中充填结核型水合物且下伏地层不含游离气时,3种介质模拟水合物底界反射能量相近㊂声波介质模拟水合物底界依然表现出与海底反射极性相同的弱反射,而弹性介质和双相介质模拟水合物底界呈负极性㊂2)当地层中充填颗粒包裹型水合物且下伏地层含游离气时,在低频条件下,3种介质模拟含水合物层底。

地震波震源加载方式对波场特征的影响全红娟;朱光明;王晋国;李桂花【摘要】目的地震波场数值模拟中,研究胀缩源、集中力源、剪切源在各向异性介质中引起的波场特征.方法采用MPML吸收边界,详细地分析了3种震源的加载机理,基于交错网格高阶有限差分法的波动方程进行数值模拟.结果不同震源在介质中对波场类型、偏振方向、各分量的能量分配等特征有较大的影响.结论结果对横波分裂的研究和野外观测系统的震源研究有一定的理论指导意义.%Aim In the numerical simulation of seismic wave, studying the wave field characteristics originated from expansion source, concentrated force source and shear source. Methods Using MPML absorbing boundary, analyzing in detail the loading mechanism of three kinds of source based on wave equation of straggered-grid high-order finite difference method. Results Different souce in the medium will have influence in wave type, polarization direction and energy distribution of every component, etc. Conclusion The result is constructive to study shear wave splitting and source research of field observation system.【期刊名称】《西北大学学报（自然科学版）》【年(卷),期】2012(042)006【总页数】5页(P902-906)【关键词】各向异性;交错网格;胀缩源;集中力源;剪切源【作者】全红娟;朱光明;王晋国;李桂花【作者单位】长安大学地质与测绘工程学院,陕西西安710054;长安大学理学院,陕西西安710064;长安大学地质与测绘工程学院,陕西西安710054;长安大学理学院,陕西西安710064;山东科技大学地质科学与工程学院,山东青岛266590【正文语种】中文【中图分类】P315.3在弹性波场模拟中，震源模拟是一个非常重要的方面，震源加载方式及震源模拟的好坏直接影响数值模拟结果。

基于有限差分法的地震波数值模拟研究综述姚铭;高刚;周游;蔡伟祥;周永娇【摘要】针对有限差分算法在进行地震波数值模拟过程中的模拟精度问题,在进行广泛文献调研后,分析了国内外关于提高有限差分法模拟精度的进展研究,主要包括提高有限差分的阶数及改进离散方式2方面。

论述了有限差分法与其他数值模拟方法的优缺点,指出了有限差分法具有占用内存小、计算量小且易于实现等诸多优点,并总结了有限差分法在进行数值模拟过程中存在的主要问题,包括震源函数、边界条件、频散问题以及稳定性分析。

同时指出了有限差分法在地震波场模拟应用中的发展趋势。

【期刊名称】《能源与环保》【年(卷),期】2017(039)010【总页数】6页(P75-79,85)【关键词】有限差分;数值模拟;震源函数;边界条件;频散问题;稳定性分析【作者】姚铭;高刚;周游;蔡伟祥;周永娇【作者单位】[1]长江大学油气资源与勘探技术教育部重点实验室,湖北武汉430100;[2]长江大学地球物理与石油资源学院,湖北武汉430100;[3]中石化重庆涪陵页岩气勘探开发有限公司,重庆408014【正文语种】中文【中图分类】P631.4地震波数值模拟技术作为波动方程正演、逆时偏移和全波形反演等问题的关键，一直都是理论研究的重要内容［1］。

其理论基础主要包括波动理论和射线理论［2］。

其中射线理论的基本原理是直接求出质点运动方程的近似解，这样虽然能够大大提高计算的速度，但其精度较低且缺少波动力学信息。

而波动理论则是通过数值计算方法求出波动方程的解，因此其模拟结果中包含了丰富的波动信息，便于了解地下地质情况，模拟结果也较精确。

基于波动理论的数值模拟方法主要包括边界元法、有限元法、有限差分法以及伪谱法等，这些方法在数值模拟过程中有其各自的优点，但问题也较突出［3-6］。

如边界元法虽然易于处理边界问题，但其计算量较大;有限元法虽然网格剖分灵活，且能够精确模拟复杂地质形态，但存在占用内存较大的问题，无法普及;伪谱法能够精确模拟地震波的传播规律，但其计算速度较慢，在进行并行计算时也会出现诸多问题，因此适用性不广。

一种基于全局优化的交错网格有限差分法印兴耀;刘博;杨凤英【摘要】在地震波场数值模拟中，交错网格有限差分技术得到了广泛的应用，但是在弹性模量变化较大时，通常会因插值而导致模拟误差增大。

旋转交错网格可以很好地克服这个缺点，因而适合于各向异性介质正演模拟。

但是对于同样大小的网格单元，旋转交错网格需要的步长比常规交错网格要大，这会使梯度和散度算子的误差增大因而更易产生空间数值频散。

针对这些问题，本文提出了旋转交错网格与紧致有限差分相结合的方法，并基于模拟退火算法进行全局优化，压制数值频散，拓宽波数范围。

数值模拟结果表明，此方法可以有效地压制数值频散，且具有较高的模拟精度。

%Staggered grid finite difference techniques have been widely used in numerical simulation of seismic wave field,but the interpolation might enlarge the simulation errors when the elastic moduli change dramatically.Rotated staggered grid can overcome this drawback thus be more applicable for modeling the anisotropic media.However,for the grid cells with the same size,rotated staggered grid needs larger step than conventional staggered grid does,conse-quently,the gradient operator and divergence operator generate larger errors which are more likely to cause spatial numerical dispersion.Aiming at the above problems,a combination of rotated staggered grid and compact finite difference method is proposed,simultaneously,global optimization is preformed based on simulated annealing algorithm to suppress numerical dispersion and to broaden the range of wavenumber.Numerical simulation results show thatthis method can effectively suppress the numerical dispersion and exhibits high simulation accuracy at the same time.【期刊名称】《地震学报》【年(卷),期】2015(000)002【总页数】11页(P278-288)【关键词】全局优化;旋转交错网格;紧致有限差分;数值频散;正演模拟【作者】印兴耀;刘博;杨凤英【作者单位】中国山东青岛 266580 中国石油大学华东地球科学与技术学院;中国山东青岛 266580 中国石油大学华东地球科学与技术学院;中国山东青岛 266580 中国石油大学华东地球科学与技术学院【正文语种】中文【中图分类】P315.3+1引言地震波数值模拟是在地震波传播理论的基础上，通过数值计算来模拟地震波在地下介质中的传播（董良国，2003），是研究地震波传播特性与地球介质参数关系的重要手段.交错网格有限差分技术已经广泛应用于各向异性介质（Igel et al，1995；Collino，Tsogka，2001；裴正林，王尚旭，2005；殷文等，2006；何燕，2008）、孔隙介质（Dai et al，1995；Zeng et al，2001；孙林洁，2012）、黏弹性介质（白晓寅，2008；孙成禹等，2010）的波动方程模拟中，但是常规交错网格在模拟非均匀性较强的介质和各向异性介质时需对介质参数进行平均或内插（陈浩等，2006），低阶插值可能会导致精度降低，高阶插值则会导致计算量增大.Gold等（1997）提出了旋转交错网格有限差分算法，Saenger等进一步研究了旋转交错网格的稳定性与频散的关系（Saenger et al，2000；Saenger，Bohlen，2004）.在旋转交错网格中采用沿网格对角线差分的方式，可以避免弹性模量的平均与插值，因而更准确.随后一些研究人员通过旋转交错网格的完全匹配层（perfectly matched layer，简写为PML）边界条件的实现，利用旋转交错网格技术对各向异性介质（李敏，刘洋，2012）和黏弹介质（严红勇，刘洋，2012）等进行了数值模拟.但由于旋转交错网格中差分是沿对角线进行，步长是同等条件下常规交错网格的■2倍（网格是正方形的情况下），若用常规的差分格式会导致离散的梯度算子和散度算子在计算时误差增大，因而容易出现数值频散.而如果使用高阶有限差分则会因需要更多的网格点而增大计算量与内存.紧致有限差分是一种隐式差分格式，在同等网格数的条件下，具有比常规差分格式更高的精度以及更低的数值频散，故成为目前关注的热点.Lele（1992）分析了紧致差分格式的分辨率特性；王书强等（2002）对弹性波方程的紧致差分格式及其与中心差分的误差进行了研究；Du等（2009）基于交错网格紧致有限差分进行了横向各向同性介质的正演模拟；Boersma（2011）利用6阶紧致有限差分求解了Navier-Stokes方程.实现有限差分时，差分系数的确定会影响数值模拟的精度.有限差分系数的求取一般包括泰勒展开和最优化算法两种方法（Liu，2013）.基于泰勒展开得到的差分系数在波数较小时精度高，但随波数增大其精度会降低；基于最优化算法的差分系数求取方法是在给定精度与差分算子长度的情况下尽可能拓宽波数范围（Holberg，1987）.Shan（2009），Kosloff等（2010），Zhou和Zhang（2011）以及Liu（2013）等对基于最小二乘算法的差分算子优化方法进行了深入研究；Zhang 和Yao（2012）则基于模拟退火算法对中心差分算子进行了优化.本文拟结合旋转交错网格与紧致有限差分技术，基于模拟退火全局优化算法对紧致差分算子进行优化，拓宽数值模拟的波数范围，在此基础上进行频散分析，最后通过数值模拟验证该方法的可行性和有效性.1 旋转交错网格紧致有限差分方法1.1 一阶速度－应力方程在弹性波正演模拟中，除了使用二阶方程外，还常常采用一阶速度－应力弹性波方程，其主要优点是无需对弹性常数进行空间差分（Virieux，1984）.这里根据运动平衡方程和本构关系给出二维情况下体力为零时各向异性介质的一阶速度－应力方程：式中：v，σ分别代表速度和应力，1表示x方向，3表示z方向；cij为介质弹性张量矩阵的元素，当极端各向异性介质时，该矩阵有21个独立参数；当二维时，不考虑y分量，故张量矩阵中与y分量有关参数为0，只有5个独立参数（c11，c13＝c31，c15＝c51，c33，c35＝c53，c55）.1.2 紧致有限差分在地震波正演模拟中，采用常规的中心有限差分往往需要较多的网格点才能比较有效地压制数值频散，不利于边界的处理（王书强等，2002）.紧致有限差分是一种隐式差分格式，只需要较少的网格点即可有效地压制数值频散，且同等阶次的精度要比常规中心有限差分格式高，因此本文采用紧致有限差分格式，弥补旋转交错网格因差分步长大导致梯度和散度算子在计算时误差增大因而更容易产生数值频散的不足.常规的2（M－1）点2 M阶紧致有限差分格式（Kim，Lee，1996）为式中，Δx为网格的空间步长，a和bn分别为差分的权系数，可以通过在第i个网格点的泰勒展开，然后对比对应系数来求解即可得到2（M－1）点2 M阶紧致有限差分格式的差分系数（Liu，Sen，2009），例如当M＝5时，可以得到8点10阶的紧致差分格式系数为a＝0.257 894 74，b1＝0.889 871 16，b2＝0.216 121 16，b3＝－0.004 701 2，b4＝0.000 151 55.1.3 旋转交错网格及完全匹配层吸收边界条件在二维旋转交错网格上，密度和速度的各个分量定义在相同位置，弹性模量和应力的各个分量定义在另外的相同位置上，如图1所示.旋转交错网格在计算时分为两步（陈浩等，2006）：首先，计算沿着对角线方向即˜x，˜z方向进行差分，得到相关物理量的对角线方向的空间一阶导数；然后，通过坐标旋转的换算关系，将得到的两个对角线方向的差分进行线性组合从而获得坐标轴方向即x，z方向的一阶空间导数，其换算关系为（Saenger et al，2000）：图1 旋转交错网格及完全匹配层吸收边界示意图x，z为坐标轴方向，x˜，z˜为对角线方向；灰色区域表示沿着x和z方向均进行衰减，白色区域表示只沿着x方向或者z方向进行衰减Fig.1 Schematic diagram of rotated staggeredgridand perfectly matched layer（PML）absorbing boundary xand zare the coordinate directions，andx˜andz˜ are the diagonal directions.Gray areas indicate that waves are absorbed along both xand z directions and white areas indicate that waves are absorbed only along the xor zdirection式中，∂／∂x˜和∂／∂z˜是沿对角线方向的导数；∂／∂x和∂／∂z是沿坐标轴方向的导数.利用波动方程进行数值模拟，一个关键问题就是边界条件.为了解决边界反射问题，引入弹性波方程的PML边界条件.陈浩等（2006）指出，在旋转交错网格中，PML吸收边界条件的处理方式及吸收效果与常规交错网格中几乎是一样的.为此，引入图1所示的PML吸收层.PML边界的基本做法是在研究区域四周引入PML，波在PML中传播时不会产生反射，并且随传播距离按一定规律衰减.当波传播到PML边界时，波场近似为零，也不会产生反射（王守东，2003）.以式（1）中的vx为例来说明旋转交错网格紧致有限差分算法PML的实现方法：式中表示n＋1／2时刻vx 在网格点（i，j）处沿坐标轴x方向的分量，表示n－1／2时刻的x方向分量，di为衰减因子，表示n 时刻σxx在网格点（i，j）处沿坐标轴x方向的导数.根据式（5）可推导出式中沿着所在对角线˜x与斜对角线˜z方向的导数值，根据求出.其中类似地可以推导出其它分量的旋转交错网格紧致有限差分下的PML控制方程. 1.4 紧致有限差分的全局优化根据平面波理论，令式中，β＝kΔx／2，0≤β≤π／2.为了使式（2）左右两边的误差尽可能小，须满足β－β＊≈0的β取值范围应尽可能地大.本文中，采用模拟退火全局最优算法（Kirkpatrick et al，1983）来优化旋转交错网格紧致有限差分算子.为此，建立以下目标函数：Kim和Lee（1996）指出，当θ接近最大值π／2时，会出现很多难以控制的误差，从而导致优化效果不佳，因此本文选择θ＝rπ／2，其中0＜r＜1.表1列出了优化的8点10阶，10点12阶，12点14阶以及14点16阶的紧致有限差分系数.表1 优化的10—16阶紧致有限差分系数Table 1 Optimized coefficients of 10—16-order compact finite-difference阶次 r a b1 b2 b3 b4／10－2 b5／10－5 b6／10－3 b7／10－5 10 0.85 0.384 0.728 0.369 －0.015 0.152 12 0.87 0.390 0.720 0.376 －0.016 0.171 －6.49 14 0.89 0.434 0.658 0.434 －0.022 0.361 －74 0.138 16 0.85 0.452 0.633 0.458 －0.0252 0.452 －114 0.322 －7.792 频散分析由式（12）可以看出，β接近于β＊的程度表征了优化的紧致有限差分算子的精度，故定义使α≈1的β取值范围越大，数值频散越小.图2给出了常规的基于泰勒展开得到的差分算子与本文提出的全局优化算子的精度对比.可以看出，对于常规的基于泰勒展开得到的紧致有限差分算子，随着空间阶次的增大，使α≈1的波数范围变得越宽，即数值模拟的精度随阶次的增大而增大.显然，本文得到的全局优化的紧致有限差分算子在相同最大误差下使α≈1的波数范围比常规紧致有限差分算子还要大.从图2中还可以看出，优化的8点10阶紧致有限差分算子的频散比常规14点16阶紧致有限差分算子小，即优化的10阶精度要高于常规的16阶精度，从而可以使用较少的网格点来达到较高的精度，节省计算内存.图2 10阶常规有限差分，10—16阶常规紧致有限差分及对应的优化紧致有限差分的频散对比（k代表波数，Δx表示空间网格间距）Fig.2 Dispersion comparison by 10-order conventional finite-difference，10—16-order conventional compact finite-difference and the corresponding optimized compact finitedifference.k denotes the wavenumber，andΔxis spatial grid spacing由图2中的10阶常规有限差分、10阶常规紧致有限差分与10阶优化紧致有限差分的频散曲线可以看出，10阶常规紧致有限差分格式的波数范围较10阶常规有限差分要宽，而经过优化的10阶紧致有限差分格式的波数范围则更进一步拓宽.设最大允许误差为εmax＝｜α－1｜0≤β≤βmax，取εmax＝0.2%，分别计算得到10阶常规有限差分、10阶常规紧致有限差分及10阶优化紧致有限差分的最大波数范围βmax分别为0.847 5，1.037 5，1.458 0，即10阶常规紧致有限差分的波数范围比常规差分格式拓宽了1.22倍，而经过全局优化后，波数范围提升至1.72倍.考虑到正演模拟中最常采用Δx＝Δz的情况，则对于旋转交错网格，Δr比常规交错网格差分步长Δx及Δz增大了2倍，因而可以采用优化系数的紧致有限差分来弥补这一不足，从而有效地压制数值频散.需要指出的是，优化系数的紧致有限差分格式也可以很方便地应用于PML边界条件，而不需要特殊的处理.3 模型试算3.1 模型1为了验证全局优化的旋转交错网格紧致有限差分的精度比常规旋转交错网格有限差分方法有所提升，首先采用一个简单的各向同性模型进行测试（模型大小为3 000m×3 000m，在测试中并未添加任何边界条件），震源采用主频30Hz的雷克子波，z方向集中力源激发，网格大小为Δx＝Δz＝10m，时间采样为Δt＝1ms.纵横波速度分别为4 000m／s和2 600m／s，模型密度为2 300kg／m3.图3给出了3种有限差分条件下400ms时刻x分量波场快照，均采用旋转交错网格系统.图3a采用常规10阶常规有限差分格式；图3b采用8点10阶，10点12阶，12点14阶，14点16阶常规紧致有限差分格式；图3c采用优化的8点10阶紧致有限差分格式.可以看出，10阶常规有限差分与10阶紧致有限差分对比，后者频散要小；而对于图3b中的紧致有限差分格式，随着空间阶次的增大，其精度越高，频散越小；对比图3b与图3c可以看出，优化的8点10阶紧致有限差分的数值频散甚至比常规14点16阶紧致有限差分还要小，这与图2所示是一致的.图4给出了从图3中提取的z＝500m处（即50个网格点处）的波形曲线，以此来对比以上3种差分格式的数值频散.可以看出，10阶常规有限差分算法在波形外出现了数值抖动，说明数值频散比较高，模拟精度相对较低，而10阶紧致有限差分和优化的10阶紧致差分精度较高，从虚线框中可以进一步看出优化的10阶紧致有限差分算法的精度更高一些.图3 10阶常规有限差分（a），10—16阶常规紧致有限差分（b）与优化的10阶紧致有限差分（c）得到的波场快照Fig.3 Snapshots of wave field by using 10-order conventional finite-difference（a），10－16-order conventional compact finite-difference（b）and optimized 10-order compact finite-difference（c）3.2 模型2图4 10阶常规有限差分，10阶常规紧致有限差分及10阶优化紧致有限差分波形对比图Fig.4 Comparison of waveforms by using 10-order conventional finite-difference（FD）（dark line），10-order conventional compact FD （blue line）and optimized 10-order compact FD （red line）模型2的几何构造如图5所示，其中凹陷构造上部及下部均为各向同性介质，最下层为横向各向同性（VTI）介质.模型纵向和横向长为3 000m，具体参数见表2.采用优化的旋转交错网格紧致有限差分方法进行正演模拟，网格大小为Δx＝Δz＝10m，时间采样为Δt＝1ms.震源采用30Hz雷克子波，在模型的（1 150m，1 500m）处以纵波源形式激发.模拟过程中，采用PML边界条件（PML层部分同样采用优化的旋转交错网格紧致有限差分算法），厚度层为30个网格点.为体现吸收效果包含了PML层，图6给出了600ms时刻x及z分量的波场快照.可以看出，震源激发的纵波传播至分界面时，发生了反射与透射，可以见到PP波、PS波等，并可在凹陷处看到断面波，波场复杂，但波前面清晰，数值频散小.从图7所示的单炮记录中也可以看出优化的旋转交错网格紧致有限差分算法具有较高的模拟精度，各类波均可得到清晰体现，并且能够证明所加PML边界条件的吸收效果很好，无明显人为边界反射产生.图5 凹陷模型示意图Fig.5 Schematic diagram of the depression model图6 600ms时刻优化的旋转交错网格紧致有限差分正演x分量（左）和z分量（右）的波场快照（虚线框外围代表PML层）Fig.6 Wavefield snapshots of x-（left）and z-component（right）at 600ms using optimized rotated staggered-grid compact finite-difference method（PML layers are outside the dashed box）表2 凹陷模型参数Table 2 Parameters of the depression model序号介质类型 c11／GPa c13／GPa c33／GPa c55／GPa 密度／（kg·m－3）31.77 42.68 31.77 13.75 2200第二层各向同性 42.33 47.04 42.34 18.82 2400第三层横向各向同性第一层各向同性87.88 21.22 67.60 22.56 2500图7 凹陷模型优化的旋转交错网格紧致有限差分正演x分量（左）和z分量（右）单炮记录Fig.7 Records of x-（left）and z-component（right）of the depression model using optimized rotated staggered-grid compact finite-difference method4 讨论与结论旋转交错网格在地震波正演模拟中由于其可以避免在模拟中因弹性模量插值或平均而产生的误差和其对边界条件处理的便捷性而被广泛应用.但由于其差分方向是沿网格的对角线方向，因此相对于常规的交错网格而言，其差分步长较大，梯度算子、散度算子在计算时更容易产生误差，这意味着旋转交错网格更容易产生数值频散，从而影响模拟的精度.若采用小网格，虽然可以压制数值频散，却又大大增加了运算量及内存.因此本文采用旋转交错网格与紧致有限差分技术相结合的方法来提高模拟精度，并在此基础上，基于模拟退火算法对紧致有限差分进行全局优化，进一步提高计算精度，压制数值频散.数值频散分析结果表明，优化的旋转交错网格紧致有限差分算法相对于普通旋转交错网格有限差分算法具有更宽的波数范围，这意味着优化的旋转交错网格有限差分算法可以在采用更大的空间采样间隔或者更高的震源主频时依然有较高的精度和较低的数值频散；且经过全局优化的10阶紧致有限差分算子比常规的16阶紧致有限差分算子具有更高的精度，即可以采用较少的网格点来达到较高的模拟精度，因此可以节省计算内存.数值模拟实验进一步证实了本文提出的优化的旋转交错网格紧致有限差分算法的正确性与可行性.凹陷模型中，各类波均具有比较清晰的响应，揭示了波场的传播规律.同时波场快照与单炮记录中无明显的人为边界反射产生，也证明了模拟采用的PML吸收边界条件的有效性.但是本文提出的方法亦存在不足之处.在求解紧致差分格式时需要解三角矩阵线性方程，较常规显式有限差分格式增加额外的计算负担，且不利于在GPU等平台直接并行，例如图3a和图3c中的400ms时长的计算时间分别为46.5s和85.6s，后者计算耗时约为前者的1.8倍.在科学计算中，效率与计算精度常常不可兼得，随着计算机科学计算能力的提高，可以采用GPU加速的LU分解算法或者使用LAPACK等优化程序包来克服求解大型矩阵的计算效率问题.在实际应用中需权衡利弊，选择适合的算法.衷心感谢审稿专家提出的宝贵意见和建议.参考文献白晓寅.2008.基于地震波衰减理论的地层吸收参数提取方法研究［D］.东营：中国石油大学（华东）地球科学与技术学院：24-54.Bai X Y.2008.The Study on Method of Layer Absorbing Parameters Extraction Based on the Theory of Seismic Wave Attenuation［D］.Dongying：School of Geosciences，China University of Petroleum：24-54（in Chinese）.陈浩，王秀明，赵海波.2006.旋转交错网格有限差分及其完全匹配层吸收边界条件［J］.科学通报，51（17）：1985-1994.Chen H，Wang X M，Zhao H B.2006.Rotated staggered grid finite-difference and the PML boundary［J］.Chinese Science Bulletin，51（17）：1985-1994（in Chinese）.董良国.2003.地震波数值模拟与反演中几个关键问题研究［D］.上海：同济大学海洋与地球科学学院：1-2.Dong L G.2003.Several Key Problems in Seismic Wave Numerical Simulation and Inversion［D］.Shanghai：School of Ocean and Earth Science，Tongji University：1-2（in Chinese）.何燕.2008.正交各向异性弹性波高阶有限差分正演模拟研究［D］.东营：中国石油大学（华东）地球科学与技术学院：45-80.He Y.2008.High-Order Finite-Difference Forward Modeling of Elastic-Wavein Orthorhombic Anisotropic Media［D］.Dongying：School of Geosciences，China University of Petroleum：45-80（in Chinese）.李敏，刘洋.2012.高阶旋转交错网格有限差分方法模拟TTI介质中横波分裂［J］.物探与化探，36（6）：934-940.Li M，Liu Y.2012.Modeling of the S-wave splitting in TTI media using high-order rotated staggered grid scheme［J］.Geophysical and Geochemical Exploration，36（6）：934-940（in Chinese）.裴正林，王尚旭.2005.任意倾斜各向异性介质中弹性波波场交错网格高阶有限差分法模拟［J］.地震学报，27（4）：441-451.Pei Z L，Wang S X.2005.A staggered-grid high-order finite-difference modeling for elastic wave field in arbitrary tilt anisotropic media［J］.Acta Seismologica Sinica，27（4）：441-451（in Chinese）.孙成禹，肖云飞，印兴耀，彭洪超.2010.黏弹介质波动方程有限差分解的稳定性研究［J］.地震学报，32（2）：147-156.Sun C Y，Xiao Y F，Yin X Y，Peng H C.2010.Study on the stability of finite difference solution of visco-elastic wave equations［J］.Acta Seismologica Sinica，32（2）：147-156（in Chinese）.孙林洁.2012.非均匀孔隙介质有限差分正演模拟方法研究［D］.青岛：中国石油大学（华东）地球科学与技术学院：64-108.Sun L J.2012.Finite-Difference Modeling Research in Heterogeneous Porous Media［D］.Qingdao：School of Geosciences，China University of Petroleum：64-108（in Chinese）.王守东.2003.声波方程完全匹配层吸收边界［J］.石油地球物理勘探，38（1）：31-34.Wang S D.2003.Absorbing boundary condition for acoustic wave equation by perfectly matched layer［J］.Oil Geophysi-cal Prospecting，38（1）：31-34（in Chinese）.王书强，杨顶辉，杨宽德.2002.弹性波方程的紧致差分方法［J］.清华大学学报：自然科学版，42（8）：1128-1131.Wang S Q，Yang D H，Yang K pact finite difference scheme for elastic equations［J］.Journal of Tsinghua University：Science and Technology，42（8）：1128-1131（in Chinese）.严红勇，刘洋.2012.黏弹TTI介质中旋转交错网格高阶有限差分数值模拟［J］.地球物理学报，55（4）：1354-1365.Yan H Y，Liu Y.2012.Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media ［J］.Chinese Journal of Geophysics，55（4）：1354-1365（in Chinese）. 殷文，印兴耀，吴国忱，梁锴.2006.高精度频率域弹性波方程有限差分方法及波场模拟［J］.地球物理学报，49（2）：561-568.Yin W，Yin X Y，Wu G C，Liang K.2006.The method of finite difference of high precision elastic wave equations in the frequency domain and wave field simulation［J］.Chinese Journal of Geophysics，49（2）：561-568（in Chinese）.Boersma B J.2011.A 6th order staggered compact finite difference method for the incompressible Navier-Stokes and scalar transport equations［J］.J Comput Phys，230（12）：4940-4954.Collino F，Tsogka C.2001.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropicheterogeneous media［J］.Geophysics，66（1）：294-307.Dai N，Vafidis A，Kanasewich E R.1995.Wave propagation in heterogeneous，porous media：A velocity-stress，finitedifference method ［J］.Geophysics，60（2）：327-340.Du Q Z，Li B，Hou B.2009.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered-grid finite difference scheme［J］.Appl Geophys，6（1）：42-49.Gold N，Shapiro S A，Burr E.1997.Modelling of high contrasts in elastic media using a modified finite difference scheme［C］∥SEG Technical Program Expanded Abstracts.Dallas，Texas：SEG：1850-1853.Holberg putational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena［J］.Geophys Prospect，35（6）：629-655.Igel H，Mora P，Riollet B.1995.Anisotropic wave propagation through finite-difference grids［J］.Geophysics，60（4）：1203-1216.Kim J W，Lee D J.1996.Optimized compact finite difference schemes with maximum resolution［J］.AIAA J，34（5）：887-893.Kirkpatrick S，Gelatt C D Jr，Vecchi M P.1983.Optimization by simulated annealing［J］.Science，220（4598）：671-680.Kosloff D，Pestana R C，Tal-Ezer H.2010.Acoustic and elastic numerical wave simulations by recursive spatial derivative operators［J］.Geophysics，75（6）：T167-T174.Lele S pact finite difference schemes with spectral-like resolution［J］.J Comput Phys，103（1）：16-42.Liu Y，Sen M K.2009.An implicit staggered-grid finite-difference method for seismic modelling［J］.Geophys J Int，179（1）：459-474.Liu Y.2013.Globally optimal finite-difference schemes based on least squares［J］.Geophysics，78（4）：T113-T132.Saenger E H，Bohlen T.2004.Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid［J］.Geophysics，69（2）：583-591.Saenger E H，Gold N，Shapiro S A.2000.Modeling the propagation of elastic waves using a modified finite-difference grid［J］.Wave Motion，31（1）：77-92.Shan G.2009.Optimized implicit finite-difference and Fourier finite-difference migration for VTI media［J］.Geophysics，74（6）：WCA189-WCA197.Virieux J.1984.SH-wave propagation in heterogeneous media：Velocity-stress finite-difference method［J］.Geophysics，49（11）：1933-1957. Zeng Y Q，He J Q，Liu Q H.2001.The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media ［J］.Geophysics，66（4）：1258-1266.Zhang J H，Yao Z X.2012.Optimized finite-difference operator for broadband seismic wave modeling［J］.Geophysics，78（1）：A13-A18. Zhou H B，Zhang G Q.2011.Prefactored optimized compact finite-difference schemes for second spatial derivatives［J］.Geophysics，76（5）：WB87-WB95.。

油气勘探化 工 设 计 通 讯Petroleum ExplorationChemical Engineering Design Communications·171·第47卷第4期2021年4月弹性波正演是认识地下复杂介质的重要手段，在弹性波有限差分正演模拟中，除了二阶质点位移方程外，常用一阶速度-应力方程，其优点是不需要对弹性参数进行空间微分，就可以得到完整的地震波场响应[1]。

除了常规网格外，Madariaga 提出的交错网格，它在不增加计算量的前提下，和常规网格相比局部精度提高了4倍，收敛速度也很快，很好地压制了数值频散[2]。

弹性波纵横波场分离是进行弹性波逆时偏移的必要条件，如若不能正确分离纵横波场，则在偏移成像过程中由于纵横波的互相干扰，导致偏移成像失败。

目前大多数学者都采用传统的Helmholtz 原理分离纵横波场。

Helmholtz 分解基于弹性波混合波场是由一个无旋场和无散场组成的，这个无旋场就是纵波的势函数，无散场就是横波的势函数。

因此可以通过对混合波场做散度运算提取纵波，做旋度运算。

但是由于旋度运算和散度运算对空间做了一阶偏导数，分离后的波场的振幅和相位都会改变，且分离后的纵横波波场与原波场的物理意义不一致。

假设原波场为位移，那么经过Helmholtz 分离后的波场为质点速度相关的一个量。

而且对其做相位矫正和振幅恢复时计算量过大，对后期的逆时偏移来算不易实现。

本文采用纵横波解耦方程有限差分交错网格形式下的正演模拟，吸收边界采用CPML 。

在提高计算效率上利用MPI 和OpenMP 进行并行加速。

1 基本原理与方法1.1 弹性波各向同性一阶速度-应力方程密度非均匀，二维弹性波各向同性一阶速度-应力方程：（1）式（1）中x z 垂直分量（在多波多分量地震勘探中，也被称之为伪横波和伪纵波）；τxx 、τzz 、τxz 为质点振动x 轴向，z 轴向的正应力，和y 平面法向的剪切应力；λ、μ则为弹性波的模量，其中μ是剪切模量；ρ表示为密度参数。

利用高阶交错网格有限差分法数值模拟VTI 介质井孔声场岳崇旺;王飞【摘要】Transversely isotropic(TI) media is a common petrophysical media. It is important to study the propa-gation characteristics of the acoutic field in the well for sonic logging theory, and it can provide the basis for the sonic log interpretation. This paper derived velocity–stress staggered-grid finite-difference equations of the elas-tic wave propagation in cylindrical coordinates for vertical transversely isotropic(VTI) media. Furthermore, it numeri-cally simulated acoustic propagation in the VTI media using finite–difference technique with two orders in time and ten orders in space. It gave the snapshots of borehole acoustic wave field in the homogeneous media at different times and calculated the full wave trains with acoustic sources located at the well axis. The calculated results show that if the coefficient of anisotropy of VTI media increases, the change of shear wave propagation has little effect, but the velocity of longitudinal wave propagation has been reduced relatively in the longitudinal direction, and has little change in the radial direction. And if the coefficient of anisotropy of VTI media increases, the first wave slowness of sonic logging will increase, and the acoustic amplitude will be slightly reduced.%横向各向同性(TI)介质是岩石地球物理中常见的一种现象，研究其井孔声场传播特征对声波测井理论以及为声波测井解释提供依据具有重要意义。

地震波交错网格高阶差分数值模拟研究摘要: 地震波数值模拟技术是勘探地球物理学中的重要组成部分，研究通过弹性波一阶速度——应力方程，采用交错网格高阶有限差分法实现了地震波在各向同性介质中的高精度的数值模拟，并采用完全匹配层( PML) 吸收边界来消除边界反射，可取得较好的效果。

通过模型的正演计算和复杂模型的处理结果表明，交错网格高阶有限差分法数值模拟是一种快速有效的地震波数值模拟方法。

关键词: 地震勘探; 交错网格; 有限差分; 数值模拟引言地震数值模拟是模拟地震波在介质中传播的一种数值模拟技术，随着地震波理论在天然地震和地震勘探中的应用，地震模拟技术便应运而生，并随着地震波理论和计算机技术的发展，地震数值模拟技术自20世纪60年代以来也得到了飞速发展，形成了目前具有有限差分法、有限元法、虚谱法和积分方程法等各种数值模拟方法的现代地震数值模拟技术。

有限差分法是偏微分方程的主要数值解法之一。

在各种地震数值模拟方法中，最早出现的数值模拟方法是有限差分法。

Ａlterman和Karal(1968)首先将有限差分法应用于层状介质弹性波传播的数值模拟中。

此后，Boore(1972)又将有限差分法用于非均匀介质地震波传播的模拟。

Alford等（1974）研究了声波方程有限差分法模拟的精确性。

Kelly等（1976）研究了用有限差分法制作人工合成地震记录的方法。

Virieux(1986)提出了应用速度——应力一阶方程交错网格有限差分法模拟P——SV波在非均匀介质中的传播。

交错网格方法提高了地震模拟的精度和稳定性，并消除了部分假想。

有限元法也是偏微分方程的数值解法之一。

Lysmer和Drake（1972）最早将有限元法应用于地震数值模拟。

Marfurt(1984)研究对比了模拟弹性波传播的有限差分法和有限元法的精度。

Seron等（1990,1996）给出了弹性波传播有限元模拟方法。

Padovani等（1994）研究了地震波模拟的低阶和高阶有限元法。