基于低降速带复杂地形的瑞雷波数值模拟研究

王道远，王鑫，李 勤

西安石油大学地球科学与工程学院陕西省油气成藏地质学重点实验室，陕西西安710065

0 引言

在地震勘探过程中，随着丘陵、盆地等近地表地区的深层次研究，发现近地表地震-地质条件下接收到的地震波场成分十分复杂，地震资料的精细处理面临着巨大挑战，数值模拟研究越发必要。在反射地震勘探中经常会接收到一类低频率、强能量的面波干扰，通常作为地震勘探中的干扰波被压制掉。但是面波中含有丰富的信息，与近地表低降速带有直接的关系。瑞雷面波勘探是一种有效获取近地表地层结构的途径，可以有效利用反射地震资料中的瑞雷面波，使其变废为宝，节省一笔非常可观的勘探费用[1]。近年来，瑞雷面波勘探研究正在进行广泛开展。

地震数值模拟是研究地震波在复杂介质中传播的有效方法之一[2]。Rayleigh首次提出了瑞雷面波理论，并于1887年证明了瑞雷波的存在[3]。Haskell在最早时提出瑞雷面波频散曲线正演Haskell算法，此算法易出现高频数值不稳定现象，而且通常计算速度偏慢[4]。为解决所面临的问题，众多学者从数值角度开展研究，简化频散方程的表现形式，解决了频散方程的不稳定性导出了新的频散方程，有效避免了高频数值的不稳定，同时计算速度提升[5-8]。

为了通过地震波场更好地分析瑞雷面波在近地表弹性介质中的表现特征、响应机理及传播规律，国内学者进行了一系列研究工作。周竹生等进行了基于弹性介质中瑞利波有限差分数值模拟，采用交错网格有限差分法对多种弹性介质模型进行了正演模拟，得到相应的波场快照和地震记录[9]。熊章强等进行了吸收边界及自由边界条件研究，正演模拟得到震源在自由界面时的瑞利面波波场，并解释了波场特征[10]。张大洲等通过对三层介质数值模拟得到波场记录，揭示了瑞雷波传播特征[11]。邹延延等模拟了简单夹层和岩溶介质模型地震波动力特征[12]。众多学者进行了复杂近地表的地震波波场模拟，得到了复杂近地表波场响应特征[13-15]。邵广周等对比、分析瑞利波的频散特征和各模式的实际情况，同时也了解了瑞利波各模式在波场中的能量分布情况[16]。苏传行等进行了起伏地形瑞雷波模拟研究[17]。

虽然上述方法对瑞雷面波进行了有效的模拟，但对于近地表低降速带中的波场特征，尤其是瑞雷面波的特征分析不够全面。基于高阶交错网格有限差分原理，对不同复杂近地表地质模型进行了数值模拟，分别从合成地震记录和波场快照两个方面分析了不同条件下瑞雷面波的响应特征及其传播规律。在此基础上，对面波频散及数值模拟边界条件等问题进行了综合、全面的分析，从而为实际工程勘探中近地表低降速带区域含有软弱夹层和异常体的情况提供了参考。瑞雷面波正演模拟有助于更准确、更全面地了解瑞雷波的传播特征，为实际工程勘探提供理论依据。

1 波动方程数值模拟基本原理

1.1 各向同性介质一阶速度-应力弹性波方程

在二维情况下，一阶速度-应力弹性波方程如下所示[18]：

式中，θxx、θzz、θxz为应力的三个分量；vx表示速度的水平分量，vz表示速度的垂直分量；λ、μ表示拉梅常数，λ=ρ（vp2-2vs2），μ=ρvs2，vp表示介质中的纵波速度，vs表示横波速度；ρ表示介质密度。

1.2 交错网格高阶有限差分格式

1.2.1 时间上的2M阶有限差分近似

通过利用交错网格求取一阶应力—速度方程时，速度分量在t+Δt/2时刻进行计算，应力分量在t时刻进行计算。将各分量进行Taylor公式展开化简得：

（2）式和（3）式相减可得2阶时间差分近似公式：

式中，Δt为时间步长。

1.2.2 空间上的2N阶有限差分近似

空间上采取2N阶有限差分近似，使用Taylor公式展开计算得：

式中：为待定的差分系数，Δx为空间网格间距。对于待定系数，根据Taylor公式的2N阶展开，有：

可以求得高阶空间差分系数。

1.3 震源设置

震源函数可以使用高斯一阶导函数子波、雷克子波等表示，在进行全波场数值模拟时选取不同的震源函数对结果影响很大[19]。

选取雷克子波，频率fm=30Hz为震源函数。

1.4 稳定性条件

稳定性条件在进行数值计算时较为关键，进行合理的参数选择，可以避免出现严重频散、结果溢出等现象，得出较好的模拟结果。采用的稳定性条件为：

1.5 边界条件

在面波的数值模拟过程中，首先要考虑自由边界情况，面波沿自由表面传播，自由边界条件必然会对面波的传播有着一定影响。其次是对吸收边界条件的探讨，在地震勘探当中，地震波在一个非均匀各项异性半空间进行传播，而在计算机模拟时，不可给予如此大的计算区域，计算是有限的，这时需要引进吸收边界条件，从而减小吸收边界条件对正演模拟的影响。

采用理想匹配层（PML）吸收边界条件处理人工边界反射，通过在待计算区域的左右两侧分别加一个吸收边界层，将所有的入射波包括纵波和横波都进行吸收，不产生反射，这样在边界处有较好的吸收效果。在近地表研究各种波时，不得不面临浅层对地震波的影响，地震波激发后，地震波迅速到达浅层的反射界面，随后反弹至自由表面，然而在自由表面的处理相对不易，产生向下的反射波的同时，也将会产生不少的转换波。对此本文采用声学与弹性介质边界方法（AEA）。AEA方法操作简便，正演模拟效果良好，其二维处理如下[20]：

λ、θzz、ρ表示自由边界以上的拉梅常数、正应力及密度。图1为不同吸收边界厚度的波场快照对比。

图1 不同吸收边界厚度的波场快照对比，0层PML（a），10层PML（b），30层PML（c）Fig.1 Comparison of wave field snapshot with different absorption boundary thickness: 0 layer of PML （a）, 10 layers of PML （b）, 30 layers of PML （c）

2 模型试算

为了模拟近地表低降速带区域的波场特征，本文设计了几种地质模型，通过增加模型的复杂程度，模拟近似真实条件下近地表及低降速带区域的瑞雷面波波场特征。

2.1 近地表低降速层介质模型

图2为简化的近地表低降速层介质模型，模拟真实情况下地下低速层对地震勘探的影响。设计模型大小为200 m×200 m，主频采用30 Hz雷克子波，各层介质的参数如表1所示，模型x方向网格步长为0.5 m，z方向网格步长为0.5 m，炮点坐标选取为（100，0），模拟时间为100 ms，得到的全波场快照如图2所示。

表1 近地表低降速层介质模型参数Table 1 Model parameters of near surface low velocity layer

图2（a）为模拟得到的近地表低降速层介质模型波场快照，图2（b）为模拟得到的地震记录。在图2中，R表示面波，S表示横波，P表示纵波，PP表示反射纵波，SS表示反射横波。从图2（b）中可以看到各波同相轴，P波斜率较小，传播速度较快，S波次之，而R波传播速度较慢。从图2（a）中可以看出，由于每层介质参数不同，地震波在不同层交界面形成反射，各层的分层情况清晰可见。入射纵波和入射横波在穿过第一层界面后形成反射波后，能量随之减弱，形成的透射纵波和透射横波继续穿过第二层界面，形成二次反射波和透射波，由于第二层和第三层分界面波阻抗较强，以致深层地震波能量偏弱。面波在自由表面由纵横波互相间叠加干涉形成，在地表传播速度较低，能量极强且远大于P波。

2.2 倾斜断层介质模型

图3为简化的近地表斜面断层模型，在工程地质调查中常见，其特点是由于受到地质作用，两侧岩块沿断裂面发生位移而导致上下面介质参数不同。设计模型大小为200 m×200 m，主频采用30 Hz雷克子波，上层介质S1参数采用：Vp1=830 m/s，Vs1=490 m/s，ρ1=1.4 g/cm3，下层介质S2参数采用：Vp2=1 400 m/s，Vs2=770 m/s，ρ2=2.0 g/cm3，模型x方向网格步长为0.5 m，z方向网格步长为0.5 m，炮点坐标选取为（100，0），模拟时间为110 ms，得到的全波场快照如图4所示。

图3 倾斜断层介质模型示意图Fig.3 Schematic diagram of inclined fault medium model

图4（a）为模拟得到的倾斜断层介质模型波场快照，图4（b）为模拟得到的地震记录。在图4中，R表示面波，S表示横波，P表示纵波，PP表示反射纵波，SS表示反射横波。从图4（a）中可以明显看出倾斜断层面的存在，地震波传播到断裂面处时方向发生偏折，产生反射波和透射波，由于两侧速度差较大，形成一强波阻抗界面，浅层面波能量最强，深层地 震波能量较少，瑞雷面波在遇到断裂面后部分能量被反弹，部分能量发生透射。从图4（b）中可以看出同相轴不以炮点为对称轴左右对称，由于地震波在通过截面时速度减小，同相轴产生一定的角度。

图4 倾斜断层介质模型波场快照（a）及地震记录（b）Fig.4 Wave field snapshot （a） and seismic record （b） of inclined fault medium model

2.3 垂直地质异常体介质模型

在地质调查中，常见较大尺度的地质异常构造现象。图5为简化的垂直地质异常体介质模型，设计模型大小为200 m×200 m，主频采用30Hz雷克子波，上层介质S1参数采用：Vp1=760 m/s，Vs1=450 m/s，ρ1=1.5 g/cm3，下层介质S2参数采用：Vp2=1 600 m/s，Vs2=900 m/s，ρ2=2.1 g/cm3，模型x方向网格步长为0.5 m，z方向网格步长为0.5 m，炮点坐标选取（100，0），模拟时间为100 ms，得到的全波场快照如图6。

图5 垂直地质异常体介质模型示意图Fig.5 Schematic diagram of vertical geological anomaly mode

图6 垂直地质异常体介质模型波场快照（a）及地震记录（b）Fig.6 Snapshot of wave field （a） and seismic record （b） of vertical geological anomaly model

图6（a）为模拟得到的垂直地质异常体介质模型波场快照，图6（b）为模拟得到的地震记录。在图6中，R表示面波，S表示横波，P表示纵波，R’表示反射面波。从图6（a）中可以看出，在地震波传播过程中，左右波形未与炮点为中心左右对称，左侧出现较多能量团聚集，并局限在一个区域内，可以确认为一处地质异常体。在异常体附近分布反射波以及拐点处产生的绕射波，各种波互相叠加干涉，震相较难准确识别。瑞雷面波在通过地质异常体时，部分能量发生反射，形成反射瑞雷面波，瑞雷面波在异常体内震荡反射，在遇到接触面后小部分能量发生透射，大部分能量继续进行反射和向外围散发，瑞雷面波在通过地质异常体后，能量发生较大衰减，传播距离落后于均匀介质。在图6（b）中可以看出，同相轴不以炮点为对称轴左右对称，发生折叠现象，能量在垂直异常体内震荡反射聚集，在两侧形成反射波和透射波。

3 结论

通过交错网格有限差分法对近地表低降速层介质模型、倾斜断层介质模型及垂直地质异常体介质模型进行了正演模拟研究，高阶差分格式可以明显削弱数值频散的影响。通过对模拟结果分析发现，近地表低降速带区域相较工程层面尺度较大，地质情况复杂，地震波场丰富。在对近地表低降速层介质模型模拟时，分层情况清晰可见，波形图完整，而真实近地表地质复杂，倾斜断层及地质异常体广泛分布，对波场能量分布影响较大。瑞雷面波在通过复杂地形时，会发生反射、透射现象，能量重新分配，能量团聚集，在地震记录上出现同相轴的多次偏折。因此，利用对低降速带复杂地形的数值模拟，分析了近地表构造对波场的影响以及瑞雷面波的响应机理和传播规律，结果对实际地质勘探具有指导意义。