有限差分法模拟地震波场
地震波场的高阶交错网格有限差分模拟
地震波场的高阶交错网格有限差分模拟霍凤斌;李振鹏;徐发;张涛【摘要】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随着地震波动理论在天然地震和油气地震中的应用，以及计算机技术的飞速发展，在现代地震数值模拟领域逐渐形成了有限差分法、有限元法、虚谱法和积分方程法等求解波动方程的方法。
复杂地表地质条件下地震波数值模拟综述
复杂地表地质条件下地震波数值模拟综述
地震波数值模拟是利用计算机技术模拟地震波传播的一种方法。
它可以提供地震波在复杂地表地质条件下传播的精确信息，是地震波分析的重要手段之
地震波数值模拟是建立在地震波方程基础上，通过求解地震波方程，计算地震波沿着地表和地下传播的精确路径和传播特性，实现地震波数值模拟的基本思路。
目前地震波数值模拟的主要方法有有限差分法、有限元法和有限体积法等。
有限差分法是最早应用于地震波数值模拟的方法，通过对地震波方程采用差分近似，构建离散的波动方程组，求解出地震波传播过程中地声速度、动压以及比推力的变化。
有限差分法的求解精度取决于网格步长，但随着网格步长的增加，计算量会急剧增加。
有限元法是近些年来被广泛应用的方法，它将地震波方程的空间域划分为多个有限元单元，并在每个单元内采用局部坐标系，把地震波方程的计算过程分解为局部有限元问题，然后采用数值积分的方法求解，极大地提高了求解效率。
有限体积法是地震波数值模拟领域近年来发展较快的一种方法，它将地震波方程抽象为一种分片计算方法，通过局部更新的方法，更有效地提高了求解的效率。
地震波数值模拟方法已经在复杂地表地质条件下得到了广泛的应用，如地表振动评估、地震波传播特性分析、城市地震波传播模拟等。
未来，地震波数值模拟将会得到更多的应用，从而为科学家们提供更多的可靠信息。
总之，地震波数值模拟是一种研究地震波传播特性的重要手段，在复杂地表地质条件下，它可以提供精确的地震波传播信息，为科学家们提供更多的可靠信息。
未来，地震波数值模拟将会得到更多的应用，从而为地震学及其它领域的研究提供更多的参考。
复杂介质地震波正演模拟方法及优化
复杂介质地震波正演模拟方法及优化摘要本文旨在探讨复杂介质地震波正演模拟方法及其优化。
我们将介绍地震波正演模拟的基本原理，同时介绍目前常用的模拟方法，并针对复杂介质中的挑战提出了一些优化措施。
通过本文的学习，读者将能够更好地理解复杂介质中地震波的正演模拟，并了解如何优化模拟结果。
1.引言地震波正演模拟是地震学中的重要研究方法，通过模拟地震波在地下介质中的传播过程，可以帮助我们解决很多实际问题，如地震勘探、地震灾害预测等。
然而，由于地下介质的复杂性，正演模拟在复杂介质中存在着一些挑战，如速度模型不准确、界面反射等问题。
因此，本文将介绍一些常用的地震波正演模拟方法，并提出一些优化措施，以改善正演模拟结果的准确性和可靠性。
2.地震波正演模拟方法地震波正演模拟方法可以分为有限差分法（F DM）、有限元法（F EM）和谱元法（S EM）等。
下面将逐一介绍它们的基本原理和适用范围。
2.1有限差分法（FD M）有限差分法是一种常用的地震波正演模拟方法，它将介质离散化为网格，通过有限差分的方式，近似求解地震波动方程。
有限差分法简单易行，适用范围广，但在复杂介质中存在一些限制，如对较大的速度变化不敏感。
2.2有限元法（F E M）有限元法是一种基于变分原理的地震波正演模拟方法。
它将介质离散化为小单元，并利用插值函数表示波场的变化。
有限元法相对于有限差分法更加灵活，适用于处理复杂介质中的问题。
然而，有限元法的计算量较大，在大规模模拟中可能存在困难。
2.3谱元法（S E M）谱元法是一种将频率域方法与网格法相结合的地震波正演模拟方法。
它首先利用傅里叶变换将地震波动方程转换为频率域方程，然后在空间域上进行离散化求解。
谱元法具有较高的精度和稳定性，适用于处理复杂介质中的地震波传播问题。
3.优化方法为了改善复杂介质中地震波正演模拟的精度和可靠性，我们提出了以下优化方法：3.1速度模型优化在复杂介质中，速度模型的准确性对地震波正演模拟结果具有重要影响。
地震波传播特性的实验与模拟研究
地震波传播特性的实验与模拟研究地震是由地壳运动引起的地震波传播特性的实验和模拟研究是地震科学中一项重要的研究内容。
通过实验与模拟研究，可以深入了解地震波在地球内部的传播规律和特性，并为地震预测与防灾提供支持和指导。
本文将从实验和模拟两个方面，对地震波传播特性进行研究，以期能为地震科学研究提供一些思路与参考。
一、地震波传播特性的实验研究地震波传播特性的实验研究通常是通过在实验室中模拟地震波的传播过程，并通过仪器设备进行观测和记录来研究。
常见的地震波传播特性实验研究方法有模型实验与震源实验两种。
1. 模型实验模型实验是将地震波传播的物理过程通过模型进行缩放和模拟。
通过建立地质模型和模拟地震源，研究人员可以模拟不同地震波传播路径和地壳结构下的地震波传播特性。
模型实验通常需要借助地震仪、地震计等设备进行观测和数据记录，以获得实验数据。
例如，1989年美国加州Loma Prieta地震后的模型实验研究，研究人员通过在室内搭建地震模型，模拟Loma Prieta地震中的地震波传播过程。
他们通过在模型中注入地震波源，观测不同地震波在模型中的传播速度、幅度衰减和力学效应等特性，研究地震波在地震中的传播规律。
2. 震源实验震源实验是通过实验室中的震源设备产生地震波源，并观测地震波在实验体（如岩石样本）中的传播特性。
这种实验方法可以更好地模拟地震中的震源产生和波传播的真实情况。
例如，1995年日本兵库地震后，研究人员利用震源实验研究了地震波在岩石样本中的传播速度和振幅衰减特性。
他们使用实验室中的震源设备产生地震波源，将岩石样本放置在震源附近，并通过地震仪观测地震波传播过程中的变化。
通过这种实验研究，他们了解到岩石样本中地震波传播速度和振幅衰减与地震中观测到的地震波特性具有一定的相关性。
二、地震波传播特性的模拟研究地震波传播特性的模拟研究是利用计算机模拟方法进行的。
通过建立地震波传播的数学模型和采用数值计算方法，可以模拟地震波在地球内部的传播过程，并预测地震波在不同地震源和地壳结构下的传播特性。
复杂介质中地震多次反射波快速正演模拟
复杂介质中地震多次反射波快速正演模拟复杂介质中地震多次反射波快速正演是地震勘探中的重要技术之一，其主要目的是通过计算复杂地质结构下的地震波传播规律，从而推断地下岩层的物性和构造特征。
在实际勘探中，地震波正演模拟技术具有非常重要的应用价值，可以帮助地震科学家快速、准确地勘探地下地质构造，并为油气勘探、水资源勘探以及矿产资源开发提供重要的技术支持。
在复杂介质中，地震波会经历多次反射、折射和散射，并受到地下介质非均匀性、垂直速度梯度、下穿速度等多种因素的影响。
因此，在进行地震波正演模拟时，需要考虑以上因素的影响，建立相应的数学模型，并通过计算求解得到相应的地震波场数据。
地震波正演模拟的过程可以简单概括为：首先，根据勘探区域的地质结构和介质模型，确定模拟所需的各种参数，如介质密度、弹性模量、泊松比等；然后，建立数学模型，包括方程组及数值计算方法；最后，通过计算求解得到地震波的波形和波场分布，从而推断地下岩层的物性和构造特征。
在地震波正演模拟中，最常用的方法是有限差分法（finite difference method，FDM）。
该方法基于波动方程，使用差分近似表示波场中的各个参数，并通过迭代求解得到波场分布。
FDM方法在计算效率和计算精度方面均有很好的平衡，因此广泛应用于地震波正演模拟中。
除了FDM方法外，地震波正演模拟还常常采用有限元法、谱方法等其他数值计算方法，以便更好地反映复杂介质中地震波的传播规律。
同时，地震波正演模拟还需结合各类物理分析方法，如射线追踪、全波形反演等，以便更好地解释波场分布和岩层结构特征。
总之，复杂介质中地震多次反射波快速正演模拟是地震勘探中的重要技术之一，通过计算复杂地质结构下的地震波传播规律，为勘探工作提供了有力的技术支持。
在未来，随着数值计算方法和物理分析方法的不断发展，地震波正演模拟技术将得到更加广泛和深入的应用。
二维地震波场的组合型紧致有限差分数值模拟
二维地震波场的组合型紧致有限差分数值模拟汪勇;石好果;周成尧;桂志先【摘要】地震波场数值模拟在地球物理勘探和地震学中具有重要的支撑作用.本文将组合型紧致差分格式用于声波和弹性波方程的数值模拟中.根据泰勒级数展开和声波方程,建立了位移场时间四阶离散格式,并将组合型紧致差分格式用于位移场空间导数的求取,然后对该差分格式进行了精度分析、误差分析、频散分析和稳定性分析.理论研究结果表明:① 该差分格式为时间四阶、空间六阶精度,与常规七点六阶中心差分和五点六阶紧致差分相比,具有更小的截断误差和更高的模拟精度;② 每个波长仅需要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数值结果，绿色实线是粗网格条件下。
有限差分法地震正演模拟程序
有限差分法地震正演模拟程序一．二阶公式推导1．二维的弹性波动方程22222222x x x ZU U U U t x z x z λμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 22222222xz z z U U U U t z x x zλμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 2．对方程进行中间差分（1）首先对时间进行二阶差分()()()()()112,,,22222n n n xi k xi k xi kx x x x U U U U t t U t U t t U t t t +--++-+-∂==∂ ()()()()()112,,,22222n n n zi k zi k zi kz z z z U U U U t t U t U t t U t t t +--++-+-∂==∂ （2）对x 方向进行二阶差分()()()()()21,,1,22222n n n xi k xi k xi kx x x x U U U U x x U x U x x U x x x +--++-+-∂==∂ ()()()()()21,,1,22222n n n zi k zi k zi kz z z z U U U U x x U x U x x U x x x +--++-+-∂==∂ （3）对z 方向进行二阶差分()()()()()2,1,,122222n n nxi k xi k xi k x x x x U U U U z z U z U z z U z z z +--++-+-∂==∂ ()()()()()2,1,,122222n n nzi k zi k zi k z z z z U U U U z z U z U z z U z z z +--++-+-∂==∂ （4）对xz 进行差分2,1,11,11,11,11,111,11,11,11,122224n nzi k zi k z zn n n nzi k zi k zi k zi k n n n n zi k zi k zi k zi k U U U U x z x zx z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==2,1,11,11,11,11,111,11,11,11,122224n nxi k xi k x x n n n n xi k xi k xi k xi k n n n n xi k xi k xi k xi k U U U U x z x z x z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==（5）把（1）（2）（3）（4）得到的结果带入波动方程3．写出差分方程()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n n n nxi k xi k xi kxi k xi k xi k nnnnnnnxi k xi k xi k zi k zi k zi k zi k U U U U U U t x U U UU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n nnnzi k zi k zi kzi k zi k zi knn n nnnnzi k zi k zi k xi k xi k xi k xi k U U U U U U t x U UUU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++即得到()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nxi k xi k xi kxi k xi k xi k n n n xi k xi k xi k n n n nzi k zi k zi k zi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nzi k zi k zi k zi k zi k zi k n n n zi k zi k zi k n n n nxi k xi k xi k xi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭二．MATLAB 程序 clear; clc;Nx=200; Nz=300;fi=30;%%%主频t_step=300;%%%%时间采样点dx=10.0;%空间采样间隔——x 方向 dz=10.0;%空间采样间隔——z 方向 dt=0.001;%时间采样间隔——1mslambda=66*1e9;%砂岩拉梅常数lamdamu=44*1e9;%砂岩拉梅常数murho=2650;%砂岩密度%%%%%%模型扩展%%%%%vp=zeros(Nz+2,Nx+2);%纵波速度vs=zeros(Nz+2,Nx+2);%横波速度c=zeros(Nz+2,Nx+2);%交叉项系数for i=1:Nz+2for j=1:Nx+2vp(i,j)=sqrt((lambda+2*mu)/rho);vs(i,j)=sqrt((mu)/rho);c(i,j)=sqrt((lambda+mu)/rho);endend%%%%%% 震源%%%%%%%t_wavelet=[1:t_step]*dt-1.0/fi;source=((1-2*(pi*fi*t_wavelet).^2).*exp(-(pi*fi*t_wavelet).^2));% 雷克子波amp=sqrt(2.0);% 振幅source_x=floor(Nx/2)+1;% 震源位置——x坐标source_z=2;% 震源位置——z坐标source_amp=zeros(Nz+2,Nx+2);% 震源振幅初始化（所有点处均为0）source_amp(source_z,source_x)=amp;% 震源振幅，在位置source_z,source_x处振幅为amp，其它位置为0 %%%%%%%%%%%%%%%%%%%%%%%%%%%U=zeros(Nz,Nx);% 弹性波x分量V=zeros(Nz,Nx);% 弹性波z分量U0=zeros(Nz+2,Nx+2);% 前一时刻的UU1=zeros(Nz+2,Nx+2);% 当前时刻的UU2=zeros(Nz+2,Nx+2);% 下一时刻的UV0=zeros(Nz+2,Nx+2);% 前一时刻的VV1=zeros(Nz+2,Nx+2);% 当前时刻的VV2=zeros(Nz+2,Nx+2);% 下一时刻的Vrecord_u=zeros(t_step,Nx);% x方向接收到的地震记录——Urecord_v=zeros(t_step,Nx);% x方向接收到的地震记录——V%%%%%% 有限差分计算U V %%%%%%for it=1:t_stepfor x=2:Nx+1for z=2:Nz+1U2(z,x)=2*U1(z,x)-U0(z,x)+(dt*dt)*(vp(z,x)^2*(1.0/(dx*dx))*(U1(z,x+1)-2*U1(z,x) +U1(z,x-1))+vs(z,x)^2*(1.0/(dz*dz))*(U1(z+1,x)-2*U1(z,x)+U1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(V1(z+1,x+1)-V1(z+1,x-1)-V1(z-1,x+1)+V1(z-1,x-1)))+source(it)*sour ce_amp(z,x)*dt*dt;V2(z,x)=2*V1(z,x)-V0(z,x)+(dt*dt)*(vs(z,x)^2*(1.0/(dx*dx))*(V1(z,x+1)-2*V1(z,x) +V1(z,x-1))+vp(z,x)^2*(1.0/(dz*dz))*(V1(z+1,x)-2*V1(z,x)+V1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(U1(z+1,x+1)-U1(z+1,x-1)-U1(z-1,x+1)+U1(z-1,x-1)));endendU0=U1;U1=U2;V0=V1;V1=V2;record_u(it,:)=U2(2,[2:201]);record_v(it,:)=V2(2,[2:201]);U=U2(2:301,2:201);V=V2(2:301,2:201);end%%%%% 绘制U 和V 的波场图%%%%%%figureimagesc(U(1:300,1:200));title('U');xlabel('x');ylabel('z');figureimagesc(V(1:300,1:200));title('V');xlabel('x');ylabel('z');%%%%% 绘制U 和V 的地震记录%%%%%%figureimagesc(record_u(1:300,1:200));title('U');xlabel('x');ylabel('t_step');figureimagesc(record_v(1:300,1:200));title('V');xlabel('x');ylabel('t_step');U分量波场图V分量波场图U分量地震记录V分量地震记录。
偏微分方程的有限差分法及地球物理应用
偏微分方程的有限差分法及地球物理应用有限差分法是一种常用的数值求解偏微分方程的方法。
它将连续的偏微分方程转化为离散的差分方程，通过近似求解差分方程，得到偏微分方程的数值解。
这种方法在地球物理学中有着广泛的应用，如地震波传播模拟、电磁场分布计算等领域。
首先，假设我们要研究地震波在地下介质中的传播，可以采用波动方程来描述地震波的传播过程。
波动方程可以写成：∂^2u/∂t^2 = c^2∇^2u其中，u是地震波场，c是地下介质中的波速。
为了用有限差分法求解波动方程，我们需要将连续的空间和时间离散化。
假设我们将空间离散化为网格点（i,j,k），其中i,j,k分别代表空间的x,y,z方向，将时间离散化为时间步长Δt。
对波动方程进行近似，我们可以得到：(u(i,j,k,t+Δt) - 2u(i,j,k,t) + u(i,j,k,t-Δt))/Δt^2 = c^2(u(i+1,j,k,t) + u(i-1,j,k,t) + u(i,j+1,k,t) + u(i,j-1,k,t) +u(i,j,k+1,t) + u(i,j,k-1,t) - 6u(i,j,k,t))/Δx^2将此差分方程应用于地震波传播模拟，我们可以得到地震波场在空间和时间上的离散解。
有限差分法在地球物理中有着广泛的应用。
例如，它可以用于模拟地震波在地下介质中的传播，帮助研究地震灾害的发生机制和地下构造的特征。
通过调整网格的大小和时间步长，可以模拟不同频率的地震波传播过程，从而了解地震波在不同介质中的传播规律。
此外，有限差分法还可以应用于电磁场的计算。
例如，在电磁勘探中，可以利用有限差分法求解麦克斯韦方程，计算电磁场在地下介质中的传播和散射过程。
通过模拟电磁场的分布情况，可以帮助研究地下矿产资源的寻找和勘探。
需要注意的是，有限差分法在应用过程中还需要考虑边界条件的处理。
通常情况下，边界条件是已知的，例如地震波在地表的边界条件可以假设为自由表面，电磁场计算中的边界条件可以假设为电场和磁场的边界条件等。