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

简介面向地震学与计算地球物理方向学习者的二维有限差分模拟资料包聚焦近地表非均质介质中地震波散射这一经典问题。非均匀的岩石成分、孔隙结构与密度分布会引发波场复杂散射与能量重分配资料配套学术论文、参考文献与可运行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_typeexp, seed2024): rng np.random.default_rng(seed) noise rng.standard_normal((nx, nz)) # 频率坐标单位 rad/m kx np.fft.fftfreq(nx, ddx) * 2.0 * np.pi kz np.fft.fftfreq(nz, ddz) * 2.0 * np.pi KX, KZ np.meshgrid(kx, kz, indexingij) 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相当。这两项检查过了模型才不会在后面的波场递推里给出偏高甚至失真的散射能量。参数经验取值范围如下表。使用场景相关函数类型a (m)ε推荐主频 f0 (Hz)强风化壳指数型1.0 5.00.12 0.2030 50冲洪积砂砾层高斯型0.5 2.00.05 0.1220 40冻土与冰碛堆积指数型2.0 10.00.08 0.1515 30ε 超过 0.2 时要特别检查背景速度 v0 是否足够高避免v0*(1 - ε)接近零甚至为负。硬岩区 v0 取 3000 m/s 时ε0.2 对应最慢 2400 m/s安全余量较大浅表低速层 v0400 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(i1,j) p(i-1,j) p(i,j1) p(i,j-1) ) c2·( p(i2,j) p(i-2,j) p(i,j2) p(i,j-2) ) ] / Δx²差分系数固定如下表。系数数值c0-5/2 -2.5c14/3c2-1/12四阶格式的振幅误差和相位误差都比二阶小一个量级以上。从计算量看要达到同样的波形保真度二阶格式需要把网格间距减半网格点数变成四倍总计算量增加约一个量级所以四阶格式虽然单步运算更复杂整体仍然划算。这个结论在散射模拟里尤其重要因为散射波本身振幅弱数值频散造成的虚假尾波会直接掩盖真实散射信息。3.2 波场递推核心代码与稳定性约束在均匀网格 Δx Δz 条件下二阶时间、四阶空间格式的稳定性条件为Δt ≤ 0.612 · Δx / v_max这个 0.612 来自 λ_max(∇²) 的离散特征值推导实际工程安全值取 0.550.6。时间步长超过这个界限波场会在高频段指数发散表现是总能量随迭代步数单调增长而不是单纯的波形失真。下面是波场递推的内核函数。import numpy as np def fd2d_wavefield_step(p0, p1, v, dt, dx): 对二维声波方程做一步递推 p0: t-dt 时刻波场快照 p1: t 时刻波场快照 v : 速度场二维数组,单位 m/s 返回 p2: tdt 时刻波场快照 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, v01200.0, eps0.15, a3.0, seed7) dt 0.55 * dx / v.max() nt int(0.8 / dt) # 总时长 0.8 s nb 30 # 吸收带宽度,单位网格数 fade np.ones((nx, nz), dtypenp.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), dtypenp.float32) p1 np.zeros((nx, nz), dtypenp.float32) rec np.zeros(nt, dtypenp.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。下面是一组实际算例的参数设定。参数计算方式算例取值模型尺寸覆盖横向炮检距 两侧吸收带400 m × 400 m网格间距 dxv_min / (10·fmax)0.5 m网格点数 nx/nz模型尺寸 / dx801 × 801时间步长 dt0.55·dx / v_max0.1 ms总时长目标最深反射双程时 ×1.20.8 s时间步数 nt总时长 / dt8000算例里 v_min400 m/sf030 Hz 则 fmax75 Hzdx ≤ 0.533 m取 0.5 m 是合理折中。接着 v_max2000 m/s 时 dt0.55×0.5/20000.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, cmapseismic, aspectauto, vmin-5e-4, vmax5e-4) ax.set_title(title) ax.set_xlabel(x (m)) ax.set_ylabel(z (m)) plt.colorbar(im, axax) plt.tight_layout() plt.savefig(snapshot_compare.png, dpi150)三个子图并排观察时主要看散射场子图是否出现与背景介质结构对应的扇形散射带。如果散射场只在震源附近出现一圈对称环而没有延伸到远道说明相关长度 a 偏小或 ε 偏低散射能量没有进入有效传播路径。反之若散射场亮斑铺满整个区域且色标饱和则可能是数值频散在放大网格扰动先加密网格再判断。4.3 检波器布置与记录道里的“假散射”近地表散射模拟很容易把数值噪声当成物理散射。判断标准很简单物理散射能量在主波到达之后才开始显现且随偏移距增大走时差逐渐增大数值干扰从第一时刻起就出现在全频带尤其在高频端与主波同时到达。检波器布置要避开四个角点那里简单吸收带处理不好仍会残留反射道间距取 25 倍网格间距也就是 12.5 m既能覆盖空间波数又不至于数据量过大。单炮记录出来之后先做 1030 Hz 的带通滤波再看散射尾波是否仍明显如果带通后尾波消失说明散射能量主要在算法引入的高频伪影里需要检查吸收带宽度或网格频散。5. 检验散射模拟可信度的三个实用动作5.1 用均匀半空间走时卡时间步把随机速度场换成常数 v1200 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-ε 参数平面以尾波包络误差最小为目标自动选出一组模型参数。我建议把这三项当成每次改模型的固定前置检查走时卡对、网格收敛比值合格、尾波趋势符合尺度预期再谈散射衰减提取和参数反演。这样二维有限差分模拟提供的不是几张漂亮的波场图而是可以放进定量分析流程的可靠合成数据。本文还有配套的精品资源点击获取