简介这套基于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/dxdvy/dt -(1/rho) * dp/dydp/dt -K * (dvx/dx dvy/dy)这里vx、vy是速度分量p是压力rho是密度K是体积模量。交错网格要解决的核心问题是p的导数应该落在哪个位置vx的导数又该落在哪个位置。如果把两者放在同一个网格点上那么所有导数都要做四点插值高频分量会被严重扭曲。交错网格的变量布局如下。网格位置存储变量物理含义(i, j) 整点p压力(i1/2, j) x半网格vxx方向速度(i, j1/2) y半网格vyy方向速度(i1/2, j1/2) 对角半网格可选剪切应力弹性介质扩展用这样vx点恰好处于两个相邻p点的正中间p对x的差分结果天然就落在vx的位置。反过来vx对x的差分结果又落在p的位置。不需要插值两个相邻点的中心差分就是二阶精度。如果网格够密这个二阶精度已经能很好地描述反射波相位。2.2 时间半层与差分更新骨架空间上错开之后时间上也采用半层更新。速度在n-1/2时刻压力在n时刻一个完整时间步分两步走先用当前压力更新速度到n1/2再用新速度更新压力到n1。这就是“交错”的第二个含义。下面是我经常用的二维声波更新骨架压缩包中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(nx1, ny); % x速度比压力多一个点 vy zeros(nx, ny1); % y速度也比压力多一个点 p zeros(nx, ny); for it 1 : 500 % 更新 vx相邻压力点的差分落在 (i1/2, j) vx(2:nx, :) vx(2:nx, :) - (dt / rho) / dx * (p(2:nx, :) - p(1:nx-1, :)); % 更新 vy同理落在 (i, j1/2) vy(:, 2:ny) vy(:, 2:ny) - (dt / rho) / dy * (p(:, 2:ny) - p(:, 1:ny-1)); % 计算速度散度并更新压力 dvx (vx(2:nx1, :) - vx(1:nx, :)) / dx; % 回到整点 dvy (vy(:, 2:ny1) - vy(:, 1:ny)) / dy; p p - dt * K * (dvx dvy); end这段代码里p(2:nx,:) - p(1:nx-1,:)取的是相邻两个压力值之差除以dx后差分位置正好在(i1/2,j)因此可以直接赋给vx。vy同理。更新压力时vx(2:nx1,:) - vx(1:nx,:)又把速度差放回到p的整点位置。整个循环里没有一处插值所有导数都是半网格上的中心差分。dt的选择不是一拍脑袋定的需要满足CFL条件后面第5章会专门讲。2.3 双相介质界面上的参数切换当模型由固体和流体两种介质组成时每个整点的rho和K不同。最简单也最常用的做法是在整点直接给p分配介质参数而半网格速度点的密度则取相邻两个整点密度的算术平均。这样处理保证了速度更新时对压力的差分不会跨越突变界面时产生单侧近似。界面上不需要引入额外的等效参数公式因为交错网格天然的变量错位已经让法向速度连续条件的离散化变得很自然。实际工程里如果介质参数差异过大比如固体速度3000 m/s、流体1500 m/s时间步长必须按最大速度来限制否则界面上的高频反射项会先发散。3. 压缩包模块拆解与Matlab实现流程3.1 文件功能对照把压缩包解开后最大的困惑不是代码难而是这些文件名各自承担什么角色。我按命名习惯和常见工程组织方式过了一遍把能确定的功能列在下面。文件名推断功能使用建议AcousticWE_RGD_Simulation.m声波方程交错网格模拟核心主模拟入口AcousticWE_RGD_Simulation.asv核心程序的历史自动保存忽略或备份TotalWES.m总波场能量统计与主控流程从它开始跑TotalWES.asv主控脚本的自动保存版本出问题时对比FCTforAW.m声波通量校正传输主流程调用FCTforEW.m弹性波通量校正扩展弹性波时用DCoef.m阻尼系数或PML衰减系数边界条件wigb.m灰色刻度波场快照绘图可视化data2.xls模型分层速度、密度参数输入文件Model2XVector.avi演示视频核对预期结果.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*5640米。如果你想让界面更陡或更平缓可以直接改变切分的行号。需要说明的是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, dy5 m空间采样小于最短波长的1/8dt0.0004 s时间步长按CFL条件缩小fm20 Hz雷克子波主频越高分辨率越好但频散越大界面深度640 m双相界面位置影响反射波到达时间PML层数30 格边界吸收厚度太少会看到边界反射如果波场快照中出现明显的条状条纹优先把dx减小一半对比一下。很多时候频散不是算法问题而是网格太粗视速度被拉长了。4.4 用wigb.m绘制快照运行结束后压力场P是三维矩阵取某一时刻用wigb绘图。wigb.m是业界常用的波场显示工具默认把矩阵的行当垂直坐标绘制出灰度波形。figure; wigb(P(:, :, 300)); % 画第 300 时间步的快照 title(t0.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_max3000 m/s若dx5 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后衰减说明边界吸收正常工作。改介质参数或网格步长后重跑一次并记录能量曲线的形态比每次肉眼判断快照更可靠。这套验证方法虽然简单却能把很隐蔽的参数错误提前暴露出来。本文还有配套的精品资源点击获取