简介这是一份基于MATLAB平台实现的初步虚谱法程序面向地球物理、地震勘探及弹性波数值模拟方向的科研人员与初学者用于快速理解并上手伪谱法求解波动方程。伪谱法结合有限差分与谱方法优势在复杂介质中模拟弹性波传播、反射和折射现象较传统差分法更适合高波数问题。压缩包共2个m文件体积约3KB代码结构简洁便于阅读和二次开发。程序主体覆盖计算网格建立、材料参数设置、初始波场与边界条件配置、波动方程求解及结果可视化等关键环节用户可按需调整网格密度、时间步长和物理参数以适配不同模拟场景。目前已有179人学习适合作为伪谱法入门与弹性波模拟实验的基础工具为后续开展地震波传播或地下结构探测等深入研究提供可扩展的起点。1. 伪谱法弹性波模拟程序先搞懂波数域再谈改参数拿到一个名为「初步虚谱法程序.rar」的压缩包解压后往往是几段脚本和一个主循环文件。如果你只是把震源频率改高、网格加密就去跑多半会看到波场像撒了一把盐一样全是噪点然后开始怀疑程序写错了。实际上伪谱法也叫虚谱法的模拟思路很简单空间导数不做有限差分而是用傅里叶变换把场量切到波数域乘一个 i·k 再反变换回来。正因如此它每个波长只需要两三个网格点就能把弹性波的频散压到很小适合做长距离传播、多震相识别和吸收边界条件研究。下面从方程推导、代码实现到参数设定和验证完整走一遍适合刚接触数值模拟的工程师和研究生。2. 伪谱法弹性波模拟的方程骨架从一阶速度应力方程到 FFT 算子2.1 先写一阶速度—应力方程比二阶位移方程更适合伪谱法弹性波最常见的呈现是二阶位移波动方程。如果把位移场的二阶时间导数写出来空间偏导会出现在二级项内部在非均匀介质里还必须对速度和密度的空间梯度做链式展开稍不留神就会丢项。而一阶速度—应力方程把未知量拆成三个应力分量 σxx、σzz、σxz 和两个速度分量 vx、vz时间导数都只出现在左边右侧只含一阶空间偏导格式非常规整。伪谱法的空间求导是可交换的线性算子面对这种一阶方程系统时只需要对每个空间导数做一次 FFT实现与推导一一对应好调试也好加震源。这也是当前主流弹性波伪谱法程序普遍采用一阶形式的原因。对于二维各向同性介质引入密度 ρ 与拉梅常数 λ、μ控制方程可以写成下面这组关系ρ ∂vx/∂t ∂σxx/∂x ∂σxz/∂z fxρ ∂vz/∂t ∂σxz/∂x ∂σzz/∂z fz∂σxx/∂t (λ 2μ) ∂vx/∂x λ ∂vz/∂z Mxx∂σzz/∂t λ ∂vx/∂x (λ 2μ) ∂vz/∂z Mzz∂σxz/∂t μ (∂vx/∂z ∂vz/∂x) Mxz其中 fx、fz 是体力源项Mij 是矩张量源项。伪谱法的任务就是计算 ∂/∂x 和 ∂/∂z。需要特别说明的是σxz 的交叉导数最容易写反符号速度方程里进入 vx 的是 ∂σxz/∂z进入 vz 的是 ∂σxz/∂x而应力方程里 σxz 对时间的导数右侧则是 ∂vx/∂z 与 ∂vz/∂x 之和。程序里一旦把这两类导数对调P 波和 S 波的偏振方向就会出错波前面会明显变形。2.2 空间导数换到波数域伪谱算子的三种实现写法在均匀网格上函数 f(x) 的导数可以用傅里叶变换写成严格形式。离散傅里叶变换把 x 方向映射为整数编号网格对应的角波数 k 取在 -π/dx 到 π/dx 之间于是 ∂f/∂x 的离散表达式为 IFFT(i·k·FFT(f))。这个式子说明伪谱法的空间导数是精确的误差不来自截断形式而来自离散采样是否满足采样定理。实际写代码时角度波数可以用各语言自带的 FFT 频率接口生成。以 numpy 为例有三种写法对应不同的 FFT 布局。第一种是用 np.gradient 做近似这只是有限差分不能算伪谱法。第二种是对 x 方向单独用 rfft 计算先对矩阵每一行做 rfft乘以 i·kx再 irfft 回来这样做一次导数只扫描一个方向速度快但代码稍显啰嗦。第三种是用二维 fft2 同时构造 KX、KZ 两个波数网格让 x 导数和 z 导数共享一次正变换。我一般倾向第三种逻辑最直观代价是多算一次反变换。典型实现如下import numpy as np nx, nz 256, 256 dx dz 10.0 # fftfreq 返回 -0.5~0.5 的归一化频率乘以 2*pi 得到角波数 kx 2.0 * np.pi * np.fft.fftfreq(nx, dx) kz 2.0 * np.pi * np.fft.fftfreq(nz, dz) KX, KZ np.meshgrid(kx, kz) # 注意 KX 的 shape 为 (nz, nx) def ddx(f): return np.fft.ifft2(1j * KX * np.fft.fft2(f)).real def ddz(f): return np.fft.ifft2(1j * KZ * np.fft.fft2(f)).real因为 fft2 默认沿每个轴对全局做周期化处理所以 KX 与 KZ 的行列方向必须与数组 shape 对齐。这里把网格定义为 (nz, nx)第一维对应 z 方向、第二维对应 x 方向那么 KX 的广播要满足每一行都是 kxKZ 的每一列都是 kz。检查时运行一次 ddx(np.ones((nz, nx))) 应得到零矩阵如果给一个正弦波场ddx 与 ddz 的结果之比应接近对应方向的波数。这两个自检常常能快速暴露坐标轴反了的问题。如果波形包含剧烈的空间跳变比如自由表面的应力边界FFT 的反变换会在跳变处出现吉布斯振荡。应对方法是把应力跳变放在两个网格点之间或者用光滑化震源避免直接产生阶跃波形。2.3 FFT 周期延拓带来的三个边界效应折返、吉布斯与混叠伪谱法是对整个计算区域做全局变换所以它隐含假设计算域在空间上周期重复。波场一旦传播到网格边缘并不是消失而是从对侧再次进入计算域。若不处理在较长的模拟时间里就能看到来自四边的虚假波前反向汇聚这是弹性波伪谱法最典型的现象之一。解决思路分为两种把计算域建得足够大让波在结束前到不了边界或者按 Cerjan 方法在边界周围设置衰减带。第二个效应是吉布斯现象。波场中含有接近阶梯状的不连续时有限项傅里叶级数会在断点附近出现过冲幅度约 9% 的振荡。地震波场里常见的不连续是自由表面应力归零这是一种面约束而非体积约束所以业内做法是不在自由表面用伪谱直接强加应力条件而是改用应力镜像法或把自由表面公式化成弹性半空间模型。第三个效应是混叠。离散傅里叶变换可表示的最大空间频率是 Nyquist 波数超过它的分量会被映射到低波数上形成假信号。伪谱法的网格距一旦大于最短波长的一半波场就会同时出现高速假影和颗粒状噪声。判断是否混叠的经验标准是看波前后面是否跟着一串与主频无关的高波数羽毛状波形。这种噪声用时间滤波无法移除只能重新加密网格或降低子波主频。伪谱法虽然理论精度高但采样定理的底线比有限差分更严格。3. 用 Python 把伪谱法弹性波程序跑通最小可运行版本3.1 模型与网格参数先把 vp、vs、rho 和波数矩阵定下来要跑通一个可用的二维弹性波伪谱法程序先把模型参数固定成均匀介质。这样不仅初期好检查后面做解析解验证也比较方便。常见参数如 vp3000 m/s、vs1732 m/s、rho2200 kg/m³其中 vs 取 vp/√3 是泊松比 0.25 的典型比例这样算出的 λ 接近 μ波形中的 P 波与 S 波振幅不会悬殊到看不出来。网格尺寸上网格距取 10 m一个 25 Hz 的 Ricker 子波在 P 波介质中的主波长约 120 m每个波长约 12 个网格点已经足够宽松。这个参数组合在普通笔记本上跑 600 步大约几十秒适合反复改参数观察效果。3.2 核心代码导数算子与应力—速度右端项的写法右端项函数负责输入当前的五个场量输出它们的时间导数。时间推进不管用哪种积分器都要反复调用这个函数所以把它写得精简清晰比什么都重要。第一步先把介质参数和拉梅常数算出来然后定义滤波器导数和右端项函数。下面是完整可运行的核心片段。import numpy as np # 模型参数 nx, nz 256, 256 dx dz 10.0 vp, vs 3000.0, 1732.0 rho 2200.0 mu rho * vs**2 lam rho * vp**2 - 2.0 * mu # 角波数网格 kx 2.0 * np.pi * np.fft.fftfreq(nx, dx) kz 2.0 * np.pi * np.fft.fftfreq(nz, dz) KX, KZ np.meshgrid(kx, kz) def ddx(f): return np.fft.ifft2(1j * KX * np.fft.fft2(f)).real def ddz(f): return np.fft.ifft2(1j * KZ * np.fft.fft2(f)).real def rhs_state(state): vx, vz, sxx, szz, sxz state # 速度方程 dvx_dt (ddx(sxx) ddz(sxz)) / rho dvz_dt (ddx(sxz) ddz(szz)) / rho # 应力方程 dsxx_dt (lam 2.0*mu) * ddx(vx) lam * ddz(vz) dszz_dt lam * ddx(vx) (lam 2.0*mu) * ddz(vz) dsxz_dt mu * (ddz(vx) ddx(vz)) return dvx_dt, dvz_dt, dsxx_dt, dszz_dt, dsxz_dt代码说明rhs_state 的输入输出都是同样结构的元组方便后面 RK4 用 k1k4 做加权求和。注意 dvx 方程中 ddz(sxz) 来自 ∂σxz/∂zdvz 方程中 ddx(sxz) 来自 ∂σxz/∂x不搞混这两个交叉导数的顺序。应力方程中的 dsxz_dt 用到 ddz(vx) ddx(vz)这也和速度方程中的交叉项互为对偶速度方程取 σ 对坐标导数之和应力方程取 v 对坐标导数之和两套关系一定要同步检查。3.3 主循环RK4 时间推进、Ricker 震源与快照保存伪谱法的空间精度已经很高时间离散一般用四阶龙格-库塔保证整体精度避免时间方向成为误差主导。RK4 的四个中间状态全部复用 rhs_state 就行。每次 rhs_state 调用要做 8 次 fft2 和 8 次 ifft2256×256 网格下 600 步约 64 次 FFT 组合在普通 CPU 上依然很快所以不需要一开始就优化。下面是主循环部分。# 时间参数 dt 0.001 # 1 msCFL 3000*0.001/10 0.3 nt 600 print_interval 50 # 初始场 vx np.zeros((nz, nx)); vz np.zeros((nz, nx)) sxx np.zeros((nz, nx)); szz np.zeros((nz, nx)); sxz np.zeros((nz, nx)) # Ricker 子波 fc 25.0 t0 0.08 def wavelet(t): arg np.pi * fc * (t - t0) return (1.0 - 2.0*arg**2) * np.exp(-arg**2) def rk4_step(state, dt): k1 rhs_state(state) s2 tuple(c 0.5*dt*k for c, k in zip(state, k1)) k2 rhs_state(s2) s3 tuple(c 0.5*dt*k for c, k in zip(s2, k2)) k3 rhs_state(s3) s4 tuple(c dt*k for c, k in zip(s3, k3)) k4 rhs_state(s4) return tuple( c dt/6.0 * (k1[i] 2.0*k2[i] 2.0*k3[i] k4[i]) for i, c in enumerate(state) ) # 爆炸源位于网格中心 src_x, src_z nx//2, nz//2 for it in range(nt): t it * dt amp 1.0e5 * wavelet(t) state rk4_step((vx, vz, sxx, szz, sxz), dt) vx, vz, sxx, szz, sxz state # 爆炸源给两个正应力分量同时加压 sxx[src_z, src_x] dt * amp szz[src_z, src_x] dt * amp if it % print_interval 0: print(fstep {it:4d} time {t:.3f}s max|vx|{np.abs(vx).max():.3e}) if it in (200, 400, 600): np.save(fvz_{it:04d}.npy, vz)RK4 的 s2、s3、s4 是对场量的显式 Euler 式推进虽然中间状态本身精度不高但四步加权组合后整体达到四阶精度。震源加载放在 RK4 主推进之后量纲上相当于每步给应力场一个增量关注的是相对波形而非绝对振幅如果关心严格物理量纲应该把源项并入 rhs_state 做统一时间积分。代码中的 amp1.0e5 是让速度场量级落在 1e-2 附近方便显示不代表真实地震震级。3.4 跑起来后先看什么快照图上的 P 波和 S 波响应运行时会每 50 步打印一次 max|vx|它会先增长然后缓慢衰减。模拟结束后用 np.load 读入 vz 快照画 imshow 图。爆炸源只产生 P 波所以 t0.3s 时应看到一个圆形波前向外扩散振幅呈现红蓝相间的双瓣结构。如果改成给 vz 加载垂直力源能看到两个波前面外圈是 P 波稍后到达的大振幅内圈是 S 波速度比约 1.73对应 vp/vs 的取值。快照图上出现方形或沿对角线延伸的条纹时优先检查 KX/KZ 坐标方向是否与数组下标一致而不是去调时间步长。4. 伪谱法程序参数怎么设CFL、网格步长与吸收边界的取舍4.1 时间步长的稳定性条件与 CFL 取值区间伪谱法的稳定性条件由最高波数 π/dx 与介质最大相速度决定。将 CFL 数定义为 C vmax·dt/dx理论上经典 RK4 对纯对流方程有较宽稳定域但弹性波是耦合的双曲系统工程上取 C ∈ [0.2, 0.4] 是安全的。时间步长超过这个范围时波形不会立刻发散而会在模拟后半段出现振幅指数增长或随机点状噪声。确认稳定性最简单的方法是连续记录 max|vx| 并看它是否随时间单调衰减任何在无源时刻的阶跃增长都说明 dt 过大。另一个常见问题是把 dt 取得过小比如为了「保险」取 0.0001结果同样的模拟时间计算量放大十倍。伪谱法的时间步长受波数上限约束并不像有限差分那样受 Courant 数下限约束所以正确的做法是先按 C0.3 反推 dt即 dt ≈ 0.3·dx / vmax然后让它在 0.20.4 之间浮动观察波形变化。若波形对 C 从 0.2 变到 0.4 出现可见差异说明时间离散精度不够此时应该提高 RK 阶数而不是继续缩小 dt。4.2 空间步长与波数上限每最短波长放几个网格点才不产生混叠伪谱法对每波长采样点数的要求很低2 个点是 Nyquist 下限小于 2 就混叠。但工程经验是每最短波长至少 46 个点否则振幅和相位都会失真。最短波速用 S 波速度 vs最短波长 λ_min vs / fmax其中 fmax 近似取 2.53 倍 Ricker 主频。例如 fc25Hz 时 fmax≈80Hzvs1732m/sλ_min≈21.6mdx 应不大于约 3.6m。如果按 P 波速度算网格很可能 S 波已经出现明显频散。这个注意点值得单独强调许多人用 25Hz 子波、dx5m觉得每个 P 波长有 24 个点非常充裕但按 S 波有效最高频算只有约 5 个点已经开始逼近混叠边界。弹性能量集中在 S 波时快照上会出现跟随在波前身后的细碎振荡这不是数值误差而是空间采样不足。稳妥做法是令 dx ≤ vs / (3·fc)在此基础上再看总网格数是否满足计算域尺寸要求。4.3 边界处理Cerjan 衰减带与 PML 在伪谱法中的取舍初步程序第一步可以不加吸收边界只要在波到达边界前终止模拟。要延长模拟时间最简单的边界处理是 Cerjan 衰减带在四周宽度 nb 的条带内每个时间步把场量乘以一个从 1 平滑降到 0 的余弦函数。nb 取 2040 格衰减系数不宜小于 0.1否则边界残留反射明显。代码实现非常直接。nb 30 alpha np.ones(nz) for i in range(nb): alpha[i] np.cos(0.5 * np.pi * i / nb) ** 2 alpha[-(i 1)] alpha[i] A np.multiply.outer(alpha, alpha) # 每个时间步末尾对全部场量统一乘衰减 for field in (vx, vz, sxx, szz, sxz): field * ACerjan 衰减带实现简单但它本质上是衰减器而非吸收边界入射角大的波仍会被反射回来边界宽度不足时会产生二次波前。伪谱法与 PML 的配合比有限差分复杂得多因为空间导数是全局算子直接把坐标复伸长代入 FFT 会破坏 PML 的局部频率偏移特性。常见做法是采用辅助微分方程在 PML 内部对时间导数项额外增加记忆变量把原本的伪谱导数替换成带衰减的差分算子。这部分实现复杂度较高不建议在第一版程序中引入先用 Cerjan 或干脆扩大计算域把精力放在波形验证上更划算。4.4 快速对照表二维弹性波伪谱法的常见参数组合用下面的对照表作为起调参数再根据波形表现逐步收紧参数推荐范围说明CFL vmax·dt/dx0.20.4RK4 常用 0.3低于 2 阶积分器要减半每最短波长点数≥4用 vs 和 2.5fc 估算不用 vpdxvs / (3·fc) 为上限网格距再大风险是 S 波混叠dt0.3·dx/vmax按最高速度取不要按最小速度Cerjan 带宽2040 格余弦平滑宽度过窄反射明显时间积分阶数4 阶 RK2 阶对长时程模拟误差偏大快照保存间隔1050 步每步保存会让磁盘写入拖慢循环表格里最关键的两行是第二行和第四行多数伪谱法程序跑出异常波形都能往上追溯到「用 P 波速度算了波长」或「按 vmin 选了 dt」。这两处修正后C 值略大一点也不会立刻发散的。5. 伪谱法程序验证与排错解析解对比和 3 个高频踩坑点5.1 解析解一致性的验证方法初至时间与波前形状验证程序最简单的办法是把模拟结果与均匀介质中的解析初至时间对比。源在网格中心 (xs, zs)检波点放在 (xr, zr)距离 r 已知那么 P 波初至时间 t_p r / vp。程序里源在 t0 附近有延迟 t0模拟出的波峰到时应当是 t_p t0。从某个检波道的时间序列中挑出 P 波第一个正峰跟理论值比对误差在 12 个时间步内就说明时间推进和速度参数正确。S 波同理再用 vs 重算一遍。如果出现 10ms 以上的系统性偏移多半是 dt 把时间步长理解错了而不是坐标单位错了。更严格的做法是对比空间快照的波前半径。t0.3s 时 P 波波前半径约 0.3×vp900m在图上沿几个方向取最大正振幅点的位置拟合圆半径应落在 ±10m 误差内也就是 ±1 个网格步长。快照验证比单点波形验证更能暴露各向异性问题比如 KX/KZ 轴反向导致的椭圆波前。5.2 高频踩坑点高频噪声、S 波缺失、边界折返的排查顺序高频噪声淹没波场是第一个高频坑。主要原因是 dt 过大或 dx 过大触发不稳定性或混叠改小 dt 并让 dx ≤ vs/(3fc)问题大概率在几十步内消失。需要特别提醒的是不要用低通滤波去掩盖高频噪声因为伪谱法的混叠噪声与真实信号在频率上重叠滤波后的波场同样不可信。第二个坑是只有 P 波没有 S 波。如果用的是爆炸源这是正常现象换垂直力源立刻能看到 S 波如果力源仍然只有 P 波检查 σxz 方程的交叉导数是否写成了 ddz(vx) ddx(vz) 之外的其它组合符号错误会直接抵消偏振分量。第三个坑是波场从对侧边界折返。先估算模拟时长内的最大传播距离 rvmax×nt×dt如果 r 大于计算域边长的一半折返一定会出现。此时在界面上加 3040 格 Cerjan 衰减带并缩短输出时长而不是简单减小振幅或加边界滤波因为滤波只会把波形抹糊。排查顺序建议固定成「时间步长→空间步长→交叉导数→边界」不要上来就怀疑波数矩阵生成有问题。绝大多数伪谱法程序第一次跑飞都出在前两个参数上。先把 dt 与 dx 放小一号跑通波形确认 P/S 波初至正确再逐步放大参数观察波形什么时候开始恶化找到那个临界点对这套伪谱法程序的控制力才算真正建立起来。本文还有配套的精品资源点击获取 伪谱法弹性波模拟：从波数域导数算子到Python实现与参数调优
简介这是一份基于MATLAB平台实现的初步虚谱法程序面向地球物理、地震勘探及弹性波数值模拟方向的科研人员与初学者用于快速理解并上手伪谱法求解波动方程。伪谱法结合有限差分与谱方法优势在复杂介质中模拟弹性波传播、反射和折射现象较传统差分法更适合高波数问题。压缩包共2个m文件体积约3KB代码结构简洁便于阅读和二次开发。程序主体覆盖计算网格建立、材料参数设置、初始波场与边界条件配置、波动方程求解及结果可视化等关键环节用户可按需调整网格密度、时间步长和物理参数以适配不同模拟场景。目前已有179人学习适合作为伪谱法入门与弹性波模拟实验的基础工具为后续开展地震波传播或地下结构探测等深入研究提供可扩展的起点。1. 伪谱法弹性波模拟程序先搞懂波数域再谈改参数拿到一个名为「初步虚谱法程序.rar」的压缩包解压后往往是几段脚本和一个主循环文件。如果你只是把震源频率改高、网格加密就去跑多半会看到波场像撒了一把盐一样全是噪点然后开始怀疑程序写错了。实际上伪谱法也叫虚谱法的模拟思路很简单空间导数不做有限差分而是用傅里叶变换把场量切到波数域乘一个 i·k 再反变换回来。正因如此它每个波长只需要两三个网格点就能把弹性波的频散压到很小适合做长距离传播、多震相识别和吸收边界条件研究。下面从方程推导、代码实现到参数设定和验证完整走一遍适合刚接触数值模拟的工程师和研究生。2. 伪谱法弹性波模拟的方程骨架从一阶速度应力方程到 FFT 算子2.1 先写一阶速度—应力方程比二阶位移方程更适合伪谱法弹性波最常见的呈现是二阶位移波动方程。如果把位移场的二阶时间导数写出来空间偏导会出现在二级项内部在非均匀介质里还必须对速度和密度的空间梯度做链式展开稍不留神就会丢项。而一阶速度—应力方程把未知量拆成三个应力分量 σxx、σzz、σxz 和两个速度分量 vx、vz时间导数都只出现在左边右侧只含一阶空间偏导格式非常规整。伪谱法的空间求导是可交换的线性算子面对这种一阶方程系统时只需要对每个空间导数做一次 FFT实现与推导一一对应好调试也好加震源。这也是当前主流弹性波伪谱法程序普遍采用一阶形式的原因。对于二维各向同性介质引入密度 ρ 与拉梅常数 λ、μ控制方程可以写成下面这组关系ρ ∂vx/∂t ∂σxx/∂x ∂σxz/∂z fxρ ∂vz/∂t ∂σxz/∂x ∂σzz/∂z fz∂σxx/∂t (λ 2μ) ∂vx/∂x λ ∂vz/∂z Mxx∂σzz/∂t λ ∂vx/∂x (λ 2μ) ∂vz/∂z Mzz∂σxz/∂t μ (∂vx/∂z ∂vz/∂x) Mxz其中 fx、fz 是体力源项Mij 是矩张量源项。伪谱法的任务就是计算 ∂/∂x 和 ∂/∂z。需要特别说明的是σxz 的交叉导数最容易写反符号速度方程里进入 vx 的是 ∂σxz/∂z进入 vz 的是 ∂σxz/∂x而应力方程里 σxz 对时间的导数右侧则是 ∂vx/∂z 与 ∂vz/∂x 之和。程序里一旦把这两类导数对调P 波和 S 波的偏振方向就会出错波前面会明显变形。2.2 空间导数换到波数域伪谱算子的三种实现写法在均匀网格上函数 f(x) 的导数可以用傅里叶变换写成严格形式。离散傅里叶变换把 x 方向映射为整数编号网格对应的角波数 k 取在 -π/dx 到 π/dx 之间于是 ∂f/∂x 的离散表达式为 IFFT(i·k·FFT(f))。这个式子说明伪谱法的空间导数是精确的误差不来自截断形式而来自离散采样是否满足采样定理。实际写代码时角度波数可以用各语言自带的 FFT 频率接口生成。以 numpy 为例有三种写法对应不同的 FFT 布局。第一种是用 np.gradient 做近似这只是有限差分不能算伪谱法。第二种是对 x 方向单独用 rfft 计算先对矩阵每一行做 rfft乘以 i·kx再 irfft 回来这样做一次导数只扫描一个方向速度快但代码稍显啰嗦。第三种是用二维 fft2 同时构造 KX、KZ 两个波数网格让 x 导数和 z 导数共享一次正变换。我一般倾向第三种逻辑最直观代价是多算一次反变换。典型实现如下import numpy as np nx, nz 256, 256 dx dz 10.0 # fftfreq 返回 -0.5~0.5 的归一化频率乘以 2*pi 得到角波数 kx 2.0 * np.pi * np.fft.fftfreq(nx, dx) kz 2.0 * np.pi * np.fft.fftfreq(nz, dz) KX, KZ np.meshgrid(kx, kz) # 注意 KX 的 shape 为 (nz, nx) def ddx(f): return np.fft.ifft2(1j * KX * np.fft.fft2(f)).real def ddz(f): return np.fft.ifft2(1j * KZ * np.fft.fft2(f)).real因为 fft2 默认沿每个轴对全局做周期化处理所以 KX 与 KZ 的行列方向必须与数组 shape 对齐。这里把网格定义为 (nz, nx)第一维对应 z 方向、第二维对应 x 方向那么 KX 的广播要满足每一行都是 kxKZ 的每一列都是 kz。检查时运行一次 ddx(np.ones((nz, nx))) 应得到零矩阵如果给一个正弦波场ddx 与 ddz 的结果之比应接近对应方向的波数。这两个自检常常能快速暴露坐标轴反了的问题。如果波形包含剧烈的空间跳变比如自由表面的应力边界FFT 的反变换会在跳变处出现吉布斯振荡。应对方法是把应力跳变放在两个网格点之间或者用光滑化震源避免直接产生阶跃波形。2.3 FFT 周期延拓带来的三个边界效应折返、吉布斯与混叠伪谱法是对整个计算区域做全局变换所以它隐含假设计算域在空间上周期重复。波场一旦传播到网格边缘并不是消失而是从对侧再次进入计算域。若不处理在较长的模拟时间里就能看到来自四边的虚假波前反向汇聚这是弹性波伪谱法最典型的现象之一。解决思路分为两种把计算域建得足够大让波在结束前到不了边界或者按 Cerjan 方法在边界周围设置衰减带。第二个效应是吉布斯现象。波场中含有接近阶梯状的不连续时有限项傅里叶级数会在断点附近出现过冲幅度约 9% 的振荡。地震波场里常见的不连续是自由表面应力归零这是一种面约束而非体积约束所以业内做法是不在自由表面用伪谱直接强加应力条件而是改用应力镜像法或把自由表面公式化成弹性半空间模型。第三个效应是混叠。离散傅里叶变换可表示的最大空间频率是 Nyquist 波数超过它的分量会被映射到低波数上形成假信号。伪谱法的网格距一旦大于最短波长的一半波场就会同时出现高速假影和颗粒状噪声。判断是否混叠的经验标准是看波前后面是否跟着一串与主频无关的高波数羽毛状波形。这种噪声用时间滤波无法移除只能重新加密网格或降低子波主频。伪谱法虽然理论精度高但采样定理的底线比有限差分更严格。3. 用 Python 把伪谱法弹性波程序跑通最小可运行版本3.1 模型与网格参数先把 vp、vs、rho 和波数矩阵定下来要跑通一个可用的二维弹性波伪谱法程序先把模型参数固定成均匀介质。这样不仅初期好检查后面做解析解验证也比较方便。常见参数如 vp3000 m/s、vs1732 m/s、rho2200 kg/m³其中 vs 取 vp/√3 是泊松比 0.25 的典型比例这样算出的 λ 接近 μ波形中的 P 波与 S 波振幅不会悬殊到看不出来。网格尺寸上网格距取 10 m一个 25 Hz 的 Ricker 子波在 P 波介质中的主波长约 120 m每个波长约 12 个网格点已经足够宽松。这个参数组合在普通笔记本上跑 600 步大约几十秒适合反复改参数观察效果。3.2 核心代码导数算子与应力—速度右端项的写法右端项函数负责输入当前的五个场量输出它们的时间导数。时间推进不管用哪种积分器都要反复调用这个函数所以把它写得精简清晰比什么都重要。第一步先把介质参数和拉梅常数算出来然后定义滤波器导数和右端项函数。下面是完整可运行的核心片段。import numpy as np # 模型参数 nx, nz 256, 256 dx dz 10.0 vp, vs 3000.0, 1732.0 rho 2200.0 mu rho * vs**2 lam rho * vp**2 - 2.0 * mu # 角波数网格 kx 2.0 * np.pi * np.fft.fftfreq(nx, dx) kz 2.0 * np.pi * np.fft.fftfreq(nz, dz) KX, KZ np.meshgrid(kx, kz) def ddx(f): return np.fft.ifft2(1j * KX * np.fft.fft2(f)).real def ddz(f): return np.fft.ifft2(1j * KZ * np.fft.fft2(f)).real def rhs_state(state): vx, vz, sxx, szz, sxz state # 速度方程 dvx_dt (ddx(sxx) ddz(sxz)) / rho dvz_dt (ddx(sxz) ddz(szz)) / rho # 应力方程 dsxx_dt (lam 2.0*mu) * ddx(vx) lam * ddz(vz) dszz_dt lam * ddx(vx) (lam 2.0*mu) * ddz(vz) dsxz_dt mu * (ddz(vx) ddx(vz)) return dvx_dt, dvz_dt, dsxx_dt, dszz_dt, dsxz_dt代码说明rhs_state 的输入输出都是同样结构的元组方便后面 RK4 用 k1k4 做加权求和。注意 dvx 方程中 ddz(sxz) 来自 ∂σxz/∂zdvz 方程中 ddx(sxz) 来自 ∂σxz/∂x不搞混这两个交叉导数的顺序。应力方程中的 dsxz_dt 用到 ddz(vx) ddx(vz)这也和速度方程中的交叉项互为对偶速度方程取 σ 对坐标导数之和应力方程取 v 对坐标导数之和两套关系一定要同步检查。3.3 主循环RK4 时间推进、Ricker 震源与快照保存伪谱法的空间精度已经很高时间离散一般用四阶龙格-库塔保证整体精度避免时间方向成为误差主导。RK4 的四个中间状态全部复用 rhs_state 就行。每次 rhs_state 调用要做 8 次 fft2 和 8 次 ifft2256×256 网格下 600 步约 64 次 FFT 组合在普通 CPU 上依然很快所以不需要一开始就优化。下面是主循环部分。# 时间参数 dt 0.001 # 1 msCFL 3000*0.001/10 0.3 nt 600 print_interval 50 # 初始场 vx np.zeros((nz, nx)); vz np.zeros((nz, nx)) sxx np.zeros((nz, nx)); szz np.zeros((nz, nx)); sxz np.zeros((nz, nx)) # Ricker 子波 fc 25.0 t0 0.08 def wavelet(t): arg np.pi * fc * (t - t0) return (1.0 - 2.0*arg**2) * np.exp(-arg**2) def rk4_step(state, dt): k1 rhs_state(state) s2 tuple(c 0.5*dt*k for c, k in zip(state, k1)) k2 rhs_state(s2) s3 tuple(c 0.5*dt*k for c, k in zip(s2, k2)) k3 rhs_state(s3) s4 tuple(c dt*k for c, k in zip(s3, k3)) k4 rhs_state(s4) return tuple( c dt/6.0 * (k1[i] 2.0*k2[i] 2.0*k3[i] k4[i]) for i, c in enumerate(state) ) # 爆炸源位于网格中心 src_x, src_z nx//2, nz//2 for it in range(nt): t it * dt amp 1.0e5 * wavelet(t) state rk4_step((vx, vz, sxx, szz, sxz), dt) vx, vz, sxx, szz, sxz state # 爆炸源给两个正应力分量同时加压 sxx[src_z, src_x] dt * amp szz[src_z, src_x] dt * amp if it % print_interval 0: print(fstep {it:4d} time {t:.3f}s max|vx|{np.abs(vx).max():.3e}) if it in (200, 400, 600): np.save(fvz_{it:04d}.npy, vz)RK4 的 s2、s3、s4 是对场量的显式 Euler 式推进虽然中间状态本身精度不高但四步加权组合后整体达到四阶精度。震源加载放在 RK4 主推进之后量纲上相当于每步给应力场一个增量关注的是相对波形而非绝对振幅如果关心严格物理量纲应该把源项并入 rhs_state 做统一时间积分。代码中的 amp1.0e5 是让速度场量级落在 1e-2 附近方便显示不代表真实地震震级。3.4 跑起来后先看什么快照图上的 P 波和 S 波响应运行时会每 50 步打印一次 max|vx|它会先增长然后缓慢衰减。模拟结束后用 np.load 读入 vz 快照画 imshow 图。爆炸源只产生 P 波所以 t0.3s 时应看到一个圆形波前向外扩散振幅呈现红蓝相间的双瓣结构。如果改成给 vz 加载垂直力源能看到两个波前面外圈是 P 波稍后到达的大振幅内圈是 S 波速度比约 1.73对应 vp/vs 的取值。快照图上出现方形或沿对角线延伸的条纹时优先检查 KX/KZ 坐标方向是否与数组下标一致而不是去调时间步长。4. 伪谱法程序参数怎么设CFL、网格步长与吸收边界的取舍4.1 时间步长的稳定性条件与 CFL 取值区间伪谱法的稳定性条件由最高波数 π/dx 与介质最大相速度决定。将 CFL 数定义为 C vmax·dt/dx理论上经典 RK4 对纯对流方程有较宽稳定域但弹性波是耦合的双曲系统工程上取 C ∈ [0.2, 0.4] 是安全的。时间步长超过这个范围时波形不会立刻发散而会在模拟后半段出现振幅指数增长或随机点状噪声。确认稳定性最简单的方法是连续记录 max|vx| 并看它是否随时间单调衰减任何在无源时刻的阶跃增长都说明 dt 过大。另一个常见问题是把 dt 取得过小比如为了「保险」取 0.0001结果同样的模拟时间计算量放大十倍。伪谱法的时间步长受波数上限约束并不像有限差分那样受 Courant 数下限约束所以正确的做法是先按 C0.3 反推 dt即 dt ≈ 0.3·dx / vmax然后让它在 0.20.4 之间浮动观察波形变化。若波形对 C 从 0.2 变到 0.4 出现可见差异说明时间离散精度不够此时应该提高 RK 阶数而不是继续缩小 dt。4.2 空间步长与波数上限每最短波长放几个网格点才不产生混叠伪谱法对每波长采样点数的要求很低2 个点是 Nyquist 下限小于 2 就混叠。但工程经验是每最短波长至少 46 个点否则振幅和相位都会失真。最短波速用 S 波速度 vs最短波长 λ_min vs / fmax其中 fmax 近似取 2.53 倍 Ricker 主频。例如 fc25Hz 时 fmax≈80Hzvs1732m/sλ_min≈21.6mdx 应不大于约 3.6m。如果按 P 波速度算网格很可能 S 波已经出现明显频散。这个注意点值得单独强调许多人用 25Hz 子波、dx5m觉得每个 P 波长有 24 个点非常充裕但按 S 波有效最高频算只有约 5 个点已经开始逼近混叠边界。弹性能量集中在 S 波时快照上会出现跟随在波前身后的细碎振荡这不是数值误差而是空间采样不足。稳妥做法是令 dx ≤ vs / (3·fc)在此基础上再看总网格数是否满足计算域尺寸要求。4.3 边界处理Cerjan 衰减带与 PML 在伪谱法中的取舍初步程序第一步可以不加吸收边界只要在波到达边界前终止模拟。要延长模拟时间最简单的边界处理是 Cerjan 衰减带在四周宽度 nb 的条带内每个时间步把场量乘以一个从 1 平滑降到 0 的余弦函数。nb 取 2040 格衰减系数不宜小于 0.1否则边界残留反射明显。代码实现非常直接。nb 30 alpha np.ones(nz) for i in range(nb): alpha[i] np.cos(0.5 * np.pi * i / nb) ** 2 alpha[-(i 1)] alpha[i] A np.multiply.outer(alpha, alpha) # 每个时间步末尾对全部场量统一乘衰减 for field in (vx, vz, sxx, szz, sxz): field * ACerjan 衰减带实现简单但它本质上是衰减器而非吸收边界入射角大的波仍会被反射回来边界宽度不足时会产生二次波前。伪谱法与 PML 的配合比有限差分复杂得多因为空间导数是全局算子直接把坐标复伸长代入 FFT 会破坏 PML 的局部频率偏移特性。常见做法是采用辅助微分方程在 PML 内部对时间导数项额外增加记忆变量把原本的伪谱导数替换成带衰减的差分算子。这部分实现复杂度较高不建议在第一版程序中引入先用 Cerjan 或干脆扩大计算域把精力放在波形验证上更划算。4.4 快速对照表二维弹性波伪谱法的常见参数组合用下面的对照表作为起调参数再根据波形表现逐步收紧参数推荐范围说明CFL vmax·dt/dx0.20.4RK4 常用 0.3低于 2 阶积分器要减半每最短波长点数≥4用 vs 和 2.5fc 估算不用 vpdxvs / (3·fc) 为上限网格距再大风险是 S 波混叠dt0.3·dx/vmax按最高速度取不要按最小速度Cerjan 带宽2040 格余弦平滑宽度过窄反射明显时间积分阶数4 阶 RK2 阶对长时程模拟误差偏大快照保存间隔1050 步每步保存会让磁盘写入拖慢循环表格里最关键的两行是第二行和第四行多数伪谱法程序跑出异常波形都能往上追溯到「用 P 波速度算了波长」或「按 vmin 选了 dt」。这两处修正后C 值略大一点也不会立刻发散的。5. 伪谱法程序验证与排错解析解对比和 3 个高频踩坑点5.1 解析解一致性的验证方法初至时间与波前形状验证程序最简单的办法是把模拟结果与均匀介质中的解析初至时间对比。源在网格中心 (xs, zs)检波点放在 (xr, zr)距离 r 已知那么 P 波初至时间 t_p r / vp。程序里源在 t0 附近有延迟 t0模拟出的波峰到时应当是 t_p t0。从某个检波道的时间序列中挑出 P 波第一个正峰跟理论值比对误差在 12 个时间步内就说明时间推进和速度参数正确。S 波同理再用 vs 重算一遍。如果出现 10ms 以上的系统性偏移多半是 dt 把时间步长理解错了而不是坐标单位错了。更严格的做法是对比空间快照的波前半径。t0.3s 时 P 波波前半径约 0.3×vp900m在图上沿几个方向取最大正振幅点的位置拟合圆半径应落在 ±10m 误差内也就是 ±1 个网格步长。快照验证比单点波形验证更能暴露各向异性问题比如 KX/KZ 轴反向导致的椭圆波前。5.2 高频踩坑点高频噪声、S 波缺失、边界折返的排查顺序高频噪声淹没波场是第一个高频坑。主要原因是 dt 过大或 dx 过大触发不稳定性或混叠改小 dt 并让 dx ≤ vs/(3fc)问题大概率在几十步内消失。需要特别提醒的是不要用低通滤波去掩盖高频噪声因为伪谱法的混叠噪声与真实信号在频率上重叠滤波后的波场同样不可信。第二个坑是只有 P 波没有 S 波。如果用的是爆炸源这是正常现象换垂直力源立刻能看到 S 波如果力源仍然只有 P 波检查 σxz 方程的交叉导数是否写成了 ddz(vx) ddx(vz) 之外的其它组合符号错误会直接抵消偏振分量。第三个坑是波场从对侧边界折返。先估算模拟时长内的最大传播距离 rvmax×nt×dt如果 r 大于计算域边长的一半折返一定会出现。此时在界面上加 3040 格 Cerjan 衰减带并缩短输出时长而不是简单减小振幅或加边界滤波因为滤波只会把波形抹糊。排查顺序建议固定成「时间步长→空间步长→交叉导数→边界」不要上来就怀疑波数矩阵生成有问题。绝大多数伪谱法程序第一次跑飞都出在前两个参数上。先把 dt 与 dx 放小一号跑通波形确认 P/S 波初至正确再逐步放大参数观察波形什么时候开始恶化找到那个临界点对这套伪谱法程序的控制力才算真正建立起来。本文还有配套的精品资源点击获取