0 引言

随着我国勘探开发逐步转向西部探区，基于起伏地表条件的地震波正演模拟、速度分析和偏移成像等方法已成为勘探开发过程中不可回避的问题^[1]。起伏地形会对地震波的传播过程产生较大影响，造成体波和面波强烈的散射和绕射，使地震波场更加复杂且波形发生畸变；许多传统的基于水平层状介质理论的地震数据采集、处理和解释方法技术已经不再适用^[2]。研究地震波在起伏地表条件下的波场特征和传播规律，对于该地质条件下地震勘探观测系统的设计、数据处理、解释以及反演算法的研究都具有重要的理论和实际意义^[3-4]。

有限差分法在地震波正演模拟中应用较为普遍。常规有限差分法处理起伏地表面临两个难题：一是复杂界面的网格化离散，二是自由表面条件的实施。在模拟复杂地表及地下复杂界面时，应用笛卡尔坐标系中的规则网格剖分，易出现阶梯状的边界，进而产生虚假的绕射波^[5]。为此，国内外研究人员针对传统有限差分方法进行了大量的研究和改进。变网格算法可用于加密起伏界面^[6-7]，但该算法依然属于规则网格类型，因而不能从根本上消除角点处的虚假绕射和散射现象。为彻底消除阶梯网格近似产生的负面影响，Hestholm等^[8]、董良国等^[9]、Tarrass等^[10]采用垂向网格映射法，将起伏边界变换为水平边界，然后在变换后的矩形网格中进行数值模拟。然而，该方法只能处理地形起伏平缓的情况，在剧烈起伏条件下易出现不稳定。

Fornberg^[11]、Komatitsch等^[12]将贴体网格引入到起伏地表正演中，实现了贴体网格下的伪谱法正演模拟，证实了贴体网格对任意起伏地形的适应性。Zhang等^[13-14]实现了基于贴体同位网格的一阶速度应力方程有限差分正演模拟算法，并在实施自由边界条件时提出了牵引力镜像法，得到了较好的模拟效果。由于同位网格下的一阶中心差分格式会引起奇偶失联高频振荡现象，需要特殊的滤波处理，其采用了DRP/opt MacCormack格式^[15]计算空间导数，基本消除了格点振荡现象，但增加了实现的复杂性；且为了获得同交错网格相同的精度，需要使用更小的网格步长，增加了计算成本。随后，Appelo等^[16]、兰海强等^[17]推导并实现了贴体网格下的二阶波动方程正演模拟算法，其空间差分精度在边界和内部都是二阶。然而二阶位移方程在泊松比较大的介质下容易不稳定，且该算法不易推广到高阶。丘磊等^[18]实现了曲坐标系下的标准交错网格(standard staggered grid，SSG)正演模拟算法，但由于曲线坐标系下的波场变量不满足交错分布，其采用四阶插值算子给出相应的缺失点的波场信息，不仅降低了模拟精度还增加了计算量。

综上所述，虽然贴体网格能够很好地拟合任意起伏地形，但其与有限差分算法的结合仍需要进行进一步的研究。此外，目前对贴体网格的研究大多集中在算法的改进上，很少有对其影响因素和适用范围进行分析的^[19]。为此，本文在对常规贴体网格正演模拟算法做了几点改进后，重点讨论了网格剖分质量对模拟结果的影响，研究了网格正交性、网格间距的突变以及不同类型的网格拼接对算法的影响，以期认清算法的适用范围。

1 曲线坐标系下的波动方程

正演模拟的第一步是网格化离散。贴体网格是一种适合复杂地表介质的网格离散方法，网格生成的原则是使离散后的网格边界与地表形态吻合，以避免阶梯边界引起的虚假散射。如图 1所示，贴体网格可以通过由计算空间到物理空间的坐标变换来获得。本文采用孙建国等^[20-21]提出的方法生成网格。

贴体网格生成之后，笛卡尔坐标系中的网格点和曲线坐标系下的网格点之间也就建立了一一对应的关系：

对式(1)分别求关于x，z的偏导，引入记号f_x≡∂f/∂x，由链锁法则可得

由式(2)整理得到

笛卡尔坐标系下基于各向同性非均匀介质的2D弹性波一阶速度应力方程为

其中：ν_x,ν_z表示质点速度；τ_xx,τ_zz,τ_xz表示应力；ρ为密度；λ，μ为拉梅常数；t表示时间。

利用链锁法则和公式(3)、(4)，可推出曲线坐标系下的波动方程：

2 全交错网格

由曲线坐标系下的波动方程(公式(5))可知，曲线坐标系下每个变量都需要计算同一格点在两个方向的空间导数；常规的标准交错网格(SSG)已不满足网格分布需求，需要复杂的波场插值，这就降低了模拟精度。同位网格虽可避免插值运算，但该网格下的中心差分会引起奇偶失联高频振荡现象，需要加入人为的滤波压制，增加了实现的复杂性、降低了模拟精度。鉴于此，本文将全交错网格^[22-24]引入曲线坐标系中，全交错网格(fully staggered grid，FSG)的变量分布如图 2所示；其网格机制的主要思想就是速度(图 2中方形)和应力(图 2中圆点)的不同分量交错定义在同一网格点上，从而满足了公式(5)的交错分布特性，避免了插值和滤波运算，提高了模拟精度、降低了实现的复杂性。因同一变量定义在同一网格的两个不同位置，所以需要分别更新计算；具体差分格式与常规SSG相同，因而简单易实现。

由于FSG仍然属于交错网格的范畴，因此可以直接使用交错网格的差分格式。且FSG的稳定性条件和频散关系等都与SSG相同，同样比同位网格具有更加宽松的频散条件，要求的每个最小波长内的网格点数(grid points per shortest wavelength,PPW)更少，从而可一定程度上降低内存和计算量。

3 曲线坐标系下的自由边界条件

自由地表模拟时，自由边界条件的实施方案是关键。由于将全交错网格引入到曲线坐标系下，因而需要对Zhang等^[13]提出的自由边界条件进行修改。本文采用牵引力镜像法计算速度分量，速度自由边界条件配合紧致交错差分格式更新应力分量。

起伏地表的自由表面条件要求牵引力为0，即T=n·τ=0。其中，为自由表面的单位法向量，。将代入T=n·τ=0可得

对式(6)关于时间求偏导，并代入公式(5)可得到速度自由表面条件：

其中：

式(8)即为二维弹性介质的速度自由边界条件，可用于更新自由边界附近的应力变量。假设采用空间四阶差分格式，则具体的应力变量计算步骤为：

1) 模型内部区域(边界之下两点)采用常规空间四阶差分格式计算速度变量的空间导数，利用公式(5)可更新相应的应力变量；

2) 自由边界上，先计算速度变量的横向导数，然后利用速度自由边界条件(公式(7))可求出速度变量的纵向导数，再根据公式(5)可实现自由边界上应力变量的更新；

3) 由步骤1)、2)可知边界及边界之下两点速度变量的纵向导数，借助四阶紧致交错差分格式(公式(9))可求出边界之下一点速度变量的纵向导数，而速度变量的横向导数可直接计算得到，从而可更新相应的应力变量。

空间四阶的紧致交错差分格式为

其中：f′_i为待求变量f在格点i的导数；Δx是网格间距。

速度变量的更新流程为：

1) 模型内部区域(边界之下两点)采用常规空间四阶差分格式计算应力变量的空间导数，利用公式(5)可更新相应的速度变量；

2) 边界附近可利用牵引力镜像法计算。

牵引力镜像法为水平地表应力镜像法的推广，应用牵引力镜像法时需要保守形式的动量方程^[13]：

式(10)中关于η方向求导的项，恰好为公式(6)中牵引力的分量，从而分别对(η_xτ_xx+η_zτ_xz)、(η_xτ_xz+η_zτ_zz)关于地表应用牵引力镜像，设置虚拟层内的牵引力分量，即可实现起伏界面的应力自由表面条件。

4 影响因素分析 4.1 算法有效性测试

采用高斯山峰模型，通过与同位网格应力镜像法对比，以检验本文算法的有效性以及宽松的频散要求。其高程表达式为

其贴体网格剖分示意图如图 3a所示，图 3b为山峰处的局部放大图。从图 3可以看出，本文网格剖分的正交性非常好。双层山峰模型的参数如下：第一层纵波速度为3 000 m/s，横波速度为1 732 m/s；第二层纵波速度为4 000 m/s，横波速度为2 310 m/s；模型网格大小为401×241，网格间距约为10 m。震源为主频18 Hz的雷克子波，时间采样间隔0.4 ms，最大记录时间2 s，炮点位于(1 000 m，20 m)处。由以上参数可知，每个最小波长内的网格点数约为4。图 4a、b为分别采用同位网格应力镜像法^[13]和本文算法计算得到的0.6 s的水平分量波场快照：图中不论是P波还是Rayleigh面波经过山峰时，都会产生反射波和转换波，且P波产生的反射Rayleigh面波能量比Rayleigh面波产生的反射Rayleigh面波能量弱很多。小山峰相当于一个二次震源，所有经过山峰的波都会激发出反射波和转换波，极大丰富了波形、增加了波场的复杂性。对比图 4a、b可知，在相同的计算参数下，由于同位网格应力镜像法需要的PPW较大，其波场快照中S波和Rayleigh面波已发生明显空间频散(箭头位置)，而本文算法则依然较好。要达到同样的效果，同位网格应力镜像法需要减小网格间距，但这就增大了内存和计算耗时。

4.2 网格正交性因素分析

为了便于分析，以水平地表模型为例，分析网格正交性对贴体网格正演模拟效果的影响。模型为均匀介质，纵波速度为3 000 m/s，横波速度为1 732 m/s。保持上边界网格固定，通过倾斜左右边界，生成三种网格，如图 5a、b、c所示，网格夹角分别为90°、45°和33.7°。正演模拟所用参数为：横纵向网格数目为401×201，横向网格间距为10 m，纵向网格间距分别为10 m(图 5a)、14 m(图 5b)和18 m(图 5c)；震源子波为主频10 Hz的雷克子波，时间采样间隔为0.5 ms，震源位于地表中心(2 000 m，10 m)位置处。图 6a、b、c分别为图 5a、b、c三种网格夹角模型在600 ms时刻对应的水平分量波场快照。对比图 6a、b、c可以看出，尽管网格越来越扭曲(正交性变差)，但正演模拟结果没有受到太大影响，都能准确模拟地震波传播过程；说明本文算法对网格的正交性没有太大要求。当震源主频提高到20 Hz时，90°和33.7°网格夹角模型在600 ms时刻的水平分量波场快照如图 7a、b所示，可以看出，90°网格夹角模型仍能正确模拟波的真实情况，而33.7°网格夹角模型则出现了较为严重的频散现象，严重影响了正演模拟效果。究其原因是由于33.7°网格夹角模型的纵向网格间距较大，从而对频散关系的要求较高。综上所述，本文算法对网格的正交性没有太大要求，只要正演模拟参数的选取满足频散关系，本算法都可得到较好的效果。

4.3 网格间距突变因素测试

本节以山谷模型为例，分析网格间距突变对贴体网格正演模拟效果的影响。模型横向5 km，纵向2 km，其左侧部分为山谷地形，右侧部分为水平地表。在横向3 km处发生了网格突变，网格剖分如图 8左栏所示：从上到下分别发生了2倍、3倍和4倍的横向网格间距突变，纵向网格间距均为10 m，最大横向网格间距分别为20、30和40 m。正演模拟所用震源为主频10 Hz的雷克子波，震源位于地表中心(2 500 m，10 m)处，时间采样间隔0.5 ms。图 8右栏为600 ms时刻对应左侧网格的水平分量波场快照。从图 8可以看出：当网格间距突变3倍时，突变界面处开始出现虚假反射，且由于S波传播速度比P波低，S波的虚假反射更明显；突变4倍时，虚假反射能量增大。分析可知，网格间距的变化引起频散关系的改变，造成突变界面两侧的相速度不同，因而会产生人为的反射。综上所述，网格剖分越均衡越好，即使网格间距需要变化，也应尽量平缓过渡；为了不出现明显的虚假反射，突变倍数最好限制在2倍以内。

4.4 网格系统拼接影响分析

本节考虑不同网格类型拼接的影响。模型如图 9a所示，网格系统的网格正交性突变十分剧烈：左侧为山峰模型，网格系统为90°夹角的正交网格；右侧为向下倾斜模型，网格系统为网格夹角26.6°的倾斜网格。模型为各向同性均匀介质，纵波速度为3 000 m/s，横波速度为1 732 m/s，模型网格大小为301×201，纵向网格间距均为10 m，正交网格系统的横向网格间距为10 m，倾斜网格系统的横向网格间距约为22.4 m。震源子波为主频10 Hz的雷克子波，时间采样间距为0.5 ms，震源位于模型中心(1 500 m，1 000 m)处。图 9b为500 ms时刻的水平分量波场快照，从图 9b可以看出，虽然网格的正交性突变较大，但正演模拟结果仍然十分稳定，且没有出现明显的虚假反射；说明本算法对不同网格系统的拼接具有较好的适应性。

5 结论与认识

本文将全交错网格推广到曲线坐标系下，结合牵引力镜像法和紧致交错差分格式，实现了起伏地表下的地震波有限差分数值模拟，并对其影响因素进行了详细研究，通过理论分析及模型试算取得了如下几点认识：

1) 在曲线坐标系下采用全交错网格，避免了标准交错网格的波场插值和同位网格的高频振荡问题，提高了模拟精度，减小了实现复杂度；且全交错网格算法和交错网格具有相同的差分模板，相比同位网格需要的PPW较小，具有更宽松的频散条件，网格间距可以取得较大，因而可在一定程度上提高计算效率。

2) 首先，贴体全交错算法对网格的正交性要求较低，只要正演模拟参数的选取满足频散关系，则任意网格夹角都可得到较好的效果。其次，网格间距的突变会引起虚假反射的产生，因而需要控制网格剖分的质量；为了较好地压制人为反射，建议突变倍数尽量限制在2倍以内。最后，本算法的适应性较强，即使不同类型的网格拼接对模拟结果也不会造成明显的影响。