基于交错网格有限差分的双相介质波场模拟与Matlab实现

简介：这套基于MATLAB平台的双相介质交错网格有限差分波场模拟程序，面向地球物理、声学、光学等领域研究波动传播的工程师与学生。程序将速度与压力分配到交错网格不同位置，可提高计算精度与稳定性，并针对双相介质交界面反射折射以及相态变化做了专门处理，适用于地震勘探、噪声控制、波导设计等场景。包内共有10个文件，包含6个M源码脚本、2个ASV自动备份、1个AVI演示视频和1个XLS数据表，整体约394KB，结构紧凑清晰，便于运行和二次开发。用户可通过调节网格尺寸或时间步长来控制模拟精度，结合演示视频和数据表格直观分析波场演化规律。目前已有259人学习下载，是掌握交错网格有限差分法、开展介质波场数值模拟的实用工具。

1. 双相介质波场模拟：为什么说交错网格是最合适的差分框架

拿到WaveFieldNumericalSimulation(StaggeredGrid).zip，最直接的想法是找主脚本改参数跑通。但如果你打算把这套双相介质波场模拟用到自己的地震勘探或水声场景里，建议先把交错网格的意图吃透。双相介质难题集中在固体-流体交界面的反射与透射，常规规则网格差分在界面附近极易产生寄生振荡，时间步稍大就发散。交错网格把速度和压力在空间上错开半个网格，在时间上也错开半个步长，数值频散比常规网格低一个量级。这套Matlab实现以AcousticWE_RGD_Simulation.m为模拟核心，FCTforAW.m做通量校正，wigb.m画快照，压缩包内还包含演示视频和模型参数。适合需要快速得到双相介质波场快照、又不想从零推写差分方程的人。

2. 交错网格有限差分：变量布局与差分格式

2.1 为什么选用速度-压力一阶方程

双相介质波场模拟的关键是界面处法向应力连续和法向速度连续。二阶位移方程在界面处需要同时约束位移、应力和速度，实现起来相当绕。把波动方程改写成一阶速度-压力方程组后，界面条件可以直接落在压力p与法向速度分量上，编程时只需要在界面网格处切换介质参数。以二维声波为例，控制方程是

• dvx/dt = -(1/rho) * dp/dx
• dvy/dt = -(1/rho) * dp/dy
• dp/dt = -K * (dvx/dx + dvy/dy)

这里vx、vy是速度分量，p是压力，rho是密度，K是体积模量。交错网格要解决的核心问题是：p的导数应该落在哪个位置，vx的导数又该落在哪个位置。如果把两者放在同一个网格点上，那么所有导数都要做四点插值，高频分量会被严重扭曲。交错网格的变量布局如下。

这样，vx点恰好处于两个相邻p点的正中间，p对x的差分结果天然就落在vx的位置。反过来，vx对x的差分结果又落在p的位置。不需要插值，两个相邻点的中心差分就是二阶精度。如果网格够密，这个二阶精度已经能很好地描述反射波相位。

2.2 时间半层与差分更新骨架

空间上错开之后，时间上也采用半层更新。速度在n-1/2时刻，压力在n时刻，一个完整时间步分两步走：先用当前压力更新速度到n+1/2，再用新速度更新压力到n+1。这就是“交错”的第二个含义。下面是我经常用的二维声波更新骨架，压缩包中AcousticWE_RGD_Simulation.m的核心循环与此一致。

% 二维声波交错网格速度-压力更新骨架 nx = 200; ny = 200; dx = 5; dy = 5; dt = 0.0006; % 空间步长 5 m，时间步长 0.6 ms rho = 2000; K = 4.0e9; % 密度 2000 kg/m3，体积模量 4 GPa c = sqrt(K / rho); % 纵波速度约 1414 m/s vx = zeros(nx+1, ny); % x速度比压力多一个点 vy = zeros(nx, ny+1); % y速度也比压力多一个点 p = zeros(nx, ny); for it = 1 : 500 % 更新 vx：相邻压力点的差分落在 (i+1/2, j) vx(2:nx, :) = vx(2:nx, :) - (dt / rho) / dx * (p(2:nx, :) - p(1:nx-1, :)); % 更新 vy：同理，落在 (i, j+1/2) vy(:, 2:ny) = vy(:, 2:ny) - (dt / rho) / dy * (p(:, 2:ny) - p(:, 1:ny-1)); % 计算速度散度并更新压力 dvx = (vx(2:nx+1, :) - vx(1:nx, :)) / dx; % 回到整点 dvy = (vy(:, 2:ny+1) - vy(:, 1:ny)) / dy; p = p - dt * K * (dvx + dvy); end

这段代码里，p(2:nx,:) - p(1:nx-1,:)取的是相邻两个压力值之差，除以dx后，差分位置正好在(i+1/2,j)，因此可以直接赋给vx。vy同理。更新压力时，vx(2:nx+1,:) - vx(1:nx,:)又把速度差放回到p的整点位置。整个循环里没有一处插值，所有导数都是半网格上的中心差分。dt的选择不是一拍脑袋定的，需要满足CFL条件，后面第5章会专门讲。

2.3 双相介质界面上的参数切换

当模型由固体和流体两种介质组成时，每个整点的rho和K不同。最简单也最常用的做法是：在整点直接给p分配介质参数，而半网格速度点的密度则取相邻两个整点密度的算术平均。这样处理保证了速度更新时对压力的差分不会跨越突变界面时产生单侧近似。界面上不需要引入额外的等效参数公式，因为交错网格天然的变量错位已经让法向速度连续条件的离散化变得很自然。实际工程里，如果介质参数差异过大，比如固体速度3000 m/s、流体1500 m/s，时间步长必须按最大速度来限制，否则界面上的高频反射项会先发散。

3. 压缩包模块拆解与Matlab实现流程

3.1 文件功能对照

把压缩包解开后，最大的困惑不是代码难，而是这些文件名各自承担什么角色。我按命名习惯和常见工程组织方式过了一遍，把能确定的功能列在下面。

.asv是Matlab在编辑过程中自动保存的文件，它不等于旧版项目，但往往保留着变量名或逻辑改动前的状态。我习惯把它复制一份改名保留，等新脚本跑出问题时再回去对比，比直接删除安全得多。

3.2 TotalWES.m主控脚本与调用链

按照命名习惯，TotalWES.m是总控脚本，负责读取模型、调用核心模拟和绘图。典型的调用路径是：读入data2.xls→ 扩展为二维介质参数 → 调用AcousticWE_RGD_Simulation.m→ 用FCTforAW.m做频散校正 → 输出波场快照。下面是我重构的调用骨架。

% TotalWES.m 主控脚本（结构恢复示例） data = xlsread('data2.xls'); % 读入分层模型参数 vp = data(:, 1); % 各层纵波速度 m/s rho = data(:, 3); % 各层密度 kg/m3 nx = 256; ny = 256; dx = 10; dy = 10; dt = 0.001; srcFreq = 25; % 雷克子波主频 Hz % 将分层参数映射到二维网格：前128行为固体，后为流体 vpGrid = ones(nx, ny) * vp(1); rhoGrid = ones(nx, ny) * rho(1); vpGrid(:, 129:end) = vp(2); rhoGrid(:, 129:end) = rho(2); % 调用核心模拟函数 [Pressure, Xcoord, Ycoord] = AcousticWE_RGD_Simulation(... vpGrid, rhoGrid, nx, ny, dx, dy, dt, srcFreq); % 通量校正：消除界面高频振荡 Pressure = FCTforAW(Pressure);

逻辑说明：vpGrid和rhoGrid是二维数组，每行每列分别代表一个网格点的介质属性。vpGrid(:,129:end)=vp(2)这一步把第129行到末尾的网格全部赋成第二种介质参数，相当于在y方向第1280米处设置水平界面。核心函数返回的Pressure是三维数组，维度为[nx, ny, nt]，分别对应x、y和时间步。FCTforAW的作用是对压力场做通量校正，它需要速度场信息来构造数值通量，调用时要注意传参顺序，否则校正效果的对比会失真。

3.3 DCoef.m与边界吸收

模拟区域如果没有吸收边界，波会从边界反射回来污染波场。DCoef.m生成的应该是边界衰减系数。一种常见的PML边界处理是在区域外侧增加几十层网格，每层施加阻尼系数

sigma(x) = sigma_max * (x / L)^2

其中L是PML厚度，x为该点到边界的距离。这个曲线是抛物线形状，最外侧的衰减系数最大。在Matlab里，DCoef.m可能返回一个与区域大小相同的衰减系数矩阵，在主循环中把它反复乘到压力或速度场上。使用时要特别注意：衰减系数不能加在内部区域，否则物理波场会被额外吸收；也不能只在压力上乘而不在速度上乘，那样会破坏方程的对称性。判断PML是否生效的方法是看快照晚期的波场是否只有向外传播的波，没有明显反射圆弧。

3.4 关于“相场”命名的误会

摘要里提到“相场”，容易让人联想到材料科学中的phase field模型。这里澄清一下：这个项目名称里的“相场”是指双相介质，也就是固体和流体两种物相，而不是求解相场方程。相场模型要额外引入序参量追踪界面演化，而这套程序只是把不同介质的密度和模量赋予不同网格点，不需要相场方程。读代码时如果发现界面处只是切换rho和K，不要误以为缺少了什么核心功能。

4. 跑通一次双相介质波场模拟：参数设置与快照输出

4.1 手工建立两层模型

不依赖data2.xls，先手工构造一个简单的两层模型，可以更快验证程序逻辑。用256×256网格，上层是固体，下层是流体。

nx = 256; ny = 256; dx = 5; dy = 5; vp = zeros(nx, ny); rho = zeros(nx, ny); % 固相：上部，速度 3000 m/s，密度 2200 kg/m3 vp(:, 1:128) = 3000; rho(:, 1:128) = 2200; % 流体：下部，速度 1500 m/s，密度 1000 kg/m3 vp(:, 129:end) = 1500; rho(:, 129:end) = 1000;

这里的界面设置在y方向第128个网格点上，实际深度为128*5=640米。如果你想让界面更陡或更平缓，可以直接改变切分的行号。需要说明的是，vp和rho作为输入给核心函数时，应该把rho也二维化。有些版本会用K = rho * vp^2来算体积模量，注意读核心程序时确认它到底接收的是速度还是模量。

4.2 雷克子波与震源加载

震源时间函数常用雷克子波，主频frequ决定了波传播的分辨率。生成方法如下。

nt = 2000; dt = 0.0004; t = (0:nt-1)*dt; fm = 20; % 主频 20 Hz delay = 0.08; % 延迟时间，保证子波起始平滑 w = (1 - 2*pi^2*fm^2*(t-delay).^2) .* exp(-pi^2*fm^2*(t-delay).^2); src = zeros(nx, ny); src(128, 128) = 1; % 震源位置坐标 for it = 1 : nt % 在每个时间步把雷克子波按权值加入压力场源点 Pressure(:,:,it) = ... AcousticWE_RGD_Simulation_step(...); % 示意 Pressure(128, 128, it) = Pressure(128, 128, it) + w(it); end

雷克子波的特点是零相位，主频处的能量集中，滞后delay秒再启动可以避免初始时刻的突变。震源位置放在模型中部，即(128, 128)附近，这样波前能同时覆盖固体和流体区域。实际代码里，AcousticWE_RGD_Simulation.m很可能已经内置了子波参数，你只需要修改srcFreq即可。

4.3 关键参数与状态对照

以下是调试过程中最常调整的参数表，建议运行前逐项确认。

如果波场快照中出现明显的条状条纹，优先把dx减小一半对比一下。很多时候频散不是算法问题，而是网格太粗，视速度被拉长了。

4.4 用wigb.m绘制快照

运行结束后，压力场P是三维矩阵，取某一时刻用wigb绘图。wigb.m是业界常用的波场显示工具，默认把矩阵的行当垂直坐标，绘制出灰度波形。

figure; wigb(P(:, :, 300)); % 画第 300 时间步的快照 title('t=0.12s 双相介质波场快照'); xlabel('X/m'); ylabel('Y/m'); % 如果觉得横纵方向反了，转置后绘图 figure; wigb(P(:, :, 300).'); % 转置后以Y为纵轴

第一张图中你会看到：震源发出的圆弧波先在上层固体中传播；遇到界面时，一部分反射回固体，一部分透射进流体；透射波在流体中速度更慢，圆弧曲率更大。第二张转置图则更适合观察反射波在垂直剖面上的延续。wigb默认会做归一化，不同时间步之间振幅对比不要靠颜色深浅，要看旁边的灰度条。

5. 排错与精度验证：CFL条件、边界吸收与FCT限流

5.1 CFL条件与NaN排查

如果程序运行到一半P变成NaN，第一反应是检查时间步长。二维交错网格声波模拟的CFL条件写作

dt <= dx / (c_max * sqrt(2))

其中c_max为整个模型的最大纵波速度，sqrt(2)对应二维情况。用上一章的两层模型计算，c_max=3000 m/s，若dx=5 m，则dt <= 5/(3000*1.414) ≈ 1.18 ms，取0.0004 s是安全的。但如果data2.xls中有更高速层，比如6000 m/s的致密岩石，时间步长就要减半。出现NaN时，我还会检查rho是否为零或负值，因为rho出现在导数前分母位置，一旦有零值，直接计算出Inf。

5.2 FCT限流如何压住数值频散

即使CFL满足，双相界面附近的波场仍可能出现高频拖尾。FCTforAW.m在这里起作用。FCT分为四步：先计算低阶通量和高阶通量，再计算反扩散通量，最后用限流器调制。调用方式是在主循环后处理压力场。

% 在 AcousticWE_RGD_Simulation.m 的时间步内 p_raw = p - dt * K * (dvx + dvy); % 原始压力更新 p_corr = FCTforAW(p_raw, vx, vy, dx, dy, dt); p = p_corr;

传入的vx和vy是当前时刻的速度场，FCTforAW需要它们来构造低阶和高阶数值通量。这个函数内部会计算出每个网格点的反扩散通量，并把通量限制在一个不产生新极值的范围内，从而消除伪振荡。验证它是否生效，可以对比校正前后同一条垂直剖线上的振幅曲线：FCT后的曲线界面反射波依然陡峭，但毛刺明显变少。

5.3 用总能量曲线做回归检查

我通常会在主控脚本里加一段总能量统计，用来判断模拟是否稳定。

% 波场总能量随时间变化（忽略边界影响） totalEnergy = squeeze(sum(sum(P .* P, 1), 2)); plot(totalEnergy); xlabel('时间步'); ylabel('总能量');

在波场未接触边界之前，总能量曲线应该基本平坦。如果曲线快速上升，说明有数值不稳定；如果曲线在波到达PML后衰减，说明边界吸收正常工作。改介质参数或网格步长后，重跑一次并记录能量曲线的形态，比每次肉眼判断快照更可靠。这套验证方法虽然简单，却能把很隐蔽的参数错误提前暴露出来。

本文还有配套的精品资源，点击获取