简介面向MATLAB环境下弹性波动方程有限差分数值模拟的完整源码包由达摩老生出品并亲测校正适合地震学、勘探地球物理等领域的研究生、工程师及相关课程学习者作为正演模拟的起步模板。该套代码将有限差分法的计算流程完全展现在源码层级便于使用者从底层理解波场如何随空间网格和时间步长逐步递推。代码覆盖了模型参数初始化、网格剖分和边界条件处理、震源加载、波场迭代计算以及地震记录和波场快照的可视化等关键环节各功能以独立脚本和函数组织模块边界清晰既支持按需调用某一部分也方便在此基础上扩展不同震源类型、吸收边界或速度模型。压缩包共31个文件主体为28个.m脚本和函数文件另有2个.mat数据文件用于辅助模拟运行并附1份docx文档整体大小仅55KB轻量便携。目前已有1556人浏览学习无论用于入门研读还是作为日常正演工具这套亲测可用的源码都具备很强的参考价值。1. 从一道波开始为什么用有限差分法模拟弹性波动方程地震波在层状介质里走一条近路结果最先到的不是直达波而是一路折射、反射、透射叠出来的复合波场。想看清弹性波在地下怎么传播解析解只存在于均匀、分层或特别对称的模型里一旦界面倾斜、速度横向变化数值模拟就成了唯一能“看”的手段。有限差分法Finite Difference Method简称 FDM是其中性价比最高的一种把连续空间切成一堆规则网格把波动方程里的偏导数换成差分近似然后用显式时间递推一步步把波场推进到目标时刻。它不像有限元那样要组装大矩阵也不像伪谱法那样依赖 FFT写起来直观、调参路径清晰在二维三维弹性波动方程的时间域模拟中一直是被使用最多的方案之一。文章后面所有内容不针对某个特定版本我用 R2023b 做演示你在 R2020a 到 2026a 之间跑结果都不会有本质差别。2. 弹性波动方程的有限差分格式先立起离散骨架2.1 解位移还是解速度—应力两个出发点的差别弹性波动方程最常用的两种写法是位移二阶方程和一阶速度-应力方程。位移形式长这样[ \rho \frac{\partial^2 u_i}{\partial t^2} \frac{\partial \sigma_{ij}}{\partial x_j} f_i ]其中 (\sigma_{ij}\lambda \delta_{ij} \nabla \cdot \mathbf{u} \mu(\partial_i u_j \partial_j u_i))。这种形式物理意义清楚但在数值上要同时离散二阶时间导数和位移的空间二阶导数对网格间距要求偏严且遇到自由表面或介质分界面时天然会产生明显的虚假振荡。一阶速度-应力方程把未知量拆成速度分量和应力分量[ \rho \frac{\partial v_i}{\partial t} \frac{\partial \sigma_{ij}}{\partial x_j} f_i ][ \frac{\partial \sigma_{ij}}{\partial t} \lambda \delta_{ij} \nabla \cdot \mathbf{v} \mu(\partial_i v_j \partial_j v_i) ]这里的优势是空间导数只有一阶可以对速度与应力分别放到交错位置取样使得二阶精度的空间模板也能达到近似四阶的相速度表现。常见做法是采用 Virieux 提出的交错网格方案将正应力放在网格角点剪应力放在网格边中点速度分量各占半边。下面的表格是二维各向同性介质中变量分布的典型约定手工写代码前最好先在纸上把这个布局画一遍。变量交错网格位置常用精度阶(v_x)((i, j1/2))4(v_z)((i1/2, j))4(s_{xx})((i1/2, j1/2))4(s_{zz})((i1/2, j1/2))4(s_{xz})((i, j))42.2 空间离散四阶中心差分的模板与系数对二维情况把速度-应力方程拆开实际需要更新五个物理量。给 (v_x) 举例更新式是[ v_x^{n1/2}(i, j1/2) v_x^{n-1/2}(i, j1/2) \frac{\Delta t}{\rho} \big( D_x s_{xx} D_z s_{xz} \big) ]其中 (D_x) 和 (D_z) 是对应方向的空间差分算子。如果使用四阶精度则[ D_x s_{xx} \approx \frac{1}{12 \Delta x}\left[ -s_{xx}(i, j1) 8s_{xx}(i,j\frac12) - 8s_{xx}(i,j-\frac12) s_{xx}(i,j-1) \right] ]这个五点半边模板在交错网格里是常用配置写成 MATLAB 代码时可以直接对矩阵切片操作避免写双层 for 循环。下面的片段展示了如何用数组运算把 (D_z) 实现出来% 输入场变量 Adx 为网格间距 % 输出 A 对 z 方向的四阶中心差分 function dAdz dz_4(A, dx) dAdz zeros(size(A)); % 内部点五点半边模板 dAdz(3:end-2, :) (1/12/dx) * ( ... - A(5:end, :) 8*A(4:end-1, :) ... - 8*A(2:end-3, :) A(1:end-4, :) ); % 边界点退化为二阶精度 dAdz(2, :) (A(3,:) - A(1,:)) / (2*dx); dAdz(end-1,:) (A(end,:) - A(end-2,:)) / (2*dx); end这段代码的思路是先用整体切片构造出内点的高阶差分再把边界处单独降阶处理。边界降阶是为了避免模板越界同时还能保持整体计算的稳定代价只是边界附近精度损失。实际做深反射地震正演时通常会把吸收边界或 PML 放在外侧内部区域仍然保持四阶。2.3 时间递推与 CFL 稳定性条件时间方向用蛙跳格式即速度的 (n-1/2) 时刻和 (n1/2) 时刻前后穿插应力按整时刻更新。这样更新一个时刻只需要两次乘加不需要解线。稳定性条件是不可省略的一环[ \Delta t \le \frac{0.606 \min(\Delta x, \Delta z)}{\sqrt{v_p^2 v_z^2}} ]实际计算中习惯把 CFL 数记作 (\alpha)令分子为 (\alpha \Delta h)取 (0.3\sim 0.5) 之间。CFL 太大波场在数千步后会出现树枝状高频噪声CFL 太小模拟停滞浪费计算时间。提示判断时间步长是否合理不要只看理论 CFL要看最高速度模型的 (v_p) 值不能只看浅层低速区。2.4 人工边界从海绵吸收到 PML无限半空间只存在于理论里数值模拟必须把边界切掉。最简单的办法是在模型外围加一圈吸收层让波幅在边界附近逐步衰减。常见做法是采用一个随距离变化的衰减因子 (\gamma(x))每次更新后乘上一次[ \mathbf{w}^{n1} \exp(-\gamma \Delta t) \cdot \mathbf{w}^{n1} ](\gamma) 从边界内第一格开始缓慢增长到边界处最大值。优点是好写、几乎不占内存缺点是反射残余有几十分贝对高精度研究不太够。PML 在二维弹性介质里更常用但代码量会多出很多而且角点处理容易出错。对大多数教学和初阶研究来说海绵吸收已经足够。3. 用 MATLAB 写一个最小可运行的二维弹性波模拟器3.1 搭建模型与网格把速度和密度填进矩阵这一步最关键的不是写循环而是把坐标、网格步长、时间步长三者的关系定死。先固定横向与纵向网格点数再设定 (dx) 和 (dz)然后根据最大速度算 (\Delta t)。我这里给出了一个完整的最小代码框架% 模型两层介质界面在 z 200 m nx 601; nz 401; dx 5; dz 5; vp 2800 * ones(nz, nx); % P 波速度 m/s vp(201:end, :) 3600; vs vp / 1.73; % S 波速度按泊松比换算 rho 2200 * ones(nz, nx); % 密度 kg/m^3 rho(201:end, :) 2600; mu rho .* vs.^2; lambda rho .* (vp.^2 - 2*vs.^2); % 时间步长取空间步长的一半来获得充足余量 dt 0.5 * min(dx, dz) / max(vp(:)); nt 2000; % 波场数组 vx zeros(nz, nx); vz zeros(nz, nx); sxx zeros(nz, nx); szz zeros(nz, nx); sxz zeros(nz, nx);这里用vp/1.73意味着泊松比为 0.25也就是 (v_p/v_s\sqrt{3})现实中很多沉积岩接近这个值。lambda和mu写成矩阵而不是标量是为了后面允许横向和纵向变速不需要重新建模。计算 (dt) 时用 0.5 乘计算出的上界是把安全系数拉满避免材料参数离散误差导致的不稳定。模型中界面用整行索引赋值代表水平两层情况。如果要做透镜体或断层只需要构造一个逻辑矩阵layer vp 3000;然后把该区域内的值直接替换。3.2 震源加载与主循环雷克子波和交错网格更新震源一般用 Ricker 子波它的特点是主瓣尖锐、往两侧迅速衰减适合模拟短脉冲破裂[ s(t) \left(1 - 2\pi^2 f_0^2 (t-t_0)^2\right) \exp\left(-\pi^2 f_0^2 (t-t_0)^2\right) ]把子波物理地加在应力场上模拟爆炸源更常见的做法是同时加到正应力分量 (s_{xx}) 和 (s_{zz}) 上以产生对称的 P 波。下面给出主循环的一部分% 雷克子波参数 f0 20; % 主频 Hz t0 1 / f0; % 时延取约为 1/f0 for it 1:nt t (it-1) * dt; amp (1 - 2*(pi*f0*(t-t0))^2) * exp(-(pi*f0*(t-t0))^2); % 进行空间差分的数组操作只展开 sxx 为例 dsxx_dx dx_4(sxx, dx); dsxz_dz dz_4(sxz, dz); dsxz_dx dx_4(sxz, dx); dszz_dz dz_4(szz, dz); vx(2:end-1, 2:end-1) vx(2:end-1, 2:end-1) dt ./ rho(2:end-1, 2:end-1) ... .* (dsxx_dx(2:end-1, 2:end-1) dsxz_dz(2:end-1, 2:end-1)); vz(2:end-1, 2:end-1) vz(2:end-1, 2:end-1) dt ./ rho(2:end-1, 2:end-1) ... .* (dsxz_dx(2:end-1, 2:end-1) dszz_dz(2:end-1, 2:end-1)); % 应力更新 sxx(2:end-1, 2:end-1) sxx(2:end-1, 2:end-1) dt .* ( ... (lambda2*mu)(2:end-1, 2:end-1) .* dvx_dx(2:end-1, 2:end-1) ... lambda(2:end-1, 2:end-1) .* dvz_dz(2:end-1, 2:end-1) ); % 震源加载在 st0 sxx(round(nz/2), round(nx/2)) sxx(round(nz/2), round(nx/2)) amp; szz(round(nz/2), round(nx/2)) szz(round(nz/2), round(nx/2)) amp; end上面的代码里 (dvx/dx)、(dvz/dz) 需要先基于当前速度场计算一次再进入应力差分这里为了节省篇幅没有重复展开。核心是先更新速度再更新应力震源只在相应时刻加在应力上不参与差分。这样能避免子波强度被空间差分算子放大。子波的时延 (t01/f0) 能保证初始时刻子波能量接近 0从而不产生人为的高频冲击。3.3 震源频率与网格分辨率的第一次校准许多人模拟结果出现强烈频散问题往往不是网不够密而是主频过高。经验法则是每个最短波长至少覆盖 8 到 12 个网格点。用公式估算[ f_{max} \approx 2.5 f_0 ]若 (f_020) Hz、(v_s1600) m/s则最小 S 波波长为 (1600/5032) m。用 (dx5) m 时一个波长只覆盖 6.4 个网格点频散已经在肉眼可见的边缘。遇到这种情况要么降低 (f_0)要么加密网格。如果不想在第一次试跑就花掉太多时间可以先用 2/3 网格尺寸跑一遍看波场如果前缘出现不明抖动再把主频下调 20%。4. 从“能出图”到“能对比”算例设计、可视化与参数敏感性4.1 均匀模型的解析参考解拿一个已知答案兜底写完代码第一步不是看复杂模型而是跑一个均匀全空间模型和解析解对比。均匀无限介质中点源产生的 P 波位移场可用远场近似大致描述[ u_r(r,t) \approx \frac{1}{4\pi\rho v_p^3 r} \frac{\partial^2 s}{\partial t^2}\left(\frac{t-r}{v_p}\right) ]我的习惯是在炮点旁边和远前方分别放两个接收器提取垂直速度分量将振幅归一化后叠加在一起与解析波形对比。代码里只需要加两行rec_x 100:100:500; rec_z round(nz/2)*ones(size(rec_x)); Vz_record zeros(nt, numel(rec_x)); Vz_record(it, :) vz(sub2ind([nz nx], rec_z, rec_x));这里的Vz_record会在时间上形成一张“道记录”可以直接用wiggle风格画成地震剖面。解析参考解的作用不是验证整个模型是否正确而是用来排除算错的地方假如波形时间错位基本是速度或距离给错假如波形形态错位多半是差分模板方向写反了。4.2 波场快照与分量分离看 P 波和 S 波的第一步弹性波模拟最让人兴奋的时刻就是看到圆环状波前从震源扩散。用imagesc显示 (s_{zz}) 分量可以直观显现 P 波前而显示剪应力 (s_{xz}) 时S 波特征会更明显因为介质受纯剪切运动时正应力变化微弱剪应力分量直接把这一特征突出出来。常用的可视化代码figure; subplot(1,2,1); imagesc((0:nx-1)*dx, (0:nz-1)*dz, szz); axis equal tight; colormap(gray); caxis([-max(abs(szz(:))) max(abs(szz(:)))]); title(正应力 szz); subplot(1,2,2); imagesc((0:nx-1)*dx, (0:nz-1)*dz, sxz); axis equal tight; colormap(gray); caxis([-max(abs(sxz(:))) max(abs(sxz(:)))]); title(剪应力 sxz);caxis使用对称范围很重要因为波场有正有负颜色映射不对称时会把零相位变成亮色看起来像多了个直流分量。做波场快照时我还会把时间点选在波刚穿过最深界面时此时反射与透射都没有完全断开信息最丰富。4.3 参数敏感性CFL、网格间距和主频分别改多少参数常见区间变化对结果的影响CFL 数0.20.6偏大时高频网格噪声快速增长偏小时耗机时网格间距最短 S 波波长的 1/101/20太粗导致频散太细使 (dt) 被迫降低子波主频勘探中常用 1040 Hz过低时分辨率不足过高时频散明显模型边界吸收层厚度 2050 网格太薄反射残影进入内部波场这里有一个反直觉的结论网格间距加倍计算量并不是增加两倍而是空间方向和时间方向同时增加总计算量约增加 (2^3) 倍。所以调低分辨率省时间的同时往往也要接受频散对波形的影响。5. 让模拟更可靠用网格收敛测试与能流检验压住误差5.1 网格收敛阶验证的操作细节高精度格式的宣传归宣传真正验证代码误差阶数的办法是对同一模型用两种网格间距各跑一遍。假设网格从 (dx) 加密到 (dx/2)如果空间格式是四阶总误差应大致降为原来的 (1/16)。执行时把接收道记录做 L2 范数比较err sqrt(sum((Vz_coarse - Vz_fine(1:2:end, :)).^2, 1)) / sqrt(sum(Vz_fine(1:2:end, :).^2, 1));这里把细网格结果每隔一个点采一次然后与粗网格结果对齐。误差小于 2% 说明格式基本正确误差在 5% 到 10% 时大概率是震源位置或接收器索引赋值错误而不是格式自身的问题。5.2 能量衰减曲线边界和吸收层同时在骗人内心活动边界吸收层好不好肉眼容易骗自己。更客观的办法是持续监测模型总能量[ E \frac12 \int_\Omega \left( \rho |v|^2 \sigma:\epsilon \right) d\Omega ]数值上每次循环后做一次全网格求和然后取 log10 画曲线。如果曲线平稳说明介质内部没有异常能量源如果曲线在某个时间点开始上升说明数值不稳定已经发生通常需要对 (dt) 或网格做回退。吸收层性能看的是能量曲线尾部下降速率衰减太快说明吸收过强造成了虚假反射衰减太慢说明返回波会被看到。E(it) 0.5 * sum(rho(:) .* (vx(:).^2 vz(:).^2)) * dx * dz;这个量还能用来对比不同 PML 参数的吸收效果对比时固定内部模型不放炮只加一个边界入射波。我一般把时间步和主频率保持不变只改吸收层厚度画在同一张图上比较尾部斜率。5.3 MATLAB 调参与跨版本兼容的几个实际做法我常遇到的角度是代码在一个版本的 MATLAB 上跑得好好的换个版本波浪号~或隐式扩展.*的行为不一致。R2023b 及以后的中期版本对数组操作做了优化但旧版本对end在切片中的使用仍有限制建议少用多层end嵌套。还有一点二维弹性波模拟最耗时间的是循环内频繁的大矩阵切片可以在开头把rho和mu预先取倒数保存成inv_rho和inv_mu这样循环里只剩乘法没有除法。如果机器内存紧张可把模型裁剪到感兴趣区并给吸收层留 20 个网格的边距。再进一步优化时可以把主循环改成parfor但必须注意Vz_record的写入位置要按循环次数索引不能直接对波场变量并行赋值否则 MATLAB 会因为数据依赖报错。调试时顺手利用warning(off, MATLAB:colon:nonIntegerIndex)并不是好习惯看到非整数索引警告通常意味着时间步或索引计算中隐藏着 bug。运行速度基线大概是600×600 网格、2000 步、单精度数据在 R2023b 下约需要十几秒到半分钟具体与 CPU 和内存带宽相关如果时间超过两分钟优先检查是不是在里面偷偷种了ones或zeros大矩阵。本文还有配套的精品资源点击获取