简介一套基于Matlab的双相介质波场数值模拟程序采用交错网格有限差分法求解波动方程兼顾计算精度与稳定性适合正在学习计算波动理论、地震勘探或声学模拟的研究生与工程技术人员参考使用。压缩包共10个文件以6个m源码文件为核心配有2个asv备份文件、avi演示录像和xls数据表格整体约394KB轻量易用。已有259人学习下载表明其在相关领域具有实际参考价值。程序涵盖总波场能量记录、声波方程模拟、FCT修正等关键脚本通过修改网格尺寸、时间步长等参数可在普通电脑上灵活权衡精度与计算开销从源码、辅助数据到可视化输出结构完整便于结合录像观察波在固体、流体等不同相态介质交界处的反射、折射与透射过程理解交错网格抑制数值色散的原理标签中的“相场”也提示了跨相态界面问题的处理思路为需要开展波场正演、介质参数反演或教学演示的读者提供了良好的二次开发基础。1. 交错网格波场模拟Matlab 里先解决“伪模”再看波场“WaveFieldNumericalSimulation(StaggeredGrid).zip”这个文件名信息量其实比很多论文标题都大StaggeredGrid 是交错网格WaveFieldNumericalSimulation 是波场数值模拟再加上 matlab、介质、波场、相场四个标签基本就是一套完整的地震波/声波模拟工具链。做勘探地球物理、超声检测、声学仿真的工程师大概率都自己攒过这样的程序。交错网格Staggered Grid是速度-应力形式波动方程的主流离散方式它能天然抑制常规中心差分出现的棋盘格伪模让低频段波形干净得多。这篇内容就把这套方案的原理、Matlab 实现、介质参数设置和相位提取技巧串起来适合正要写第一个有限差分程序、或者已经写了但波场里有明显“锯齿”和“棋盘”的读者。2. 交错网格的差分原理与离散系数为什么速度要错开 1/2 网格2.1 速度-应力方程与普通中心差分的棋盘格问题很多从标量波动方程起步的读者会习惯写这样的二阶形式∂²p/∂t² c² ∇²p然后对时间和空间都做中心差分。这种做法实现简单但在非均匀介质里需要处理二阶导数的“波动混合项”而且当速度差异大时常规中心差分很容易在奇偶节点之间产生解耦形成棋盘格形状的寄生振荡。这种伪模不是数值误差大小的问题而是离散格式本身“放行”了物理上不存在的短波长模式。交错网格换了一条思路先引入速度分量 vx、vz 和应力/压力 p把二阶波动方程改写成一阶双曲系统ρ ∂vx/∂t ∂p/∂x ρ ∂vz/∂t ∂p/∂z ∂p/∂t λ (∂vx/∂x ∂vz/∂z)这里的 ρ 是密度λ 是体积模量。速度分量被放在整网格点之间的半网格位置上p 放在整网格点上。导数运算不再需要跨越两个整网格去取中心差而是在相邻的整网格和半网格点之间做差分信息传播路径变短奇偶节点之间的数值联系被重新建立起来。对比项常规中心差分交错网格变量定义位置全部在整数网格压力整数、速度半网格空间导数跨度跨 2 个网格跨 1 个网格起算棋盘格伪模容易出现天然抑制存储开销少一个速度数组需多存 vx、vz略增边界处理单变量边界较直白速度、压力边界错位处理这也是为什么在 matlab 里做波场模拟尤其是有强介质分界面的时候我更倾向直接用速度-应力交错网格而不是高阶 Laplacian 方案。2.2 交错网格差分系数同一套公式3 阶和 4 阶差一倍交错网格的空间导数是“错位采样”的。以∂p/∂x为例输出位置在半网格 i1/2 处二阶精度格式长这样∂p/∂x ≈ (p(i1) - p(i)) / Δx这是最简单的交错差分代码里几乎不用付出额外计算成本。如果要让波场更光滑、频散更小可以升级到四阶∂p/∂x ≈ (9/8)(p(i1) - p(i)) / Δx - (1/24)(p(i2) - p(i-1)) / Δx对应的系数组是 a1 9/8 ≈ 1.125a2 -1/24 ≈ -0.04167。很多第一次写交错网格的人会惯性套用常规中心差分的 -1/2、1/2 系数结果波场出现奇怪的振荡其实是系数表用错了。交错网格的系数是为了“输出点位于输入点之间”设计的和普通中心差分不是一套来源。精度阶数系数 a1系数 a2系数 a32 阶1.0004 阶1.125-0.04166706 阶1.171875-0.0651040.0046875同样要注意时间推进的 CFL 稳定性条件。二维声波交错网格里实际跑 matlab 时可以用一个偏保守的经验公式dt ≤ 0.6 * Δx / (vmax * Σ|a_k|)其中 vmax 是模型里最大波速Σ|a_k| 是空间差分系数绝对值之和。二阶格式里这个和是 1.0四阶是 1.1667CFL 限制反而比二阶更严格一点换精度同时也要缩小 dt。2.3 为什么参数顺序会影响整个模拟成功率另一个容易被忽略的点是网格间距和震源主频的匹配。交错网格虽然能抑制伪模但如果每个波长里采样点太少数值频散依然会让波场尾部拖出一串“小尾巴”。经验上每一最短波长内至少要有 8 到 10 个网格点。网格确定后最高可靠主频 fmax 受控于最小速度和网格间距fmax ≈ vmin / (8 * Δx)所以模型设计顺序应该是先确定目标频带 → 再挑选 Δx → 最后用 CFL 条件推 dt。反过来先瞎选 dt 再调网格的流程经常会卡在“波场发散”和“频散严重”之间反复试错。3. 用 Matlab 写一个速度-应力交错网格模拟循环3.1 模型参数和网格初始化的最小代码在 matlab 里实现交错网格并不复杂关键是变量排列和更新顺序要清晰。下面这段代码定义了一个 240 × 240 网格的二维模型默认介质速度 2000 m/s中央放一个半径为 20 个网格的高速透镜体% 网格与时间参数 nx 240; nz 240; dx 2; dz 2; % 空间步长单位 m x (0:nx-1) * dx; z (0:nz-1) * dz; vmax 3000; % 模型最大速度用于 CFL dt 0.6 * dx / (vmax * 1.0); % 二阶格式系数和 1.0 nt 1000; % 总时间步 % 介质参数 c zeros(nx, nz); c(:) 2000; % 背景速度 m/s [xx, zz] meshgrid(x, z); idx (xx - 120).^2 / (20^2) (zz - 80).^2 / (20^2) 1; c(idx) 2500; % 高速透镜 rho zeros(nx, nz) 2200; % 密度 kg/m^3 lambda rho .* c.^2; % 体积模量 % 定义交错网格变量 P zeros(nx, nz); % 压力定义在整网格 Vx zeros(nx, nz); % x 方向速度定义在 x 半网格 Vz zeros(nx, nz); % z 方向速度定义在 z 半网格提示这里把所有变量都声明成同样尺寸的二维数组交错效果通过索引偏移实现不是真的把数组切一半。这样做的好处是矩阵运算不需要做维度重排速度更快坏处是边界附近的半网格点在更新时会被“留空”正好可以配合吸收边界一起处理。3.2 主循环先更新速度再更新压力交错网格时间推进用的是蛙跳格式。每一步先由压力梯度更新速度分量再由速度散度更新压力。发震源直接用 Ricker 子波加到压力数组的对应点上% Ricker 子波参数 fc 25; % 主频 Hz t (0:nt-1) * dt; ricker (1 - 2*(pi*fc*(t-1/fc)).^2) .* exp(-(pi*fc*(t-1/fc)).^2); % 震源位置位于高速透镜上方 isrc 120; jsrc 40; for it 1:nt % 1. 压力梯度更新速度注意半网格索引位移 Vx(:, 1:end-1) Vx(:, 1:end-1) - (dt ./ rho(:, 1:end-1)) ... .* (P(:, 2:end) - P(:, 1:end-1)) / dx; Vz(1:end-1, :) Vz(1:end-1, :) - (dt ./ rho(1:end-1, :)) ... .* (P(2:end, :) - P(1:end-1, :)) / dz; % 2. 速度散度更新压力 div (Vx(:, 2:end) - Vx(:, 1:end-1)) / dx ... (Vz(2:end, :) - Vz(1:end-1, :)) / dz; P(2:end-1, 2:end-1) P(2:end-1, 2:end-1) ... - dt * lambda(2:end-1, 2:end-1) .* div(1:end-1, 1:end-1); % 3. 加载 Ricker 震源 P(jsrc, isrc) P(jsrc, isrc) ricker(it); % 4. 吸收边界处理下节给出 P P .* damp; end这里的核心逻辑有三处需要注意。第一Vx 的更新用的是 P 在 x 方向上的相邻差分Vx(:,1:end-1)对应的是P(:,2:end) - P(:,1:end-1)这正是“速度定义在半网格”的索引表达。第二P 的更新只用内部区域2:end-1因为最外圈的速度值是无效边界强行参与计算会引入反射。第三震源加载必须放在边界衰减之前还是之后取决于边界带离震源的距离如果吸收带离得足够远顺序无所谓但放在主更新完成后更直观。3.3 简单吸收边界衰减带比零边界好用模拟区域边缘如果直接截断波场到达边界会反射回来干扰有效记录。完整的解法是 PML 和 CPML但对多数二维测试模型衰减带已经够用。在进入时间循环前生成一个 damping 系数数组npad 20; % 吸收带宽度 alpha 0.02; % 衰减强度 damp ones(nx, nz); for i 1:nx for j 1:nz % 按到四边的距离计算衰减 d min([i-1, nx-i, j-1, nz-j]); if d npad damp(i, j) exp(-(alpha * (npad - d))^2); end end end衰减强度 alpha 的取值很关键。太小不起作用波场反射后会在窗口里造成大面积干涉条纹太大则会因为压力在边界区被快速吸收形成一个新的“人工反射面”。一般 alpha 在 0.01 到 0.05 之间且吸收带宽度至少要覆盖一个主频波长。如果想更省事也可以把衰减带函数写成一个子函数方便不同模型复用。4. 介质模型与波场记录从均匀速度到分层模型的参数设计4.1 分层介质与异常体的 matlab 快速建模地质模型不可能一直用均匀速度场。常见的做法是用速度网格配合掩膜数组生成分层结构然后再叠加上高速或低速异常体。这样一个二维模型在 matlab 里的基本写法是% 速度模型两层介质 低速透镜 c2 ones(nx, nz) * 1800; % 第一层速度 c2(:, 140:end) 2400; % 第二层速度深度 z 超过 140 网格 % 在两层交界处放一个低速透镜 [idx2] ((xx-120).^2 / 15^2 (zz-150).^2 / 15^2) 1; c2(idx2) 1500; % 密度用经验关系生成 rho2 0.31 * c2.^0.25 * 1000; lambda2 rho2 .* c2.^2;低速透镜是波场模拟里特别值得加的对象。低速体内部波速低、波长短最容易触发数值频散同时它在层状背景中会产生强反射和绕射是检验吸收边界和相位提取效果的好标的。密度这里用了一个经验幂律关系不追求地层学精度但保证了速度和密度之间没有矛盾量级。4.2 参数设多少参考值表和设计顺序模型参数并不是越大越好。对同一套 240×240 网格dx 从 2 提升到 1计算量上升四倍但物理上的新信息可能并没有增加多少。推荐先把介质最小速度和目标主频钉死再反推网格间隔目标主频最小速度推荐 dx推荐 dt二阶10 Hz1500 m/s10 m约 4 ms25 Hz1500 m/s5 m约 2 ms50 Hz1500 m/s2 m约 0.8 ms100 Hz1500 m/s1 m约 0.4 ms同样重要的是把总时长设成足够波场穿越整个模型。可以用nt 2 * sqrt((nx*dx)^2 (nz*dz)^2) / vmin / dt先估一个值再手动加 20% 余量。dt 如果设置过大模拟通常不是缓慢恶化而是直接 NaN出错时先查 CFL 条件而不是查代码逻辑。4.3 波场快照记录与合成地震道提取模拟产生的波场如果不记录就是白算。常见的做法是把每个采样步长都存下来但那样内存涨得飞快。更实用的方案是每 N 步存一张快照以及在固定接收点记录完整时间序列% 接收点坐标 rec_pos [30, 60; 30, 120; 30, 180; 180, 120]; nrec size(rec_pos, 1); % 预分配记录数组 snap_interval 20; nsnap floor(nt / snap_interval); snapshots zeros(nx, nz, nsnap); seismo zeros(nrec, nt); for it 1:nt % ……省略主循环更新…… % 保存检波点数据 for ir 1:nrec seismo(ir, it) P(rec_pos(ir, 1), rec_pos(ir, 2)); end % 保存波场快照 if mod(it, snap_interval) 0 snapshots(:, :, it / snap_interval) P; end end接收点如果正好压在交错网格的速度点或边界点上取值可能不是整数网格压力值所以布置时要留几格安全距离。检波点数据后面画 wiggle 剖面或者做相位分析都是直接用的这段矩阵。4.4 频散排查看到拖尾先检查每波长采样数模拟结果最常见的特征异常是波前面后跟着一串周期性的“波纹”这不是介质里真实存在的波而是空间频散。排查时先把波毛估估网格支持的最高频率fmax vmin / (λ_min / (Δx))当最短波长里少于 8 个网格点时频散会非常明显少于 5 个网格点时波场基本不能看。解决办法有三个方向减少 dx、降低震源主频、或者把空间差分提升到四阶。前两个是调整参数第三个是修改差分算子。四阶格式代价不大但主循环里所有P(:,2:end)-P(:,1:end-1)都要改成包含二阶修正项的形式边界索引也要同步扩大 padding不是只改一行系数。5. 把波场变成相场用瞬时相位做界面识别与精度验证5.1 用 Hilbert 变换从快照里提取瞬时相位模拟完成后压力场只是一个数值矩阵但波的“相位”却携带了界面位置和波型变化的关键信息。这里所说的相场就是指对波场快照做复数追踪后得到的瞬时相位分布。matlab 里做这一步非常直接对每一道信号做 Hilbert 变换再取相位角% 取某一张快照例如第 50 张 field squeeze(snapshots(:, :, 50)); % 沿着水平方向对逐道信号做解析延拓 analytic hilbert(field.).; phase_snap angle(analytic); % 瞬时相位 unwrapped unwrap(phase_snap, [], 2); % 按行展开这里对快照做 Hilbert 变换而不是对时间序列做原因是我们关心的是某个时刻波场沿空间方向的相位结构。hilbert(field.).是为了让 matlab 默认沿第二维做变换最终想得到的是波峰和波谷之间的关系。相位剖面图里介质界面处会出现明显的相位过零也就是从正到负或从负到正的翻转位置这些翻转点的空间坐标应该和速度模型里的界面吻合。5.2 相位剖面和模拟精度验证的三步对比用相场做验证我一般做三步。第一步是分别画出第 40、60、80 张快照的瞬时相位图看透镜体边界上是否有一条稳定的相位翻转带第二步沿着透镜左右两侧各取一条水平相位曲线用find(diff(sign(unwrapped)) ~ 0)找出过零点位置再和真实边界坐标比对第三步比较二阶和四阶格式下提取到的过零偏移量如果偏移超过 2 个网格说明空间离散精度不够此时再回查参数表调整 dx 或差分阶数。相位过零位置比压力幅值更稳定受震源强度影响小在信噪比变差时依然能保持较好的可辨识度。读完这一套流程再看波场快照的顺序应该是先看棋盘格有没有消失再检查边界反射是否被压制最后用瞬时相位剖面确认介质界面的空间位置是否准确。这三项都通过了交错网格模拟程序才算真正闭环。本文还有配套的精品资源点击获取