企业官网建设流程全解析

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 的空气管道，声速 c=340 m/s，网格间距 dx=0.01 m，那么网格点数就是1000。时间轴不是随意取的，要由库朗数反过来决定：

(dt = r·dx / c)

如果取 r=0.5，dx=0.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，存储空间基本够用，精度和稳定性相对平衡。

下面这个表总结了不同库朗数下的一维声波模拟表现，是我在不同实验里观察到的规律，可以给新手一个大致参考：

2.3 时间递推格式

一维声波方程的时间二阶导用中心差分近似后，得到的递推格式是：

p(n+1,i) = 2 p(n,i) − p(n−1,i) + r² [p(n,i+1) − 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 激励源与波形

声波模拟需要一个声源，最常用的是高斯脉冲、正弦脉冲或者由用户自定义波形。高斯脉冲的好处是频谱平滑、频带可控，激发出来的波不会带过多高频抖振。如果硬生生给一个阶跃信号，频谱里无限宽带，数值上会引起吉布斯振铃，边界越陡峭，振铃越明显。

施加声源的方式有硬源和软源两种。硬源直接强行把某个网格点的值覆盖为给定波形，实现简单但会在这个网格点产生人工散射，影响附近场的形态。软源则是在递推公式计算完新值之后，再把源项叠加到对应点上，相当于给介质加了一个外力，物理上更自然。

我在一维验证中经常用初始条件高斯脉冲法，也就是在 t=0 时刻直接把场的初始分布设成一个高斯波包，整个区域里没有任何边界源，这样最有利于先检查递推格式本身的正确性，排除源项实现引入的Bug。

2.5 边界条件对比

一阶Mur边界条件的实现公式是：

p_new(1) = p(2) + (r−1)/(r+1) * (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(n+1,i,j) = 2 p(n,i,j) − p(n−1,i,j) + r² [p(n,i+1,j) + p(n,i−1,j) + p(n,i,j+1) + p(n,i,j−1) − 4 p(n,i,j)]

这里库朗数条件变为 r = c·dt/dx ≤ 1/√2，原因是在两个方向上同时传播时，波走的实际距离更长。我当年做二维模拟时在这个问题上栽过跟头，先是按一维习惯取r=1，结果二维场的角点附近总是爆，后来查资料才意识到条件变了。

二维的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 常见问题速查表

4.6 从一维实验到实际项目的一些体会

我当初把这段程序跑通之后，最大的收获不是代码本身，而是对“数值解必须和物理机制对得上”这件事有了直觉。比如看到波碰到界面后反射波和透射波幅度关系，会自然想起声阻抗匹配的概念。如果程序中加入了不同密度的介质层，你会看到反射系数随着两侧声阻抗差变化，这比单纯看公式要生动得多。

再后来我把这个程序做成一个用来验证超声检测回波信号的辅助工具，只需要在代码里改换换能器位置和介质层参数，就能快速估算出回波信号的大致时延和幅度，帮实验同事节省了大量试错时间。这种“用数值实验替代盲目实测”的工作方式，在我看来就是这个题目最大的实用价值。

对新手来说，不用急着上手复杂界面，第一步把一维代码跑通，第二步加不同介质层看反射和透射，第三步扩展到二维看散射现象，每一步都是稳固的积累。

5. 后续扩展建议

从一维声波到二维，甚至三维，最核心的框架不变，只是递推模板变大、边界更复杂。

我建议你下一步做两个方向。第一是加多层介质板，在网格里放一块和空气声速、密度不同的材料区域，观察波在平面界面的反射和透射是否和理论公式相符。这个实验能帮助你理解阻抗失配对反射强弱的影响。第二是在网格里设置一个障碍物（比如圆形或多边形），观察波绕射形成的阴影区，这在无损检测领域是常见研究内容。

如果对PML感兴趣，可以在两个方向上各加10层的衰减吸收层，衰减系数沿层内渐变，比如:

sigma_max = 0.8 * (m + 1) / dx; % 粗略的经验公式

实际调参过程中，你会发现自己逐渐掌握一套“先稳定、再吸收、后精度”的调试顺序。这跟做实验、调设备是同一个逻辑，没有捷径，但多跑几次心里就踏实了。

我个人做下来的体会是，声波时域差分程序是一种“看起来简单、做起来充满细节”的典型数值仿真工具。代码量不大，但每个参数都有它存在的理由，每一处边界处理都可能让你学到新的数值分析方法。把它当敲门砖，后面学弹性波、电磁波、地震波，思路都是相通的。