简介近地表非均质介质中的地震波散射模拟是地震学数值研究的重要方向。资源面向地震学、勘探地球物理方向的科研人员与研究生围绕二维有限差分方法在复杂近地表模型中的波动传播计算可作为理解波动方程时空离散、散射波场分析与非均匀地质模型构建的入门参考。压缩包共7个文件包体仅16KB以markdown说明文档、Python脚本、数值建模代码为主另含JOSS论文引用条目与备份文件等结构清晰便于对照文献复现模拟流程。已有39人浏览学习在小型专题资料中具备一定参考热度。读者可借此获取完整的模型参数设定、程序实现框架及结果解读思路并通过论文与代码的相互对应掌握从理论推导到数值实验落地的关键环节为近地表波场模拟及相关科研工作提供实用支撑。1. 近地表非均质性让合成地震记录脏在哪近地表几十到几百米的风化层、低速透镜体和不规则界面是地震勘探里最难绕开的干扰源。波穿过这些随机速度扰动时会发生地震波散射在炮记录上表现为尾波串、同相轴横向抖动、走时扰动和高频能量提前衰减分层介质模型一概模拟不出来。二维有限差分2D FD是处理这个问题性价比最高的工具时间域直接给出全波场散射、转换波、面波一次算完而且不因速度模型复杂而失效。这里把从方程离散、随机介质建模、边界处理到散射能量定量分析的完整链路讲清楚每个参数都给出可直接落地的参考值。做近地表正演、散射衰减标定或观测系统论证的工程师可以照着这套方案直接搭自己的模拟流程。2. 速度-应力交错网格二维有限差分地震波模拟的离散基础2.1 为什么选速度-应力方程而不是位移二阶方程位移形式的二阶弹性波方程只含一个未知向量场节点少内存省。但它对弹性参数求空间导数在近地表强非均质模型里λ、μ、ρ 在相邻网格上可以是几倍到十几倍的跳变直接离散会产生虚假的界面反射和数值不稳定。速度-应力形式把方程写成一阶偏微分方程组差分只作用于波场分量材料参数以乘积形式参与运算界面处理干净得多。配合交错网格staggered grid空间导数在两个半网格位置取样四阶精度只需一个较短的模板频散特性远好于同位网格。实际代价是每个网格点要维护 5 个二维数组vx、vz、txx、tzz、txz而不是 2 个内存带宽开销翻倍。近地表模型动辄几万乘几万个网格点在 GPU 上通常没有压力CPU 上则要注意数组按列访问的缓存友好性避免按行遍历跨步访问造成 cache miss。常见做法是直接把 z 方向作为数组第二维、把内层循环绑在 z 上让相邻网格在内存里连续。2.2 交错网格二维数组的差分更新代码二维各向同性弹性介质的速度-应力方程组如下交错网格的核心就是把这组方程按半网格位置错开求差import numpy as np def elastic_step(vx, vz, txx, tzz, txz, lam, mu, rho, dt, dx, dz): 2D 速度-应力交错网格单步更新空间二阶、时间二阶。 交错约定数组形状均为 (nx, nz) vx[i,j] 位于 (i0.5, j) vz[i,j] 位于 (i, j0.5) txx[i,j], tzz[i,j] 位于 (i, j) txz[i,j] 位于 (i0.5, j0.5) 边界行 [0,:], [-1,:], [:,0], [:,-1] 不在这里更新 由后续的 PML / 自由表面模块接管。 # 材料参数插值到半网格位置 rho_vx 0.5 * (rho[:-1, :] rho[1:, :]) # ρ 在 (i0.5, j) rho_vz 0.5 * (rho[:, :-1] rho[:, 1:]) # ρ 在 (i, j0.5) mu_txz 0.25 * (mu[:-1, :-1] mu[1:, :-1] mu[:-1, 1:] mu[1:, 1:]) # μ 在 (i0.5, j0.5) # 速度空间导数用于应力更新 dvx_dx (vx[1:, :] - vx[:-1, :]) / dx # 在 (i1, j) dvz_dz (vz[:, 1:] - vz[:, :-1]) / dz # 在 (i, j1) dvx_dz (vx[:, 1:] - vx[:, :-1]) / dz # 在 (i0.5, j0.5) dvz_dx (vz[1:, :] - vz[:-1, :]) / dx # 在 (i0.5, j0.5) # 应力更新 diag lam 2.0 * mu txx[1:-1, 1:-1] dt * (diag[1:-1, 1:-1] * dvx_dx[:-2, 1:-1] lam[1:-1, 1:-1] * dvz_dz[1:-1, :-2]) tzz[1:-1, 1:-1] dt * (lam[1:-1, 1:-1] * dvx_dx[:-2, 1:-1] diag[1:-1, 1:-1] * dvz_dz[1:-1, :-2]) txz[1:-1, 1:-1] dt * mu_txz[1:, 1:] * (dvx_dz[1:-1, 1:] dvz_dx[1:, 1:-1]) # 应力空间导数用于速度更新 dtxx_dx (txx[1:, :] - txx[:-1, :]) / dx # 在 (i0.5, j) dtzz_dz (tzz[:, 1:] - tzz[:, :-1]) / dz # 在 (i, j0.5) dtxz_dx (txz[1:, :] - txz[:-1, :]) / dx # 在 (i0.5, j0.5) dtxz_dz (txz[:, 1:] - txz[:, :-1]) / dz # 在 (i0.5, j0.5) vx[1:-1, 1:-1] dt / rho_vx[1:, 1:-1] * (dtxx_dx[1:, 1:-1] dtxz_dz[1:-1, :-2]) vz[1:-1, 1:-1] dt / rho_vz[1:-1, 1:] * (dtxz_dx[:-2, 1:-1] dtzz_dz[1:-1, 1:]) return vx, vz, txx, tzz, txz这段代码的切片错位很容易写错核对时只看一件事差分结果落在哪个物理位置目标数组那个位置的分量方程需要哪个导数。比如 txx 定义在整数格点 (i, j)它需要 ∂vx/∂x 和 ∂vz/∂z也都必须在 (i, j) 取值vx 的半格位置决定了 ∂vx/∂x 必须取 vx[i,j] 与 vx[i-1,j] 之差体现在代码里就是dvx_dx[:-2, 1:-1]与txx[1:-1, 1:-1]的行错位。参数方面dt是时间步长dx、dz是网格间距二阶精度下差分模板就是相邻两点一阶差商没有更多系数。注意这里的 μ 在 txz 节点取的是四点算术平均。强非均质模型里两侧 μ 差一个量级时更严格的做法是取四点调和平均对应弹性界面上切向应力连续的物理要求ρ 则保持算术平均。这个细节在高扰动幅度ε 超过 15%的随机介质里会直接影响散射能量不是可以忽略的数值洁癖。2.3 网格距、CFL 条件与数值频散的约束表二维交错网格的稳定性条件习惯写成 Courant 数形式$$ C v_{max} \cdot dt \cdot \sqrt{\frac{1}{dx^2} \frac{1}{dz^2}} \le 1 $$二阶空间差分取 C≈0.9 是安全值换成四阶空间算子Levander 模板后上界明显变紧工程上我一般直接压到 0.5 以下多花的计算时间换来的稳定余量值得。频散方面要注意介质最低速度决定最短波长近地表风化层速度可以低到 300–500 m/s控制频散的不是基岩速度。参数符号参考值依据震源主频f050 Hz可控震源/锤击的典型主频有效截止频率fmax2.5 × f0 125 HzRicker 谱衰减到可忽略处最低速度vmin400 m/s未固结风化层最短波长λminvmin / fmax 3.2 m网格距dxλmin / 10 0.32 m二阶精度每波长至少 10 个点时间步长dtdx / (vmax·√2) ≈ 0.13 msCFL取 0.1 ms 留余量记录时长tmax500 ms保证散射尾波完整上表里 vmax 取基岩速度 1800 m/s 算 CFLvmin 只用于算频散这两个角色不要搞混。近地表正演的惯例是先做一次分辨率测试把 dx 减半跑同一测线对比直达波到时差小于半个采样间隔才算网格够密。这一步每个新工区都要做因为风化层速度的横向变化比纵向上更没规律表里的 vmin 只是初猜值。3. von Kármán 随机介质近地表非均质模型的二维生成与参数标定3.1 三种自相关函数与 von Kármán 谱的选择近地表低速带的非均质性尺度跨越几十厘米到几十米露头剖面和井中速度测量都显示其波数谱在多个数量级上呈幂律衰减具有自相似特征。确定性建模在这里行不通只能把介质当作随机介质用自相关函数或功率谱密度描述涨落的统计特征。工程上常用三种自相关函数高斯型谱高频成分急剧衰减生成的介质过于光滑散射波在高频端被低估指数型谱对应 von Kármán 族中赫斯特指数 κ0 的特例粗糙度固定von Kármán 谱则通过 κ 连续调节小尺度粗糙度更贴合实测的幂律谱特征是近地表随机介质建模的默认选择。von Kármán 谱的二维形式通常写成$$ P(k_x, k_z) \frac{C}{\left[1 (k_x a_x)^2 (k_z a_z)^2\right]^{\kappa 1}} $$其中 a_x、a_z 是横向和纵向自相关长度κ 是赫斯特指数。谱在相关长度对应的波数附近转折大于该波数的部分按幂律衰减这决定了散射波的能量分配相关长度越接近波长散射越强远小于波长的涨落只产生瑞利散射能量贡献很小。3.2 谱域滤波与二维傅里叶反变换生成随机扰动场生成随机介质最直接的算法是谱域滤波先产生高斯白噪声乘上 von Kármán 功率谱的平方根再做二维傅里叶反变换回到空间域。核心代码很短但零频和归一化处理不能省def von_karman_medium(nx, nz, dx, dz, ax, az, eps, kappa0.1, seed42): 生成标准差为 eps 的零均值二维随机扰动场 δv/v。 参数: ax, az: x/z 方向自相关长度米近地表取 az ax eps: 扰动幅度相对速度的标准差如 0.10 表示 ±10% kappa: 赫斯特指数0~0.3 范围内取值 rng np.random.default_rng(seed) w rng.standard_normal((nz, nx)) # 1. 白噪声 kx np.fft.fftfreq(nx, ddx) * 2.0 * np.pi # 2. 波数网格rad/m kz np.fft.fftfreq(nz, ddz) * 2.0 * np.pi KX, KZ np.meshgrid(kx, kz) P 1.0 / (1.0 (KX * ax)**2 (KZ * az)**2)**(kappa 1) P[0, 0] 0.0 # 3. 去掉零频均值成分 W np.fft.fft2(w) # 4. 频域滤波 field np.fft.ifft2(W * np.sqrt(P)).real field field / field.std() * eps # 5. 归一化到指定标准差 return field滤波和反变换那两行就是整个方法的全部物理其余步骤都是保证统计量的。第 2 步的波数网格用np.fft.fftfreq直接生成物理单位波数避免索引与波数的换算错误第 3 步把零频置零否则反变换后扰动场带非零均值等效于改变背景速度第 5 步用field.std()重新归一化保证无论 κ、a_x、a_z 取什么值输出场都有确定的标准差 ε这样不同参数组合的散射能量才有可比性。生成扰动场之后速度模型为v(x,z) v0(z) * (1 field)密度扰动可以独立生成也可以按 Gardner 关系取 δρ/ρ ≈ 0.8·δv/v。注意 field 只代表相对扰动v0(z) 里的背景梯度近地表典型的强速度梯度要保留否则模拟的是没有背景折射的均匀介质散射几何完全失真。提示相关长度必须被网格分辨。a_x 或 a_z 小于 2~3 倍网格距时频谱的衰减段超出奈奎斯特波数生成的场在网格尺度上接近白噪声模拟出的散射全是数值频散伪影。3.3 ε、相关长度与赫斯特指数的标定经验随机介质参数不是反演目标时按经验范围给初值就够了要做定量散射衰减标定就得把 (a_x, ε) 当作二维网格搜索的自由度。近地表常见取值范围如下参数参考范围物理含义与说明εδv/v5%–15%风化层可达 20%超过 20% 建议考虑参数平均与更高阶差分a_x10–50 m透镜体、河道砂体的横向尺度a_z2–10 m层状单元厚度a_z / a_x 取 0.2–0.5κ0.05–0.3越小越粗糙0 对应指数型介质δρ/ρ≈ 0.8·δv/vGardner 关系密度扰动不必单独扫描标定流程的常见做法是固定震源和观测系统在 (a_x, ε) 平面上按 3×3 或 5×5 网格生成随机介质各跑 10 次正演取散射尾波包络平均与实测炮记录的尾波能量做最小二乘匹配。这里有个物理上限要记住波数小于震源带宽下限的涨落其散射波不在记录里相关长度远大于 100 m 的慢变化基本表现为走时扰动而非散射不要把它当自由参数扫否则反演不唯一。4. 自由表面与 PML 吸收边界下的近地表地震波散射模拟4.1 自由表面的应力清零与镜像延拓陆上近地表正演不能省略自由表面否则面波和虚反射缺失散射记录的能量构成完全不对。自由表面条件是牵引力为零法向正应力和切向应力在界面上为零。交错网格的实现要点是把自由表面放在 txx/tzz/txz 的一整行节点上这一行恰好落在 z0每步更新后做两件事把该行的 tzz 和 txz 清零再对速度场做镜像延拓补足半网格位置的取值def apply_free_surface(vx, vz, tzz, txz, j_surface0): 自由表面位于 jj_surface 的应力节点行。 tzz[j_surface, :] 0.0 txz[j_surface, :] 0.0 vx[j_surface, :] vx[j_surface 1, :] # 切向速度偶延拓 vz[j_surface, :] -vz[j_surface 1, :] # 法向速度奇延拓镜像延拓的方向由交错约定决定vx 是切向分量越过表面的镜像与内部同号vz 是法向分量镜像反号。对应到代码里是vx[j_surface] vx[j_surface1]还是取负号取决于自由表面落在应力行上还是速度行上先画一张网格位置图再写代码比对着公式蒙要快。做散射研究时还要注意自由表面产生的瑞利波能量非常强尾波分析前通常要做 FK 滤波或时变窗压制否则散射尾波会被面波淹没。4.2 PML 阻尼剖面及参数选择侧面和底边用 PML 吸收散射波比固定边界干净得多。散射波来自四面八方入射角接近掠射时 PML 吸收效果会明显下降所以阻尼剖面的构造比分层介质模拟更讲究。经典 PML 的阻尼系数沿吸收层内随距离做多项式渐变代码构造如下def pml_profile(L, dx, vp_max, R1e-3, m2): 构造 PML 层内各网格点的阻尼系数 d单位 1/m。 参数: L: PML 层厚度网格数 dx: 网格间距 vp_max: 模型最大纵波速度 R: 理论反射系数取 1e-3 ~ 1e-4 m: 渐变多项式阶数2~3 d0 -(m 1) * vp_max * np.log(R) / (2.0 * L * dx) x (np.arange(L) 0.5) * dx d d0 * (x / (L * dx)) ** m return dd0 的公式来自 PML 理论反射系数与阻尼积分的关系式m2 的二次渐变在反射强度和数值稳定性之间比较均衡。参数选择上L 取 10–20 个网格点R 取 1e-3 到 1e-4m 取 2。PML 区域本身也参与差分更新只是每步更新后要对吸收层内的波场乘以衰减算子或对分裂场加衰减项实现时务必把 PML 区内的速度也纳入总模型的最大速度否则 CFL 条件按内部速度算会低估。提示经典 PML 对掠射角入射的散射波吸收不理想工程上更稳的是 CPML卷积 PML额外维护 3 到 6 个记忆变量数组。先用均匀介质做一次吸收测试确保 PML 反射比直达波低 40 dB 以上再进入随机介质模拟。4.3 雷克子波、爆炸源与一个最小正演循环震源用雷克子波主频 f0 由勘探分辨率要求决定实现时注意子波要整体平移以保证因果性def ricker(f0, dt, nt): 雷克子波主频 f0Hz采样间隔 dts长度 nt。 t np.arange(nt) * dt - 1.5 / f0 # 平移 1.5 个主频周期保证起始振幅接近零 return (1.0 - 2.0 * (np.pi * f0 * t)**2) * np.exp(-(np.pi * f0 * t)**2)爆炸源注入两个正应力分量震源在 (isz, isx) 网格点主循环的骨架如下wavelet ricker(f0, dt, nt) rec np.zeros((nt, nr)) # nr 个检波器的垂直分量记录 for it in range(nt): amp wavelet[it] txx[isz, isx] amp * dt # 源项对齐微分方程乘 dt 保持离散一致性 tzz[isz, isx] amp * dt vx, vz, txx, tzz, txz elastic_step( vx, vz, txx, tzz, txz, lam, mu, rho, dt, dx, dz) apply_free_surface(vx, vz, tzz, txz) apply_pml(vx, vz, txx, tzz, txz, pml_param) # 各方向阻尼衰减 rec[it, :] vz[iz_rec, ix_rec] # 陆上检波器记录垂直质点速度源项乘不乘 dt各套代码习惯不一线性问题里绝对幅度本来就不重要但要保证整套流程里量纲一致否则换震源强度后散射能量分析没有可比性。爆炸源主要激发纵波适合看散射尾波要研究面波与散射的耦合就得用作用在自由表面的垂直力源。检波器记录 vz 对应陆上垂直分量检波器记录 (txxtzz)/2 则对应压力分量两种记录的能量衰减规律不同后面做 Q 标定时要分开处理。5. 散射波分离与散射衰减 Q 的定量读取5.1 背景场相减与多次实现的统计分离随机介质模拟最容易犯的错误是把单次随机实现的结果当成真值。单炮记录里既有介质涨落引起的散射也有背景速度结构本身的响应两者混在一起。最干净的分离办法是跑两次正演一次用背景模型 v0(z)一次用 v0(z)·(1field)两次的震源、网格、边界参数完全一致相减得到的就是纯散射波场# u_bg: 背景模型正演记录无随机扰动 # u_pt: 含随机介质的正演记录 u_scat u_pt - u_bg # 纯散射波含散射尾波与直达波的振幅扰动这个差场的能量就是散射造成的附加能量走时扰动导致的波形错位也包含在内。要做统计意义上的衰减分析还需要对多个随机实现做系综平均相干场是 N 次实现的平均散射强度是各次实现相对平均场的方差。N 取 50 次以上统计量才稳定少于 20 次时尾波包络的起伏会淹没你要测的衰减趋势。5.2 谱比法读取散射衰减 Q散射衰减的定量指标是等效 Q测量手段用谱比法取两个偏移距上直达波的同一相位窗比较振幅谱的斜率。二维介质里要先扣掉几何扩散项振幅按 r 的 -1/2 衰减剩余衰减才归因于散射加非弹性def scattering_q(amp1, amp2, dt, r1, r2, t1, t2): 两个偏移距直达波窗的谱比反推 Q。 amp1, amp2: 已截取直达波时窗的波形 r1, r2: 对应偏移距用于扣除 2D 几何扩散 (1/sqrt(r)) t1, t2: 对应到时用于 ln(A2/A1) -pi*f*(t2-t1)/Q n len(amp1) win np.hanning(n) A1 np.fft.rfft((amp1 - amp1.mean()) * win) A2 np.fft.rfft((amp2 - amp2.mean()) * win) f np.fft.rfftfreq(n, ddt) m (f 10.0) (f 100.0) # 只取信噪比高的频带 # 扣几何扩散: 2D 中振幅比例因子 sqrt(r2/r1) ln_ratio np.log(np.abs(A2[m]) / np.abs(A1[m])) 0.5 * np.log(r2 / r1) slope np.polyfit(f[m], ln_ratio, 1)[0] Q -np.pi * (t2 - t1) / slope return Q谱比法的坑主要在时窗窗太长会混入反射波窗太短则频谱分辨率不够一般取 1.5 到 2 倍子波主周期。实测资料里还有一个更实用的路径是测尾波 QQc用单散射模型拟合尾波包络的衰减斜率不需要第二个台站对横向速度变化不敏感适合近地表这种观测系统极不对称的场景。把散射 Q 和随机介质参数 (a_x, ε) 放在同一张图上就能直接读出哪个参数组合把高频能量消耗得最多——这是后期做全波形反演时给介质参数加先验约束的重要依据。本文还有配套的精品资源点击获取 近地表随机介质与地震波散射：二维交错网格正演模拟全解析
简介近地表非均质介质中的地震波散射模拟是地震学数值研究的重要方向。资源面向地震学、勘探地球物理方向的科研人员与研究生围绕二维有限差分方法在复杂近地表模型中的波动传播计算可作为理解波动方程时空离散、散射波场分析与非均匀地质模型构建的入门参考。压缩包共7个文件包体仅16KB以markdown说明文档、Python脚本、数值建模代码为主另含JOSS论文引用条目与备份文件等结构清晰便于对照文献复现模拟流程。已有39人浏览学习在小型专题资料中具备一定参考热度。读者可借此获取完整的模型参数设定、程序实现框架及结果解读思路并通过论文与代码的相互对应掌握从理论推导到数值实验落地的关键环节为近地表波场模拟及相关科研工作提供实用支撑。1. 近地表非均质性让合成地震记录脏在哪近地表几十到几百米的风化层、低速透镜体和不规则界面是地震勘探里最难绕开的干扰源。波穿过这些随机速度扰动时会发生地震波散射在炮记录上表现为尾波串、同相轴横向抖动、走时扰动和高频能量提前衰减分层介质模型一概模拟不出来。二维有限差分2D FD是处理这个问题性价比最高的工具时间域直接给出全波场散射、转换波、面波一次算完而且不因速度模型复杂而失效。这里把从方程离散、随机介质建模、边界处理到散射能量定量分析的完整链路讲清楚每个参数都给出可直接落地的参考值。做近地表正演、散射衰减标定或观测系统论证的工程师可以照着这套方案直接搭自己的模拟流程。2. 速度-应力交错网格二维有限差分地震波模拟的离散基础2.1 为什么选速度-应力方程而不是位移二阶方程位移形式的二阶弹性波方程只含一个未知向量场节点少内存省。但它对弹性参数求空间导数在近地表强非均质模型里λ、μ、ρ 在相邻网格上可以是几倍到十几倍的跳变直接离散会产生虚假的界面反射和数值不稳定。速度-应力形式把方程写成一阶偏微分方程组差分只作用于波场分量材料参数以乘积形式参与运算界面处理干净得多。配合交错网格staggered grid空间导数在两个半网格位置取样四阶精度只需一个较短的模板频散特性远好于同位网格。实际代价是每个网格点要维护 5 个二维数组vx、vz、txx、tzz、txz而不是 2 个内存带宽开销翻倍。近地表模型动辄几万乘几万个网格点在 GPU 上通常没有压力CPU 上则要注意数组按列访问的缓存友好性避免按行遍历跨步访问造成 cache miss。常见做法是直接把 z 方向作为数组第二维、把内层循环绑在 z 上让相邻网格在内存里连续。2.2 交错网格二维数组的差分更新代码二维各向同性弹性介质的速度-应力方程组如下交错网格的核心就是把这组方程按半网格位置错开求差import numpy as np def elastic_step(vx, vz, txx, tzz, txz, lam, mu, rho, dt, dx, dz): 2D 速度-应力交错网格单步更新空间二阶、时间二阶。 交错约定数组形状均为 (nx, nz) vx[i,j] 位于 (i0.5, j) vz[i,j] 位于 (i, j0.5) txx[i,j], tzz[i,j] 位于 (i, j) txz[i,j] 位于 (i0.5, j0.5) 边界行 [0,:], [-1,:], [:,0], [:,-1] 不在这里更新 由后续的 PML / 自由表面模块接管。 # 材料参数插值到半网格位置 rho_vx 0.5 * (rho[:-1, :] rho[1:, :]) # ρ 在 (i0.5, j) rho_vz 0.5 * (rho[:, :-1] rho[:, 1:]) # ρ 在 (i, j0.5) mu_txz 0.25 * (mu[:-1, :-1] mu[1:, :-1] mu[:-1, 1:] mu[1:, 1:]) # μ 在 (i0.5, j0.5) # 速度空间导数用于应力更新 dvx_dx (vx[1:, :] - vx[:-1, :]) / dx # 在 (i1, j) dvz_dz (vz[:, 1:] - vz[:, :-1]) / dz # 在 (i, j1) dvx_dz (vx[:, 1:] - vx[:, :-1]) / dz # 在 (i0.5, j0.5) dvz_dx (vz[1:, :] - vz[:-1, :]) / dx # 在 (i0.5, j0.5) # 应力更新 diag lam 2.0 * mu txx[1:-1, 1:-1] dt * (diag[1:-1, 1:-1] * dvx_dx[:-2, 1:-1] lam[1:-1, 1:-1] * dvz_dz[1:-1, :-2]) tzz[1:-1, 1:-1] dt * (lam[1:-1, 1:-1] * dvx_dx[:-2, 1:-1] diag[1:-1, 1:-1] * dvz_dz[1:-1, :-2]) txz[1:-1, 1:-1] dt * mu_txz[1:, 1:] * (dvx_dz[1:-1, 1:] dvz_dx[1:, 1:-1]) # 应力空间导数用于速度更新 dtxx_dx (txx[1:, :] - txx[:-1, :]) / dx # 在 (i0.5, j) dtzz_dz (tzz[:, 1:] - tzz[:, :-1]) / dz # 在 (i, j0.5) dtxz_dx (txz[1:, :] - txz[:-1, :]) / dx # 在 (i0.5, j0.5) dtxz_dz (txz[:, 1:] - txz[:, :-1]) / dz # 在 (i0.5, j0.5) vx[1:-1, 1:-1] dt / rho_vx[1:, 1:-1] * (dtxx_dx[1:, 1:-1] dtxz_dz[1:-1, :-2]) vz[1:-1, 1:-1] dt / rho_vz[1:-1, 1:] * (dtxz_dx[:-2, 1:-1] dtzz_dz[1:-1, 1:]) return vx, vz, txx, tzz, txz这段代码的切片错位很容易写错核对时只看一件事差分结果落在哪个物理位置目标数组那个位置的分量方程需要哪个导数。比如 txx 定义在整数格点 (i, j)它需要 ∂vx/∂x 和 ∂vz/∂z也都必须在 (i, j) 取值vx 的半格位置决定了 ∂vx/∂x 必须取 vx[i,j] 与 vx[i-1,j] 之差体现在代码里就是dvx_dx[:-2, 1:-1]与txx[1:-1, 1:-1]的行错位。参数方面dt是时间步长dx、dz是网格间距二阶精度下差分模板就是相邻两点一阶差商没有更多系数。注意这里的 μ 在 txz 节点取的是四点算术平均。强非均质模型里两侧 μ 差一个量级时更严格的做法是取四点调和平均对应弹性界面上切向应力连续的物理要求ρ 则保持算术平均。这个细节在高扰动幅度ε 超过 15%的随机介质里会直接影响散射能量不是可以忽略的数值洁癖。2.3 网格距、CFL 条件与数值频散的约束表二维交错网格的稳定性条件习惯写成 Courant 数形式$$ C v_{max} \cdot dt \cdot \sqrt{\frac{1}{dx^2} \frac{1}{dz^2}} \le 1 $$二阶空间差分取 C≈0.9 是安全值换成四阶空间算子Levander 模板后上界明显变紧工程上我一般直接压到 0.5 以下多花的计算时间换来的稳定余量值得。频散方面要注意介质最低速度决定最短波长近地表风化层速度可以低到 300–500 m/s控制频散的不是基岩速度。参数符号参考值依据震源主频f050 Hz可控震源/锤击的典型主频有效截止频率fmax2.5 × f0 125 HzRicker 谱衰减到可忽略处最低速度vmin400 m/s未固结风化层最短波长λminvmin / fmax 3.2 m网格距dxλmin / 10 0.32 m二阶精度每波长至少 10 个点时间步长dtdx / (vmax·√2) ≈ 0.13 msCFL取 0.1 ms 留余量记录时长tmax500 ms保证散射尾波完整上表里 vmax 取基岩速度 1800 m/s 算 CFLvmin 只用于算频散这两个角色不要搞混。近地表正演的惯例是先做一次分辨率测试把 dx 减半跑同一测线对比直达波到时差小于半个采样间隔才算网格够密。这一步每个新工区都要做因为风化层速度的横向变化比纵向上更没规律表里的 vmin 只是初猜值。3. von Kármán 随机介质近地表非均质模型的二维生成与参数标定3.1 三种自相关函数与 von Kármán 谱的选择近地表低速带的非均质性尺度跨越几十厘米到几十米露头剖面和井中速度测量都显示其波数谱在多个数量级上呈幂律衰减具有自相似特征。确定性建模在这里行不通只能把介质当作随机介质用自相关函数或功率谱密度描述涨落的统计特征。工程上常用三种自相关函数高斯型谱高频成分急剧衰减生成的介质过于光滑散射波在高频端被低估指数型谱对应 von Kármán 族中赫斯特指数 κ0 的特例粗糙度固定von Kármán 谱则通过 κ 连续调节小尺度粗糙度更贴合实测的幂律谱特征是近地表随机介质建模的默认选择。von Kármán 谱的二维形式通常写成$$ P(k_x, k_z) \frac{C}{\left[1 (k_x a_x)^2 (k_z a_z)^2\right]^{\kappa 1}} $$其中 a_x、a_z 是横向和纵向自相关长度κ 是赫斯特指数。谱在相关长度对应的波数附近转折大于该波数的部分按幂律衰减这决定了散射波的能量分配相关长度越接近波长散射越强远小于波长的涨落只产生瑞利散射能量贡献很小。3.2 谱域滤波与二维傅里叶反变换生成随机扰动场生成随机介质最直接的算法是谱域滤波先产生高斯白噪声乘上 von Kármán 功率谱的平方根再做二维傅里叶反变换回到空间域。核心代码很短但零频和归一化处理不能省def von_karman_medium(nx, nz, dx, dz, ax, az, eps, kappa0.1, seed42): 生成标准差为 eps 的零均值二维随机扰动场 δv/v。 参数: ax, az: x/z 方向自相关长度米近地表取 az ax eps: 扰动幅度相对速度的标准差如 0.10 表示 ±10% kappa: 赫斯特指数0~0.3 范围内取值 rng np.random.default_rng(seed) w rng.standard_normal((nz, nx)) # 1. 白噪声 kx np.fft.fftfreq(nx, ddx) * 2.0 * np.pi # 2. 波数网格rad/m kz np.fft.fftfreq(nz, ddz) * 2.0 * np.pi KX, KZ np.meshgrid(kx, kz) P 1.0 / (1.0 (KX * ax)**2 (KZ * az)**2)**(kappa 1) P[0, 0] 0.0 # 3. 去掉零频均值成分 W np.fft.fft2(w) # 4. 频域滤波 field np.fft.ifft2(W * np.sqrt(P)).real field field / field.std() * eps # 5. 归一化到指定标准差 return field滤波和反变换那两行就是整个方法的全部物理其余步骤都是保证统计量的。第 2 步的波数网格用np.fft.fftfreq直接生成物理单位波数避免索引与波数的换算错误第 3 步把零频置零否则反变换后扰动场带非零均值等效于改变背景速度第 5 步用field.std()重新归一化保证无论 κ、a_x、a_z 取什么值输出场都有确定的标准差 ε这样不同参数组合的散射能量才有可比性。生成扰动场之后速度模型为v(x,z) v0(z) * (1 field)密度扰动可以独立生成也可以按 Gardner 关系取 δρ/ρ ≈ 0.8·δv/v。注意 field 只代表相对扰动v0(z) 里的背景梯度近地表典型的强速度梯度要保留否则模拟的是没有背景折射的均匀介质散射几何完全失真。提示相关长度必须被网格分辨。a_x 或 a_z 小于 2~3 倍网格距时频谱的衰减段超出奈奎斯特波数生成的场在网格尺度上接近白噪声模拟出的散射全是数值频散伪影。3.3 ε、相关长度与赫斯特指数的标定经验随机介质参数不是反演目标时按经验范围给初值就够了要做定量散射衰减标定就得把 (a_x, ε) 当作二维网格搜索的自由度。近地表常见取值范围如下参数参考范围物理含义与说明εδv/v5%–15%风化层可达 20%超过 20% 建议考虑参数平均与更高阶差分a_x10–50 m透镜体、河道砂体的横向尺度a_z2–10 m层状单元厚度a_z / a_x 取 0.2–0.5κ0.05–0.3越小越粗糙0 对应指数型介质δρ/ρ≈ 0.8·δv/vGardner 关系密度扰动不必单独扫描标定流程的常见做法是固定震源和观测系统在 (a_x, ε) 平面上按 3×3 或 5×5 网格生成随机介质各跑 10 次正演取散射尾波包络平均与实测炮记录的尾波能量做最小二乘匹配。这里有个物理上限要记住波数小于震源带宽下限的涨落其散射波不在记录里相关长度远大于 100 m 的慢变化基本表现为走时扰动而非散射不要把它当自由参数扫否则反演不唯一。4. 自由表面与 PML 吸收边界下的近地表地震波散射模拟4.1 自由表面的应力清零与镜像延拓陆上近地表正演不能省略自由表面否则面波和虚反射缺失散射记录的能量构成完全不对。自由表面条件是牵引力为零法向正应力和切向应力在界面上为零。交错网格的实现要点是把自由表面放在 txx/tzz/txz 的一整行节点上这一行恰好落在 z0每步更新后做两件事把该行的 tzz 和 txz 清零再对速度场做镜像延拓补足半网格位置的取值def apply_free_surface(vx, vz, tzz, txz, j_surface0): 自由表面位于 jj_surface 的应力节点行。 tzz[j_surface, :] 0.0 txz[j_surface, :] 0.0 vx[j_surface, :] vx[j_surface 1, :] # 切向速度偶延拓 vz[j_surface, :] -vz[j_surface 1, :] # 法向速度奇延拓镜像延拓的方向由交错约定决定vx 是切向分量越过表面的镜像与内部同号vz 是法向分量镜像反号。对应到代码里是vx[j_surface] vx[j_surface1]还是取负号取决于自由表面落在应力行上还是速度行上先画一张网格位置图再写代码比对着公式蒙要快。做散射研究时还要注意自由表面产生的瑞利波能量非常强尾波分析前通常要做 FK 滤波或时变窗压制否则散射尾波会被面波淹没。4.2 PML 阻尼剖面及参数选择侧面和底边用 PML 吸收散射波比固定边界干净得多。散射波来自四面八方入射角接近掠射时 PML 吸收效果会明显下降所以阻尼剖面的构造比分层介质模拟更讲究。经典 PML 的阻尼系数沿吸收层内随距离做多项式渐变代码构造如下def pml_profile(L, dx, vp_max, R1e-3, m2): 构造 PML 层内各网格点的阻尼系数 d单位 1/m。 参数: L: PML 层厚度网格数 dx: 网格间距 vp_max: 模型最大纵波速度 R: 理论反射系数取 1e-3 ~ 1e-4 m: 渐变多项式阶数2~3 d0 -(m 1) * vp_max * np.log(R) / (2.0 * L * dx) x (np.arange(L) 0.5) * dx d d0 * (x / (L * dx)) ** m return dd0 的公式来自 PML 理论反射系数与阻尼积分的关系式m2 的二次渐变在反射强度和数值稳定性之间比较均衡。参数选择上L 取 10–20 个网格点R 取 1e-3 到 1e-4m 取 2。PML 区域本身也参与差分更新只是每步更新后要对吸收层内的波场乘以衰减算子或对分裂场加衰减项实现时务必把 PML 区内的速度也纳入总模型的最大速度否则 CFL 条件按内部速度算会低估。提示经典 PML 对掠射角入射的散射波吸收不理想工程上更稳的是 CPML卷积 PML额外维护 3 到 6 个记忆变量数组。先用均匀介质做一次吸收测试确保 PML 反射比直达波低 40 dB 以上再进入随机介质模拟。4.3 雷克子波、爆炸源与一个最小正演循环震源用雷克子波主频 f0 由勘探分辨率要求决定实现时注意子波要整体平移以保证因果性def ricker(f0, dt, nt): 雷克子波主频 f0Hz采样间隔 dts长度 nt。 t np.arange(nt) * dt - 1.5 / f0 # 平移 1.5 个主频周期保证起始振幅接近零 return (1.0 - 2.0 * (np.pi * f0 * t)**2) * np.exp(-(np.pi * f0 * t)**2)爆炸源注入两个正应力分量震源在 (isz, isx) 网格点主循环的骨架如下wavelet ricker(f0, dt, nt) rec np.zeros((nt, nr)) # nr 个检波器的垂直分量记录 for it in range(nt): amp wavelet[it] txx[isz, isx] amp * dt # 源项对齐微分方程乘 dt 保持离散一致性 tzz[isz, isx] amp * dt vx, vz, txx, tzz, txz elastic_step( vx, vz, txx, tzz, txz, lam, mu, rho, dt, dx, dz) apply_free_surface(vx, vz, tzz, txz) apply_pml(vx, vz, txx, tzz, txz, pml_param) # 各方向阻尼衰减 rec[it, :] vz[iz_rec, ix_rec] # 陆上检波器记录垂直质点速度源项乘不乘 dt各套代码习惯不一线性问题里绝对幅度本来就不重要但要保证整套流程里量纲一致否则换震源强度后散射能量分析没有可比性。爆炸源主要激发纵波适合看散射尾波要研究面波与散射的耦合就得用作用在自由表面的垂直力源。检波器记录 vz 对应陆上垂直分量检波器记录 (txxtzz)/2 则对应压力分量两种记录的能量衰减规律不同后面做 Q 标定时要分开处理。5. 散射波分离与散射衰减 Q 的定量读取5.1 背景场相减与多次实现的统计分离随机介质模拟最容易犯的错误是把单次随机实现的结果当成真值。单炮记录里既有介质涨落引起的散射也有背景速度结构本身的响应两者混在一起。最干净的分离办法是跑两次正演一次用背景模型 v0(z)一次用 v0(z)·(1field)两次的震源、网格、边界参数完全一致相减得到的就是纯散射波场# u_bg: 背景模型正演记录无随机扰动 # u_pt: 含随机介质的正演记录 u_scat u_pt - u_bg # 纯散射波含散射尾波与直达波的振幅扰动这个差场的能量就是散射造成的附加能量走时扰动导致的波形错位也包含在内。要做统计意义上的衰减分析还需要对多个随机实现做系综平均相干场是 N 次实现的平均散射强度是各次实现相对平均场的方差。N 取 50 次以上统计量才稳定少于 20 次时尾波包络的起伏会淹没你要测的衰减趋势。5.2 谱比法读取散射衰减 Q散射衰减的定量指标是等效 Q测量手段用谱比法取两个偏移距上直达波的同一相位窗比较振幅谱的斜率。二维介质里要先扣掉几何扩散项振幅按 r 的 -1/2 衰减剩余衰减才归因于散射加非弹性def scattering_q(amp1, amp2, dt, r1, r2, t1, t2): 两个偏移距直达波窗的谱比反推 Q。 amp1, amp2: 已截取直达波时窗的波形 r1, r2: 对应偏移距用于扣除 2D 几何扩散 (1/sqrt(r)) t1, t2: 对应到时用于 ln(A2/A1) -pi*f*(t2-t1)/Q n len(amp1) win np.hanning(n) A1 np.fft.rfft((amp1 - amp1.mean()) * win) A2 np.fft.rfft((amp2 - amp2.mean()) * win) f np.fft.rfftfreq(n, ddt) m (f 10.0) (f 100.0) # 只取信噪比高的频带 # 扣几何扩散: 2D 中振幅比例因子 sqrt(r2/r1) ln_ratio np.log(np.abs(A2[m]) / np.abs(A1[m])) 0.5 * np.log(r2 / r1) slope np.polyfit(f[m], ln_ratio, 1)[0] Q -np.pi * (t2 - t1) / slope return Q谱比法的坑主要在时窗窗太长会混入反射波窗太短则频谱分辨率不够一般取 1.5 到 2 倍子波主周期。实测资料里还有一个更实用的路径是测尾波 QQc用单散射模型拟合尾波包络的衰减斜率不需要第二个台站对横向速度变化不敏感适合近地表这种观测系统极不对称的场景。把散射 Q 和随机介质参数 (a_x, ε) 放在同一张图上就能直接读出哪个参数组合把高频能量消耗得最多——这是后期做全波形反演时给介质参数加先验约束的重要依据。本文还有配套的精品资源点击获取