弹性介质中二维波动方程数值模拟
二维TTI介质的纯P波波动方程数值模拟
二维TTI介质的纯P波波动方程数值模拟张千祥;王德利;周进举【摘要】声波各向异性数值模拟对地震数据处理和解释起着重要的作用.基于Tsvankin提出的精确色散关系,通过平方根近似,在时间-波数域中推导出二维TTI 介质纯P波声波波动方程,并利用快速展开法(Rapid Expansion Method,REM)进行了数值模拟.与传统的有限差分法求解二维TTI介质耦合方程和傅里叶有限差分法在时间上进行波场外推相比,该方法的模拟结果精度更高,计算速度更快,并且成功去除横波分量.【期刊名称】《石油物探》【年(卷),期】2015(054)005【总页数】8页(P485-492)【关键词】声波各向异性数值模拟;纯P波声波方程;快速展开法;有限差分法;傅里叶有限差分法【作者】张千祥;王德利;周进举【作者单位】吉林大学地球探测科学与技术学院,吉林长春130026;吉林大学地球探测科学与技术学院,吉林长春130026;吉林大学地球探测科学与技术学院,吉林长春130026【正文语种】中文【中图分类】P631地震数值模拟是地震勘探方法研究的前提和基础,在地震勘探和地震学的各项研究及生产工作中都扮演着重要的角色[1]。
常用的地震波场数值模拟方法主要有几何射线法、波动方程法和积分方程法[2]。
波动方程模拟方法中的有限差分法由于计算速度快、占用计算机内存小而被广泛应用。
很早时候各国地球物理学家就对各向同性介质和各向异性介质弹性波地震数值模拟进行了深入研究。
近年来,周进举等[3]利用高阶旋转网格有限差分法研究了复杂介质下弹性波数值模拟。
在对地下的各向异性介质进行弹性波数值模拟时,由于弹性波方程复杂,各向异性参数多,导致模拟计算量大,耗时长,增加了弹性波偏移和干涉的难度。
为了解决这些问题,我们采用声波各向同性近似理论,通过设定弹性波中的横波速度为零来简化计算参量,在保证模拟精度的条件下提高计算效率。
然而,由于地下介质的不均匀性,这种各向同性的声学假设常常是不恰当的。
二维地震波场的组合型紧致有限差分数值模拟
二维地震波场的组合型紧致有限差分数值模拟汪勇;石好果;周成尧;桂志先【摘要】地震波场数值模拟在地球物理勘探和地震学中具有重要的支撑作用.本文将组合型紧致差分格式用于声波和弹性波方程的数值模拟中.根据泰勒级数展开和声波方程,建立了位移场时间四阶离散格式,并将组合型紧致差分格式用于位移场空间导数的求取,然后对该差分格式进行了精度分析、误差分析、频散分析和稳定性分析.理论研究结果表明:① 该差分格式为时间四阶、空间六阶精度,与常规七点六阶中心差分和五点六阶紧致差分相比,具有更小的截断误差和更高的模拟精度;② 每个波长仅需要5.6个采样点,且满足稳定性条件的库郎数为0.792,可以使用粗网格和较大时间步长进行计算.所以该方法具有占用内存少、计算效率高和低数值频散等优势.最后,本文进行了二维各向同性完全弹性介质的声波和弹性波方程的数值模拟,实验结果表明本文提出的方法具有更高的计算精度,能够大幅度的节约计算量和内存需求,对于三维大尺度模型问题具有更好的适应性.【期刊名称】《地球物理学报》【年(卷),期】2018(061)011【总页数】16页(P4568-4583)【关键词】组合型紧致差分;声波方程;弹性波方程;数值模拟;数值频散;稳定性条件【作者】汪勇;石好果;周成尧;桂志先【作者单位】油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100;胜利油田勘探开发研究院西部分院,山东东营 257001;油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100;油气资源与勘探技术教育部重点实验室(长江大学),武汉 430100;长江大学地球物理与石油资源学院,武汉 430100【正文语种】中文【中图分类】P6310 引言地震波场数值模拟是应用数值方法模拟地震波在地下介质中的传播过程，计算在地面或地下各观测点地震记录的一种数值模拟方法，是地球物理勘探和地震学的重要基础.目前，常用的数值模拟方法主要有射线追踪法和波动方程法，其中波动方程法有伪谱法、有限元法、边界元法、谱元法和有限差分法(Alterman and Karal，1968；Graves，1996；董良国等，2000a；Saenger et al.，2000，2004).随着数值模拟技术的发展和生产实践的要求，围绕着提高有限差分计算效率(黄超和董良国，2009；唐佳等，2016)、模拟精度(Liu，2013；李振春等，2016)、算法稳定性(Virieux, 1986；董良国等，2000；杜启振等，2015)、压制数值频散、处理复杂介质(Fletcher et al.，2008；王璐琛等，2016；程玖兵等，2013)和吸收边界条件(Berenger，1994；Pan et al.，2012)等方面，研究人员已经发展了多种方法.特别是在如何压制数值频散和提高计算效率方面，杨顶辉等在这方面做了大量相关工作(Yang et al.，2009，2012a,b，2014；Yang and Wang, 2010，Ma and Yang，2017；Tong et al.，2013；Zhou et al.，2015；He et al.，2015；Liu et al.，2017；Jing et al.，2017)，取得了许多有意义研究成果.提高差分格式数值计算精度最直接的方法就是在差分计算时增加网格节点数，但这也增加了计算量和所需的存储空间.紧致差分方法恰好能够较好地解决这个矛盾，同时紧致差分是一种隐式差分格式，具有较好的稳定性，这些优势也使得它成为目前研究较多的有限差分方法之一.1989年，Dennis和Hudson(1989)针对Navier-Stokes方程提出了空间四阶的紧致格式，能够使用粗网格获得较高的精度.1992年，Lele(1992)构造了求解一阶导数和二阶导数的紧致差分格式.Adams 和Shariff(1996)提出了紧致ENO格式，用于求解激波湍流相互作用问题.Chu和Fan(1998)提出了三点组合型紧致差分格式，并将其用于求解对流扩散方程.随后人们针对不同的问题和方程，提出了许多种不同的紧致差分格式，广泛用于空气动力学、流体力学和电磁波方程等方面.在地震波场数值模拟方面，王书强等(2002)和Zhou和Greenhalgh(2003)先后将Lele(1992)提出的紧致差分格式用于弹性波方程的数值模拟中.Yang等(2003)年给出各向异性情况下的紧致差分方法以及声波和弹性波数值模拟结果，并给出相应的吸收边界条件.Du等(2009a,b)利用紧致交错网格差分方法对横向各向同性介质地震波场进行了数值模拟，Liu和Sen(2009)研究了任意阶空间导数的紧致格式和一阶空间导数的交错紧致格式，分析了它们的精度和计算效率，并将其用于声波和弹性波数值模拟中.Kosloff等(2010)提出可以利用当前点两边各任意多个点计算导数值的一般紧致格式，并基于Remes方法求取差分系数.杨宽德等(2011)研究了Biot流动和喷射流动耦合作用下波传播的FCT紧致差分模拟方法.此外，Chu和Stoffa(2010)、Liu(2014)还研究了频率域地震波数值模拟中的紧致差分方法.而组合型紧致差分格式在地震波场模拟中的应用，尚未见到文献报道.本文在声波和弹性波方程时间四阶离散格式的基础上，将组合型紧致格式应用到位移场空间偏导数的求取，实现了二维各向同性介质地震波场的数值模拟.1 组合型紧致有限差分方法原理早期的紧致差分格式是基于Hermite多项式构造而来的，在Hermite多项式中，相邻三个节点的函数值、一阶和二阶导数值关系可以表示为其中fi表示节点函数值，和分别表示一阶和二阶导数值.1992年，Lele(1992)对Hermite公式进行了扩展，构造了求解一阶导数和二阶导数的紧致差分格式，表示为(2)其中a,b,c,α,β为差分系数，h表示网格间距.对上述Lele差分格式进行泰勒公式展开和求解差分系数，可以得到一阶导数五点中心六阶精度差分格式：(3)二阶导数五点中心六阶精度差分格式：(4)由公式(3)和(4)可以看出，紧致差分(Compact Difference，以下简称CD)格式的特点就是用相邻节点的函数值求解若干节点上的导数值.而常规中心差分(Central Finite Difference，以下简称CFD)仅求解中心节点的导数值.另一方面，也可以认为常规中心差分是紧致差分的特例，即公式(2)中的α=β=0，例如一阶和二阶导数的六阶精度差分格式表示为(5)其中b1=0.75, b2=-0.15, b3=0.0167, c0=-2.7222, c1=1.5, c2=-0.15,c3=0.0111.常规中心差分如果要得到六阶精度的差分格式，需要用到七个节点，而紧致差分只需要五个点，这使得紧致差分格式比常规差分处理边界节点上更为方便.且一般情况下，在相同网格间距时，紧致差分格式比常规差分具有更高的精度，具有更小的数值频散(王书强等，2002).Chu和Fan(1998)等构造了精度更高的三点六阶组合型紧致差分格式(Combined Compact Difference，以下简称CCD)：(6)与公式(3)和(4)对比可以看出，上述组合型紧致差分格式只需要相邻的三个节点就可以同时求得一阶和二阶导数的六阶精度近似值，比常规紧致差分的节点数更少.公式(6)中的一阶导数和二阶导数是耦合的，既可以同时求出也有利于波形的保真性.2 声波方程组合型紧致有限差分方法分析二维声波方程可以表示为(7)式中u(x,z)为地震波位移，v(x,z)为地震波速度.利用截断的泰勒公式表示n+1和n-1时刻的位移场可以得到：(9)两式相加，略去高次项，并代入声波方程，得到位移场时间四阶精度的差分格式：(10)公式(10)为位移场的三层显式差分格式，利用它就可以进行地震波场时间层的推进.在公式中含有位移对x和z的二阶和四阶导数，对这些导数将利用组合型紧致差分格式(6)进行求取.假设模型纵向和横向节点数为m和n，h为空间网格大小.利用公式(6)求公式(10)中的空间偏导数∂2u/∂z2和∂2u/∂x2，矩阵表示为AE=BU， FC=UD,(11)其中A和C为公式(6)左端的差分系数矩阵，大小分别为2m×2m和2n×2n；E和F为待求位移场空间一阶和二阶导数矩阵，大小分别为2m×n和m×2n；B和D为公式(6)右端的差分系数矩阵，大小分别为2m×m和n×2n；U为位移场矩阵，大小为m×n.这些矩阵分别表示为(12)(13)(15)(16)(17)(18)公式(11)可以用追赶法进行求解.由公式得出位移场ui,j可同时求得它的空间一阶和二阶导数∂ui,j/∂x、∂ui,j/∂z、∂2ui,j/∂x2和∂2ui,j/∂z2，表示为E=A-1BU， F=UDC-1,(19)差分公式(10)中的空间四阶导数可以在求得二阶导数后，将二阶导数值当作U，再次使用公式(19)进行求取.2.1 精度分析不论是利用CCD格式还是利用常规的七点六阶CFD或五点六阶CD格式，进行声波方程数值模拟时，它们在时间层推进方式上是相同的，即都是利用差分格式(10)，时间差分为四阶精度.CD、CFD和CCD不同之处在于它们分别利用公式(4)、(5)和(6)计算位移空间高阶导数.虽然在声波方程的差分格式中并没有空间一阶偏导数，但在下文中的第4节进行弹性波场数值模拟中，存在空间的混合偏导数，需要通过一阶导数的迭代求取，所以这里对这三种格式计算的一阶和二阶导数的截断误差进行对比，结果见表1，α、β、a、b和c为公式(2)中的差分系数.从表1可以看出，虽然三种方法都能达到空间六阶精度，但截断误差有较大的差别.CFD和CD差分格式计算一阶偏导数的截断误差约是CCD的40倍和4.4倍，计算二阶偏导数的截断误差约是CCD的36倍和8.5倍，这说明CCD方法的具有更高的差分精度.表1 一阶和二阶导数不同差分方法的截断误差比较Table 1 Truncation errors in various difference schemes for the first and second derivative calculationsαβabc截断误差uxCFD003/2-3/51/1036/(7!)×f(7)h6CD1/3014/91/904/(7!)×f(7)h6CCD/≈-0.9/(7!)×f(7)h62ux2CFD003/2-3/51/1072/(8!)×f(8)h6CD2/11012/113/110≈-17/(8!)×f(8)h6CCD/-2/(8!)×f(8)h62.2 误差分析通过模拟二维平面谐波初值问题并计算模拟误差，定量分析CCD方法的数值模拟精度.二维平面波初值问题可以表示为(20)其中v是平面波的波速，θ 0是初始时刻波阵面法线方向(即传播方向)与x轴的夹角，f0是平面简谐波的频率.上述初值问题的解析解为(21)在二维波场数值模拟中，假设f0=20 Hz，θ0=π/4，均匀介质的波速为3600 m·s-1，模型长度和深度均为2000 m，纵横向网格长度相同，采样时间为1 s.在不同的空间网格长度和时间步长条件下，计算数值模拟的相对误差.相对误差定义为(22)式中表示数值解，u(tn,xi,zj)表示解析解.研究表明(Yang et al.，2014)，优化的近似解析离散化方法(Optimal Nearly-Analytic Discrete，以下简称ONAD)也具有较高的数值模拟精度和计算效率，为了比较，在此计算CCD、CD、ONAD和CFD四种方法在不同模拟参数情况下的相对误差曲线见图1.从图1可以看出，随着空间网格长度和时间步长的增加，四种方法的相对误差均会逐渐增加，并且随模拟时间的增加而增加.便于比较，将实验结果列于表2.从表中可以看出，四种情况下，误差相对关系基本一致，例如当空间步长为20 m，时间步长1 ms时，CCD算法的相对误差仅为0.066%，CD次之，为0.089%，CFD相对误差为0.943%，这与精度分析的结论是一致的.需要说明的是，本次实验的ONAD算法相对误差比CFD方法大，这是由于ONAD只有空间四阶精度，而本文的CFD空间精度为六阶.当采用较细网格时，如空间步长10 m，时间步长0.5 ms时，四阶ONAD方法的误差要小于六阶CFD方法，这也说明了ONAD方法在提高模拟精度方面也有积极的意义.用较小的空间步长(10 m)和时间步长(0.5 ms)时，CCD算法的精度显著提高，相对误差仅有0.0008%，并且相对误差增长缓慢，这说明了采用CCD格式的声波模拟结果精度较高，能进行较长时间的地震波场数值模拟.图1 四种数值模拟算法的相对误差曲线(a) CCD; (b) CD; (c) ONAD; (d) CFD.Fig.1 Relative error-time curves of the numerical simulation using four different methods and models(a) CCD; (b) CD; (c) ONAD; (d) CFD.表2 四种方法在不同情况下的最大相对误差(%)比较Table 2 Relative errors in four different schemes and parameters模拟参数CCDCDONADCFDdx=10 m,Δt=0.5ms0.00080.00160.1190.596dx=15 m,Δt =0.5 ms0.0120.0150.690.62dx =15m,Δt =1 ms0.0080.0190.5560.629dx =20 m,Δt =1 ms0.0660.0892.160.943 2.3 频散分析在数值模拟时，如果空间网格长度使用过大，会造成较大的求解误差，并会产生所谓的数值频散现象，所以频散关系分析既是判断数值模拟方法优劣的重要方法，也是确定空间网格大小的重要依据.通过对声波方程的组合型紧致差分格式进行数值频散分析，进一步分析该方法的适用条件.首先考虑x方向，令简谐波解=Aexp×[I(ihkx+jhkz-kv′nΔt)]，且：(23)如果数值模拟时不存在误差，则A,B,C应该满足以下关系：B=(Ikx)A, C=-(kx)2A,(24)但在差分的实际计算中，所产生的数值波数(又称修正波数)与解析波数不同.令：(25)将公式(23)和(25)代入差分格式(6)，可以得到方程组：(26)其中：解上述方程组可得：(27)同理z方向可得：(28)式中φx=φcosθ，φz=φsinθ，且φ=kh.θ是波的传播方向与x轴的夹角，k为波数，h为空间步长.将公式(23)、(27)和(28)代入声波方程的差分格式(10)中，并令x=kv′Δt，α=vΔt/h.其中v′为数值波速，v是实际波速，Δt是时间步长，α称为Courant数，利用欧拉公式，化简可得CCD频散关系：]]/12.(29)通过上述频散关系确定φ后可以解得对应的x.定义数值波速与真速度的比值为(30)在压制数值频散方面，优化的近似解析离散化方法(ONAD)方法也具有较好的效果，所以选取CFD、CD和ONAD 三种方法与之进行比较.在理想情况下，如果不存在数值频散则速度比γ恒等于1.γ偏离1越大，则说明该方法的数值频散越严重，反之则说明该方法能更好地压制数值频散.图2显示了以上四种方法在不同的α和θ的时候的速度比曲线.取φ∈[0,π]作为横坐标，它是波数与空间步长的乘积，单位波长内采样点数N=2π/φ，所以横坐标也可以看作N由逐渐减小至2.从图中的速度比曲线可以看出：①随着空间采样点数的减小，四种方法的频散现象逐步加剧.CCD、CD和ONAD方法的数值频散比CFD方法要小，其频散曲线更趋近于1，这也说明了紧致差分(CCD和CD)和近似解析离散类(ONAD)方法在压制数值频散方面的都具有较大的优势，都适合使用较大的空间步长，进而提高计算效率；②假设数值速度在理论速度的99.9%～100.1%以内表示不存在频散，则CCD、CD、ONAD和CFD 对应的φ的极小值分别为1.12、0.98、1.06和0.91，所以它们在最小主波长内需要使用的网格点数分别为5.6、6.4、5.9和6.9个，即CCD和ONAD方法对于空间网格长度的要求最低，其次是CD方法，CFD方法最为严格；③CCD方法的空间网格长度过大时，速度比曲线上翘，数值速度大于真速度，这与另外三种方法相反.这一结论也能从图3c中可以看出，当φ=π，即一个周期内取2个样点，蓝色的速度比曲线大部分要大于1.图2 不同差分方法的速度比曲线(a) θ=0,α=0.25; (b) θ=π/8,α=0.25; (c)θ=π/4,α=0.25; (d) θ=0,α=0.45; (e) θ=π/8,α=0.45; (f) θ=π/4,α=0.45.Fig.2 Velocity ratio curves of the numerical simulation in different methods图3 不同差分方法的方位角-速度比曲线(a) CFD; (b) CD; (c) CCD; (d)ONAD.Fig.3 Azimuth-speed ratio curves of the numerical simulation in different methods图4 空间网格32 m，时间步长1 ms，420 ms时刻波场快照(a) CFD ; (b) CD; (c) CCD; (d) ONAD.Fig.4 Snapshots at 440ms: Δx=32 m，Δt=1 ms为研究该差分方法的数值各向异性特征，即速度比与传播方向的关系，绘制极坐标速度比曲线，如图3所示.图中的极角表示地震波传播方向与x轴的夹角，极径表示速度比.四种方法在空间采样(φ=π/3)满足频散要求时，都没有各向异性特征.在空间采样不足时(φ=7π/5，φ=9π/5)均会出现速度各向异性，0°(或90°)与45°方向上的速度差别最大.利用二维模型来说明CCD方法在压制数值频散方面的优势.模型长度和深度均为4 km，纵波速度3600 m·s-1，震源子波为f(t)=sin(2πf0t)exp×主频20 Hz.根据频散分析结论，CCD方法在最小主波长内只需要使用约5.6个点，所以设置空间网格间距32 m，时间步长1 ms.图4为四种方法数值模拟的波场快照，可以看出CCD和ONAD方法的结果清晰，没有频散现象，CD方法在0°和90°方向上存在微弱的频散现象，而CFD的频散非常严重，与前面理论分析结果一致.2.4 稳定性分析稳定性条件也是有限差分数值模拟中一个非常重要的问题，是影响差分方法计算效率的重要因素.采用Fourier方法(张文生，2006)对声波方程的组合型紧致差分格式进行稳定性分析.定义：(31)其中为空间网格大小，kx和kz为x和z方向的视波数，i、j和n分别为空间和时间网格坐标，代入公式(6)中可得：(32)令θ1=kxh，-π/2≤θ1≤π/2，并利用欧拉公式，将上式化简可得：(33)解方程可得：=Aξn/h2,(34)其中(35)同理，定义：=eIh(kxi+kzj), =eIh[kxi+kz(j-1)],=eIh[kxi+kz(j+1)], =eIh(kxi+kzj),=eIh[kxi+kz(j-1)],=eIh[kxi+kz(j+1)],且令θ2=kzh， -π/2≤θ2≤π/2，代入公式(6)可以得到：=B ξn/h2,(37)其中(38)公式(10)是一个三层显式差分格式，为了分析其增长矩阵，将其改写为：(39)式中=φneIh(kxi+kzj).将公式(34、35、37和38)代入(∂2u/∂x2=exp[Ih(kxi+kzj)]和(∂2u/∂z2=exp[Ih(kxi+kzj)]，并代入声波差分格式(39)中，化简并写成矩阵格式为(40)其中(41)为增长矩阵，其中α=vmaxΔt/h为该差分格式的Courant数.为使差分格式(40)满足稳定性条件，需要其增长矩阵的谱半径小于等于1.直接求解该增长矩阵的特征值较为复杂，用二分法求该问题的解(宋国杰，2008)，计算得到的近似稳定性条件为α=vmaxΔt/h<0.792.2.5 计算效率分析采用2.3节所用模型进行数值模拟，采样时间1 s，分析CCD、CFD、CD和ONAD四种方法的计算效率，各方法使用的数组大小、个数和算法时间复杂度等列于表3.表3 不同方法计算效率比较Table 3 Comparison of computational efficiency between different methodsCCDCFDCDONAD网格半径3753每个波长采样点5.66.96.45.9Δx=Δz/(m)32262831数组/个1281225Courant0.790.980.930.527Δt/ms7774时间复杂度O(n3)O(n3)O(n2)O(n2)时间/s0.6670.8681.041.79从表3中可以看出：CCD和ONAD差分方法的网格半径最少，仅需要3个点，这对于边界处理最有利，但ONAD空间精度仅为四阶，要低于CCD的六阶精度.CFD和CD差分方法虽然空间精度能达到六阶，但网格半径也随之增加.CCD和ONAD均为低数值频散的差分方法，满足无数值频散要求的空间网格长度大于CFD和CD方法.但是ONAD方法需要计算位移场的梯度，这就增加了数组的个数和计算时间.同样的声波方程，CCD仅需要12个数组，而ONAD需要25个，在网格长度相近的情况下，CCD方法占用内存最小.从稳定性来看，要达到同样的时间四阶精度，CFD和CD方法的Courant数最大，CCD次之，ONAD最小.考虑空间步长的不同，CCD、CFD和CD能允许较大的时间步长(7 ms).从算法的时间复杂度来看，CCD和CD算法的复杂度要大于CFD和ONAD，其中n表示网格点数.值得注意的是，在模型大小相同的时候，由于网格长度的差异，会造成四种算法的空间网格节点数n不同，也就造成了实际所需运算时间的不同，表中运算时间是在同一台计算机和相同编程环境下模拟得到的.综合模拟精度、空间网格长度、时间步长和占用内存四个因素来看，CCD方法进行声波方程数值模拟时，它具有高精度(时间四阶、空间六阶)、高计算效率(较大Courant数、占用内存少)和低数值频散(粗网格)的优势.需要说明的是，表3中四种方法各参数的选取时，均不考虑计算精度，空间网格和时间步长仅满足无数值频散和稳定性的要求.但在实际计算时，为了精度考虑，二者均需要按要求适当减小，这势必会进一步增加计算量和存储量.3 弹性波方程组合型紧致有限差分方法各向同性完全弹性介质中的二维波动方程可以表示为(43)其中，u(x,z)和w(x,z)表示x和z方向的质点位移，VP(x,z)和VS(x,z)表示纵横波波速.同声波方程一样，利用公式(8)和(9)，并代入波动方程，将u对时间的导数转换为u和w对空间的导数，可以得到时间四阶精度差分格式：(44)同理可以得到w的时间四阶精度差分格式为(45)弹性波差分格式与声波方程不同，公式(44)和(45)中除了位移u的空间二阶导数∂2ui,j/∂x2和∂2ui,j/∂z2外，还有二阶混合偏导数∂2ui,j/∂x∂z和其他四阶偏导数.假设偏导数∂2/∂x2，∂2/∂z2，∂/∂x和∂/∂z的差分算子为和则二阶混合偏导数∂2/∂x∂z的差分算子为四阶偏导数∂4/∂x4的差分算子表示为即对于二阶混合偏导数，可以在求出∂ui,j/∂x和∂ui,j/∂z之后，再次使用公式(19)，对另一个方向求导，从而求出空间二阶混合偏导数∂2ui,j/∂x∂z.同样，对差分公式中出现的u和w的其他高阶偏导数，也可以按照上述同样的方法进行迭代求取，这里不再赘述.4 模型试算为了验证方法的数值性能，我们采用组合型紧致方法对公式(43)和(44)表示的二维各向同性完全弹性介质的弹性波方程进行数值模拟.4.1 均匀介质模型模型长度和深度均为10 km，介质为泊松体，纵波速度3460 m·s-1，横波速度2000 m·s-1，在模型中心处激发10 Hz主频的雷克子波.由于最小速度为2000 m·s-1，按频散分析的结果，每个波长内需采样5.6个点，所以设置网格长度为35 m.同时，由于最大速度为3460 m·s-1，按稳定性条件，取时间步长为8 ms.长时程数值模拟是数值模拟中的一个重要内容(Chen,2006; Li et al.,2012; Gao et al.,2016).对于大尺度模型和长时程传播，即便是十分微弱的误差也会随着传播时间的增加而不断累积，最终可能导致不容忽视的数值假象.为了分析CCD方法的长时程数值模拟结果，数值模拟时没有使用吸收边界条件，图5(a—d)显示的是1 s、2 s、5 s和10 s时刻的位移水平分量波场快照，地震波场特征清晰，无可见的数值频散现象.为了比较，对该模型采用相同的参数，使用时间四阶空间六阶精度CFD格式进行数值模拟，图5(e—h)为相同时刻的波场快照.可以看出，在满足频散要求和稳定性条件的情况下，CCD方法具有计算效率较高、长时程稳定的特点，能适应较长时间的数值模拟，计算结果较为理想.在相同参数下，CFD结果存在很明显的频散现象，并且随时间的增加，频散更加明显，严重地影响了数值模拟的效率和精度.图5 均匀介质模型弹性波数值模拟波场快照(a) CCD(1 s); (b) CCD(2 s); (c)CCD(5 s); (d) CCD(10 s); (e) CFD(1 s); (f) CFD(2 s); (g) CFD(5 s); (h) CFD(10 s).Fig.5 Wave field snapshots using elastic wave equation in homogeneousmedium(a—d) are the snapshots at 1,2,5,10 s for the displacement component in X direction using CCD; (e—h) are the snapshots at 1,2,5,10 s for the displacement component in X direction using CFD.4.2 水平层状介质模型设置水平层状介质模型，共20层，每层厚度400 m，模型总厚度8 km，模型长度同样是8 km.第一层纵波速度为3500 m·s-1，往下各层纵波速度比上一层增加100 m·s-1.模型为泊松体，即纵横波波速比为1.732，所以模型中最小和最大波速为2020 m·s-1和5400 m·s-1.在地面中间处(X=4 km，Z=0 m)激发10 Hz雷克子波震源，人工边界采用PML边界条件(刘有山等，2012)处理.按前文分析结果，设置纵横向网格长度为36 m，时间步长为5 ms.图6显示的是地面接收到的地震记录，记录长度8 s，由于直达波能量明显强于反射波，所以图中对地震记录进行了增益处理.多层模型中地面模拟接收到的地震记录非常清晰，没有数值频散和不稳定数值结果，地震记录中的直达波、反射波和转换波显示清楚，说明组合型紧致有限差分算法可以有效地模拟弹性波在多层各向同性介质模型中的传播.为了检验本文算法长时程数值模拟的结果，图7给出在较粗网格条件下(网格长度36 m，时间步长5 ms)，CCD和CFD方法分别得到的检波点(X=4 km，Z=0 m)处的波形图，并与利用CFD方法在较细的网格(网格长度10 m，时间步长1 ms)情况下得到的近似理论解的参考波形进行对比.图7a显示三种情况下得到的反射波到达时间一致，从图7(b，c)的局部放大结果来看，在粗网格和长时程模拟条件下，CCD方法得到的地震记录平滑无数值频散，与之前的分析结果一致，且CCD方法和近似解析解的波形曲线基本重合，这说明了CCD方法的正确性.而CFD方法在相同条件下存在较严重的数值频散，且模拟时间越长，模拟结果与近似解偏差越大.这些结果和对比说明，CCD方法能够长时间计算稳定性，而且也能很好地压制数值频散.5 结论与建议有限差分法是地震波场数值模拟的重要方法之一，在目前常规有限差分数值模拟中为了获得更高的模拟精度，需要使用更多的节点，这样会增加计算量和所需的存储空间，还会增加人工边界处理的难度.本文从组合型紧致差分格式出发，建立了时间四阶、空间六阶的声波方程的差分格式，对该差分格式进行了模拟精度、频散曲线和稳定性分析，并与常用的几种格式进行了综合比较，该格式可同时计算空间一阶和二阶导数，且具有使用节点少(3个)，低数值频散(每个波长采样5.6个点)，较高稳定性(Courant数为0.792)和模拟精度(空间6阶)的特点，适用于大尺度和较长时程地震波场数值模拟.文中还建立了各向同性介质弹性波方程的差分格式，并进行了数值模拟，模拟结果显示地震波场特征清晰，说明了该方法的适用性.该方法不仅可以用于弹性介质的波场模拟，还可以进一步推广到二维或三维的各向异性介质、黏弹介质和双相介质等复杂介质的数值模拟中，为今后研究地震波传播规律、逆时偏移、全波形反演等工作提供了一种有效的方法.组合型紧致差分格式仅需使用三个节点，就能使一阶和二阶空间偏导数的差分精度达到传统差分方法的六阶精度，因而能够大幅度节约内存需求和计算量，对于三维大尺度模型的正演和反演等问题均具有重要意义.图6 水平层状介质模型弹性波方程数值模拟地面地震记录(a) X方向位移分量； (b) Z方向位移分量.Fig.6 Seismic records of the model′s surface using elastic wave equation in horizontal layers media(a) Displacement component in X direction； (b) Displacement component in Z direction.图7 接收点(4 km，0 km)处接收到的反射波水平位移波形图(a) 整道反射波记录对比； (b) 4000～4500 ms局部放大对比； (c) 6000～6500 ms局部放大对比.其中，红色加点实线是粗网格条件下的CCD数值结果，绿色实线是粗网格条件下。
二维地震正演模拟方法技术研究
二维地震正演模拟方法技术研究一、本文概述随着地球物理学的深入发展和油气勘探的不断推进，二维地震正演模拟方法技术在地震勘探领域的应用越来越广泛。
该技术通过模拟地震波在地下介质中的传播过程，为地震资料解释、储层预测和油气勘探提供重要的理论支撑和实践指导。
本文旨在深入研究二维地震正演模拟方法技术，探讨其基本原理、发展历程以及当前的研究热点和难点，为进一步提高地震勘探的精度和效率提供理论支持和技术保障。
本文将对二维地震正演模拟方法技术的基本概念进行阐述，包括其定义、特点以及应用领域等。
接着，回顾二维地震正演模拟方法技术的发展历程，分析其在不同阶段的主要特点和优缺点。
在此基础上，重点探讨当前二维地震正演模拟方法技术面临的主要挑战和难点，如复杂地质条件下的模拟精度问题、大规模计算的效率问题等。
针对这些挑战和难点，本文将进一步分析现有的解决方案和发展趋势，如基于高性能计算的并行化技术、基于人工智能的反演方法等。
同时，结合具体的应用案例，分析二维地震正演模拟方法技术在油气勘探、矿产资源调查等领域的实际应用效果，以验证其有效性和可靠性。
本文将对二维地震正演模拟方法技术的未来发展进行展望，提出可能的研究方向和应用前景。
通过本文的研究，旨在为推动二维地震正演模拟方法技术的发展和应用提供有益的参考和借鉴。
二、二维地震正演模拟理论基础二维地震正演模拟是地球物理学中一种重要的方法，其理论基础主要基于波动方程和地震波的传播原理。
在二维空间中，地震波的传播受到介质速度、密度、弹性等因素的影响，这些因素决定了波场的空间分布和时间变化。
理解和应用波动方程是二维地震正演模拟的关键。
波动方程是描述波在介质中传播的基本方程，对于地震波而言，常用的波动方程有弹性波方程和声波方程。
在二维正演模拟中，我们通常采用声波方程，因为它相对简单且能够较好地模拟地震波的主要特征。
声波方程描述了声波在弹性介质中的传播规律，包括波速、振幅、相位等参数的变化。
二维波动方程 python
二维波动方程 python
二维波动方程是描述二维空间中波动现象的数学模型。
在数值
计算中，我们可以使用Python来模拟和求解二维波动方程。
在Python中，我们可以使用科学计算库如NumPy和SciPy来进行数值
计算，也可以使用可视化库如Matplotlib来展示结果。
首先，二维波动方程的一般形式可以写成：
∂^2u/∂t^2 = c^2(∂^2u/∂x^2 + ∂^2u/∂y^2)。
其中，u是波函数，t是时间，x和y分别是空间中的两个维度，c是波速。
为了求解这个方程，我们可以使用有限差分法（finite difference method）来离散化偏微分方程。
具体来说，我们可以将
空间和时间分别离散化，然后通过迭代求解来模拟波动的演化过程。
在Python中，我们可以使用NumPy来创建表示波函数的二维数组，并使用SciPy来进行数值计算。
例如，我们可以使用SciPy中
的solve_ivp函数来求解偏微分方程的数值解。
另外，我们还可以
使用Matplotlib来可视化波函数随时间和空间的演化过程。
总的来说，使用Python来求解二维波动方程涉及到数值计算和可视化两个方面，需要结合使用NumPy、SciPy和Matplotlib等库来完成。
当然，具体的实现方式会根据具体的问题和求解方法而有所不同。
希望这个回答能够帮助到你理解如何在Python中处理二维波动方程。
各向异性介质中二阶弹性波方程模拟PML吸收边界条件
弹性波散射现象的数值模拟与分析
弹性波散射现象的数值模拟与分析引言弹性波散射是指当弹性波在遇到不同介质或物体边界时发生的反射、折射和散射现象。
研究弹性波散射现象对于地震勘探、地质灾害预测以及材料科学等领域具有重要意义。
本文将介绍弹性波散射现象的数值模拟与分析方法，以及其在实际应用中的意义和挑战。
一、数值模拟方法1. 有限差分法有限差分法是一种常用的数值模拟方法，适用于求解弹性波方程。
该方法将连续的空间和时间离散化，通过有限差分近似来求解偏微分方程。
有限差分法简单易行，适用于各种边界条件和复杂介质情况。
然而，由于网格剖分的限制，有限差分法对于大尺度、高频率的问题计算量较大。
2. 有限元法有限元法是一种基于离散化方法的数值模拟方法，适用于求解各种复杂边界条件和非均匀介质情况下的弹性波散射问题。
该方法将连续的物理域分割成有限个小单元，通过插值函数和基函数来逼近解的形式。
有限元法具有较高的计算精度和灵活性，但对于大规模问题的计算量较大。
3. 边界元法边界元法是一种基于边界积分方程的数值模拟方法，适用于求解边界上的弹性波散射问题。
该方法通过将边界上的积分方程离散化，将问题转化为求解线性方程组的形式。
边界元法适用于各种复杂边界条件和介质情况，具有高效的计算速度和较小的存储需求。
二、数值模拟与实际应用1. 地震勘探地震勘探是一种通过观测地震波在地下传播和散射的信息来获取地下结构和物性的方法。
数值模拟可以帮助预测地震波在地下的传播路径和散射特性，从而指导地震勘探的设计和解释。
通过模拟不同介质和地下结构的散射现象，可以提高地震勘探的效率和准确性。
2. 地质灾害预测地质灾害预测是一种通过分析地下介质和构造的变化来预测地质灾害风险的方法。
数值模拟可以模拟地震波在地下的传播和散射过程，从而帮助预测地质灾害的发生概率和影响范围。
通过模拟不同地质条件下的散射现象，可以提高地质灾害预测的准确性和可靠性。
3. 材料科学材料科学是一门研究材料性质和结构的学科，对于材料的弹性波散射现象的研究具有重要意义。
PML边界条件下二维粘弹性介质波场模拟
的粘 弹性 体 ， 种粘 滞 性 表 现 为 基 质 内部存 在 “ 这 蠕
变 ” “ 弛” “ 、松 、 迟滞 ” 等现 象 。地震 波在 大地 中传 播
时， 质点 的振 动能 量 转 化 为 其 它各 种 形 式 的能 量 ，
目前 较 为普 遍 的说 法 即介 质 内部 的 热 效 应 现 象 。 这种 能量 的转 化使 地震 波在 粘弹 性介 质 中传播 时 ， 其 高频成 分 很 容 易 被 吸 收 ， 幅 近 似 指 数 规 律 衰 振
位 置发 生偏 离 。
的二 阶有 限差 分法 的一种 推 广 ［ 。此 外 ， ５ ］ 同有 限差 分法 和有 限元 法相 比 ， 伪谱 法不仅 克 服 了它们 对 高
频成 分 的限制 ， 以实 现 全 频 带 地震 波 模 拟 ， 且 可 而 节省 了前 两种 方 法所需 要 的大量 内存 。 我们 使 用伪谱 法模 拟粘 弹性 介 质 中的波 场 ， 并 引入 最 佳匹 配 层 （ ＭＬ） 理 边 界 ， 得 原 本 对 空 Ｐ 处 使 间做 的二维傅 里 叶变换 变成 了一维 傅里 叶变 换 ， 大 大 提高 了伪谱 法 的计 算效 率 。通过 模 型试算 ， 析 分 了粘 性 系数 对地 震波 的 吸收 和衰减 。
伪 谱 法 是 对空 间坐 标 通 过快 速 傅 里 叶 变换 在 时间 域作 差分运 算 ， 避免 了求 偏 导数 。由于伪谱 法 可 以看 成 是高 阶 网格 差 分 法 当其 空 间 差分 精 度 的
阶数 达到无 穷 时 的极 限 隋况 ， 故伪 谱 法可视 为 传统
减， 波形 和相位 失 真 ， 得分 辨率 降低 ， 使 反射 界 面 的
Ｐ Ｌ边 界 条 件 下 二 维 粘 弹 性 介 质 波 场 模 拟 Ｍ
大学物理-波动方程
2
谱方法的优点是精度高，适用于大规模问题求解， 且能够处理复杂的边界条件和初值条件。
3
谱方法的缺点是计算量大，需要较高的编程技巧 和计算资源，且对非线性问题的处理较为困难。
06 波动方程在物理中的应用
声波传播
声波传播
波动方程可以描述声波在介质中的传播规律 。通过求解波动方程，可以得到声波的传播 速度、振幅和相位等信息。
有限差分法的优点是简单直观，易于编程实现，适用于规则区域的问题求解。
有限差分法的缺点是对不规则区域和边界条件的处理较为复杂，且精度相对较低。
有限元法
01
有限元法是一种将连续的波动问题离散化为有限个相互连接的子域（即有限元 ）的方法，通过将波动方程转化为有限元方程组，然后求解该方程组得到波动 问题的数值解。
大学物理-波动方程
contents
目录
• 波动方程概述 • 一维波动方程 • 二维波动方程 • 三维波动方程 • 波动方程的数值解法 • 波动方程在物理中的应用
01 波动方程概述
波动方程的定义
波动方程是描述波动现象的基本数学 模型，它描述了波动在空间和时间上 的变化规律。
波动方程通常表示为偏微分方程，其 中包含未知函数（如波动位移或速度 ）及其偏导数。
地震定位与测深
利用地震波的传播规律，可以进行地震定位和测深，以了解地球内 部结构和构造。
地震灾害评估
地震波的传播特性可以为地震灾害评估提供重要信息，如地震烈度、 震源深度和地表破裂带等。
THANKS FOR WATCHING
感谢您的观看
偏微分方程的形式
三维波动方程通常采用偏微分方程的形式，包含了波动传播的空间 和时间信息。
三维波动方程的解法
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
弹性介质中二维波动方程数值模拟
作者：赵英杰王智亮孙雄伟田继东何丹
来源：《西部资源》2016年第06期
摘要：本文对于弹性介质中二维波动方程，介绍了求解地震波方程的ONADM方法，然
后用该方法对二维波动方程在均匀介质和水平层状速度模型进行数值模拟，结果表明，在均匀介质中得到的了清晰、无频散的波场快照，水平层状速度模型的地震记录清晰，直达波和反射波区分明显。
可见，ONADM方法是一种有效的地震波模拟算法，对提高地震勘探精度，提高地震防震减灾具有一定的指导意义。
关键词：波动方程；ONADM；波场快照
4. 结论
本文在采用NAD类算法，在粗网格的条件下，对均匀介质和水平层状介质进行模拟，获得了清晰的波场快照，表明该方法能够较好地压制数值频散，得到清晰的地震波场快照和地震记录。
该方法对我们认识地下介质构造，提高地震勘探精度获得更加准确的地震成像结果具有重要意义。
同时，运用该方法进行地震波场模拟对防震减灾也具有一定的指导意义。
参考文献：
[1] 宋国杰.三维弹性波方程的改进近似解析离散化方法及波场模拟.硕士毕业论文.
[2] 杨桂通，张善元.弹性动力学.中国铁道出版社， 1988.
[3] Yang D H， Teng J M， Zhang Z J， et al. A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media. Bull.Seis.Soc.Am.， 2003，93（2）：882—890.
[4] 牛滨华. 半空间介质与地震波传播[M]. 石油工业出版社， 2002.
[5] 宋亚辉，孙琪皓，聂振邦，等. 裂隙介质下天然气水合物CO2-CH4相平衡实验研究[J]. 西部资源， 2016（1）：25-26.
[6] 孔然，唐喜，谢梦谦. 地震勘探中的低信噪比处理分析[J]. 西部资源， 2016（5）.
致谢：
1.大北区块构造数值模拟研究，中石油集团股份有限公司科技公关重大专项，
2011zx05046—3.
2.41674141，分数阶粘性地震波场的高精度有限差分算法研究，国家自然科学基金面上项目.
3.西南石油大学地震波传播与成像（2015CXTD07）青年科技创新团队.
4.X151514KLX0项目名称：各向异性弹性介质中的拟声波控制方程逼近研究.。