1. 整体设计与思路拆解1.1 声波数值模拟的核心需求这阵子好几个朋友问我声波时域差分怎么做我干脆把自己一套经验整理出来顺便把关键坑位标出来。所谓“利用matlab编程实现的声波的时域差分方程”说白了就是用Matlab把描述声波传播的偏微分方程离散成差分方程然后一步步往前推得到波场随时间演化的过程。这套东西在室内声学仿真、超声无损检测成像、声呐信号模拟以及地震波传播研究里都有直接应用也是很多做计算声学方向的人最早接触的数值仿真实验之一。你把声波在空气或水里的传播想象成水波丢石头后的涟漪只是声波是压缩波介质中的压力在局部偏高和偏低之间振荡能量以一定速度向外推进。要捕捉这个过程如果介质均匀、边界规则可以用解析解但实际工程遇到的是复杂形状、分层结构、障碍物绕射、非均匀介质解析解根本写不出来。这时只能把空间和时间都切碎直接在网格上用差分近似代替微分一步一步模拟波前怎么走、碰到界面怎么折射反射、进到障碍物后面怎么绕射衰减。Matlab做这件事的优势在于矩阵运算本身就是向量化的编程模型声波场每个时刻无非就是一个数组前后时刻之间做线性组合就能得到新数组天然契合。再加上内置的绘图函数你可以边算边看波场演化调试起来非常直观。1.2 时域差分方法为什么适用很多人上来就问频域方法做声学不是更成熟吗为什么非要用时域差分我的看法是这两种路线解决的问题有重合但不完全一样。频域方法通常假设场量随时间按简谐规律变化把波动方程化简为亥姆霍兹方程然后求解稳态响应。这对于单一频率或者窄带信号很有效比如某个固定频率的设备噪声对房间声压级分布的影响。但如果是宽带信号比如脉冲声源、爆炸波、超声换能器发出的短脉冲频域方法要么逐个频率点算效率很痛苦要么用快速算法做逆变换数学复杂度又上来了。时域差分方法不一样。你只需要把激励源在时域上写成函数施加到网格上后面的演化就真实地按物理过程走了。脉冲源天然覆盖宽频带一次模拟就能同时提取多个频点的响应。更重要的是时域方法处理非线性效应——比如大振幅声波下声速不再是常数——只需要在每个时间步上修改递推公式里的参数即可不需要重新推导整套频域理论。回到标题里“时域差分方程”这几个字本质上就是把二阶波动方程(p_tt c^2 p_xx)改成时间上的二阶中心差分和空间上的二阶中心差分组合成一个显式递推式。其中 p 是声压c 是声速t 下标表示对时间求二阶导x 下标表示对空间求二阶导。这个显式格式有一个著名限制就是库朗数 r c·dt/dx 必须不超过1否则时间步长太大波的因果性就会被破坏数值解直接爆掉。具体原因后面我详细说这是整个方法里最需要心领神会的一点。1.3 边界截断与有限计算域实际做仿真时你不可能一直把所有空间的网格都铺满计算域总是有限的。问题在于真实声波在自由空间里应该无反射地一直往外传播但有限计算域的边界如果处理不当波到了边缘就会像撞到墙一样反弹回来干扰你真正关心的内部观测区域。所以在设计阶段就要决定边界方案。最简单的办法是让网格足够大在关心的区域之外留出足够厚的“牺牲层”在数值模拟停止之前边界反射波还来不及传回观测区。这个办法实现成本最低但网格多、算得慢只适合验证程序逻辑。实际项目里更常用的是吸收边界条件其中对入门者最友好的是Mur一阶吸收边界只需要在计算域边界点额外套一条单向波近似公式实现代码也就几行。更高吸收率的是完美匹配层PML吸收效果确实好但实现复杂、内存开销大。我第一次做波动仿真时也是一步一步从零反射边界切到PML的所以这篇文章后面会有专门的篇幅展开实验技巧。1.4 Matlab给这个项目带来的效率与调试便利用Matlab写声波时域差分还有一个隐性收益调试门槛特别低。你可以设置断点看任意网格点在任意时刻的数值也可以用多条figure窗口同时看波场快照和时域波形还能用surf(x, t, p)把整个时空演变画成瀑布图。这种可视化的直观程度对理解波的传播机理帮助巨大。但Matlab也有它的软肋。如果代码里到处是for循环遍历所有网格点内部运行时会被慢速的循环解释机制拖慢。正确做法是把空间上的递推操作整段向量化用一次数组切片计算替代循环效果往往快上十倍百倍。这部分优化经验在接下来的实操章节里会具体展开。2. 核心细节解析与实操要点2.1 数值网格设计物理问题的几何尺寸决定计算域计算域里再把空间离散成网格点每个网格点存一个声压值。以我常用的一个一维模型为例假设模拟一段长 10 m 的空气管道声速 c340 m/s网格间距 dx0.01 m那么网格点数就是1000。时间轴不是随意取的要由库朗数反过来决定(dt r·dx / c)如果取 r0.5dx0.01 m那么 dt≈1.47×10⁻⁵ s。这个 dt 意味着每秒大约有68000个时间步。仿真100 ms就是6800步Matlab里跑起来是很快的。有些初学者误以为 dx 越小越精确就把 dx 设得特别小结果网格点数太多内存爆掉或者算半天没反应。要注意dx 的最小值取决于你关心的最高频率。一个经验法则是一个波长内至少要保证10个网格点否则空间色散误差会明显到让波形产生虚假振荡。比如关心的最高频率是 10 kHz对应的波长大约 0.034 m那么 dx 取 0.001 m 也能覆盖但再加细到 0.0001 m 就属于浪费计算资源了。2.2 CFL条件的物理直觉CFL条件全称Courant–Friedrichs–Lewy条件当时是三位数学家在偏微分方程数值解的论文里提出的。道理简单得像赶公交车每推进一步的时间 dt 里波只能走 c·dt 那么远而你的差分模板在空间上只“看得见”相邻网格点信息从一格传到下一格至少要需要一个时间步 dt。如果 c·dt dx那就意味着在这一个时间步里波已经跨过了几格距离但当前差分模板还没来得及把它传过去数值解就会出现不稳定也就是通常所说的“爆了”。用数学式写就是(r c·dt / dx ≤ 1)实际操作里我不建议把 r 设成1.0虽然理论上稳定但离散色散误差会让波形略微滞后于真实波前。更稳妥的是取0.4到0.9之间。我常用的值是0.5存储空间基本够用精度和稳定性相对平衡。下面这个表总结了不同库朗数下的一维声波模拟表现是我在不同实验里观察到的规律可以给新手一个大致参考r 值稳定性风险数值色散倾向时间步推进效率适用场景0.1极低纵向扩散明显慢网格验证、教学演示0.5低轻微色散中等通用仿真最常用0.9存在一定风险靠近波前处可见拖尾快追求效率且熟悉调参1.0临界稳定易受边界干扰边界处可能出现假反射最快纯理想情况2.3 时间递推格式一维声波方程的时间二阶导用中心差分近似后得到的递推格式是p(n1,i) 2 p(n,i) − p(n−1,i) r² [p(n,i1) − 2 p(n,i) p(n,i−1)]其中 n 标记时间步i 标记空间网格点。这个格式是三层的要算新时刻的值需要上一时刻和上上时刻的场值。所以在主循环开始前需要准备两个历史数组初始时刻通常是零场再加上初始扰动或源项。我之所以用中心差分而不是前向或后向差分是因为中心差分对二阶导的误差是二阶小量也就是说你每把网格细化一半误差会缩小到原来的四分之一。如果为了省事用前向差分误差会变成一阶小量同样细化一半网格误差只缩到一半收敛太慢了波形形状也会明显畸变。还有一点值得提醒时间步进里p(n−1,i) 的历史值一定不要覆盖掉。我见过不少初学者在循环里把 p_old p; p p_new; 写成 p p_new; p_old p; 结果新旧时刻的值互相污染波场很快变成一团噪声。2.4 激励源与波形声波模拟需要一个声源最常用的是高斯脉冲、正弦脉冲或者由用户自定义波形。高斯脉冲的好处是频谱平滑、频带可控激发出来的波不会带过多高频抖振。如果硬生生给一个阶跃信号频谱里无限宽带数值上会引起吉布斯振铃边界越陡峭振铃越明显。施加声源的方式有硬源和软源两种。硬源直接强行把某个网格点的值覆盖为给定波形实现简单但会在这个网格点产生人工散射影响附近场的形态。软源则是在递推公式计算完新值之后再把源项叠加到对应点上相当于给介质加了一个外力物理上更自然。我在一维验证中经常用初始条件高斯脉冲法也就是在 t0 时刻直接把场的初始分布设成一个高斯波包整个区域里没有任何边界源这样最有利于先检查递推格式本身的正确性排除源项实现引入的Bug。2.5 边界条件对比边界类型递推方程物理含义实现难度刚性壁面p0 或 ∂p/∂n0完全反射极低Mur一阶单向波近似弱吸收低PML坐标拉伸/衰减层强吸收高一阶Mur边界条件的实现公式是p_new(1) p(2) (r−1)/(r1) * (p_new(2) − p(1))这个公式的思路是把边界上的波动当作沿法线方向外传的单向波来处理公式里同时用了当前时刻相邻点的值和上一步的值形成一种“朝向边界”的传播关系。实际效果是垂直入射的波吸收得不错但斜入射波会有部分反射所以吸收效果有方向依赖。如果只需要短时间模拟且观测区域离边界足够远Mur边界性价比很高。如果是二维或三维仿真斜入射波的比例大幅上升Mur边界效果下降这种时候就值得上PML。PML的思路是在最外层设置一个衰减介质层让进入里面的波逐步衰减而不是反弹回来。PML参数调得好时反射系数可以做到 −60 dB 以下但要注意媒质参数在层内连续渐变否则界面处还会产生数值反射。3. 实操过程与核心环节实现3.1 一维声波FDTD完整Matlab代码下面给出一套可运行的一维声波时域差分程序。这个版本包含了初始高斯脉冲源、Mur吸收边界、以及简单的可视化。我用的是Matlab R2020b版本语法比较保守绝大多数版本都能直接跑。% 一维声波时域有限差分FDTD仿真 % 物理参数 c 340; % 声速m/s dx 0.01; % 空间步长m L 10; % 计算域长度m Nx round(L / dx) 1; % 网格点数目 x (0:Nx-1) * dx; % 离散参数 r 0.5; % 库朗数 dt r * dx / c; % 时间步长 T 0.03; % 仿真时长s Nt round(T / dt); % 时间步数 % 初始高斯波包参数 x0 3.0; % 波包中心位置 sigma 0.08; % 波包空间宽度 % 场初始化 p_now zeros(1, Nx); p_prev zeros(1, Nx); % 设置初始条件 p_now exp(-((x - x0).^2) / (2 * sigma^2)); % 吸收边界系数Mur 一阶 coefMur (r - 1) / (r 1); % 绘图准备 figure(Position, [100 100 800 300]); % 时间递推 for n 1:Nt % 内部点差分递推 p_next zeros(1, Nx); p_next(2:Nx-1) 2 * p_now(2:Nx-1) - p_prev(2:Nx-1) ... r^2 * (p_now(3:Nx) - 2 * p_now(2:Nx-1) p_now(1:Nx-2)); % 边界Mur 一阶吸收 p_next(1) p_now(2) coefMur * (p_next(2) - p_now(1)); p_next(Nx) p_now(Nx-1) coefMur * (p_next(Nx-1) - p_now(Nx)); % 更新历史场 p_prev p_now; p_now p_next; % 每隔若干步绘图 if mod(n, 200) 0 plot(x, p_now, b-, LineWidth, 1.2); xlim([0 L]); ylim([-1.2 1.2]); grid on; title(sprintf(t %.5f s, n * dt)); xlabel(x (m)); ylabel(声压); drawnow; end end这段代码的第一步是初始化场第二步进入主循环第三步是边界处理第四步是可视化。整个流程不到四十行却把一维声波传播的数值模拟完整跑通非常适合作为后续加PML、加源项、加非均匀介质的模板。3.2 关键参数选择与验证思路上面的程序里有几个参数需要特别关注。首先是 dx 的选择0.01 m 对应最高频率大致在 3.4 kHz 附近因为一个波长对应34个网格点精度已经足够。如果关心50 kHz的超声波dx 至少得压到0.0005 m网格数和计算量都会大幅上升这就是为什么高频超声模拟的时间成本如此之高。库朗数 r 取0.5用上面公式算出 dt。很多新手会犯的错误是直接从任意给定的 dt 出发而不根据 dx 和 c 反推结果算出来库朗数大于1程序一秒钟就发散。判断发散有一个简单特征场值出现数千倍于声源量级的高频振荡那基本就是CFL超限了。验证递推精度的第一步是拿程序算一个最简单场景比如去掉边界条件下让脉冲在管内来回反射记录某点的波形和解析镜像法计算的结果对比。如果两者对准时间差在几个时间步以内核心递推格式基本正确。第二步再打开Mur边界确认波到达边界之后反射波幅度明显下降。3.3 可视化方法一维声波仿真的可视化最简单直观的就是波场-时间瀑布图。你可以把每个时间步的 p_now 存下来最后用:imagesc(x, (0:Nt-1)*dt, p_history); xlabel(位置 (m)); ylabel(时间 (s)); colorbar;来看波动的完整传播路径。这种图特别适合观察反射波、透射波和绕射波的时距关系。比如入射波到达某个界面后你会看到一条弯曲的散射双曲线从界面点向外扩散这就是物理上的散射信号。Matlab的动画绘制虽然方便但每次都 plot 会拖慢速度所以我一般是每50个时间步才更新一次绘图。如果嫌手写循环图太占空间也可以用 VideoWriter 录制导出成视频方便后续整理到PPT里。3.4 二维扩展方向一维程序跑通之后二维扩展是自然的下一步。二维声波方程形式如下(p_tt c² (p_xx p_yy))差分格式变成五点的空间模板p(n1,i,j) 2 p(n,i,j) − p(n−1,i,j) r² [p(n,i1,j) p(n,i−1,j) p(n,i,j1) p(n,i,j−1) − 4 p(n,i,j)]这里库朗数条件变为 r c·dt/dx ≤ 1/√2原因是在两个方向上同时传播时波走的实际距离更长。我当年做二维模拟时在这个问题上栽过跟头先是按一维习惯取r1结果二维场的角点附近总是爆后来查资料才意识到条件变了。二维的Mur边界也有对应的公式但代码复杂度明显上升。如果只做验证实验更推荐直接在计算区域外侧加一圈吸收层简易PML通过连续加几次衰减系数降低反射。4. 常见问题与排查技巧实录4.1 程序发散先查库朗数发散是时域差分程序最常见的故障。观察到的现象是运行几步后场值出现断崖式增大图形上像打碎的玻璃没有规律可言。多数情况下就是库朗数大于1了。排查步骤很简单在代码开头加一行r c * dt / dx; if r 1 error(CFL 条件不满足); end这行代码能挡住大多数低级错误。还有一种隐蔽的情况是你在程序某处改了 dx但 dt 还按旧值算导致库朗数悄悄越限。建议把所有网格参数集中放在脚本顶部所有后续参数都从它们推导不要散落在各个段落里改。4.2 边界反射干扰观测结果用Mur边界时如果反射波仍然明显先确认你的波是否接近垂直入射。如果入射角大于45度Mur一阶边界的吸收能力就明显下降。两个补救办法一是把计算域边界到观测区的距离加大让反射波在仿真时间内来不及回到观测区二是直接升级到PML。对刚入门的人来说先加大距离往往最省事因为PML本身就有一大堆参数需要调。怎么判断看到的反射是真实的还是数值反射呢一个办法是把边界改回刚性反射边界跑一遍对比两者差异。如果Mur边界解和刚性边界解在某些时刻差异不大说明吸收边界没有起到作用如果差异明显但Mur边界解里仍然有小幅度反射那是吸收不完全的数值反射。4.3 波形出现拖尾或畸变当你看到脉冲传播一段时间后波前后面跟着一串轻微振荡这就是数值色散。它来自差分格式对物理色散关系的近似误差。高频分量在离散网格中传播速度略微慢于低频分量所以本来一个光整的脉冲会逐渐展宽并拖尾。如果我关心的是波前到时和振幅峰值少量色散是可以接受的。如果色散严重到无法辨认波形该考虑减小dx和dt维持库朗数不变提高空间分辨率或者用更高阶的空间差分格式。以我经验dx取最高频率波长的1/20到1/30色散基本可控。如果必须用粗网格还要减少色散四阶空间差分比二阶格式的色散误差低很多但递推公式会写成五个网格点参与而不是三个编码量和边界处理难度也随之上升。4.4 计算效率太低Matlab里for循环一层套一层是大忌。例如三维声波模拟如果写三层for嵌套去遍历每个网格点更新公式一个时间步可能要几十毫秒一万步就是几百秒很容易让人崩溃。正确做法是所有内场更新尽量用冒号切片表达p_next(2:Nx-1) ... p_now(2:Nx-1) ...这种写法在Matlab内部调用的是编译后的循环效率远高于解释执行的for循环。二维也可以完全用矩阵切片p_next(2:Nx-1,2:Ny-1) ...如果向量化之后仍然慢可能是网格规模太大了这时候要重新审视dx是否选得过于保守。曾经有个朋友把空气中声波模拟的dx设成0.0001 m导致网格数高达十万个点一个时间步的计算格外吃力。事实上对于可听声0.005 m的网格已经能满足精度他这个设置纯粹是做无用功。4.5 常见问题速查表现象可能原因解决方式数值爆掉库朗数超过临界值减小dt或增大dx保证r≤1后期出现疑问反射边界吸收不完全加缓冲层或改用PML波形拖尾空间网格过粗降低dx并同步调整dt波速滞后差分格式色散误差增大r不超过限制并细化网格计算慢循环遍历或网格过细向量化内场递推检查dx必要性初值边界互相干扰历史数组更新顺序错误p_new计算完成后再依次更新p_prev和p_now4.6 从一维实验到实际项目的一些体会我当初把这段程序跑通之后最大的收获不是代码本身而是对“数值解必须和物理机制对得上”这件事有了直觉。比如看到波碰到界面后反射波和透射波幅度关系会自然想起声阻抗匹配的概念。如果程序中加入了不同密度的介质层你会看到反射系数随着两侧声阻抗差变化这比单纯看公式要生动得多。再后来我把这个程序做成一个用来验证超声检测回波信号的辅助工具只需要在代码里改换换能器位置和介质层参数就能快速估算出回波信号的大致时延和幅度帮实验同事节省了大量试错时间。这种“用数值实验替代盲目实测”的工作方式在我看来就是这个题目最大的实用价值。对新手来说不用急着上手复杂界面第一步把一维代码跑通第二步加不同介质层看反射和透射第三步扩展到二维看散射现象每一步都是稳固的积累。5. 后续扩展建议从一维声波到二维甚至三维最核心的框架不变只是递推模板变大、边界更复杂。我建议你下一步做两个方向。第一是加多层介质板在网格里放一块和空气声速、密度不同的材料区域观察波在平面界面的反射和透射是否和理论公式相符。这个实验能帮助你理解阻抗失配对反射强弱的影响。第二是在网格里设置一个障碍物比如圆形或多边形观察波绕射形成的阴影区这在无损检测领域是常见研究内容。如果对PML感兴趣可以在两个方向上各加10层的衰减吸收层衰减系数沿层内渐变比如:sigma_max 0.8 * (m 1) / dx; % 粗略的经验公式实际调参过程中你会发现自己逐渐掌握一套“先稳定、再吸收、后精度”的调试顺序。这跟做实验、调设备是同一个逻辑没有捷径但多跑几次心里就踏实了。我个人做下来的体会是声波时域差分程序是一种“看起来简单、做起来充满细节”的典型数值仿真工具。代码量不大但每个参数都有它存在的理由每一处边界处理都可能让你学到新的数值分析方法。把它当敲门砖后面学弹性波、电磁波、地震波思路都是相通的。