企业官网建设流程全解析

简介：超声声场FDTD计算程序是一款基于MATLAB的有限差分时域数值模拟工具，面向从事超声传播研究、HIFU高强度聚焦超声治疗规划、超声成像与无损检测的科研和工程人员。资源压缩包共10个文件，约331KB，其中MATLAB脚本（.m）为核心算法实现，.mat数据文件用于存放介质参数与仿真结果，.dll与.mexglx为编译后的加速接口文件，.c源码可供深入修改底层代码，便于在MATLAB环境中直接运行、调试和二次开发。程序通常包含网格初始化、时间步进循环、边界处理和后处理等模块，可计算声压、声强在复杂介质中的分布，预测HIFU治疗剂量，并通过调整频率、功率、声源形状及介质属性，适应不同实际场景。已有407人浏览学习，研究者可借此快速搭建超声声场仿真环境，理解FDTD算法的实现逻辑与优化方法，并结合示例数据开展进一步的对比实验与参数分析。 上个月我在调一套单阵元超声换能器的声场参数，瑞利积分和角谱法都跑了一遍，结果近场干涉结构始终和实测水听器数据对不上。后来把主程序换成超声声场FDTD计算程序，把透声层、边界反射和近场旁瓣一起纳进来，才算把问题解释通。这篇就把我从波动方程离散到写出一套能跑、能验证的FDTD程序的完整过程整理出来。重点不在把公式推导一遍，而在哪些参数会毁掉结果、边界怎么收、怎么确认你算出来的声场确实是对的。如果你正在做超声换能器设计、无损检测仿真或者声场时域分析，这篇可以省掉你至少一周的试错时间。

1. 为什么超声声场仿真我最后选了FDTD

1.1 瑞利积分和角谱法解决不了的问题

绝大多数人一开始做超声声场，第一反应都是瑞利积分。它确实快，把换能器表面法向振速分布积分一下，空间任意一点的声压就出来了，远场近似和轴线响应都挺准。但瑞利积分的前提是“无限大均匀介质”，换能器周围的边界、透镜、仿体层只要一出现，它的积分核就没有办法处理多次散射和层间折射。我做组织仿体实验时中间夹了两层不同厚度的PVA薄膜，瑞利积分对这个局部扰动完全无感，算出来和实测差了快30%。

角谱法本质上是平面波分解，比瑞利积分灵活一些，但也卡在“计算域两侧介质均匀”这个隐含假设上。多层介质不是不能做，每层界面都要重算透射和反射系数，斜入射时还会冒出数值泄漏。再加上它天然适合频域，想拿时域波形做脉冲回波分析，就得每帧都做一遍正向和逆变换，效率反而下去了。

1.2 有限元方法在时域问题上的代价

你可能问，COMSOL这类有限元软件不是也能算吗？能算，而且频域下非常好用，算一张单频谐波声场图很漂亮。但如果激励是短脉冲或宽带信号，它每个频率点都要重新组装一次刚度矩阵，时域显式求解又要求网格尺寸小到离谱，自由度数量比同等精度的FDTD大一到两个数量级。我试过一次三维时域压电换能器仿真，网格剖分完直接吃了30多G内存，机箱风扇转得跟起飞一样，最后算完一个case花了整晚。

FDTD的好处在于三点：第一，它直接给时域波形，脉冲传播过程一步一帧看得清清楚楚；第二，异质介质只是给每个网格点赋不同的声速和密度，改数组就行，界面处天然按法向连续性处理；第三，时间显式迭代，代码简洁，后期做参数扫描和反演很方便。

1.3 三种方法的选型对比

如果你要做的是换能器近场结构、多层介质中的传播、脉冲回波，或者后续想加非线性项，FDTD是综合成本最低的一条路。

2. 声学FDTD的核心离散逻辑：拆方程、摆网格、踩住两条硬约束

2.1 把二阶波动方程换成一阶速度-压力方程组

直接离散二阶声压波动方程不是不行，但遇到介质阻抗突变时，二阶中心差分会在界面处产生非物理振荡。更好做的做法是把方程降阶，写成速度-压力耦合的一阶方程组：

$$\rho \frac{\partial \mathbf{v}}{\partial t} = -\nabla p$$

$$\frac{\partial p}{\partial t} + K \nabla \cdot \mathbf{v} = 0$$

其中 $K = \rho c^2$ 是体积模量。这组方程对应的物理含义很直观：压力梯度驱动质点运动，质点运动散度反过来改变局部体积和压力。界面处只需要保证法向振速连续和压力连续，离散时非常自然。

2.2 交错网格与蛙跳更新

声学FDTD沿用电磁学的Yee网格思路。在二维笛卡尔网格里，压力p定义在网格中心(i,j)，水平振速u定义在水平半格点(i+1/2,j)，垂直振速w定义在垂直半格点(i,j+1/2)。时间上也交错，p在nΔt时刻，速度在(n+1/2)Δt时刻。这样空间和时间都做到中心差分，保证二阶精度。核心更新伪代码大致是：

# p 尺寸 Nx×Nz, u 尺寸 (Nx+1)×Nz, w 尺寸 Nx×(Nz+1) # K = rho * c^2，网格步长 dh，时间步长 dt # 第一步：由速度散度更新压力 p[:, :] -= dt * K * ( (u[1:, :] - u[:-1, :]) / dh + (w[:, 1:] - w[:, :-1]) / dh ) # 第二步：由压力梯度更新速度 u[1:-1, :] -= dt / rho * (p[1:, :] - p[:-1, :]) / dh w[:, 1:-1] -= dt / rho * (p[:, 1:] - p[:, :-1]) / dh

实际代码里需要小心数组边界不同长的问题，我上面刻意保留了交错偏移，就是为了强调相邻半格点不要少算一格。很多新手的第一个bug就出在这里。

2.3 CFL条件与网格分辨率两条硬约束

两条硬约束决定了你的计算不会发散的底线和结果可信度的上限。

第一是CFL稳定性条件。二维声学FDTD要求 $c_{max}\Delta t \le \frac{\Delta h}{\sqrt{2}}$，其中 $c_{max}$ 是计算域内最大声速。原因很直观：一个时间步内，波最多跨越一个网格对角线长度。我一般取理论上限的0.9倍左右，但实际操作中发现0.5到0.8更稳，尤其是加了PML吸收层之后。

第二是网格分辨率PPW，即每个波长内至少多少个网格点。经验值是至少10个，建议15到20个。低于8个时数值色散会非常明显，声束会变糊，焦点位置会漂移。举个例子：水中1MHz超声，波长1.5mm，取75μm的网格就是每波长20个点；取CFL=0.9，时间步长约为31.8ns；模拟50μs的传播过程大约需要1570步。在这个尺度下二维算起来非常轻松，几秒到几十秒就能跑完。

3. 写程序时最容易被忽视的五个部件

3.1 计算域与介质参数数组

不要一开始就想把所有几何细节都做进去。先画清楚计算域里有什么：换能器在哪个位置，介质层在哪，厚度多少，声速密度是多少。每个网格点存两个基本参数——密度ρ和声速c，然后实时算出K=ρc²。非均匀介质就是把不同区域的网格值填成不同数组，界面处可以用阶梯近似，也可以取体积加权平均，后者在斜入射时数值散射更小。

有一个细节要特别注意：换能器压电层本身的阻抗比水高很多，如果直接用压力源激励，界面处会看到明显的反射波。实际换能器表面是用背衬层吸声的，仿真里要记得在换能器背面加吸收层，否则这部分反射会混进声场里。

3.2 激励源设定

超声仿真里最常用的激励是高斯包络正弦脉冲，也叫tone burst：

$$s(t) = A \exp\left(-\frac{(t-t_0)^2}{\sigma^2}\right) \sin(2\pi f_0 t)$$

$t_0$ 一般取3到5个周期，$\sigma$ 决定带宽。要注意的是$t_0$不能太小，否则起始时刻脉冲已经被截断出一段高频毛刺；$\sigma$也不能太大，不然频带太窄，所谓“脉冲”会退化成连续波。我习惯让脉冲包含4到8个完整周期。

为什么不直接加连续波？因为连续波会让计算域里的PML反射和多次混响持续叠加，你根本分不清看到的干涉条纹是真实的稳态场还是边界反射造成的驻波。脉冲激励可以按到达时窗把直达波和后续反射波分开，对验证程序太方便了。

3.3 边界处理心得

这是整个程序里最容易翻车的地方。最简单的是把外边界设成硬边界（振速为0），但这样算出的声场会被刚性边界反射污染。稍微好一点用一阶Mur吸收边界，垂直入射效果不错，斜入射就会有一小部分反射回来。我最终用的是PML吸收层，严格成品做法是分裂场PML。核心思想是在计算域外圈加一层有损耗的介质，让波在层内指数衰减，再通过坐标伸缩让阻抗匹配，消除界面反射。

一个能直接参考的经验参数：

• PML层厚L取10到15个网格点；
• 层内吸收系数按幂分布 $\sigma(x) = \sigma_{max} \cdot (x/L)^m$，$m$取2到3；
• 最大吸收系数 $\sigma_{max} = -\frac{(m+1) c \ln(R)}{2L}$，反射系数R取$10^{-3}$到$10^{-5}$之间。

这个公式不是拍脑袋写的，它是从坐标伸缩PML的反射系数推导结果整理出来的。照着设置，实测残余反射比Mur边界低一个数量级以上。

3.4 探针与快照

超声声场FDTD的价值不只是最后出一张彩色云图，更重要的是可以任意位置放探针点，记录随时间变化的压力波形，模拟水听器测量。探针数据可以存在内存里，最后统一写盘。全场的压力快照则每隔若干时间步存一次，文件格式建议用npy、zarr或HDF5，别用文本CSV堆几百个G。

另外强烈建议：先别急着上三维。二维程序跑通、验证正确之后，改成三维的代码工作量其实不大，就是把数组从二维扩成三维，但内存和计算时间会涨几十倍。在二维下把所有物理验证做完，再上三维，节约的调试时间非常可观。

4. 调试记录：四个把结果毁掉的坑和排查链路

4.1 PML层内参数填错，边界“闷响”不吸收

症状：声场快照里，波不但没从边界消失，还在边界附近形成一圈明亮的振荡带，持续很长时间不衰减。

排查链路：我先怀疑PML层不够厚，加到30层，没用；又怀疑吸收系数太小，把σ_max上调两个数量级，结果反而更糟，界面处出现更强的反射。最后检查PML区域的ρ和c数组才发现，PML层里的介质参数填成了空气值，阻抗完全不匹配，波在PML入口直接被反射了。

修复：让PML层内每个网格的ρ和c与相邻介质完全一致，只改σ吸收系数。改完再跑，边界反射残影基本消失。这个坑的原因很简单：PML只有在阻抗匹配时才“透明”，任何参数不一致都会让PML本身变成一面镜子。

4.2 网格分辨率不足，焦点位置漂移

症状：算出来的轴线声压最后一个极大值位置比理论值偏前了好几个毫米，旁瓣结构也不对称。

排查链路：第一反应是源设置有问题，换了不同激励时长，焦点位置不变；又怀疑边界反射，把计算域扩大了一倍，问题依旧。最后想到可能是网格太粗导致数值色散。一检查，我的PPW只有6个网格/波长，远低于经验值15。粗网格下相速度偏慢，波前传播距离被低估，焦点自然往前偏。

修复：把网格从波长/6加密到波长/20，焦点距离与理论值的偏差从约10%降到1%以内。此后我给自己定了一个规矩：任何模型先检查PPW，低于12的网格结果不看。

4.3 时间步长取太满，数值发散成NaN

症状：程序跑到几百步后，压力场出现高频振荡，再跑一二十步直接变NaN。

排查链路：我当时CFL取了理论极限的0.99，自认为卡得很准。结果发现PML层内由于吸收系数梯度的存在，局部等效波速并不完全等于介质声速，导致局部CFL刚好突破稳定上限。这其实是个很典型的“局部稳定性”问题。

修复：把全局时间步长从CFL=0.99降到0.6，其余代码没动，程序稳定跑到10万步。之后我建议把CFL设置为0.5到0.8之间，为了让计算加快而卡上限，最后往往花更多时间找NaN。

4.4 一帧都没算完，硬盘先爆掉

症状：跑了一次三维程序，设置每步都存全场快照，结果算了三分之一的步数后磁盘空间告警。

排查链路：三维网格240万点，每帧float32就要约10MB，一万步就是100GB。我当时完全没做步长抽取和浮点压缩。

修复：一是隔N步存一帧，N按输出需求来选；二是把全场快照降成float16或单精度numpy数组；三是只针对观测面保存，比如x-y剖面或者特定深度的二维切片。改完之后，同样一个仿真数据量下降了两个数量级。这是工程问题，不是物理问题，但处理不好照样会中断整个任务。

5. 验证声场：怎么确认算出来的图不是噪声

5.1 先做一维平面波验证

任何FDTD程序写完，第一件事不是去看二维云图，而是先简化为纯一维网格，从源发出一个平面波，观察压力和速度沿传播方向的波形。对比理论解，距离每增加一个波长，相位延迟应该是2π，幅值保持不变。如果一维都出现幅值衰减或相位偏移，说明离散格式或者介质参数数组有问题。这个测试只需要几十行代码、跑几百步，比在二维里猜bug快得多。

5.2 轴线声压曲线与解析解对比

二维验证通过后，再算一个圆盘活塞换能器的轴线声压分布。均匀激励圆盘在近场会出现一系列极大极小值，最后一个极大值大约落在近场长度 $N = a^2/\lambda$ 附近，远场区域则按1/z衰减。

我常用的验证案例：圆盘半径a=5mm，频率1MHz，水中声速1500m/s，波长1.5mm，近场长度约16.7mm。FDTD算出最后一个轴上极大值在17mm左右、偏差1mm以内，-6dB声束宽度也符合菲涅尔近似，我就认为这个程序的物理正确性过关了。

5.3 看声场快照的物理直觉

数值没问题，还要看物理上的合理性。近场区域应该看到明显的干涉条纹，换能器边缘会产生衍射旁瓣，远场主瓣逐渐变宽并平滑。我建议画图时用声压包络或dB图，而不是直接画瞬时压力。瞬时压力场是正负交替的，容易让人误以为声场断裂，包络图才能清晰展现声束走向和旁瓣结构。

另外，换能器正前方轴线幅值应该随距离有起有伏，但整体趋势不会出现突然跳变；任何明显的“断崖”或“突刺”，都值得回头检查网格和边界。

这个程序现在已经成为我的一套“虚拟声场实验台”。后续我又往里面加了随频率变化的衰减模型和层状组织参数，很多实验前预判的问题都能先在仿真里过一遍。如果你正准备从零写超声声场FDTD，我的建议是：先别急着上三维，把一维、二维的平面波和圆盘验证跑通，物理正确性确认以后再扩展。我第一次把边界从反射改成PML、看到波“穿墙而过”毫无反射时，才觉得这个程序真正能当仿真工具用，而不是一个会画彩图的玩具。祝你好运，搞出能放心交付的声场结果。

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