简介一套面向弹性波数值模拟的伪谱法虚谱法入门程序基于MATLAB平台实现适合地球物理、地震学、地质勘探及波动问题研究的初学者快速上手。伪谱法通过快速傅里叶变换FFT求解波动方程兼顾有限差分法的直接性与谱方法的高精度尤其擅长处理高波数、复杂介质中的波传播、反射与折射模拟。程序压缩包仅3KB包含2个M脚本文件分别承担模型参数设置和主计算流程结构精简便于逐行阅读、修改与二次开发目前已有179人学习下载。代码覆盖计算网格建立、材料参数设定、初始波场初始化、边界条件处理、波动方程求解及结果可视化等核心流程可演示弹性波在非均匀介质中的传播特征。使用者可通过调整网格密度、时间步长、物性参数和边界类型快速适配不同模拟场景是掌握伪谱法实现细节并开展拓展研究的基础参考。1. 伪谱法模拟弹性波为什么这个程序只用几条 FFT刚拿到一版初步虚谱法程序的人最容易犯的误判是以为它和有限差分一样要面对一个巨大的稀疏矩阵。实际上伪谱法把空间导数搬进波数域对整个波场做 FFT乘上 ik再做 IFFT一次偏导数就完成。弹性波的正演程序在一个时间循环里只会反复出现这几条 FFT代码量比有限差分法少一大截换来的是网格可以更粗、数值频散更小代价是边界处理和稳定性条件比想象中挑剔。下面我按自己落地时的顺序把这类伪谱法弹性波程序从导数算子到参数设计完整拆开适合刚接触伪谱法模拟的读者也适合正在调试旧代码、想让正演结果更可靠的人。我的原则是先用二维均匀模型跑通再考虑分层与三维。2. 从弹性波动方程到谱域求导导数算子与时间递推怎么搭伪谱法的原理不复杂关键是把导数算子和时间递推的骨架搭对。这一章把数学上要用的东西压到最少直接对应后面程序里的函数先把“为什么用谱导数”和“代码怎么写”绑在一起。2.1 为什么空间导数要走波数域把偏微分方程变成乘除法经典有限差分用少数邻近点拟合局部斜率伪谱法用的是全网格信息。对任意离散场 u傅里叶变换把空间域转到波数域空间导数变成复数乘法∂u/∂x F⁻¹[ i·kx · F[u] ]。这里 kx 是沿 x 方向的波数负频率部分按 FFT 的固有排列处理。只要场是带限的、网格足够细这个导数算子的精度几乎到机器精度不会像差分模板那样自带数值频散。直白地说伪谱法在同样网格密度下波前保真度比二阶差分高得多把网格放宽到每波长 4 个点仍能维持可用精度这是很多人把伪谱法拿来做大网格快速试探的原因。弹性波和声波的最大区别在于它有横波和纵波两种传播速度。均匀各向同性介质里速度由拉梅参数 λ、μ 和密度 ρ 决定vp sqrt((λ2μ)/ρ)vs sqrt(μ/ρ)。我们解的不是一个标量压力场而是两个速度分量 vx、vz 和三个应力分量 σxx、σzz、σxz。一阶速度-应力方程组可以写成ρ·∂vx/∂t ∂σxx/∂x ∂σxz/∂z fxρ·∂vz/∂t ∂σxz/∂x ∂σzz/∂z fz以及应力的三个方程 ∂σxx/∂t (λ2μ)·∂vx/∂x λ·∂vz/∂z、∂σzz/∂t λ·∂vx/∂x (λ2μ)·∂vz/∂z、∂σxz/∂t μ·(∂vx/∂z ∂vz/∂x)。这组方程里的空间导数全部可以交给同一个谱求导函数处理。所以教学版伪谱法程序的核心不是解大型方程组而是一遍遍重复“FFT → 乘 ik → IFFT”这个动作。看懂这个算子整个主循环就清楚了。2.2 波数向量怎么排决定了导数算子对不对FFT 的波数排列有固定约定这也是最先出错的地方。MATLAB 的 fft 输出索引 0 对应零频正频率在 1 到 N/2负频率在 N/21 到 N-1。二维网格的 x 方向波数向量要写成function kvec wavenumber(n, dx) % n: 该方向网格数, dx: 网格间距(m) % 返回与 fft 输出顺序一致的波数向量, 单位 rad/m if mod(n, 2) 0 kvec (2 * pi / (n * dx)) * [0:n/2-1, -n/2:-1]; else kvec (2 * pi / (n * dx)) * [0:(n-1)/2, -(n-1)/2:-1]; end end偶数网格是最常见的情况这里把正半波数和负半波数直接拼起来对应 fft 的输出顺序。注意乘的是 2π 而不是 1因为 FFT 对离散序列的定义含 2π 因子空间导数对应的连续波数是 2π/(N·dx) 的整数倍单位要写成 rad/m后面乘 i·k 才量纲正确。谱求导函数本身只有几行。x 方向和 z 方向我习惯分开写避免二维广播维度出错function du ps_deriv_x(u, kx) % 沿 x 方向(第2维)做傅里叶导数 U fft(u, [], 2); du real(ifft(bsxfun(times, U, reshape(kx, 1, [])), [], 2)); end function du ps_deriv_z(u, kz) % 沿 z 方向(第1维)做傅里叶导数 U fft(u, [], 1); du real(ifft(bsxfun(times, U, reshape(kz, [], 1)), [], 1)); end这段代码里fft 沿指定维度做变换然后乘上对应波数向量ifft 之后取实部。理论上 i·k 乘完应该得到纯虚部与实部相加的结果由于数值舍入会带极小虚部直接取实部是常规操作。dx、dz 哪怕相等也不要混用同一个波数向量否则斜向传播的波会出现方向性失真。2.3 一阶弹性波方程的最小时间递推骨架空间导数有谱算子后就剩时间递推。伪谱法程序里最常见的时间方案是蛙跳式二阶中心差分即速度场和应力场交替更新。为了后面方便套阻尼带通常写成一个函数在一个时间步内先更新应力再更新速度function [vx, vz, sxx, szz, sxz] pseudo_elastic_step(... vx, vz, sxx, szz, sxz, lam, mu, rho, kx, kz, dt) % 1) 由速度场求应变率更新应力 dvxdx ps_deriv_x(vx, kx); dvxdz ps_deriv_z(vx, kz); dvzdx ps_deriv_x(vz, kx); dvzdz ps_deriv_z(vz, kz); sxx sxx dt * ((lam 2 * mu) .* dvxdx lam .* dvzdz); szz szz dt * (lam .* dvxdx (lam 2 * mu) .* dvzdz); sxz sxz dt * mu .* (dvxdz dvzdx); % 2) 由应力场求加速度更新速度 dsxxdx ps_deriv_x(sxx, kx); dsxzdz ps_deriv_z(sxz, kz); dsxzdx ps_deriv_x(sxz, kx); dszzdz ps_deriv_z(szz, kz); vx vx dt ./ rho .* (dsxxdx dsxzdz); vz vz dt ./ rho .* (dsxzdx dszzdz); % 阻尼带与震源由调用方处理保持函数只做物理更新 end这个骨架里所有空间导数调用同一套谱算子应力与速度的耦合方式来自前面的一阶方程组。注意更新应力时用的是速度的空间导数更新速度时用的是应力的空间导数方向不能颠倒。如果写反程序表现为波场能量缓慢发散而不是立刻崩掉很容易被误判成格式不稳定。另外 dt 的单位是秒速度单位 m/s应力单位 Parho 单位 kg/m³代入数值前先统一单位。伪谱法的主要误差来源在时间离散空间导数本身几乎不引入数值频散。因此 dt 的选择比有限差分更敏感第 4 章会给经验范围。这里先记住一个判断方法骨架跑通后把 dt 减半看结果是否变化如果波形随 dt 明显变化说明时间步长还没进入收敛区。这一点和有限差分的习惯差异很大是伪谱法程序常被忽略的特征。3. 把程序跑起来目录结构、震源加载和最小可运行脚本这一章直接给一个能跑的二维均匀模型。先别急着加分层、加地形、加吸收边界把最小闭环搭出来确认递推和震源没问题再往上面堆功能。3.1 按五个文件拆开不要写成一个巨型脚本我经手过的教学版伪谱法程序十有八九是单个脚本从头写到尾参数、模型、递推、画图全部塞在一起。这样跑通可以但排错和改参数很痛苦。我一般会拆成五个文件职责清晰也方便逐步验证。最小目录结构是文件职责输出main_pseudo_elastic.m读参数、建模、调用主循环、保存结果波场快照build_model.m生成 vp、vs、rho、lambda、mu 和阻尼带模型数组ricker_source.m生成震源时间函数离散时间序列pseudo_elastic_step.m执行一次谱导数加应力/速度更新更新后的波场ps_deriv_x.m / ps_deriv_z.m沿 x / z 的傅里叶导数导数场这个拆分的好处是模型不对只看 build_model导数不对只测 ps_deriv递推发散只查 pseudo_elastic_step。如果你已经习惯用 Fortran 或 Python 组织代码职责边界也可以照搬不是非得用 MATLAB。3.2 最小主程序从读速度模型到逐帧输出波场快照先写一个均匀全空间模型。vp、vs、rho 全数组赋值一次lambda、mu 也直接按公式生成避免一开始就被分层界面反射干扰。% main_pseudo_elastic.m clear; clc; close all; % ---------- 模型参数 ---------- nx 512; nz 512; % 网格数 dx 5.0; dz 5.0; % 网格间距单位 m vp_true 2000; vs_true 1150; rho0 2200; vp vp_true * ones(nz, nx); vs vs_true * ones(nz, nx); rho rho0 * ones(nz, nx); lam rho .* (vp.^2 - 2 * vs.^2); mu rho .* vs.^2; % ---------- 时间参数 ---------- nt 800; dt 0.0005; % 总步数和步长步长 0.5 ms % ---------- 震源 ---------- f0 15; t0 1.2 / f0; % 主频 15 Hz延迟约 80 ms src_t ricker_source(f0, t0, dt, nt); sx nx / 2; sz nz / 2; % 震源在模型中心 sigma 2 * dx; % 空间平滑半径 [zz, xx] ndgrid(1:nz, 1:nx); spatial exp(-((xx - sx).^2 (zz - sz).^2) / (2 * sigma^2)); % ---------- 主循环 ---------- vx zeros(nz, nx); vz zeros(nz, nx); sxx zeros(nz, nx); szz zeros(nz, nx); sxz zeros(nz, nx); kx wavenumber(nx, dx); kz wavenumber(nz, dz); snap vx; % 先存最后一帧初测用 for it 1:nt [vx, vz, sxx, szz, sxz] ... pseudo_elastic_step(vx, vz, sxx, szz, sxz, lam, mu, rho, kx, kz, dt); % 暂只加水平分量便于观察P/S波前 vx vx dt * 1e12 * src_t(it) .* spatial; snap vx; end % ---------- 保存 ---------- save(snap_vx_final.mat, snap, dx, nz, nx, dt);这段程序里几个关键点值得解释。震源幅度 1e12 是经验值目的是让波场振幅落在 1e-6 到 1e-4 的可观察量级实际建模时应该按真实物理量标定但做正演合理性检查时不需要追求绝对振幅。空间平滑半径 sigma 取两个网格间距比单点源温和很多可以明显抑制伪谱法里点源造成的吉布斯振荡。lam、mu、rho 用数组而不是标量后面换成分层模型时主循环一行都不用改。ricker_source 的生成方式也是初学容易写错的地方。常见写法是function w ricker_source(f0, t0, dt, nt) % f0: 主频(Hz), t0: 延迟(s), dt: 时间步长(s), nt: 总步数 t (0:nt-1) * dt - t0; w (1 - 2 * pi^2 * f0^2 .* t.^2) ... .* exp(-pi^2 * f0^2 .* t.^2); end延迟 t0 必须大于零并覆盖子波主瓣否则震源从第一个时间步就注入突变脉冲相当于激励一个极宽频带高频成分立刻超过网格分辨率波场很快颗粒化。t0 取值至少 1/f0我习惯用 1.2/f0确保子波起始段接近零。3.3 运行后的第一张图怎么判断程序正常跑完 nt 步后画出 vx 快照。均匀全空间模型里的预期是一个近似圆形的波前从中心向外扩散外圈是 P 波内圈或紧随其后的是 S 波。伪谱法的波前边缘应该连续光滑不应该有可见的“锯齿”。同时观察最大振幅如果 max(abs(vx(:))) 在 1e-8 到 1e-4 之间说明震源量级基本合理如果出现 NaN 或振幅超过 1e2多半是 dt 太大或震源被重复叠加。如果波前形状明显不是圆形优先检查 dx 与 dz 是否同时进入 wavenumber。常见翻车是把 x 方向的波数向量误用于 z 方向导数导致垂直导数是错的波前面变成椭圆甚至沿对角线发散。初测阶段只保存最后一帧 vx 就够了不要把所有时刻都存成三维数组二维模型全时刻保存也能轻松到几个 GB三维模型直接内存爆炸。4. 网格、时间步长与吸收边界四个必调参数和取值依据伪谱法程序的参数比有限差分少但每个参数都和稳定性、分辨率直接相关。这一章给出四个必调参数的估算方法和初始值按这个顺序调能少走很多弯路。4.1 网格间距 dx先算最短波长再乘安全系数决定 dx 的是模型里最慢速度和期待的最高有效频率。伪谱法的理论极限是每个波长两个网格点也就是 Nyquist 极限但实际正演中两个点会造成波形走样和方向性误差。弹性波模拟我按最短波长至少 4 个点来取稳妥时取 5 个点。最短波长 λmin vs_min / fmax其中 vs_min 是模型里的最小横波速度fmax 不是震源主频而是震源频谱里还有有效能量的最高频率约为主频的 2 到 3 倍。于是dx ≤ vs_min / (4 * fmax)例如 vs_min 1150 m/sf0 15 Hz取 fmax 40 Hzλmin 28.75 mdx 应不超过 7.2 m。为了留余量用 dx 5 m 合适。真正限制 dx 的是最低速度不是 P 波速度这个顺序不要搞反。4.2 时间步长 dtCourant 数按 0.2 起步不是 0.8伪谱法的稳定性是新手最容易踩的坑。有限差分程序里许多人习惯把 Courant 数放到 0.5 甚至 0.8伪谱法不行。原因是谱导数算子看到的最高波数接近 Nyquist这些高频成分在显式时间递推里的稳定域很窄。我按 vmax·dt/dx 来估计vmax 取模型里的最大 P 波速度初始取 Courant 数 C ≤ 0.2再减半验证一次。以 vp 2000 m/s、dx 5 m、C 0.2 为例dt C·dx/vmax 0.0005 s正好是 0.5 ms。如果 dt 太大现象不会是立刻发散而是先出现高频噪声波前边缘开始毛糙随后能量骤增变成 NaN。这种“先携带噪声、再崩溃”的过程很容易被误判成模型缺陷实际只是稳定条件没满足。也可以用四阶 Runge-Kutta 替换时间递推稳定域更大但每个时间步要算四次谱导数成本翻倍。初学阶段建议先用二阶跑通后再考虑换 RK4。4.3 吸收边界在模型四周压一圈阻尼带伪谱法自带的周期性边界意味着波从右边出去就会从左边绕回来。处理它的常见手法是在模型四周压一圈阻尼带每个时间步末乘一个衰减系数让波在到达边界前按指数衰减。我常用的参数是阻尼带宽 n_pml 30 个网格点衰减系数从内边界 0 渐变到外边界 α_maxα_max 的经验公式取 2.5·vmax/(n_pml·dx)单位 1/s渐变函数用二次型。代码片段% 在 build_model.m 中生成阻尼系数数组 vmax max(vp(:)); n_pml 30; % 阻尼带宽单位网格数 alpha_max 2.5 * vmax / (n_pml * dx); alpha zeros(nz, nx); for i 1:n_pml % 上、下边界 ratio ((i - 0.5) / n_pml)^2; alpha(i, :) alpha_max * ratio; alpha(end-i1, :) alpha_max * ratio; end for j 1:n_pml % 左、右边界 ratio ((j - 0.5) / n_pml)^2; alpha(:, j) max(alpha(:, j), alpha_max * ratio); alpha(:, end-j1) max(alpha(:, end-j1), alpha_max * ratio); end damp exp(-alpha * dt); % 每步的衰减系数数组主循环里每个波场更新前都要乘 dampvx vx .* damp; vz vz .* damp; sxx sxx .* damp; szz szz .* damp; sxz sxz .* damp;乘阻尼带的位置很关键。如果每步乘一次等效于一个平滑吸收不会引入额外频散如果只在少数几步乘波会在带内产生内部反射。出现边缘反射时先加宽带宽再微调 α_max不要只动一个。α_max 太大时阻尼带本身会像一道硬界面产生新反射这处参数调和起来像玄学但按顺序来就很快。4.4 Ricker 主频 f0 和空间平滑半径f0 决定震源频带也直接参与最短波长计算。f0 提高 2 倍最短波长减半dx 要跟着减半计算量按二维是 4 倍、按三维是 8 倍。因此在试探性模拟里f0 习惯取偏低数值比如 10-20 Hz先把波前面形态、P/S 分离和层位响应看明白再提高 f0 做精细计算。空间平滑半径 sigma 取 1 到 2 个 dx 即可太大会把震源等效成大面积气枪信号改变波场的高频特征。另一点常被忽略的是提高 f0 时一定要同步检查 dx。很多人只改了震源频率觉得波形会变细致结果却看到网格噪声原因就是最短波长已经掉到 4 个网格点以下。四个参数汇总成初始建议表参数默认值取值范围建议影响后果dx, dz5 m按 vmin/(4*fmax) 估算网格太粗出现伪频散dt0.5 msvmax*dt/dx 0.2dt 太大先噪声后发散n_pml30 格20~50 格太窄边界反射α_max2.5vmax/(n_pmldx)按带宽联动微调太大带内反射f015 Hz模型最小尺度/20~40太高计算量骤增提示这些默认值来自二维均匀模型。换到分层模型或三维前先按本节公式重算 dx 和 dt不要沿用上一个模型的数值。5. 伪谱法程序常见问题排查五条踩坑记录这五条是我在调试伪谱法弹性波程序时遇到的最典型的坑按出现频率排序。无论哪一条先把模型退化到均匀介质、震源放中心、只输出一个波场分量再排查效率最高。5.1 边缘出现一圈强反射甚至波形绕回另一侧现象快照放到中后期模型四角出现明亮弧线随后在另一侧出现对称的“新波前”。原因伪谱法 FFT 自带周期性边界程序如果没有阻尼带或者阻尼带只覆盖了上边界波从左边界出去就会从右边界绕进来。解决检查四个边界阻尼带必须覆盖四边。调参时先固定 n_pml 50α_max 按经验公式确认反射完全消失后再逐步减小带宽到 30。常见错误是阻尼带只作用于上边界这种程序在下方震源时看似没事一旦震源靠近上边界就会暴露。5.2 波场出现密集颗粒状噪声而不是光滑波前现象波前周围出现大量短波长抖动或整个快照看起来像“砂纸”。原因空间分辨率不足或震源频谱太宽高频成分超过网格能表示的波数。还有一种隐蔽情况是震源没做空间平滑单点源在波数域是常数谱相当于所有波数等幅注入。解决先按 4.1 估算 fmax 和 dx确认震源主频不是过高再把震源换成高斯平滑空间分布sigma 至少 1 到 2 个 dx。如果噪声只在局部区域出现检查是不是速度界面太硬波数域对阶梯模型的吉布斯效应需要模型平滑过渡来抑制。5.3 命令行报错“不是内部或外部命令”或“无法定位程序输入点”现象换到新电脑或双击脚本运行时报某某命令不认识、动态库找不到入口点和算法代码本身无关。原因伪谱法程序很少是单一可执行文件一般依赖 MATLAB、Python 或 Fortran 编译环境。这类报错绝大多数是系统 PATH 没配置好或运行库版本混用而不是代码逻辑错误。解决先确认执行环境。MATLAB 脚本就在 MATLAB 编辑器里运行不要在 cmd 里直接敲文件名Python 端的 numpy/scipy 建议用 conda 统一管理装完重启终端如果仍然找不到 conda 命令先看系统 PATH 是否把 conda 目录加进去。遇到“无法定位程序输入点”这类报错多半是动态链接库版本混用重装对应运行时或把依赖库一并拷贝到程序目录通常能解决。这类环境问题占新手问询的比例相当大浪费的时间很可惜。5.4 内存不足问题出在保存全部时刻的波场现象nt 步后内存不够程序被杀或 MATLAB 无限转圈开始用硬盘交换页面。原因二维模型 512×512 网格、3000 步只保存 vx就是 512×512×3000×8 字节约 6.3 GB。如果五个波场全保存30 GB 起。三维模型更是灾难。解决两步走。先改单精度 float32波场内存直接减半再只保留关键时层每 10 步采样一次或把快照落到磁盘而不是保留在内存。要做合成记录时每个接收点数据单独按 1×nt 抽取不要全波场都留在内存。初测阶段只保存最后一帧足够判断对错。5.5 运行很久后能量不降反升最后变成 NaN现象前几百步波形正常某个时刻开始局部出现极大值随后整个数组变成 NaN 或 Inf。原因时间步长已经进入不稳定区或时间递推顺序写错。先更新应力再更新速度的同步格式误差会随步长累积如果 dt 一开始就偏大崩溃点通常出现在波形传过最快速度区域之后而不是第一步。解决把 dt 减半再跑一次。如果减半后稳定说明是稳定条件问题按 Courant 0.2 重新取 dt如果仍然发散检查阻尼带是否被应用了两次或符号写反尤其注意 damp 数组是否出现负值。另一个诊断技巧是单独把震源幅度提高一个量级线性方程在幅度增大时本应等比例放大如果小幅度稳定、大幅度崩溃说明时间格式对高频成分没有稳定抑制问题仍然在 dt 或阻尼带。注意上面五条是按出现频率排的。处理任何一条之前先用均匀介质、中心震源、单分量输出的配置把问题从模型复杂性里剥离出来。6. 进阶验证把弹性波场拆成 P 波和 S 波检查程序到底对不对6.1 用散度和旋度分离两类波均匀介质里弹性波场是两类波的叠加P 波是无旋场对应速度场的散度S 波是无散场对应速度场的旋度。用谱导数可以直接分离div_p ps_deriv_x(vx, kx) ps_deriv_z(vz, kz); % P 波 s_curl ps_deriv_x(vz, kx) - ps_deriv_z(vx, kz); % S 波如果程序实现正确div_p 快照里只能看到外圈 P 波波前s_curl 快照里只能看到 S 波波前互不串扰。这一步比直接看 vx 更严格因为 vx 图上两种波混在一起只能凭速度差猜。div_p 图里如果出现明显的 S 波波前说明速度-应力方程的耦合项写错或 lambda、mu 的分配有误。6.2 用理论走时做定量检查均匀模型里点源走时可以直接算接收点距震源 rP 波理论到时是 r/vpS 波是 r/vs。从某一帧快照数出两个波前半径再对照 dt 累加后的传播时间误差在两三个网格以内说明从导数算子、时间递推、震源加载到阻尼带是完整的链路。超过这个范围先查单位再查速度模型里 vp、vs 与 lam、mu 的换算是否一致。这个错误在代码里很难一眼发现却会让所有波速整体变慢。我现在的习惯是每次改模型参数后先跑一遍均匀模型的走时对照再做实际的正演或反演计算。这一步十几分钟成本能把大量“看起来能跑但数值不对”的问题拦在正式计算之前。修伪谱法程序这些年我最大的体会是别过度相信“图好看”波前形状对、振幅量级对、走时对三个条件同时满足这个程序才值得继续投入。希望这个验证流程对你也有用希望帮到你。本文还有配套的精品资源点击获取