复杂地表具有地表起伏剧烈、近地表速度低和速度变化范围大等特征，为保证数值计算模拟的精度，利用波动方程有限差分计算时若采用均匀网格有限差分法，则在低速介质区域网格步长必须取得比较小，空间网格步长的选取必须要满足采样定理，小于半个最小波长，最小波长与模型中最小速度(低速层的速度)成正比，若在一个波长上的采样点数过少会产生数值频散现象.而网格空间小步长采样会造成对高速介质区域出现过采样问题，增加计算量和计算机内存的消耗，同时小步长均匀网格要求有较小的时间步长，计算时间会增加，高速层区域的稳定性条件约束着时间步长的选择.受限于空间步长与模型最大速度的比值，空间步长越小，或者模型速度分布范围越大，相应的时间步长也要变得更小才能满足稳定性条件，否则会引起数值计算的不稳定.

为了能够有效解决这一问题，国内外学者做了许多研究，比如在模拟声波的传播和起伏地表条件下弹性波的传播过程中采用网格纵向步长变换法(Jastram and Tessmer, 1994).但该方法只考虑了地下介质速度随深度增加而增大的情况，对于横向介质速度的变化并不能处理好，同时此方法对于在介质区域内存在局部的高速或低速介质时也不适用.后来学者提出相似的不规则网格的方法.采用变网格方法一个不可避免的问题就是在网格相邻区域网格结点的计算(Opršal and Zahradník, 1999; Pitarka, 1999; Wang et al., 2001).早期的研究对这一问题的解决方法主要是采用插值法进行处理，但插值法容易加强算法的不稳定性，同时在实际应用时不能够较好的选取到合适的步长，插值方法的选取也没有一个完整的理论分析依据.后来有学者从波动方程出发，利用高阶有限差分方法，通过改进高阶有限差分系数的求解方法，消除在步长变化处由于步长的不一致而产生的计算不稳定，同时相对的减少计算量.有学者曾利用高阶差分方法模拟地震波在复杂介质条件下的传播.他们采用变网格方法依然只是纵向空间步长变化的单向变网格方法，而且这种通过改进高阶差分系数的方法需要在每个计算结点确定差分格式中每个相邻结点的步长，同时，这种方法虽然在低速层可以减小网格步长，但不能避免下层高速介质层水平方向出现过采样，该方法同样的对于近地表起伏及速度横向变化的模型依然不能很好的处理(孙成禹等，2008；张慧和李振春，2011；孙成禹和丁玉才，2012).亦有学者采用自适应网格，利用小波函数的多尺度分解的特性模拟地震波的传播(裴正林，2005；宓铁良等，2009), 该方法通过判定在初始网格点处的小波系数与阈值之间的大小决定是否对网格细化，但是利用小波函数会产生较大的计算量，同时该方法中对于粗细网格相邻处依然使用的是插值法来处理.非规则网格有限差分法将是提高复杂地震波模拟精度的主要方法(王庆等，2015).

在地震波有限差分数值模拟中不可避免的会出现数值频散现象，究其原因是用有限差分法求解波动方程时，偏导数用泰勒级数展开得到，空间导数算子的计算精度和准确性不仅依赖于近似的阶数还取决于离散化网格间距的大小，即使在各向同性的介质中，波动方程的有限差分解也有数值相速度和群速度，它们不同于波在真实介质中传播的相速度和群速度，随着波的传播，波前形状会发生变化并逐渐扩散开导致波形失真，称为数值频散现象(Sun and Trueman, 2005).后来将通量校正传输法，应用到声波和弹性波有限差分数值模拟和逆时偏移中来消弱数值频散，但这种方法会使不需要校正的地方分辨率降低，同时还会增加一定的计算量(李文杰等，2011；周竹生，2012)；对有限差分的差分系数和波动方程进行优化来压制数值频散，取得了良好的效果，然而压制数值频散的优化参数需要凭借经验才能选取合适的值，给数值计算带来了很大难度(Etgen, 2007; Liu et al., 2009).

通过上面的分析，当介质区域的速度变化范围较大时，在不同速度的介质区域采用不同的空间网格步长进行处理就显得比较合理.本文在早期学者研究的基础上，提出一种步长自适应的矩形网格生成法，在给定介质速度模型的条件下，随着速度的不同在网格划分时空间步长自动增加或减少，实现了对介质速度模型的精细化网格划分.利用有限差分在不均匀网格上模拟计算的一个关键问题就是网格过渡处的结点如何计算，本文从二元函数的泰勒展开式出发，结合相邻的几个结点，直接计算过渡处结点的拉普拉斯算子，通过分析求解得到一组相邻结点系数的计算公式，这些公式在四个方向上具有通用性，避免在相邻区域进行插值计算的同时，保持了算法的稳定性.对于频散现象，本文提出了一种新的方法，在使用有限差分法求解波动方程时，采用中值定理与有限差分法相结合构造新的离散化波动方程来处理数值频散问题，并经过推导给出了确切的频散校正参数值，我们称这种方法为抗频散有限差分法.

本文中的算法基于二阶声波方程，公式为

(1)

其中U=U(x, z, t)为标量波场，V为介质速度.

本文中对模型中不同速度区域自适应网格步长的划分是由模型中速度的分布范围决定.设模拟区域为Ω，区域中速度V∈(V_min, V_max)，V_min和V_max分别为区域中介质速度的最小值和最大值，在模拟区域Ω中进行网格划分时，可以划分的网格类型有种.相应的在网格划分时不同步长取值的速度阈值V_treshold可以有种取值，各阈值的选取应满足公式为

(2)

(3)

图 1为速度阈值的选取示意图.

图 1

Figure 1

图 1 速度阈值选取示意图 Figure 1 A schematic view of the speed threshold

假设在V_treshold¹内的速度区域网格步长为Δh, 那么其他的速度阈值对应的速度区域内，网格步长是以Δh的二进格式取值，即：

(4)

各个速度区域步长的取值与速度阈值之间的对应关系为

(5)

其中在速度的分布范围内，取Δh时，，取Δh_R时，V_treshold^R=V_max.图 2为采用本文自适应网格划分方法得到的不均匀模型的结果.

图 2

Figure 2

图 2 自适应网格划分示意图 Figure 2 Adaptive meshing schematic

图 3

Figure 3

图 3 五点差分格式 Figure 3 Five point finite difference

一般的用有限差分法计算波场的空间导数项时，在低波数区域的有限差分法可以近似的符合，在高波数区域会产生频散(朱生旺等，2007).为了在数值计算时能够尽可能地减少数值频散的产生，在空间网格步长的选取时应满足一定的条件.设地震波场的最大波数为k_max，它由子波的最高频率f_max和介质的最低速度V_min决定，即：

(6)

设k_α为可以忽略频散的最高波数，即有k < k_α.一般的取k_α=1/(2Δh)(朱生旺等，2007)，若使得数值计算的结果不受频散的影响，则Δh的选择必须满足：k_α>k_max，即：

(7)

上式(7) 就是细网格步长Δh在选取时应该满足的条件.同理，对于其他速度区域的介质，其最小速度就是速度阈值，此时网格步长应该满足：

(8)

因此在网格步长的选取时满足式(7) 和式(8) 即可.

假设在计算区域Ω中，各个空间步长对应的介质区域在Ω中所占比例为，那么可以计算得到采用本文自适应网格划分方法需要的网格点比全域采用细网格进行划分节省的网格点比例为

(9)

如图 2所示的模型，三种不同步长网格所占比例分别是0.158、0.158、0.2、0.484，由(9) 式可以计算得到，利用自适应步长法可以节省的计算网格结点比例为96.27%.

同时，在利用有限差分法计算时，时间步长及空间步长的选择应该满足一定的稳定性条件.由Courant-Friedrichs-Lewy(CFL)条件有：

(10)

由上(10) 式可以推知，利用自适应步长法在不同网格区域，计算时满足的稳定性条件为

(11)

将(2) 和(4) 式代入到(11) 式中，不同步长区域对应的选择满足公式为

(12)

由(12) 式可见，本文提出的自适应步长方法，对时间步长的要求降低，只要低速区域时间步长满足稳定性条件则其在其他区域依然能够满足稳定性条件.若是整个计算区域均采用细网格划分，为了满足在高速层计算的稳定性，时间步长Δt只能选择比较小的值.同样的以图 2所示的模型为例，若在低速层的网格步长取4 m, 在满足低速层稳定性条件时，时间步长Δt最大可以取2 ms, 若是采用细网格划分，则要满足高速层的稳定条件，则时间步长Δt最大可以取到0.1 ms.而采用自适应网格法，时间步长的选取则只要满足在低速层的稳定性要求即可.相比而言对于图 2所示模型，自适应网格法理论上对相同的模拟时间可以最大节省约90.25%的计算时间.

设计算区域中步长为Δh_R-1的网格结点有I×J个，波场可以表示为：U_{i, j}^k，采用其他步长的网格点可以表示为：

(1) 考虑到在低速层处的网格采用细网格步长进行划分，这种精细的划分在模拟中已然能够准确描述地质的特征，因此在细网格的内部结点处，采用五点有限差分格式.

步长为的细网格点，与其相邻的其余四个点分别为.处的五点差分格式为

(13)

(2) 对于在不同网格类型相邻处的结点波场值，其中细网格结点与其相邻的粗网格，在某一方向上会缺失其相应的一些细网格结点的波场值，这对计算这些结点造成困难.对于这种情形，本文利用二元函数的泰勒展开式，通过与计算结点相邻的几个结点在计算结点处的二元泰勒展开式的组合，构建一个方程组，由这个方程组直接求解在相邻处细网格结点波场值的拉普拉斯算子，然后使用时间二阶差分迭代计算相邻处结点在每一时刻的波场值，这种方法不需要求解那些缺失的结点波场值，通过方程的推导，这种方法在不同方向上的相邻结点处，其构成的方程组中的方程只是在方程的排列顺序和正负号上有着轻微的改变，对于相邻的结点，按照一定的顺序排列则与其相乘的系数的计算公式具有统一的形式.

假设步长为Δh的细网格与其相邻的粗网格步长为Δh_R-1.以左边界处的相邻网格结点计算为例，设在一个细网格的左边是一个粗网格如图 4所示.

图 4

Figure 4

图 4 相邻网格左边界处 Figure 4 Left border of different grid step

当前计算结点为边界结点U_{i, j, 0, m}^k，若使用传统的方法求解则缺失U_{i, j, 0, m}^k左边相邻结点的波场值，为此利用U_{i, j, 0, m}^k周围的五个点求解，在左边界时用到的五个结点分别为：左下粗网格U_{i-1, j+1}^k，下边细网格U_{i, j, 0, m+1}^k，左上粗网格U_{i－1, j}^k，上边细网格U_{i, j, 0, m+1}^k，右边细网格U_{i, j, 1, m}^k，各点的位置如图 4中所示.上述五个点在U_{i, j, 0, m}^k处的二元泰勒展开式，舍去三阶及以上的高阶偏导数项有：

(14)

(15)

(16)

(17)

(18)

公式(14) 至(18) 分别乘以系数a_i(i=1, …5)，使得它们满足公式为

(19)

将式(19) 写为方程组形式有：

(20)

解方程得各系数的计算公式为

(21)

(22)

(23)

(24)

然后在时间上采用二阶中心差分，则U_{i, j, 0, m}^k的迭代计算公式为

(25)

其余三个方向上的边界结点计算式各结点的分布如图 5所示，图 5a至图 5c分别为右、上、下三个方向，相应的也可以得到这三个方向上系数的求解方程组，式(25) 至(27) 分别为右、上、下三个边界系数的方程组.公式为

图 5

Figure 5

图 5 粗细网格右、上、下边界处 Figure 5 Right, top and bottom border of different grid step

(26)

(27)

(28)

(26) 和(20) 两个方程组的区别只是在于前两个方程的符号发生变化，但由于右边项为零，因此这两个方程组是相同的，和式(20) 相比较，式(27) 和(28) 是第一个和第二个方程的位置进行了调换，并且正负号发生变化，但是由于右边项为零，因此前两个方程实际并没有变化，仅仅改变了位置，同时方程位置改变的还有第四个和第五个方程，所以整体的上面的两个式子和式(20) 是相同的，因此它们求解系数a_i(i=1, …5) 的计算公式不变.

(3) 对于网格过渡处的细网格四个角点，在差分计算时需要单独计算，对于角点处结点本文采用与计算细网格内部结点相同的差分格式，即式(12).但不同的是在计算每个角点时需要判断角点周围的网格类型，从而确定x₁、x₂、z₁及z₂的值.

如图 6所示若要计算其左上角节点U_{i, j}^k=U_{i, j, 0, 0}^k处波场值则需要判断U_{i－1, j－1}^k、U_{i, j－1}^k以及U_{i-1, j}^k点的速度所属区间.假设步长为的网格与其相邻的网格步长可能取值为.

图 6

Figure 6

图 6 细网格角点示意图 Figure 6 Fine grid corner schematic diagram

利用五点差分格式求解左上角点U_{i, j, 0, 0}^k时，x₂与z₂为细网格步长，即x₂=z₂=.一般的对于间距不相等的网格点利用五点差分格式计算时，假设中心点为U_{i, j}^k，其余四个点分别为：U_{i-1, j}^k，U_{i+1, j}^k，U_{i, j-1}^k，U_{i, j+1}^k，那么其x方向上的二阶导数可表示为

(29)

同理，z方向上的二阶导数可表示为

(30)

最终一般的五点差分格式为

(31)

在已知x₂与z₂的情况下，利用(31) 式进行计算需要得到x₁与z₁的值.而x₁的取值需要通过判断U_{i-1, j－1}^k、U_{i-1, j}^k结点速度所属区间确定：

(32)

而z₁的值需要通过判断U_{i-1, j－1}^k、U_{i, j-1}^k结点的类型确定，其判断方法与x₁的方法相似，如下：

(33)

假设，式(32) 对应的x₁取值的网格分布情况如图 7所示，图中绿色小圆圈表示左上角点U_{i, j, 0, 0}^k在计算时用到的左边结点.

图 7

Figure 7

图 7 左边结点分布情形 Figure 7 The left node distribution

式(33) 对应的z₁取值的网格分布情况如图 8所示，图中绿色小圆圈表示左上角点U_{i, j, 0, 0}^k在计算时用到的上边结点.

图 8

Figure 8

图 8 上边结点分布情形 Figure 8 The top node distribution

(4) 对于步长为Δh_R网格结点，本文采用偶数阶精度的高阶有限差分格式，公式为

(34)

其中为差分系数(Liu and Sen, 2009).

一般的在进行地表低速层进行模拟时，会产生比较严重的数值频散现象，这对地震波模拟的记录产生严重的干扰，为此如何在近地表模拟时很好的消除数值频散现象也是地震波数值模拟的一个重点.为本文利用中值定理的原理，提出一种能够很好的消除数值频散现象的有限差分格式.

首先来考虑下面这个初值问题，公式为

(35)

(36)

现用显式欧拉方法将(34) 式和(35) 式在时间上用中值定理加权近似表示为

(37)

(38)

式中Δt为时间步长，0 < α < 1，0 < β < 1分别为加权系数.将(37) 式代入(36) 式，经整理可得：

(39)

对(38) 式中Δtf(u_n)和q(u_n+1)分别用Taylor展开法近似，公式为

(40)

(41)

将(39) 式和(40) 式分别带入到(38) 式中经整理忽略高阶导数项可得：

(42)

现将q(u_n－1)运用中心差分法近似为

(43)

将(42) 式带入到(41) 式中得：

(44)

假设u=u(x, y, z, t)表示波场，那么(43) 式中的q(u_n)= ，现在对声波方程(1) 在二维情况下变形，令v=v(x, z)，，则有：

(45)

对波场进行正演外推时，将(44) 式代入(43) 式中对时间做离散处理可得：

(46)

现在令(45) 式可化简为

(47)

数值频散主要是因相速度与群速度不一致造成的，那么当相速度和群速度近似相等时，数值频散就会明显地消弱，这样可以解出抗频散的优化校正参数，在空间二阶近似下：

(48)

其中，在实验算例中p(θ)取算术平均值3/4.而在空间四阶近似时，可以求出频散优化校正参数b=1/6.

在波动方程的数值模拟中，一个关键的问题就是边界条件，因为在实际的地震波的传播是在地下无限区域中传播的，而在模拟计算时计算区域是有限的，这样就人为的增加了边界.这种人为增加的边界会产生边界反射波，干扰波场中的有效地震波信息，因此在对地震波进行数值模拟时需要消除这些人工边界反射波.

本文中采用二阶傍轴近似边界和衰减边界条件相结合的吸收边界条件(Liu and Sen, 2010)，使用二阶傍轴近似能够吸收部分的边界反射，其余的反射能量通过衰减边界条件吸收，这样可以在增加少量计算量的同时很好的消除人工边界反射.

傍轴近似边界条件是在边界网格点处用单向波动方程代替原波动方程，从而达到边界吸收的效果.一般的，在某一传播方向上较小的范围内用单向波动方程可以正确的表示波在该方向上的吸收，而在其他方向上的波仍会被反射.实际上，由于对单向波动方程的有限差分离散化方法以及其误差项等，可能使得在这些边界网格点处计算的波场值是不正确的.而且，傍轴近似条件在边界处采用单向波动方程，使得边界与内部网格点由于采用的波动方程不同，引起比较尖锐的变化.为此本文通过将傍轴近似边界条件和衰减边界条件相结合，在边界与内部网格点之间增加一些过渡层，减小由于波动方程不同而引起的波场的变化.

将扩充边界后的计算区域Ω划分为如图 9所示三部分：(Ⅰ)内部区域Ω^iner；(Ⅱ)过渡区域Ω^trans；(Ⅲ)边界Γ.在三个区域，分别采用不同的波动方程计算，区域Ω^iner内各点的波场值采用原二阶波动方程计算，其具体的差分格式由区域内网格的类型决定.

图 9

Figure 9

图 9 计算区域示意图 Figure 9 The calculation area schematic diagram

边界Γ上各点的波场值采用单向波动方程B₂计算，以下边的区域为例, 公式为

(49)

区域Ω^trans内各点的波场值是两种波场值的加权平均，其他方向的过渡区域方程可以通过将式(1) 与其他方向上的B₂方程相结合得到.

区域Ω^trans减小了内部区域Ω^iner与边界Γ之间由于波动方程的不同产生的突变，随着区域Ω^trans层数的增加，边界的反射波会逐渐的减小.当吸收层数N=1时，就退化为传统的傍轴近似吸收边界.

对某计算时刻，整个波场的计算分为三个步骤：

(1) 利用原波动方程(1) 式计算区域Ω^iner和Ω^trans的波场值U_i⁽¹⁾；

(2) 利用单向波动方程(49) 式计算Ω^trans和Γ的波场值U_i⁽²⁾；

(3) 对区域Ω^trans处的波场值进行加权平均，公式为

(50)

其中U_i⁽¹⁾和U_i⁽²⁾分别为在层B_i处原波动方程和单向波动方程的波场值，β_i为在层B_i处的加权系数，其取值范围为0, 1，U_i为层B_i处的最终波场值.

权值β_i的取值方式可有多种形式：指数型、抛物线型、线性型等，在本文中选取β_i指数变化，即：

(51)

其中，d_i为B_i距离区域(I)Ω^iner的边界的距离，L为增加的层的厚度.

为验证上述算法的正确性、有效性，设计如图 10所示的低速夹层模型，震源位于水平地表的中心处.

图 10

Figure 10

图 10 低速夹层模型 Figure 10 Velocity Layers Model

对于上述模型分别采用常规网格和本文方法进行模拟计算.模型大小为2500 m×2500 m，采用常规网格计算时步长为10 m，利用本文方法时网格步长的取值有两种5 m和10 m，采集时间均为1.5 s.震源为主频30 Hz的雷克子波，差分长度M=4，在计算区域外层增加10层的衰减层.

上图 11、图 12和图 13分别为采用常规粗网格方法、本文方法以及常规细网格法得到的单炮道集记录.由图 11可以看出，采用常规10 m步长网格得到的道集记录中，在低速夹层的界面放射波存在严重的数值频散现象，而且夹层下界面的频散程度要强于上界面的频散，究其原因是在低速夹层的空间采样过低造成.与图 11相比较，图 12中采用本文方法得到的道集记录中各界面的反射波清晰整洁，基本不存在频散现象，因此本文方法能够提高模拟的精度，有效压制数值频散.图 12和图 13相比较可以看出，本文方法得到的模拟结果与细网格方法得到的结果具有相近的精度和效果，可见本文方法是可以几乎达到常规细网格方法的结果的.本文采用的方法其内存需求量仅为常规5 m步长小网格的75%，计算时间节省约33%，相比较采用10 m步长的常规网格，内存的需求增加约15%，计算时间增加了约5%.

图 11

Figure 11

图 11 常规网格10 m步长单炮记录 Figure 11 Conventional grid step length of 10 m single shot records

图 12

Figure 12

图 12 自适应网格步长为5 m和10 m单炮记录 Figure 12 Adaptive grid spacing of 5 m and 10 m single shot records

图 13

Figure 13

图 13 常规细网格网格步长为5 m单炮记录 Figure 13 Conventional fine mesh grid spacing of 5 m single shot records

为验证本文方法对复杂地表的处理能力，构建图 14为复杂地表两层介质模型.速度模型网格结点大小为：1001×465，地表上层速度填充以声速340 m/s，地下介质速度从上至下分别为：500 m/s、1000 m/s、2000 m/s、3000 m/s和4500 m/s, 空间网格步长分别为1 m和4 m，时间步长为0.1 ms，共采集1.5 s时长.震源采用主频30 Hz的雷克子波，位于网格(254, 36) 处，采用的差分长度为M=4，在计算区域外增加10层的衰减层.

图 14

Figure 14

图 14 复杂地表模型 Figure 14 Complex surface model

图 15为常规4 m粗网格得到的单炮道集记录.从图 15中可以看出，由于近地表存在低速层，对模型采用大网格计算时，在低速层产生的了非常严重的数值频散现象，以至于掩盖了下层介质的反射波信息.图 16为本文方法得到的道集记录，从图 16中可以看出，相比较于常规粗网格方法，本文方法得到的道集记录中能够正确的记录地下各个界面的反射波.从图 16可以看出，由于地表起伏，下层界面的反射波及地表直达波的同相轴发生剧烈的扭曲.而图 17为常规1m细网格方法得到的道集记录，由图可见采用细网格划分能够很好的抑制近地表低速层产生的频散现象.比较图 16和图 17可知，本文方法得到的道集记录与常规细网格方法得到的道集记录有着相似的精度，但比之常规粗网格方法有着明显的优势.

图 15

Figure 15

图 15 常规4 m粗网格道集记录 Figure 15 Conventional 4 m coarse grid gather records

图 16

Figure 16

图 16 本文方法得到的道集记录 Figure 16 Gather records obtained by this method

图 17

Figure 17

图 17 常规1m细网格道集记录 Figure 17 Conventional 1m fine grid gather records

由图 15至图 17可以看出，利用本文方法模拟复杂地表条件下的地震波传播，也可以得到比较精确的模拟结果，基本不存在数值频散.在复杂地表条件下一般的模拟时均采用小步长网格进行划分，利用本文方法计算时，内存的需求减少了约67%，而计算时间减少约45%.上面两个模型的计算环境均为：Windows 7操作系统，CPU为Intel Core2 Extreme X7800 2.6GHz，内存：2 GB，IDE为Code::Block 12.11.由上面的两个实验可以验证采用本文方法能够节省计算时间，提高计算效率同时计算机内存也有很好的节省.

本文通过构建一种自适应网格划分方法，提升了对模型的划分精度.同时在不同网格类型相邻区域采用二元函数泰勒展开式的组合，直接求解相邻区域结点的空间拉普拉斯算子，避免使用插值法，减少了计算量.同时对采用有限差分模拟易产生数值频散的缺点，本文提出一种能够抑制频散的差分格式，并给出在不同差分阶数时，抗频散校正参数的最优取值.通过不同的模型实验结果分析可得，本文提出的方法用于地震波的模拟是正确有效的，起伏地表模拟说明本文方法能够很好的处理地表起伏条件.实验结果表明，相比较于细网格方法，本文方法在提升模拟精度的同时，能够有效的减少内存需求以及减少计算时间.