基于二维有限差分模拟的非均质近地表地震波散射分析

简介：面向地震学与计算地球物理方向学习者的二维有限差分模拟资料包，聚焦近地表非均质介质中地震波散射这一经典问题。非均匀的岩石成分、孔隙结构与密度分布会引发波场复杂散射与能量重分配，资料配套学术论文、参考文献与可运行Python脚本，适合入门数值模拟并理解波场传播机制。压缩包共7个文件，整体仅16KB，以md说明文档、bib文献、py脚本及建模代码目录为主，另含备份zip，覆盖从理论推导到程序实现的核心链条。目前已有39人浏览学习，适合个人自学或课程对照。借助论文与代码可掌握网格离散、震源设置、散射波频谱及能量衰减分析思路；配套说明文档与论文有助于快速理清模型假设，为地震危险性评价与勘探资料解释提供模拟参考。

1. 近地表非均质性下地震波散射几乎无处不在，但常规模拟常常假装它是均匀的

风化壳、冲洪积扇、沙漠沙丘和冻土表层并不讲理，几十米的碎石层里，速度扰动往往达到背景值的 10%～20%，尺度从厘米级到数十米级都有。入射波撞上这些随机分布的速度异常，会产生前向散射和后向散射为主的强尾波，地震记录上的表现就是高频衰减、波形畸变和道间不相关噪声突然增大。二维有限差分模拟是理解这个现象最直接的数值实验手段：把近地表速度场用一张二维数组存下来，给定震源子波，按时间步递推更新每个网格点的波场值，就能同时得到散射波场的传播快照和合成炮集。它的价值不在于替代三维计算，而在于降维之后仍保留非均质体的空间相干特征，参数扫描效率和可解释性都更可控。下面从怎么描述非均质介质开始，把整套流程拆开讲。

2. 近地表非均质性的随机介质建模与地震波散射分类

2.1 判断散射机制先看波长与扰动尺度的相对关系

入射地震波遇到速度扰动体时，散射强弱并不只由扰动幅度决定，更关键的是扰动尺度与波长的比值。设背景波速为 v0、震源主频为 f0，有效波长近似为 λ = v0 / f0。把非均质体的自相关长度记作 a，工程中粗略按 a/λ 分档：比值远小于 0.1 时，单个扰动体散射很弱，但数量多、累积效应明显，高频成分衰减快，容易在资料处理时被误当成吸收衰减；比值落在 0.1 到 1 之间时，单个扰动体已经把能量重新辐射到各个方向，并与后续波叠加产生长尾波，这是近地表数据里最难压制的“不可聚焦能量”；比值超过 1 以后，波前面基本按几何路径走，射线近似已经够用。

真实的近地表极少是单一尺度。风化壳里既有毫米级矿物颗粒，也有米级风化碎块和数十米级剥蚀面，实际波场是多种散射机制叠加的结果。用随机介质模型描述时，不需要把每种尺度都显式建模，只要用自相关长度和扰动标准差把能量谱的集中区间框住，再选对指数型或高斯型自相关函数，就能复现大部分散射特征。

2.2 用自相关函数在波数域构造二维速度扰动场

速度模型写为：

v(x, z) = v0 * (1 + ε * r(x, z))

其中 r(x,z) 是均值为 0、标准差为 1 的随机场。扰动场的空间相关性由自相关函数决定：指数型 C(r) = exp(-|r|/a) 保留了更丰富的小尺度高频成分，适合强风化壳、破碎带；高斯型 C(r) = exp(-r²/a²) 能量集中在相关长度 a 附近，适合冰碛物、充填通道这类块状堆积体。

在二维有限差分模拟里，随机场必须先按网格点生成，再直接写进速度二维数组。用谱方法构造最快：先产生高斯白噪声，经过与目标自相关函数匹配的波数域滤波器，再反变换回空间域。这样生成的扰动场自动满足统计特征，不会因逐点插值产生人工纹理。下面这段代码可以生成 200 m × 200 m、网格间距 1 m 的指数型随机介质模型。

import numpy as np from numpy.fft import fft2, ifft2 def random_velocity(nx, nz, dx, dz, v0, eps, a, corr_type='exp', seed=2024): rng = np.random.default_rng(seed) noise = rng.standard_normal((nx, nz)) # 频率坐标，单位 rad/m kx = np.fft.fftfreq(nx, d=dx) * 2.0 * np.pi kz = np.fft.fftfreq(nz, d=dz) * 2.0 * np.pi KX, KZ = np.meshgrid(kx, kz, indexing='ij') k2 = KX**2 + KZ**2 if corr_type == 'exp': # 指数型：一阶低通，高频衰减较慢 filt = 1.0 / (1.0 + (a * np.sqrt(k2)) ** 2) ** 0.75 else: # 高斯型：高频快速衰减 filt = np.exp(-0.25 * (a ** 2) * k2) spec = fft2(noise) * filt r = np.real(ifft2(spec)) r = (r - r.mean()) / r.std() return v0 * (1.0 + eps * r)

代码里kx和kz分别对应模型 x 方向和 z 方向的波数，meshgrid展开后与fft2输出的频率索引一一对应，滤波后反变换得到的二维数组与模型网格坐标对齐，不会出现错位。a以米为单位，乘到波数上形成无量纲量，波数越高、滤波衰减越强，等效于把原始白噪声的小尺度分量压制下来。eps是相对扰动强度，a是相关长度，二者是后续有限差分模拟中最需要反复试的两个参数。

模型生成后要立刻验证两点：一是扰动幅度与v0*(1 ± eps)是否吻合，二是沿任意水平剖面算自相关，看相关系数降到 1/e 的距离是否与a相当。这两项检查过了，模型才不会在后面的波场递推里给出偏高甚至失真的散射能量。参数经验取值范围如下表。

ε 超过 0.2 时要特别检查背景速度 v0 是否足够高，避免v0*(1 - ε)接近零甚至为负。硬岩区 v0 取 3000 m/s 时，ε=0.2 对应最慢 2400 m/s，安全余量较大；浅表低速层 v0=400 m/s 时同样取 0.2，最慢只有 320 m/s，空间采样稍微不够就会出现明显网格频散，这时宁可降低 ε 或提高 v0。

3. 二维有限差分的交错网格格式与波场递推实现

3.1 四阶空间差分比二阶格式更省计算量

近地表介质速度横向变化剧烈，低阶差分会产生比物理散射更明显的网格散射：同一个波前在不同网格间距下算出的走时不一致，这种数值各向异性会污染散射尾波的形态，且无法通过提高震源信噪比消除。工程实现中最常用的组合是二阶时间精度加四阶空间精度，即在时间方向做中心差分，空间上用五个点的四阶中心差分算子。

对二维标量声波方程：

∂²p/∂t² = v(x,z)² · ∇²p

拉普拉斯项离散为：

∇²p(i,j) ≈ [ c0·p(i,j) + c1·( p(i+1,j) + p(i-1,j) + p(i,j+1) + p(i,j-1) ) + c2·( p(i+2,j) + p(i-2,j) + p(i,j+2) + p(i,j-2) ) ] / Δx²

差分系数固定如下表。

四阶格式的振幅误差和相位误差都比二阶小一个量级以上。从计算量看，要达到同样的波形保真度，二阶格式需要把网格间距减半，网格点数变成四倍，总计算量增加约一个量级；所以四阶格式虽然单步运算更复杂，整体仍然划算。这个结论在散射模拟里尤其重要，因为散射波本身振幅弱，数值频散造成的虚假尾波会直接掩盖真实散射信息。

3.2 波场递推核心代码与稳定性约束

在均匀网格 Δx = Δz 条件下，二阶时间、四阶空间格式的稳定性条件为：

Δt ≤ 0.612 · Δx / v_max

这个 0.612 来自 λ_max(∇²) 的离散特征值推导，实际工程安全值取 0.55～0.6。时间步长超过这个界限，波场会在高频段指数发散，表现是总能量随迭代步数单调增长，而不是单纯的波形失真。下面是波场递推的内核函数。

import numpy as np def fd2d_wavefield_step(p0, p1, v, dt, dx): """ 对二维声波方程做一步递推 p0: t-dt 时刻波场快照 p1: t 时刻波场快照 v : 速度场二维数组,单位 m/s 返回 p2: t+dt 时刻波场快照 """ nx, nz = p1.shape p2 = np.zeros_like(p1) c0, c1, c2 = -2.5, 4.0/3.0, -1.0/12.0 # 四阶拉普拉斯算子,只计算去掉边界带后的内部区域 lap = (c0 * p1[2:-2, 2:-2] + c1 * (p1[3:-1, 2:-2] + p1[1:-3, 2:-2] + p1[2:-2, 3:-1] + p1[2:-2, 1:-3]) + c2 * (p1[4:, 2:-2] + p1[:-4, 2:-2] + p1[2:-2, 4:] + p1[2:-2, :-4])) # 变系数: (v*dt/dx)^2 是逐点不同的 k = (v * dt / dx) ** 2 p2[2:-2, 2:-2] = (2.0 * p1[2:-2, 2:-2] - p0[2:-2, 2:-2] + k[2:-2, 2:-2] * lap) return p2

这段代码里所有切片都精确对齐了同样的大小：p1[3:-1]、p1[1:-3]、p1[4:]、p1[:-4]与内部区域p1[2:-2]的长度一致，因此lap不需要额外裁剪。k不是标量而是与v同形的二维数组，乘到lap上时自然实现了速度场的逐点变系数作用，这正是近地表非均质性进入方程的唯一路径。

边界处理在上面的函数里故意空了出来。对于散射研究，边界反射对记录前 200 ms 影响不大，但长尾波会与边界反射混在一起。最简单的工程做法是沿模型四周设置 20 到 40 个网格点的衰减带，每一步把波场乘以衰减因子；更好的方案是卷积完全匹配层，它对大角度入射波吸收更干净，不会像简单衰减带那样反射残余能量。

3.3 震源加载与时间循环的最小实现

震源子波常用 Ricker 子波：

s(t) = (1 - 2π²f0²(t-t0)²) · exp(-π²f0²(t-t0)²)

这个子波无零频成分，主频恰好为 f0，工程参数好控制。主循环里把子波振幅直接加到震源所在网格点的波场上，然后调波场递推、做边界吸收、把检波点位置的波场值写入记录道。

dx = 1.0 v = random_velocity(nx, nz, dx, dx, v0=1200.0, eps=0.15, a=3.0, seed=7) dt = 0.55 * dx / v.max() nt = int(0.8 / dt) # 总时长 0.8 s nb = 30 # 吸收带宽度,单位网格数 fade = np.ones((nx, nz), dtype=np.float32) for i in range(nb): w = (i + 1.0) / nb fade[i, :] *= w fade[-1 - i, :] *= w fade[:, i] *= w fade[:, -1 - i] *= w p0 = np.zeros((nx, nz), dtype=np.float32) p1 = np.zeros((nx, nz), dtype=np.float32) rec = np.zeros(nt, dtype=np.float32) f0 = 30.0 t0 = 1.0 / f0 isrc, jsrc = nx // 2, nz // 4 irx, jrx = isrc + 40, jsrc for it in range(nt): t = it * dt wav = (1.0 - 2.0 * (np.pi * f0 * (t - t0)) ** 2) * \ np.exp(-(np.pi * f0 * (t - t0)) ** 2) p1[isrc, jsrc] += wav p2 = fd2d_wavefield_step(p0, p1, v, dt, dx) p0 = p1.copy() p1 = p2.copy() # 简单吸收带,每步乘一次衰减因子 p1 *= fade rec[it] = p1[irx, jrx]

代码里的t0 = 1.0 / f0把子波峰值对齐到第一个主周期上，避免子波在零时刻就有明显非零值；wav的单位与波场振幅一致，实际工程中还要按源强度标定，但相对振幅分析不需要这一步。rec数组每一行记的是当前时刻检波点位置(irx, jrx)的波场值，等价于单道地震记录。

4. 二维有限差分模拟工作流里的参数联动与炮集处理

4.1 空间步长与时间步长要先满足两个不等式

网格步长不是拍脑袋定的。先看频带：Ricker 子波主频 f0 的有效频带上限按 fmax ≈ 2.5·f0 估算，空间步长要保证最慢速度对应的最小波长上至少有 10 个网格点：

dx ≤ v_min / (10 · fmax)

再看稳定性：时间步长必须满足 CFL 条件，而 CFL 条件里的 v 取最大值 v_max，否则高速层先发散：

dt ≤ 0.612 · dx / v_max

两个不等式是串行的，先由速度和主频定 dx，再由最大速度定 dt。下面是一组实际算例的参数设定。

算例里 v_min=400 m/s，f0=30 Hz 则 fmax=75 Hz，dx ≤ 0.533 m，取 0.5 m 是合理折中。接着 v_max=2000 m/s 时 dt=0.55×0.5/2000=0.1375 ms，取 0.1 ms 留出足够安全余量。801×801 网格的二维模拟在单张消费级显卡上跑毫秒级时间步，总时间约几分钟，这个量级才轮得到做参数扫描。

4.2 用波场快照和差分炮集分离散射能量

递推完成后，除了检波器位置的标量记录，还能在特定时间步把整个波场快照存成二维数组并可视化。判断散射强度最直观的做法是同一套采参数下，分别跑均匀背景模型和随机介质模型，再把两组波场快照或炮集相减。差值场就是速度扰动引起的散射响应，这样能把直达波和规则反射从画面里约掉，散射尾波的传播路径会非常清晰。

import matplotlib.pyplot as plt # 时间索引切换到 0.3 s 左右的快照 it_snap = int(0.30 / dt) fig, axes = plt.subplots(1, 3, figsize=(15, 4.5)) titles = ['total field', 'backgroud field', 'scattered field'] fields = [p_hetero, p_homo, p_hetero - p_homo] for ax, title, field in zip(axes, titles, fields): im = ax.imshow(field.T, cmap='seismic', aspect='auto', vmin=-5e-4, vmax=5e-4) ax.set_title(title) ax.set_xlabel('x (m)') ax.set_ylabel('z (m)') plt.colorbar(im, ax=ax) plt.tight_layout() plt.savefig('snapshot_compare.png', dpi=150)

三个子图并排观察时，主要看散射场子图是否出现与背景介质结构对应的扇形散射带。如果散射场只在震源附近出现一圈对称环，而没有延伸到远道，说明相关长度 a 偏小或 ε 偏低，散射能量没有进入有效传播路径。反之，若散射场亮斑铺满整个区域且色标饱和，则可能是数值频散在放大网格扰动，先加密网格再判断。

4.3 检波器布置与记录道里的“假散射”

近地表散射模拟很容易把数值噪声当成物理散射。判断标准很简单：物理散射能量在主波到达之后才开始显现，且随偏移距增大走时差逐渐增大；数值干扰从第一时刻起就出现在全频带，尤其在高频端与主波同时到达。

检波器布置要避开四个角点，那里简单吸收带处理不好仍会残留反射；道间距取 2～5 倍网格间距，也就是 1～2.5 m，既能覆盖空间波数又不至于数据量过大。单炮记录出来之后，先做 10～30 Hz 的带通滤波，再看散射尾波是否仍明显；如果带通后尾波消失，说明散射能量主要在算法引入的高频伪影里，需要检查吸收带宽度或网格频散。

5. 检验散射模拟可信度的三个实用动作

5.1 用均匀半空间走时卡时间步

把随机速度场换成常数 v=1200 m/s，其他参数不变，跑一组直达波记录。理论走时 t_theory = 源检距 / v0，与数值记录波峰到时误差应小于一个时间采样步。若误差偏大，优先怀疑dt取值接近稳定性上限，或 Ricker 子波峰值对齐方式有偏差。这一步能在 10 分钟内完成，值回票价。

5.2 用 2 倍和 4 倍加密网格做收敛性检查

把 dx 从 1.0 m 减到 0.5 m，其他物理参数不变，看同一检波器记录。四阶空间格式的差分误差应为 O(Δx⁴)，所以网格减半后同一地震道的残差能量应下降约 16 倍。实际操作时算一个比值即可：残差能量 = sum((rec_coarse - rec_fine)²) / sum(rec_fine²)。比值低于 0.05，说明当前网格已收敛；高于 0.2，说明散射模拟仍被网格频散主导，继续加密网格直到比值稳定。工程上不必追求绝对收敛，控制在 0.05 以内即可。

5.3 用尾波包络斜率核对非均质尺度

散射尾波的包络衰减斜率与自相关长度 a 直接相关。a 小则尾波衰减快，a 大则尾波持续时间长。把不同 a 值代入模型，画出尾波包络对数衰减曲线，如果 a 从 1 m 增加到 5 m 而包络几乎不变，多半是模型里最大扰动量没坐落在目标频段，需要增加小尺度分量或减小 a 对应的谱截止。也可以用二维网格搜索扫过 a-ε 参数平面，以尾波包络误差最小为目标，自动选出一组模型参数。

我建议把这三项当成每次改模型的固定前置检查：走时卡对、网格收敛比值合格、尾波趋势符合尺度预期，再谈散射衰减提取和参数反演。这样二维有限差分模拟提供的不是几张漂亮的波场图，而是可以放进定量分析流程的可靠合成数据。

本文还有配套的精品资源，点击获取