不均匀风化花岗岩地层模型被动源地震模拟

高 启 朱 伟 郝 龙 邓小虎 许 可 杨俊杰

（1.长江大学地球物理与石油资源学院，湖北 武汉 430100；2.武汉拓盟能源科技有限公司，湖北武汉 430074；3.武汉地震工程研究院有限公司，湖北 武汉 430071）

0 引言

瑞雷面波勘探是工程物探中的常用方法。在城市、铁路和公路附近，地表振动强烈，主动源面波勘探可能难以应用，而被动源面波勘探利用场地背景噪声中的面波，数据采集的环境要求低，是当前工程物探领域重点研究的方法［1-4］。

20 世纪50 年代，国外学者开始研究被动源面波记录探测地层结构［5］。在国内，被动源勘探研究始于20 世纪60 年代，2000 年以后应用于探测地下横波速度结构［6-7］。李翀等［8］利用被动源面波勘探方法划分南昌某地区地层结构，确定了古河道的大致位置。刘国峰［9］利用被动源面波勘探对内蒙古浅覆盖区的矿区厚度进行划分。邵广周［10］采用主、被动源面波联合勘探，准确探测黄土覆盖区的分层结构。基于模型的正、反演研究对被动源数据的采集、处理、反演和解释具有重要的指导作用，赵东［11］证明空间自相关方法提取频散曲线与实际频散曲线拟合度更高。白帅［12］对比分析不同速度模型的模拟结果，为被动源地震成像提供了有力的理论依据。

花岗岩在我国东南地区广泛出露,与原岩相比，风化花岗岩的力学性质变化较大，增加了市政和道路工程施工的难度。花岗岩的不均匀风化表现出不同的形式,孤石是一种常见的不均匀风化产物，对地铁盾构施工的影响较大［13-15］。对山区道路桥梁、隧道和边坡施工而言，探测不同风化程度的花岗岩的分布具有重要意义。

在道路改扩建工程中，场地情况复杂，背景噪声大，花岗岩风化程度的探测与评价是工程物探的难题。目前，被动源勘探在这方面的应用相对较少。本研究根据道路工程勘探和建设施工中的实际情况，设计花岗岩不均匀风化模型，模拟生成被动源地震记录，分析不同模型的地震记录的差异。

1 波动方程正演原理

在具有线状分布的随机强震动情况下，被动源勘探有可能采用与主动源勘探相似的线性观测系统，可用二维模型模拟被动源地震记录。假设二维模型的上表面为地面，设置自由表面边界条件［16］；在其他边界应用完全匹配层边界条件［17］。波动方程采用一阶速度—应力格式的弹性波方法，数值模拟方法为交错网格有限差分法。数量众多的随机震源设置在地表附近。

1.1 弹性波方程

一阶速度—应力格式的弹性波方程如公式（1）［18-19］。

式中：ρ为密度；vx、vz表示速度；σxx、σxz、σzz表示应力；f x、f z表示体力量；λ、μ为拉梅系数。

1.2 边界条件

模型边界条件的设置如图1 所示。假设模型的上表面为地表，采用自由表面边界条件。其他三个边界设置完全匹配层（PML）边界条件。

图1 边界条件示意

1.2.1 吸收边界条件。PML 边界条件是当前弹性波数值模拟中应用最多、效果最好的吸收边界条件［20-21］。通过设置合理的参数，它可以衰减绝大部分入射波的能量。应用PML 边界条件时，对波场变量沿坐标轴方向分解，在边界层内利用阻尼衰减相应方向的波场变量的分量。阻尼系数的表达式为式（2）［22］。

式中：R为理论反射系数；δ为吸收层厚度，取值越大吸收效果越好；x为吸收层网格点到边界的距离。

1.2.2 自由表面边界条件。自由表面边界条件的方法有多种，目前模拟精度较好的是声学弹性边界近似法。其基本思想是利用声学—弹性边界近似代替自由表面，令正应力σzz在自由表面处直接为零，同时还考虑自由表面上下横向应力保持连续的条件。在二维数值模拟时的处理方式为式（3）［21］。

式中：σzz、ρ分别为自由表面上的垂向正应力、密度；λ、μ为拉梅常数；ρ1和μ1分别为自由表面下介质的密度和拉梅常数。

1.3 交错网格有限差分法

交错网格有限差分法由Vireux 在进行SH 和PSV 波正演时建立，能稳定地应用在非均匀模型中，对流固边界的适应性也较好［23-24］。交错网格有限差分法可以与PML 边界条件和自由表面边界条件相结合，实现在有限大小的模型中模拟瑞雷面波的传播。交错网格有限差分示意如图2 所示，其主要特点是应力和位移分布在网格的不同位置，且能满足波场变量差分计算的要求。

图2 交错网格有限差分示意

1.4 震源

被动源面波勘探常采用傅里叶变换法求取地震记录的频波谱，提取频散曲线，反演横波速度。这说明时间域的波形不是重要影响因素。随机噪声在地震记录中的波形可能没有特定的形态。因此，震源子波可以用雷克子波表示为式（4）。

式中：A为振幅；f m为主频；td为延迟时间。

模拟被动源记录，震源设置较为关键。震源是一组随机震源的组合。各个震源的加载位置、振幅、主频和延迟时间都应具有一定的随机性。

1.5 均匀模型测试

本研究建立了一个大小为1 000 m×600 m 的均匀模型测试模拟方法。地层的纵波速度为1 000 m/s，横波速度为400 m/s，密度为2 000 kg/m3。网格大小0.6 m，时间步长0.08 ms。单一震源设置在模型中部地表附近。50 ms时刻垂直速度分量的波场快照如图3所示，可以清楚地观测到瑞雷面波、横波和纵波。

图3 50 ms时刻垂直速度分量的波场快照

2 不均匀风化花岗岩地层模型的被动源面波正演与分析

2.1 地层模型

本研究构建了三个不均匀风化花岗岩地层模型，如图4 所示。三个模型的大小均为3 000 m×600 m，网格大小均为0.6 m。花岗岩的风化程度分为全风化、强风化和中风化三类。模型一［图4（a）］的第一层为全风化层，厚度为50 m；第二层为强风化层，厚度100 m；第三层为中风化层，厚度为450 m。模型二［图4（b）］的第一层为全风化层，厚度为20 m；第二层为强风化层，厚度为80 m；第三层为中风化层，厚度为500 m。模型三［图4（c）］的第一层为全风化层，厚度为20 m；第二层为强风化层，厚度为480 m；在强风化层上部残留一个球形的中风化花岗岩体，垂向厚度约为200 m，顶部与全风化层相连。各层花岗岩的纵波速度、横波速度和密度见表1。在工程地震勘探中，检波器排列长度一般较短。本研究模型的长度远远大于排列长度，主要是为了布置数量众多的震源。

表1 不同风化程度花岗岩的纵波速度、横波速度和密度

图4 三个不均匀风化花岗岩模型

2.2 随机震源

被动源地震模拟需要大量震源。这些震源激发的空间位置，频率、出现时间和振幅须满足一定的随机性。在有限的空间和时间内模拟地震波传播并记录，必须对这些参数的分布范围进行控制。5 000 个随机震源的参数统计结果如图5所示。由图5可知，震源的激发位置在水平方向分布较均匀，在垂直方向分布在地表以下50 m 范围内；震源主频相对较低；延迟时间分布范围约为250～750 ms。这5 000 个随机震源在模拟之前确定，因此在三个模型的模拟中，震源参数是一致的。

图5 震源参数直方图

2.3 波场快照

波场快照能够展示地震波在模型中的传播特征，对观测系统的布置具有重要的指导意义。三个花岗岩模型在50 ms 时刻的垂直速度分量的波场快照如图6所示。由图6可知，在不同风化程度的花岗岩地层和岩体中地震波场具有明显不同。这说明地层速度的变化，对波传播具有重要影响。

图6 三个模型在50 ms时刻的垂直速度分量的波场快照

2.4 时间记录

在地表附近记录速度垂直分量形成时间记录。在数值模拟中，地表每一个网格点都可以是接收点，从而形成高密度的地震记录。三个模型垂直速度分量的时间记录如图7所示，由于模型长度达3 000 m，记录道数多，图中按一定的间隔显示了模型中部约2 000 m 范围内的记录。由图7 可知，不同的花岗岩地层模型，时间记录存在差异，但在波形图中并不明显。三个模型记录两两相减结果的波形如图8 所示，表明三个时间记录存在明显差异。

图7 三个模型垂直速度分量的时间记录

图8 三个模型时间记录的两两差值

3 结语

在我国东南部山地丘陵地区，浅层风化花岗岩地质体的勘探是道路工程建设中的重要工作内容。在道路改扩建工程等背景噪声强烈的地区，被动源面波勘探可能是一种有效的方法。本研究构建了三个不均匀风化的花岗岩地层二维模型，设计随机震源，采用弹性波方程正演，获得了被动源地震记录和波场快照。当模型结构发生变化时，波场快照和地震记录均存在明显的变化，表明被动源记录对模型的响应特征不同，能够反映不同模型的特征。