基于混合网格的有限差分地震波场正演模拟

1、 基于混合网格的有限差分地震波场正演模拟 摘要：为克服单一的交错网格或者规则网格在波场计算中存在的不足，本文提出了一种基于混合网格的时域有限差分算法。该算法在有效计算区域内采用四阶中心差分算子的规则网格有限差分法，与混合网格方法相比较具有较高的计算效率并能节省存储空间；在边界处理中则引入基于交错网格的完全匹配层吸收边界条件，能够有效地抑制波场模拟正演算法中存在的频散问题以及人工边界造成的反射干扰。在理论分析的基础上，我们解决了规则网格与交错网格计算区域边界网格匹配与变量计算的问题，给出了算法实现的详细步骤，最后通过Marmousi模型等模拟算例进行了正演波场模拟，从模拟结果中验证了该算法的有效

2、性。关键字：声波方程；混合网格；有限差分；完全匹配层1 引言现在永久磁体在许多技术中都是重要的组成部分。经过近40年的研究,稀土永磁体经历了第一代稀土永磁体(SmCo5) 、第二代稀土永磁体(Sm2(Co,Fe)17)和第三代稀土永磁体(NdFeB) 。目前,第三代永磁体(NdFeB)的磁能积已经接近其理论极限值。现代高技术(信息与通讯技术、航空航天技术、磁悬浮交通技术、核磁共振技术、纳米技术等)的发展对永磁体的性能与质量提出了更高的要求。2 原理地震波的传播过程是一个复杂的问题。波动方程的组成包括波动函数相对于空间坐标和时间坐标的偏导数计算。虽然弹性波波动方程包含的信息量更充分，使用弹性波方

3、程计算的数据更精确，但是这种方法的计算量非常大，计算消耗昂贵。并且决定场数据的主要因素是P波，可以近似的不考虑S波的影响.因此在实际的应用中，往往应用声波方程来近似的描述地震波的传播规律目前数值模拟中常用的波动方程可分为二阶波动方程和一阶位移-应力波动方程。2.1 规则网格有限差分法无损流体媒质中声波的二阶偏微分支配方程如下18： (1)其中表示在某一时刻上计算区域内二维空间任一点处的位移，表示纵波的传播速度。对方程(1)进行有限差分离散化。2.1.1 波动方程二阶差分近似先取二阶差分精度格式，也即忽略小量和的二阶及其高阶指数项。从而得到离散化的声波方程的二阶差分格式如下： (2)2.1.2

4、波动方程四阶差分近似由于有限差分法的缺点是数值频散问题，因离散化求解波动方程而产生的伪波动，是差分方程所固有的本质特征。有限差分法中的数值频散包括空间频散和时间频散两部分，而空间频散是主要方面19。且数值分析表明20-21，网格频散随网格间距的增大而加剧，而减小网格间距必然会增大计算量。但只要空间网格间距和时间网格间距均满足采样定理，不产生时间和空间假频，随着差分方程阶数的提高，对物性参数改变就越小，进而数值频散越少。所以，可以通过提高差分方程阶数的方法来提高波场模拟精度，压制数值频散。本文用四阶差分逼近声波波动方程(1)式中的空间偏导数，时间仍采用二阶差分精度格式。从而得到声波方程的四阶差分

5、精度离散化结果如下： (3)使用上式(3)递推波场时，在靠近边界的一层使用二阶差分公式(2)计算，从第二层开始用式(3)计算。2.1.3 波动方程六阶差分近似本文用六阶差分逼近声波波动方程(1)式中的空间偏导数，时间仍采用二阶差分精度格式。从而得到声波方程的六阶差分精度离散化结果如下： （4）使用上式(4)递推波场时，在靠近边界的一层使用二阶差分公式(2)计算，从第二层开始用式(3)计算，第三层开始使用式（4）计算。2.2 交错网格有限差分法 考虑均匀，无损流体煤质中的声波传播，其直角坐标系下的一阶位移-应力波动方程如下22-23： （4） （5）其中c为声速，为媒质密度，p为压力，为质点速度

6、矢量。将方程(4)(5)式展开为一组一阶二维偏微分方程 ，形式如下24-26 ： （6） （7） （8）图(1) 二维各向同性介质声波方程交错网格示意图 如图(1)所示，建立了声波的时域交错网格有限差分的计算元胞模型，由图可见压力和速度分量在空间各节点上的分布特点：二维情况时，网格分为一个个的小矩形，两个速度分量分别固定沿不同的坐标方向排布在矩形的边上，每个速度分量被两个压力分量环绕，而压力分量则位于矩形的四个顶角，被各个速度分量所环绕。本文在具体的实现过程中，将两个速度分量设置为中间变量，而压力分量设置为所需求解的位移分量。同样对方程(6)-(8)式进行差分离散化。2.2.1 波动方程二阶差

7、分近似声波偏微分方程(6)-(8)式的二阶精度中心有限差分离散化形式如下： (9) (10) (11)2.2.2 波动方程四阶差分近似声波偏微分方程(6)-(8)式的四阶精度中心有限差分离散化形式如下： (12) (13) (14) 上式中， = =在本文的计算中，媒质密度的值取为1,位移和中间变量的初始值均设为零，利用该初始值按上述公式得到新的变量值来替代前一时刻的值，而新的变量值又是由前一时刻的值得到，用此方式来计算不需要存储波场的整个时刻的值。而对于震源点的处理，因为压力源是一个附加项27，故它可以直接加在(6)式上，对方程式的差分离散化方法同上。2.3 解的稳定性条件2.3.1 规则网

8、格有限差分格式的稳定性条件二阶差分格式的稳定条件为： 四阶差分格式的稳定条件为： 六阶差分格式的稳定条件为： 为波速，为时间步长，为x和z方向的空间网格间隔20 ，此处假定两种空间网格间隔相等。2.3.2 交错网格有限差分格式的稳定性条件二阶差分格式的稳定条件：四阶差分格式的稳定条件： ，分别为x和z方向的空间网格间隔。当网格间隔相等时，仍可同2.3.1中所示，表示为。通过上面几种差分精度的稳定性条件比较可知，随着交错网格高阶差分精度的提高，其稳定性要求仅有稍许提高，故没有必要为提高差分精度而加密差分网格，影响计算的效率28。2.4 边界处理由于计算机容量和速度等运行环境的限制，波场数值模拟计

9、算只能在有限的区域内进行，而这种有限区域的边界除了自由地表外就是人为的边界，称为人工边界。为了模拟开域的波场传播过程，本文在计算区域的截断边界处分别采用了比较典型的Clayton-Engquist边界条件29，Reynolds边界条件30和完全匹配层吸收边界条件13。2.4.1 Clayton-Engquist边界条件Clayton-Engquist边界条件形式如下： (17)其中n为边界的外法线方向，在区域的左边界，右边界，上边界和下边界处的差分格式为： (18)其中，l和d分别表示横向网格数和纵向网格数。2.4.2 Reynolds边界条件Reynolds边界条件形式如下：， (15)可得

10、区域的左边界，右边界和下边界处的差分形式如下： (16)其中，l和d分别表示横向网格数和纵向网格数。2.4.3 完全匹配层吸收边界条件在PML媒质中，假设压力分量P分解为两个子分量Px和Pz,其声波支配方程形式如下31 25： (19)上式中，为PML媒质中的衰减因子，在上文中已经指明将取为1，故在此可直接忽略该参数。使用交错网格的有限差分法对支配方程(19)式进行离散化，结果如下：(20) 在本文中，PML衰减因子变化采用以下函数形式32：，i=0,1,2, 上式中，B为衰减幅度因子，也即最大衰减值，是PML外侧边界处的值，在本文中取为800；为PML层的厚度。3 算法具体实现在本论文中提出

11、采用混和网格时域有限差分法进行地震波场数值模拟，在具体的实现过程中存在规则网格和交错网格重叠边界处的网格匹配问题和变量求取问题，下面将对这些问题给出具体的解决原理与步骤。 3.1 原理如图1所示，在具体的计算中，将整个区域分为五个部分，分别进行计算。区域1由图1中所示的区域LU，L,LD三个部分组成，在区域1中，吸收边界层中的衰减因子参数是不同的，在区域LU和区域LD中，衰减因子，均不为零，而在区域L中，；区域2由图1中所示的区域RU，R，RD三个部分组成，同样，在区域RU和区域RD中，衰减因子，均不为零，而在区域R中，；区域3为图1中所示的区域U，该边界层中，；区域4为图1中所示的区域D，同

12、样，该边界层中，；区域5为图1中所示的计算区域，在该区域中衰减因子全为0。如图2所示，黑色加粗的边界线处即是计算区域的最外层，也是吸收边界层的最内层，黑色矩形的上边界为地表，黑色矩形区域以内是时域有限差分的计算区域，矩形区域以外是PML吸收边界层。红色线边界层设置为炮点所在层。蓝色边界线所构成的矩形边界为计算区域的次外层边界。图2中，边界区域使用的是交错网格，图中显示了位移变量和中间变量的具体交错位置，在计算区域使用规则网格，该区域内只有位移变量。3.2 具体实现步骤先计算时域有限差分计算区域，也即区域5内的下一时刻的位移变量的值,但计算区域的最外层边界，也即图2中黑色矩形边界上的位移变量值并

13、不在此步骤中计算，而是在下面计算吸收边界区域中的位移变量值时求得，也即在该步骤中，只需要计算图2中蓝色边界线上以及边界线以内的下一时刻位移变量值。除炮点所在网格点处的位移变量外，蓝色边界线上的下一时刻所有位移变量的值均由二阶有限差分公式(2)式求得，边界线区域内的下一时刻位移变量的值均由四阶有限差分公式(3)式求得。接着计算图2中蓝色矩形区域以外的下一时刻的所有变量的值，包括下一时刻位移变量和中间变量的值。本文在算法实现时，区域1,2,3,4均是从整个区域的最外层向内层逐层计算出下一时刻的各个变量值。从图2中可以看到，整个区域的最外层边界上只有位移变量，在编程实现中，该层上同一时刻的位移变量值

14、与其次外层上的位移变量值相等，故只要求得其次外层上的位移变量值，便可计算出最外层上的位移变量值。从区域1开始计算，由公式(20)即可求得区域1中所有交错网格点上的下一时刻位移变量值和中间变量值，要注意的是，中间计算区域与PML吸收边界层的交叠边界上的下一时刻位移变量值也是在该步骤中计算求得，也即计算区域的最外层边界上的位移变量是位于交错网格上的。本文是按LU,L,LD的顺序分别计算下一时刻各个变量值的，先分别计算出LU区域中的各个变量值，先利用(20)式中的后两式将交错网格上的下一时刻所有中间变量，的值计算求得，再利用(20)式中的前两式将交错网格上的下一时刻所有位移变量，的值计算求得，再同理

15、计算L，LD区域中的下一时刻中间变量值和位移变量值。然后按同样的方法来分别计算区域2，区域3和区域4。由于本论文中，在地表接收波形记录值，故在吸收边界的上边界层，也即在区域3的计算中，实现波形记录值的获取，并且在区域1,2,3,4的计算过程中，还要利用方程(19)式中的第一个式子计算出计算区域的最外层边界上的下一时刻位移变量，也即的值。用步骤中所求的下一时刻的位移变量，的值和中间变量，的值来更新当前时刻的位移变量，的值和中间变量，的值。交换步骤和步骤中计算求得的时域有限差分计算区域内(包括计算区域最外层边界线上)的不同时刻的位移变量值，也即用当前时刻的位移变量的值来更新上一时刻的位移变量的值，

16、下一时刻的位移变量的值来更新当前时刻的位移变量的值。不断的循环步骤，直至达到所希望的时间步数。 图1 整个区域的划分图图2 整个区域的混合网格示意图4 结果与分析4.1 Marmousi模型Marmousi速度模型横向范围为0-9200米，横坐标网格间距为12.5米，共737网格；纵向深度0-3000米，纵坐标网格间距4米，共750网格，速度分布如图(2)所示。炮点分布于地表，横坐标范围从3000米到8975米，炮间距25米，共240个炮点。检波点分布于各炮点左侧，最大炮检距2575米，最小炮检距200米，道间距25米，共96个检波点，即96道。采样率为4毫秒，记录长度3秒，共750个记录点。

17、以Marmousi模型第一炮的地震数据为例，炮点位于地表3000米处，接收点分布在425-2800米间，接收道的第1个检波点位于425米处，第96个检波点位于2800米，中间94个检波点平均分布于此区域内，如图(3)所示。 图(2) Marmousi速度模型 图(3) Marmousi模型第一炮地震数据 震源选用雷克子波，主频f=20Hz,完全匹配吸收层厚度为20。为满足差分法稳定性条件，时间步长取为0.5毫秒，计算长度为6000个时刻，每8个时刻采样一次。4.1.1 数值模拟结果比较实验结果图如下：图(4) 不加边界条件的地震记录 (a) Clayton-Engquist边界条件 (b) R

18、eynolds边界条件图(5) 加不同边界条件的地震记录 (a) 交错网格二阶有限差分法 (b) 交错网格四阶有限差分法图(6) 交错网格有限差分法 (a) 混合网格二阶有限差分法 (b) 混合网格四阶有限差分法图(7) 混合网格有限差分法图(4)为不加入边界条件，只是将左右边界各外推20道，即左右各延拓40个网格点的结果，图(5)至图(7)均是延拓后的结果。图(5)分别为在计算区域边界处使用Clayton-Engquist边界条件和Reynolds边界条件后的波形记录结果，图(6)-(7)为加入相同完全匹配吸收层边界，但在有效计算区域内采用不同有限差分法的结果。图(6)分别为交错网格二阶差分

19、格式，四阶差分格式的正演结果。图(7)分别为混合网格二阶差分格式，四阶差分格式的正演结果。分析图(4)至图(7)可知，在边界处理中，使用Clayton-Engquist边界条件和Reynolds边界条件都不能有效的减少人工边界造成的反射干扰，但使用PML吸收边界条件，就可以有效的减小人工边界造成的反射干扰，得到较清晰的波纹信息。比较图(6)中的(a)和(b)，图(7)中的(a)和(b)可知，差分精度较低时，会出现明显的频散现象，波纹变得模糊，提高差分精度能有效的抑制频散，使波纹信息更加明显。比较图(6)和图(7)中的(a)，图(6)和图(7)中的(b)可知，虽然在有效的计算区域内，随着差分阶数

20、的增加，交错网格差分的局部精度能迅速提高，但当差分阶数还未达到一定高度的情况下，使用规则网格得到的结果更优于交错网格，从结果图即可看出，图(7)中的频散明显弱于图(6)中相同差分阶数的频散，波纹信息更加清晰。4.1.2 两种网格类型的计算效率比较从计算效率的角度进行分析，对文中的混合网格及交错网格有限差分法的计算时间进行对比，结果如表(1)所示，由表中数据对比可知，在同样的计算条件下（64位的win7系统），使用混合网格有限差分法较之交错网格有限差分法的计算效率更高。有限差分法PML层厚度(网格层数)25101520混合网格(s)6567697175交错网格(s)141159178194223

21、 表(1) 使用混合网格和交错网格计算效率结果4.2 匀速介质模型 采用匀速介质模型的横向和纵向网格数均为50，横向空间步长为10m,纵向空间步长为5m，采用频率为20Hz的Ricker子波激发，模型介质的速度为2000m/s，震源置于模型(30m,0m)处。为了保证整个计算过程稳定，时间步长设定为0.5ms，总计算时间样点数为600。采用不同的吸收层厚度，分别利用本文中的混合网格以及交错网格法进行计算，吸收层厚度所占网格点数依次为1，5，10，15。得到如图(8)至图(9)所示的波场快照结果图，波场快照记录时间为0.12s。 a b c d 图(8) 采用交错网格时不同PML对应的波场快照a

22、 PML=1; b PML=5; c PML=10; d PML=15 a b c d 图(9) 采用混合网格时不同PML对应的波场快照a PML=1; b PML=5; c PML=10; d PML=15 由于文中使用的交错网格有限差分法和混合网格有限差分法的两个工程对数值处理的算法差异，以致两个工程得到的波形记录数值不在同一量级上，故两组的波形记录图会有些差异，但并不影响结论的比较。可分别从各组的比较中得知，采用交错网格有限差分法，如图(8)b所示，在PML=5时，仍得不到较好的波形记录值，PML=10时，如图(8)c所示，波形记录值达到较好的结果；而采用混合网格有限差分法，如图(9)所

23、示，在PML=5时，已经能得到较好的波形记录值。由此可见，我们采用的混合网格有限差分法较交错网格有限差分法能更好的适用于数值模拟。5 结论本文采用了一种基于混合网格的时域有限差分法。由于必须在有限的计算区域内进行波场模拟计算，对于边界处非介质因素产生的反射现象，本文引入完全匹配层吸收边界条件，并达到了较好的消除边界反射的效果。考虑到提高差分阶数虽然能使交错网格差分的局部精度得到迅速提高，但同时计算的复杂性也随之递增，本文在不要求差分阶数很高的前提下，采用了混合网格有限差分法，也即在吸收边界层中使用交错网格法，而在有效计算区域内仍使用高阶空间差分算子的规则网格有限差分法，最后通过模拟算例结果验证

24、，采用该混合网格有限差分法能在节省内存和提高计算效率的同时得到更理想的结果，得到的波形信息比整个计算区域使用交错网格法时更加清晰。参考文献：1 Dmith W.D., The Application of Finite-element Analysis to body wave propagation problemsJ.Geophysics.Roy.Astr.Soc.,1975,42:747-768.2 GAZDAG J. Modeling of the acoustic wave equation with transform methodsJ.Geophysics,1981,46：854

25、-859. 3 KOSLOFF D,BAYSAL E. Forward modeling by the Fourier methodJ.Geophysics,1982，47:1402-1412.3 Seybert,A. F.,Cheng, C.Y.R.,Wu,T.W., The solution of coupled interior/exterior acoustic problem using the boundary element methods, J.Acoust.Soc.Am., 1990,88(3) 1612-1618.4 Kelly K R et al. Synthetic s

26、eismograms: A finite-difference approachJ. Geophysics, 1976,41:2-27.5 Alterman Z.and Karak F.C. Propagation of elastic wave in layered media by finite-difference methodsJ.Bull.,Seism.Soc.Am.,1968,58:367-398.6 Alford R M,Kelly K R,Boore D M. Accuracy of finite-difference modeling of the acoustic wabe

27、 equationJ.Geophysics,1974;39:838842.7 Mdariga R.Dynamics of an expanding circular faultJ. Bull.,Seism.Soc.Am., 1976,65,163-182.8 Virieux Jean. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference methodJ.Geophysics,1984,49(11):1933-1957.9 董国良,马在田等.一阶弹性波方程交错网格高阶差分解法J.地球物理学

28、报,2000,43(3):411-419.10 董良国,马在田等.一阶弹性波方程交错网格高阶差分解法稳定性问题J.地球物理学报,2000,43(6):856-864.11 裴正林.任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟J.石油地球物理勘探,2004,39(6):629-634. 12 L.Gilles, S.C.Hagness,and L.Vazquez. Comparison between Staggered and Unstaggered Finite-Difference Time-Domain Grids for Few-Cycle Temporal Optical S

29、oliton Propagation.Journal of Computational Physics 161, 379400 (2000).13 Berenger J. A perfectly matched layer for the absorbing of eletromagnetics waves J. Journal Computation Physics,1994,114:185-200.14 Qing-Huo Liu and Jianping Tao,The perfectly matched layer for acoustic waves in absorptive med

30、iaJ. Acoust. Soc. Am.102 (4), October 1997.15 X. Yuan, D. Borup, M. Berggeren, R. Eidens, and S. A. Johnson. Formulation and validation of Berengers PML absorbing boundary for the FDTD simulation of acoustic scattering. IEEE Pans. Ultrason., Ferroelect., Freq. Contr., vol. 44, no. 4, 1997, pp. 816-8

31、22.16 X. Yuan, D. Borup, J. Wiskin, M. Berggeren, and S. A. Johnson. Simulation of acoustic wave propagation in dispersive media with relaxation losses by using FDTD method with PML absorbing boundary condition.IEEE Pans. Ultrason., Ferroelect., Freq. Contr., vol. 46,no. 1,1999, pp.14-23.17 Guido Kn

32、eib and Claudia Kerner. Accurate and efficient seismic modeling in random mediaJ. Geophysics 1993;58, 576-588.18 Bjorn Engquist. Computational high frequency wave propagation J. Acta Numerica, 2003: 181266.19 吴国忱,王华忠.波场模拟中的数值频散分析与校正策略J .地球物理学进展,2005,20: 58-65.20 李锐坚.VPS高阶有限差分正演模型.J.西安石油学院学报,1996,11(2):38-41.21 张文生.波动方程成像方法及其计算.M.北京:科学出版社,2009.22 黄文武.弹性波和声波的时域仿真方法研究D.天津大学,2005. 23 Schneider J.B.，Wagner C.L. and Broschat S.L.，Implementation of Transparent Sources Embedded in Acoustic Finite-Difference Time-Domain GridsJ. Journal of the Acoustical Society of America,