1 引言

复杂地表条件地区地震资料的采集和处理效果在很大程度上取决于地形起伏、地面地质风化、地表岩性变化和表层地质结构的复杂程度.其中，地形起伏和表层地质结构导致严重的静校正问题，同时还引起强散射噪音，常常干涉并掩盖来自深层的有效地震反射信号，严重降低地震资料的信噪比；地面地质风化和非均匀地表岩性变化波散化近地表散射波和深部反射信号，波形破碎带严重，半随机半相干的近地表散射噪音弥漫整个剖面.因此，精确模拟复杂近地表条件下地震波场特征和传播的规律对复杂地区地震资料采集、处理和解释具有重要的理论价值和实际意义.

对勘探地震而言，复杂近地表主要由两种介质构成.第一种是边界介质，包括起伏地形和表层地质结构，往往具有很强的速度反差，导致强边界散射(反射)，形成炮集上的主要能量分布；第二种是体积介质，包括地表下的非均匀低降速带介质和近地表结构内部的非均匀介质，形成体散射(反射)波，干涉和波散化边界反射波，导致波形破碎和衰减.精确模拟这两种介质构成对地震波传播的影响是本文的研究内容.复杂近地表地震波传播数值模拟常用的有限差分法、伪谱法和有限元法^[1]，将边界散射和体散射融为一体进行模拟.有限差分法和伪谱法中自由表面需要过于精细的网格化处理^[2].自由表面坐标变换的伪谱法模拟^{[3, 4]}和有限差分法模拟^[5]稳定性受到地形起伏程度的制约.有限元方法可以通过复杂的几何关系和任意大小的单元对起伏自由边界做灵活的处理^[6].另外，上述常用模拟方法的共同特点是通过网格化隐式地使用边界连续条件，对不规则界面上波的反射/透射进行全数值模拟时难以达到足够高的精度，也难以分离来自不同类型介质的波场.

精确描述起伏地形、显式应用边界连续条件的起伏地表数值模拟方法得到了广泛的研究，包括离散波数边界积分方程法^{[7, 8]}、边界积分方程广义反射/透射矩阵方法^{[9, 10]}、直接边界元地震模拟方法^{[11, 12]}、间接边界元地震模拟方法^{[13, 14]}.这些边界离散方法在几何上精确描述起伏地形和近地表构造的不规则边界，能够精确模拟复杂近地表中的边界介质；显式应用近地表地层界面的连续边界条件，实现半解析的地震波数值模拟；分区处理近地表复杂结构，可分别模拟来自复杂近地表不同部位的波.但这些只对近地表结构边界进行离散的方法不能模拟体积介质，即地表下的非均匀低降速带和近地表结构内部非均匀介质.边界元、有限差分和有限元混合方法^[15]可以模拟复杂近地表边界+体积介质模型，但存在三种方法之间的数值衔接问题.

本文将复杂近地表中边界介质与体积介质分离，分别模拟地震波在复杂近地表中传播的边界散射和体散射，这需要考虑复杂近地表中边界-体积之间的相互耦合作用.Fu等^[16]提出了一种求解Lipmann- Schwinger波动积分方程的体积元数值方法，可以精确模拟非均匀体积介质产生的强散射波.通过将边界积分方程与Lipmann-Schwinger体积分方程相结合，本文研究一种模拟复杂近地表地震波传播的边界元-体积元波动方程数值模拟方法.其中，导致强边界散射(反射)的起伏地形和表层地质结构等边界介质用边界元地震模拟方法；而产生体散射(反射)波的地表下非均匀低降速带介质和近地表结构内部的非均匀介质用体积元法模拟.在边界元-体积元数值模拟方法中，与体积分有关的散射项采用了高精度的高斯积分数值计算方法，确保对随机强非均匀体介质散射的地震模拟精度.这种方法除了具有上述边界元地震模拟方法的优点外，还能够模拟带起伏地表的分区非均匀地质结构，有效模拟复杂地表下非均匀介质对地震波场的体散射.通过几种不同复杂地表模型的数值模拟计算，本文研究了复杂近地表条件对地震波传播的影响.

2 方法原理与基本方程

广义Lipmann-Schwinger积分方程由边界波动积分方程与Lipmann-Schwinger体积分方程两部分构成，描述地震波在分区非均匀介质(一种具有不规则边界、内部为随机非均匀介质的多层结构)中的传播.该方程的解由源波场、边界散射波场和体散射波场三项构成.图 1所示为起伏地形及近地表非均质低降速层，该散射介质模型Ω 由起伏的自由表面Γ₀、不规则的底界面Γ₁ 和人工截断边界Γ_∞ 形成闭合区域，即Γ=Γ₀+Γ₁+Γ_∞ ，背景速度为v₀，内部速度异常分布为v(r).假设震源位于r₀ 处，震源项在频率域可表示为:

(1)

式中S(ω)为震源谱，δ(r_-r₀)为脉冲函数.

该散射介质中地震波的传播满足如下的Helmholtz稳态标量波动方程:

(2)

式中k=ω/v(r)为波数.地震波场u(r)可表示为入射波场u⁰(r)、边界散射u₁^s(r)和体积散射u₂^s(r)三项之和，即

(3)

其中第一项u₀(r)为背景介质中的入射波场，考虑点源激发，则u₀(r)可表示为

(4)

第二项u₁^s(r)为起伏地形及近地表结构中的不规则边界贡献的边界散射波场，满足如下均匀背景介质中的边界积分方程:

(5)

式中，∂/∂n为边界法向导数.第三项u₂^s(r)为非均质低降速层产生的体散射波场，满足如下的Lipmann- Schwinger体积分方程:

(6)

式中k₀ =ω/v₀ 为背景波数，O(r′)=v₀²/v²(r′)-1为速度扰动，O(r′)=0时则退化成均匀背景介质.G(r，r′)为背景介质中的格林函数，对于二维介质，G(r，r′)=iH₀⁽¹⁾ (k₀|r_-r′|)/4，其中H₀⁽¹⁾ 为汉克尔函数；对于三维介质，G(r，r′)=e^{ik₀|r_-r′|}/(4π|r_-r′|).

将式(4)~(6)代入方程(3)，即可得到广义 Lipmann-Schwinger积分方程:

(7)

式中Ω=Ω +Γ，C(r)=$\frac{\theta }{2\pi }$，θ为角点P(r)对区域Ω 的张角(如图 2所示).

广义Lippmann-Schwinger 积分方程(7)是边界元-体积元地震模拟方法的核心方程，属于第二类 Fredholm 积分方程.根据Fredholm 积分方程理论，数学上可以证明，在Neumann 和Dirichlet边界条件下，上式对于内边值问题(r′，r₀ ∈Ω)，其解存在，且惟一.由于格林函数在r￫r′ 时出现奇异积分，虽然数学上该奇异点为可去奇点^[16]，但为了确保数值计算的稳定性，需要进行奇异积分处理，可采用格林函数的离散波数展开^[7]或奇异积分在奇异点展开的解析处理方法^[12].广义Lippmann-Schwinger 积分方程自然满足无穷远处边界条件，即位移场，由于基本解为Hankel函数或e^ikr/r的形式，而在无穷远处趋于0.对于模型两端的人工截断，可采用无限元法处理人工截断边界的吸收问题^[17].

3 积分方程的数值离散

广义Lippmann-Schwinger 积分方程中，震源项u⁰(r)为已知，由计算频率和源点r₀ 确定，u₁^s(r)和u₂^s(r)分别为散射体边界和内部各点上的未知量.积分方程离散常用的方法为配置法，以二维问题为例，首先把边界Γ 划分为L个互不重叠的线性单元Γ_e(e=1，2，…，L)，域Ω 剖分为M个三角形或四边形单元Ω_e(e=1，2，…，M).假定边界加上内部共有N个节点，并统一节点号.对于各个变量(r，u，∂u/∂n)在单元内任意点的值，利用插值型函数Φ ，可以近似地表示成相邻节点I₁ 至I₂ 的线性和，例如位移函数可表示为:

(8)

式中ξ 为一个单元的局部坐标.这样，方程(7)中的各项积分可写为:

(9)

(10)

(11)

其中，δ_lj和δ_ij为Kronecker脉冲函数.上述离散单元积分可利用高精度的高斯积分数值算法计算.方程(7)可改写为算子的形式:

(12)

其中，u⁽¹⁾ 为边界Γ 上的地震响应，t为u⁽¹⁾ 的法向导数，u⁽²⁾ 为Ω 域内的地震响应，f为入射波场，H、G为边界积分算子，K为体积分算子.

对于如图 1所示表层含有非均匀介质的复杂近地表结构，设Γ₀ 的离散点数为N₀，不规则底界面Γ₁ 的离散点数为N₁，非均质区域Γ 内的离散点数为NP.对于自由边界Γ₀ 上的离散点，由自由边界条件可知t(r)=0.设定地层边界Γ₁ 上面介质的密度为ρ₀，速度为v₀，下面介质的密度为ρ₁，速度为v₁.根据边界连续条件，边界Γ₁ 上面与下面的位移、应力连续，即u₊ (r)=u_- (r)，μ₁u₊ (r)= -μ₀u_- (r).其中μ₀ =ρ₀v₀²，μ₁ =ρ₁v₁².利用方程(12)，对所有节点建立联立方程组，得到如下总体矩阵方程组:

(13)

式中左边的各个子矩阵分别为利用公式(9)、(10)、(11)计算得到的边界Γ₀ 上的各个离散点和边界Γ₁上的各个离散点以及非均质区域内的各个散射点的相互作用积分算子矩阵；式中右边的各个子向量分别为利用公式(4)计算得到的点源激发对边界和边界Γ₁ 上的各个离散点以及非均质区域内的各个散射点的作用产生的震源项向量.

例如H₀₀为利用公式(9)计算得到的边界Γ₀ 上的各个离散点对其上的其他离散点的相互作用积分算子矩阵，即

(14)

F₀ 为利用公式(4)计算得到的点源激发入射到边界Γ₀ 上的各个离散点作用而产生的震源项向量，即

(15)

求解公式(13)的总体系数矩阵，就可以得到公式(13)左边未知数向量对应的边界和内部各点的波场值.对于多层结构的复杂近地表结构，需要对各层计算的系数矩阵H，G，K按照应力和位移边界连续条件进行总体系数矩阵组装，最后求解总体系数矩阵方程，得到所有层边界和内部各点的波场值.

方程(13)包含了穿过整个模型的地震波的信息.对方程(13)进行求解，回代入公式(7)就可求得任一接收点处的地震响应波场值，这个值是在频率域给出的.设接收点的个数为N_R，则接收点单个频率的波场值向量U_R 可由下面的公式给出:

(16)

式中H_R0、H_RP、H_R1、G_R1 分别为利用公式(9)、(10)、(11)计算得到的各个接收点与自由边界Γ₀ 上的各个离散点和边界Γ₁ 上的各个离散点以及非均质区域内的各个散射点的相互作用积分算子矩阵.公式(16)中去掉[F_R]项就可以得到没有直达波的波场.通过对接收点的不同频率的波场值进行逆FFT 变换，即可得到所要的合成地震记录.

利用边界-体积分波动方程进行复杂地表地震数值模拟的具体步骤为:

(1) 针对具体模型，建立包含震源项的边界-体积分方程.

(2) 选用不同的离散单元与形函数，对边界积分方程进行边界元离散，用配置法对体积分方程进行离散.

(3) 对由(2)式得到的单元系数矩阵组装成总体系数矩阵.

(4) 求解总体系数矩阵并回代入对应的边界-体积分方程，得到频率域的波场值.

(5) 对不同频率的波场值进行逆FFT，变换到时间域，得到接收点处的合成地震记录.

在实际的计算过程中，由于地震勘探中震源位于近地表且仅需要计算自由表面的地震响应，根据这些特点，采用由深至浅逐层进行式(13)中的系数矩阵H、G和K的数值计算、组装和消元，消去两个相邻子域公共边界上的耦合数据，如此计算直到地表层，得到最后的矩阵方程，其维数仅由地表层的节点数确定.此矩阵隐含着穿过整个模型的地震波的信息.将此矩阵与地表多测线多震源相结合，可计算出地表任意处的观测波场.为提高算法的效率，我们采用了改进的分块高斯消去法来求解矩阵方程.在整个计算过程中，全局系数矩阵需要的最大内存限于最大一个子域的总节点数.这样既节省了内存又缩短了计算时间.

4 数值解与解析解的对比分析

本文对如图 3所示的半圆峡谷和沉积谷地两个半空间模型进行了数值解与解析解的比较.由于这两个模型形态在自然界中的普遍性，常常被作为两个经典例子来研究局部地形和不规则近地表结构对地震波传播的影响.两个模型的半径均为a，半空间背景介质的SH-波速度为β，密度为ρ，由于模型的对称性，其SH-波入射的解析解^{[18, 19]}可表示为Bessel函数的线性和.模型两端(x=±a)90°的拐角点是测试各种地震数值模拟方法和算法精度的最佳靶点.

定义无量纲化频率η 为半圆峡谷或沉积谷直径与入射波长的比值，即

(17)

式中，λ 为入射波波长.在半径a不变的情况下，η越小，表示入射波的频率越小，反之亦然.

图 4为半圆峡谷模型的地震响应，图 4a为60° SH-波入射地表各点对不同η 值的振幅响应.可见地形的突变会引起地震波响应的巨大变化，数值解(虚线)与解析解(实线)之间吻合很好，微小的误差源于半圆弧的单元逼近误差.图 4b为SH-波垂直入射地表各点的地震合成记录，震源为Ricker 子波，计算频率为0~4Hz, 可见半圆地形两端垂直拐角处波场的强烈干涉作用.

在均匀沉积谷模型中，假定谷内介质SH-波速度β和密度为ρ 满足:ρ/ρ=1.5和β/β=2，图 5a为 SH-波垂直入射地表各点对不同η 值的振幅响应，数值解(虚线)与解析解(实线)基本一致，说明方法的正确性.图 5b为SH-波垂直入射地表各点的地震合成记录，震源为Ricker子波，计算频率为0~4Hz, 半圆谷底的反射和两端垂直拐角点散射波清晰可见.

离散单元的尺寸大小是影响数值模拟精度和效率的主要因素，下面以沉积谷模型为例分析离散单元长度的选择.采用不同离散单元长度对均匀沉积谷模型进行数值模拟，取η = 1 时地表各点的地震振幅响应曲线如图 6a所示，实线为解析解，虚线为数值解.可见，在SH-波垂直入射情况下，每个波长用两个离散单元的结果与解析解在半圆的垂直拐角处稍有偏差，用三个离散单元就能达到很好吻合.增大入射角度，地面地震波振幅响应变化剧烈，图 6b模拟结果表明每个波长需要用五个离散单元才能达到与解析解较好的吻合.时间域的合成地震记录对离散单元长度的变化相对于上述频率域的振幅响应而言显示不敏感，因此，合成地震记录对离散单元尺寸的要求要弱些，通过大量的地震数值模拟表明，每个波长用三个离散单元基本上能满足实际需要.

我们用边界元-体积元法模拟沉积谷中填充随机非均匀介质的地震响应频率响应.图 7 为非均匀沉积谷模型，谷内随机非均匀介质由三个参数确定: SH-波平均速度β、平均密度ρ 和速度扰动量δ.

对大部分近地表而言，地表下低降速带的速度扰动量变化一般为5%~15%之间.例如，对谷内平均速度β=3000m/s进行δ=10% 的扰动，实际谷内随机速度分布在2700~3300m/s之间.假定谷内介质平均速度β和平均密度为ρ 满足:ρ/ρ=1.5和β/β=2，取速度扰动量δ 为0%、5%、10%和15%，对随机非均质沉积谷模型进行数值模拟，图 8a为 SH-波垂直入射地表各点对η =1 时振幅响应.可见，对于δ>10%后，谷内随机非均匀体积介质对地面地震响应产生较大的影响，波场振幅扰动非常强烈，特别是在半圆两端垂直拐角点的散射波场变化剧烈.图 8b为取δ=10%时，SH-波垂直入射地表各点对不同η 值的振幅响应，可见，谷内随机非均匀体介质对高频地震响应产生更大的影响.

上述频率域地震振幅响应的对比研究表明，本文的方法对近地表带自由表面的复杂结构进行地震波传播模拟具有很高的模拟精度，对含有地层尖灭、断层、陡峭拐角等复杂地质构造具有很好的适应能力.下面我们用边界元-体积元地震数值模拟方法来模拟地震波在复杂近地表中的传播.

5 数值算例

对于实际地震勘探来讲，震源一般放置在自由界面下10~30m, 检波器在自由界面下十几厘米至几十厘米，炮检距一般为30 m~3km.这样的地震数据采集观测设置导致了起伏地表的影响很大程度上与震源和检波器附近的地表几何形状、近地表结构和低降速带介质的非均质性密切相关，特别是近地表散射体尺度与地震主波长相当时会引起较强的散射.

我们首先对简单的、起伏程度不同的复杂近地表模型进行数值模拟，以便直观了解复杂近地表对地震响应的影响模式和影响程度，如图 9所示.地震模拟计算采用Ricker子波，频率范围0~40 Hz, 震源位于自由界面以下20m.为了直观展示起伏地表的散射影响，模拟计算时滤除了直达波.由图可见，平地表合成记录地震反射规则清晰，两个山包形状的地表起伏引起很强的散射波分布，特别是山根下波场干涉严重，但地表散射和地层反射波都具有较好的相干性.相比之下，当近地表层介质为随机非均匀时，非均匀体介质产生的漫散射及其对边界波的干涉使整个记录波散化，波前破碎.增加近地表层随机非均匀介质的速度扰动到δ=15% 时，相干散射波场进一步随机化，变为半随机半相干的散射噪音，弥漫整个剖面，干涉和掩盖地层反射波.

图 10a为一个包含复杂地表的盐丘模型，速度分布如图所示.模拟该模型的主要目的是对来自盐丘的强反射能量与复杂地表强散射能量进行比较.震源位于水平距离1350m, 自由表面下25m 处，数值模拟得到的自由表面接收合成地震记录如图 10b所示.由图可见，地表强散射噪音能量很大，弥漫整个剖面，具有一定的相干性，干涉和淹没来自地表下各地层和盐丘的反射信号.图 10c为把起伏地表变平进行数值模拟得到的合成记录，来自盐丘的强反射在此炮集上很突出.由此可见，剧烈的起伏地表是降低地震资料品质的主要原因之一.

图 11是针对中国西部某地区起伏地表模型进行地震数值模拟得到的合成地震炮集记录.本区地表起伏剧烈、高差近400 m, 低降速带横向变化较大，对地震采集数据品质造成严重影响，主要表现为采集的地震资料信噪比低.地震模拟计算采用 Ricker子波，频率范围0~40 Hz, 3 个震源位置如图所示，位于自由表面以下25m, 图 11(b~d)为自由表面接收合成的3个炮集记录.由图可见，近地表散射分布受地表起伏和低降速带横向变化的影响，形成很强的、延续至中深层的散射噪音环境，破坏来自地下深部的反射.另外，崎岖地表和相关的近地表低速异常带为来自地下深部的反射波提供一个再次散射的平台，进一步增强剖面上中深层的散射噪声.

6 结论

复杂近地表边界元-体积元地震模拟方法将边界积分方程与Lipmann-Schwinger体积分方程相结合，边界元法模拟起伏地形和表层地质结构边界对地震波传播的影响，体积元法模拟起伏地表下非均质低降速层对地震波的体散射影响.该方法的主要优点为几何上精确描述含崎岖地表在内的近地表结构中各种不规则界面，实现精确模拟复杂近地表对地震波传播的边界散射；显式应用近地表地层界面的连续边界条件，实施半解析-半数值的地震模拟；有效模拟复杂地表下随机非均匀介质对地震波场的体散射，实现边界散射与体散射的有机结合；可以分别模拟来自近地表不同地质体的反射(散射)波，实现直观表达地质元素与地震波场构成的相互关系.

利用边界元-体积元地震模拟方法对复杂近地表模型进行的地震数值模拟结果表明，起伏地表和低降速带的横向变化是影响地震资料品质的主要原因，主要表现为复杂近地表引起下伏构造反射走时畸变；崎岖自由表面散射形成很强的相干性散射噪音，干涉并掩盖来自下伏地层的反射波；低降速带随机非均匀介质引起的体散射波散化整个波场，形成很强的半随机半相干近地表散射噪音环境，严重降低采集地震资料的信噪比.