基于二次场的二维大地电磁有限元法数值模拟

二维地震波场的组合型紧致有限差分数值模拟汪勇;石好果;周成尧;桂志先【摘要】地震波场数值模拟在地球物理勘探和地震学中具有重要的支撑作用.本文将组合型紧致差分格式用于声波和弹性波方程的数值模拟中.根据泰勒级数展开和声波方程,建立了位移场时间四阶离散格式,并将组合型紧致差分格式用于位移场空间导数的求取,然后对该差分格式进行了精度分析、误差分析、频散分析和稳定性分析.理论研究结果表明:① 该差分格式为时间四阶、空间六阶精度,与常规七点六阶中心差分和五点六阶紧致差分相比,具有更小的截断误差和更高的模拟精度;② 每个波长仅需要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数值结果，绿色实线是粗网格条件下。

二次插值的时间域激电2.5维有限元数值模拟张志勇;周峰;李泽林【期刊名称】《石油地球物理勘探》【年(卷),期】2015(050)005【摘要】通过研究基于三角单元二次插值的有限元2.5维时间域激电正演,引入Cole-Cole模型推导了2.5维复电阻率满足的边值问题,采用Guptasarma提出的滤波算法,实现了频率域激电到时域的转换;从刚度矩阵的特点出发、结合基于符号分析的线性方程组直接解法,设计实现了适合时间域2.5维激电正演的快速算法.该算法通过对计算区边界的近似处理,使处理后的刚度矩阵只与供电频率及波数有关,而与供电点位置无关,从而有效减少线性方程组直接解法中矩阵分解的工作量;设计了基于图论理论的矩阵重排与填入元分析算法,实现了高效的矩阵LDLT分解;优化了计算流程,在正演计算过程中对于同一剖分结构只需进行一次符号分析,相同频率和波数条件下所有供电点只需一次LDLT分解.最后,利用快速算法,研究了异常体的视电阻率、视极化率、视频散率等参数的异常特征.【总页数】8页(P999-1006)【作者】张志勇;周峰;李泽林【作者单位】中南大学地球科学与信息物理学院,湖南长沙410083;东华理工大学核工程与地球物理学院,江西南昌330013;东华理工大学核工程与地球物理学院,江西南昌330013;东华理工大学核工程与地球物理学院,江西南昌330013;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074【正文语种】中文【中图分类】P631【相关文献】1.基于二次场算法的时间域激电2.5D有限元正演 [J], 周峰;张志勇;蓝泽鸾;林文东2.频率域与时间域激电法的应用对比 [J], 王淦;薛江;黄镇3.频率域-时间域激电对比试验分析 [J], 沈礼锋;单京豪4.基于双二次插值的探地雷达有限元数值模拟 [J], 戴前伟;王洪华;冯德山;陈德鹏5.频率域激电有限元数值模拟 [J], 杨晓弘;何继善;童孝忠因版权原因，仅展示原文概要，查看原文内容请购买。

大地电磁场三维地形影响的矢量有限元数值模拟的报告，600
字
本报告旨在研究大地电磁场分部数值模拟的三维空间影响。

空间中的电磁波是按其位置、强度和传播方向变化的，因此进行空间建模必须要考虑三维影响。

本报告运用矢量有限元数值模拟来建模大地电磁场的三维空间影响，以便了解电磁场的特性。

首先，在已知的三维地形情况下，使用自定义的算法来生成三维矢量网格。

三维网格由节点和单元构成，其中每个节点代表一个位置，而单元由节点连接而成。

然后，使用有限元的技术来求解每个节点的大地电磁场。

该模型在求解过程中充分考虑到了三维地形的影响，模拟出的地形图能够清楚地展示出电磁场在该地形上的变化，进而帮助我们更好地理解电磁场的数学特性。

本研究的结果表明，使用矢量有限元数值模拟来建模大地电磁场在三维空间中的影响是可行和有效的。

在建模过程中，通过对三维地形的准确建模，使得模型能够更为精确地反映电磁场在实际空间的变化，从而更好地加深了人们对电磁场的理解，也促进了电磁场的研究。

综上所述，通过矢量有限元数值模拟来模拟大地电磁场在三维空间中的影响，可以更加准确地模拟出电磁场的变化规律，从而为我们深入了解电磁场提供了新的研究方向。

二维大地电磁反演matlab程序二维大地电磁反演是地球物理学中的一种重要方法，它利用地下电磁场的测量数据来推断地下介质的电导率分布。

在本文中，我们将介绍如何使用MATLAB编写二维大地电磁反演程序，并详细解释程序的实现原理和步骤。

在开始编写程序之前，首先需要明确二维大地电磁反演的基本原理。

在地球物理勘探中，使用了一种称为大地电磁法的方法来探测地下的电导率分布。

这种方法是通过在地面上放置电磁感应线圈，然后通过测量感应线圈接收到的地下电磁信号来推断地下介质的性质。

在这个过程中，我们需要解决一个反问题，即从测量数据推断地下电导率分布。

二维大地电磁反演程序的实现基于有限元方法，它将地下介质划分为许多小的网格单元，并利用有限元法的原理来求解电磁场的分布。

具体而言，程序分为以下几个步骤：1. 网格划分：首先，需要将地下区域划分成一系列的小网格单元。

这些网格单元可以是正方形、矩形或任意形状，根据具体情况选择合适的网格类型。

网格的大小和密度可以根据需要进行调整。

2. 电磁场模拟：在每个网格单元内，根据电磁场的基本方程，利用有限元法来模拟电磁场的分布。

这涉及到求解一个关于电场和磁场的偏微分方程组，可以使用数值方法如有限差分法或有限元法求解。

3. 数据处理：通过模拟计算得到的电磁场分布，可以计算出在感应线圈上的电磁信号。

这些信号可以与实际测量数据进行比较，从而得到测量数据与模拟数据之间的差异。

常用的数据处理方法包括正演计算和反演算法。

4. 反演算法：通过调整地下介质的电导率分布，使模拟数据与实际测量数据之间的差异最小化。

常用的反演算法有最小二乘法、共轭梯度法等。

这些算法可以通过迭代的方式来求解，每次迭代都会更新电导率分布，直到满足反演精度要求为止。

5. 结果展示：最后，将得到的电导率分布结果进行可视化展示。

可以使用各种图形绘制函数，如contour、pcolor等，将电导率分布以图像的形式呈现出来。

这样可以直观地观察到地下电导率的空间分布特征。

海洋可控源电磁法二维有限元正演及反演海洋可控源电磁法（Marine Controlled-Source Electromagnetic，MCSEM）是一种新兴的海洋地球物理勘探技术。

它通过在海洋中设置控制电源，利用电磁场对海底地质进行探测，可以获取地底结构的信息。

本文将从MCSEM的原理入手，介绍MCSEM的二维有限元正演和反演方法，并讨论其应用及优缺点。

一、MCSEM原理MCSEM利用控制电源发出高频电信号，该信号在海洋中传播时，会激发海底地质物体中的电流。

这些电流在海水中会产生电磁场，通过检测电磁场的变化，可以解析出地底物质的电导率、磁导率等物理参数，从而获取地底结构信息。

MCSEM主要有两种控制电源：直流电源和交流电源。

直流电源具有较大的侵入深度和较好的低频响应，适用于大区域的浅层勘探；交流电源具有较好的高频响应，适用于小区域的深部勘探。

控制电源的设置可以根据勘探的需求进行调整。

二、二维有限元正演MCSEM二维有限元正演是指将地下介质分布描述为二维平面的电导率分布，采用有限元理论计算出该模型产生的电磁响应。

MCSEM的二维有限元正演主要包括以下步骤：1.建立数学模型建立地下介质的二维平面电导率模型，将所要研究的地质构造物体分为多个区域。

2.建立有限元网格将地下介质划分为若干小块，每个小块中的电导率均为常数。

在每个小块内部建立一个节点，并通过有限元网格连接所有的节点。

3.设置边界条件和时间离散在模型的边缘处设置边界条件，确定控制电源和检测电极的位置。

对时间进行离散化，并设置时间步长。

4.求解定态矩阵根据有限元法原理，求解定态矩阵，包括系统矩阵、电势向量和电流向量。

5.计算电场和磁场响应根据电场和磁场的计算公式，使用有限元法计算出电场和磁场的响应。

三、二维有限元反演MCSEM二维有限元反演是指利用已知的电磁响应数据，反推出地下介质的电导率分布。

MCSEM的二维有限元反演主要包括以下步骤：1.误差函数的定义确定反演的误差函数，通常采用观测数据与模拟数据之间的二次差值作为误差函数。

电磁场有限元方法
电磁场有限元方法（finite element Method,FEM）是电磁场分析和设计中一种新兴的解析方法，它将电磁场问题看作是一个数学方程组，然后用”有限元”的数值求解方法进行求解。

可以简单的理解电磁有限元方法的原理就是，先将物理场先用几何拼装的对象表示，用有限个节点(Node)和有限个单元(Element)来组合起来，并对每一个单元内的所有量(如场、势等)的作量线性拟合，这样就将复杂的电磁场问题拆分成几何元素相互连接在一起的小片状，甚至可以定义为0维，1维，2维，3维电磁场问题，可以作出相应的对应有限元元素，比如三维空间就有单元四面体和单元六面体，这样子就可以将这些有限元元素拼成一个完整的电磁场，并且在每个单元内使用坐标系，用均匀格点的方法将微分方程数值插值，以达到计算的目的。

因此求解此式的核心就是有限元的概念，它的基本思想就是对一个复杂的模型分割成若干小几何实体，在这些小几何实体上需要求解的量的取值用某种连续的样条函数的插值来表示，给定一族几何实体上的及其边界条件，可以求出各个点上的量的值。

二维地震波场的五点八阶超紧致有限差分数值模拟周诚尧;汪勇;蔡伟祥;桂志先【摘要】首先将迎风机制引入五点八阶超紧致有限差分(CCD8)格式,得到五点七阶迎风超紧致(UCCD7)格式,并对两种格式进行数值频散分析和精度分析;其次建立了二阶声波方程的位移场时间四阶离散格式,将五点CCD8格式和五点UCCD7格式分别应用于位移场空间导数的求取,并分析这两种格式的稳定性条件;最后基于PML 边界条件,将上述两种格式分别应用于声波方程的均匀介质、水平层状介质及Marmousi模型的数值模拟和波场特征分析及对比.研究结果表明:相较于普通紧致差分,五点CCD8格式具有小截断误差、高模拟精度、低数值频散、高稳定性、所需网格点数少的优点;引入迎风机制后,声波方程的五点UCCD7格式稳定性得到进一步提高;模型试算的结果验证了五点CCD8格式适用于复杂介质的数值模拟,模拟精度和计算效率都高.【期刊名称】《石油物探》【年(卷),期】2019(058)002【总页数】11页(P176-186)【关键词】五点八阶超紧致差分;五点七阶迎风超紧致差分;数值模拟;数值频散;稳定性条件【作者】周诚尧;汪勇;蔡伟祥;桂志先【作者单位】油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100;油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100;中石化重庆涪陵页岩气勘探开发有限公司,重庆408014;油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100【正文语种】中文【中图分类】P631地震正演模拟是探索地震波在介质中传播规律的重要手段,基于波动方程的数值模拟方法是一种地震正演模拟方法[1]。

目前波动方程法主要包括伪谱法[2]、有限元法[3]、边界元法[4]、谱元法[5]以及有限差分法。

有限元法和有限差分方法相比,二者精度相同时,有限元法耗时长[6]。

有限差分法因其具有简单灵活、计算效率高以及占用内存小等优势被广泛应用于地震波场数值模拟[7]。

高次插值有限元在大地电磁正演中的应用王培杰【摘要】大地电磁测深正演问题是反演和解释的基础,正演的精度决定了反演的准确性,在常规的有限单元计算中常用二次或三次插值函数来进行计算,较少使用高次插值,其主要原因是高次插值的形函数计算更复杂,而且一般情况下低次插值计算精度能满足要求,且计算量少.推导出了任意阶次插值有限元的形函数计算公式,并编程实现了基于高次插值有限元法的大地电磁测深一维正演,对2个模型进行了计算,并和解析法进行对比,验证了程序的正确性.由该算例应用得知,高次插值有限元能提升计算的精度.【期刊名称】《中州煤炭》【年(卷),期】2018(040)005【总页数】4页(P136-139)【关键词】大地电磁测深;正演;数值模拟;有限元;高次插值【作者】王培杰【作者单位】长江大学油气资源与勘探技术教育部重点实验室,湖北武汉 430100;长江大学地球物理与石油资源学院,湖北武汉 430100【正文语种】中文【中图分类】P631.30 引言大地电磁测深法(Magnetotelluric Sounding，简称MT)是地球物理学中重要的勘探方法，由原苏联学者Tikhonov[1]和Cagniard[2]于20世纪50年代提出，它以岩石的电性差异为基础和前提，利用天然交变电磁场来研究地球电性结构。

它具有低成本、勘探深度大、不受高阻屏蔽、对低阻有较高分辨率的优点，广泛应用于地质构造和资源勘查等方面[3]。

当前，国内外很多学者对MT的正演问题做了大量的研究工作，Coggon[4]首先实现了有限单元法的电磁法正演模拟，William等[5]采用矩形网格剖分实现了MT 二维正演模拟；Rijo[6]引入通用性网格剖分方法，提升了计算精度和速度，使有限单元法成为MT正演的有效工具；Wannamaker等[7]采用矩形单元与三角单元混合剖分，开发了经典的MT二维正演程序PW2D；Franke等[8]采用了非结构化的三角剖分，进一步研究了MT的二维正演模拟。

基于Petri网的二维电磁场有限元参数化建模方法
王曙鸿;邱捷;励庆孚
【期刊名称】《中国电机工程学报》
【年(卷),期】2002(22)2
【摘要】将参数化绘图技术引入电磁场有限元建模中 ,设计人员只需建立一张参数化的图形 ,可方便快速地生成系列化有限元场域图形 ,进行不同设计方案的有限元分析。

采用Petri网描述参数化的有限元场域模型 ,提出了基于Petri网的参数化推理方法 ,该方法具有基于图的参数化模型的优点 ,且与参数化实体建立顺序无关。

基于Petri网的参数化模型解决了参数化设计中的循环约束判断问题。

经实例表明 ,基于Petri网的参数化建模方法方便实用 ,它将成为电磁场有限元分析软件的一个重要功能。

【总页数】5页(P107-111)
【关键词】三维电磁场;有限元;参数化建模;Petri网
【作者】王曙鸿;邱捷;励庆孚
【作者单位】西安交通大学电气工程学院
【正文语种】中文
【中图分类】TM151.1
【相关文献】
1.基于Excel-VBA与CATIA的拱坝有限元参数化建模方法 [J], 王飞
2.基于Petri网IP网QoS策略控制方法建模 [J], 宋婷禹
3.基于层次着色Petri网的网构软件性能建模与仿真分析方法 [J], 徐倩;应时;贾向阳;耿江屹;李琳
4.基于ANSYS的金属软管参数化有限元建模方法 [J], 韩淑洁
5.基于Isight的冲压驱动桥壳参数化有限元建模方法 [J], 刘博林;谢里阳;张娜;罗义建
因版权原因，仅展示原文概要，查看原文内容请购买。