简介超声声场有限差分时域FDTD计算程序是一套基于MATLAB的仿真工具面向超声成像、无损检测与高强度聚焦超声治疗领域的科研人员和工程技术人员用于求解超声波在非均匀介质中的传播问题获得声压、声强及声场分布为换能器设计和治疗剂量预测提供参考。压缩包共10个文件其中4个M脚本文件为核心计算与流程控制涵盖初始化、时间步进、边界处理、声场显示和数据保存等模块3个MAT数据文件存放介质参数及计算结果另有DLL、MEXGLX和C源码用于底层加速与跨平台兼容整体大小331KB便于下载与二次修改。已有407人浏览学习程序注释友好、模块清晰使用者可修改网格尺寸、声源频率、吸收边界及组织参数复现HIFU焦点区的高强度声压分布也可扩展至超声成像、无损检测或多阵元相控阵声场研究适合科研与教学深入参考。 超声声场FDTD计算程序说白了就是拿时域有限差分法求解波动方程把超声波在空间里怎么传播、怎么反射、怎么聚焦一步一步在电脑里算出来。我在做超声换能器设计方案验证和工业无损检测可行性评估时自己用Python写过一版二维的整体代码量不大跑起来也挺快但真正把稳定性、边界条件、声源加载这些细节处理好还是有几个绕不过去的坎。这篇就把我从零搭这个程序的完整思路、关键参数、核心代码以及踩坑记录整理出来给准备用数值手段做声场分析的朋友做个参考。1. 为什么自己写FDTD而不是用现成软件1.1 这个方法适合解决什么问题超声声场的数值仿真最常碰到的需求有三类一是换能器设计阶段想看某种阵元排列下声束的形状、焦点的位置、旁瓣的高低二是无损检测方案评估想预判一个已知缺陷在超声回波里长什么样三是医学超声里评估聚焦超声的治疗区域或者超声成像里声场在组织中的分布。这三类问题共同点在于声源几何往往是自定义的介质可能是多层结构、甚至声速空间分布不均匀而且我们往往不只关心稳态声场还关心声波从发射到碰到障碍物再到被接收的完整瞬态过程。FDTD这类时域方法天然贴合这个需求因为它每一步算的都是真实时间上的波场快照。相比之下商业软件如COMSOL、PZFlex功能很全但在做参数扫掠、把声场计算函数嵌入到自动优化循环时操作成本高、授权限制也多。自己写一版灵活度完全在自己手里。1.2 和其他声场计算方法放在一起比很多刚接触声场仿真的人会问为什么不用更简单的解析公式或者用有限元、边界元我把几种常见方案放在一张表里差别就很清楚了。方法物理假设优势局限解析法/瑞利积分无限大均匀介质、刚性障板速度快、适合远场纯预测难以处理反射、折射、非均匀介质射线法高频近似、声线追踪计算量小、适合大尺度聚焦区和衍射效应算不准确有限元 FEM频域或时域离散复杂几何建模能力强大场域网格量大、瞬态模拟成本高边界元 BEM均匀介质、边界积分降低一维维度处理非均匀介质困难、矩阵稠密FDTD时域全波求解瞬态过程直观、支持非均匀介质需要满足网格与稳定条件、内存随场域增长FDTD在超声这个领域最大的优势是它对介质的描述非常“朴素”。在每个网格点上单独设置声速、密度界面上的反射和透射自动就被波动方程计算出来了不需要像有限元那样显式处理边界条件。再加上显式时间步进天生适合并行加速所以它一直是我做超声声场快速验证的首选。1.3 一个明确的程序定位我给自己写的程序定了一个明确目标二维声场、均匀或分层介质、PML吸收边界、支持点源和阵列源激励能输出每一个时间步的声压场快照并且可以导出指定位置的时域波形。不上非线性不做热粘性损耗先保证声传播的物理过程正确。事实也证明把它跑通并验证稳定之后再往里头加障碍物散射、加声速梯度、加换能器聚焦延时都是很自然的扩展。2. 核心算法与关键参数2.1 控制方程与显式差分格式二维线性声波方程用声压可以写成[ \frac{\partial^2 p}{\partial t^2} c^2 \left( \frac{\partial^2 p}{\partial x^2} \frac{\partial^2 p}{\partial y^2} \right) ]其中c是介质声速。FDTD的做法是把这个偏微分方程在时间和空间上做中心差分。时间上保留二阶精度空间上对x和y方向分别做二阶中心差分整理后得到显式递推格式[ p^{n1}{i,j} 2p^n{i,j} - p^{n-1}{i,j} R^2 \left( p^n{i1,j} p^n_{i-1,j} p^n_{i,j1} p^n_{i,j-1} - 4p^n_{i,j} \right) ]其中 ( R c\Delta t / \Delta x )就是常说的Courant数。这个格式写起来极简每步只用到上一时刻和上上时刻的相邻网格值所以计算效率高、内存占用低。[!NOTE] 初学者最容易踩的坑是想当然地把时间步长设得很小以为越小越稳定。实际上这个格式的稳定性受Courant数约束和时间步长的绝对值无关关键是 ( R\frac{c\Delta t}{\Delta x} ) 不能超过一个上限。这也是后面数值发散最常见的根因。如果介质中有声速突变界面只需要在每个网格点单独赋不同的c值差分格式本身不用改反射和透射就自动被算出来了。这一点非常适合做分层介质中的声传播模拟。2.2 网格尺寸和时间步长怎么定网格尺寸的选择基本决定了一个FDTD程序能不能算出可信的结果。核心原则是每个波长内至少要放10到20个网格点。超声工程里一般取最关心频率对应的波长按 ( \Delta x \le \lambda_{min} / 15 ) 来设。我实际使用下来15个点是个不错的平衡点少于10个点会出现明显的数值色散算出来的波前面会“发毛”聚焦焦点也会变钝。举例来说水中的声速约1500 m/s如果中心频率是1 MHz那么波长 ( \lambda 1.5\ mm )网格取 ( \Delta x 0.1\ mm ) 就够精细。但如果频率提到5 MHz波长缩到0.3 mm网格就要取0.02 mm。同样是50 mm×50 mm的计算域1 MHz只需要500×500个网格5 MHz则要2500×2500个网格计算量差了25倍。这个数量级概念在项目初期一定要有不然很容易把算例规模设计得过大跑一次等半天还不知道是哪里卡住。时间步长则直接由网格尺寸和Courant数决定。二维情况的理论稳定上限是 ( R \le 1/\sqrt{2} )实际我习惯取 ( R0.4\sim0.5 )留出充足的稳定裕度。对应的 ( \Delta t R \cdot \Delta x / c )。我曾经试过卡着上限取 ( R0.7 )短期内没发散但算到几千步后能量开始异常增长最后查了很久才发现是稳定余量不足。频率声速波长推荐网格参考Δt1 MHz1500 m/s1.5 mm0.1 mm约3.3 ns2 MHz1500 m/s0.75 mm0.05 mm约1.7 ns5 MHz1500 m/s0.3 mm0.02 mm约0.7 ns2.3 声源激励与边界条件声源的加载是最影响结果物理意义的一个环节。最简单的点源直接在每个时间步往源网格点上叠加一个时域信号就行。我常用Ricker子波作为激励信号因为它没有直流分量频谱峰位于设计频率附近适合宽频超声激励。表达式如下[ s(t) \left(1 - 2\pi^2 f^2 (t - t_0)^2\right) \exp\left(-\pi^2 f^2 (t - t_0)^2\right) ]其中f是峰值频率( t_0 )一般取1.2到1.5个周期。用Python实现就是def ricker_wavelet(t, f_peak): t0 1.5 / f_peak arg np.pi * f_peak * (t - t0) return (1 - 2 * arg**2) * np.exp(-arg**2)如果是阵列换能器可以对每个阵元独立加载带延时的激励信号这样就能形成波束偏转或聚焦效果。平面波源也很简单在计算域一侧整行网格同时加权激励即可。边界处理是整个程序里最容易被低估的部分。如果计算域边界直接设为零声压很快就能看到强反射波从边界传回内部把正常声场完全淹没。我的做法是加PML吸收边界。PML的思想是在计算域外围构造一层阻抗匹配的损耗介质让入射波进去后被指数衰减吸收掉而不是反射回来。实现时通常会引入复坐标伸缩也就是把空间导数替换为带衰减因子的形式。实际工程里厚度取10到20个网格衰减系数从内向外按二次函数渐变最大值用 ( \sigma_{max} 3c / (2\Delta x) \cdot \ln(1/R) ) 来估算R取 ( 1\times10^{-3} ) 到 ( 1\times10^{-5} ) 的量级。按这个参数配出来的PML反射水平完全够用。[!TIP] 如果你用的是二阶声压方程PML实现会比较绕。我更建议在后续做PML时切换到一阶速度-压力方程在交错网格上更新质点振速和声压PML的引入会自然很多。这个转换只增加一个中间数组代码复杂度并不高。3. 从零到跑通完整实现与验证3.1 程序主循环的骨架代码整个程序的核心就是时间递推循环我用numpy向量化实现代码非常短。初始化阶段先分配三个二维数组分别保存下一时刻、当前时刻和上一时刻的声压场。import numpy as np nx, ny 500, 500 dx 0.1e-3 # 网格尺寸单位m c 1500.0 # 声速单位m/s dt 0.45 * dx / c # 时间步长Courant数取0.45 nt 2000 # 总时间步 # 声压场p_next, p_current, p_prev pp np.zeros((nx, ny)) pn np.zeros((nx, ny)) p_next np.zeros((nx, ny)) C2 (c * dt / dx) ** 2 # 激励源设置 src_x, src_y nx // 2, ny // 2 for n in range(nt): # 递推更新内部网格点 p_next[1:-1, 1:-1] ( 2 * pn[1:-1, 1:-1] - pp[1:-1, 1:-1] C2 * (pn[2:, 1:-1] pn[:-2, 1:-1] pn[1:-1, 2:] pn[1:-1, :-2] - 4 * pn[1:-1, 1:-1]) ) # 加载源信号 p_next[src_x, src_y] ricker_wavelet(n * dt, 1e6) # 边界处理PML或吸收层代码可替换为简单零边界 p_next[0, :] p_next[-1, :] 0 p_next[:, 0] p_next[:, -1] 0 # 推进时间层 pp, pn, p_next pn, p_next, pp if n % 100 0: # 每100步保存一帧云图 np.save(fframe_{n:04d}.npy, pn)这段代码里边界处用了简单的零边界启动时内侧区域不会立刻受边界影响所以也常用来做前几百步的快速自检。真正跑完整仿真务必换成PML。3.2 一个可复现的测试案例这里给一个适合新手的标准测试在一个大小为 ( 40 \lambda \times 40 \lambda ) 的均匀水介质计算域中从中心放置一个Ricker点源频率1 MHz网格按 ( \lambda/20 ) 划分固定边界时刻设为零。由于是各向同性介质中的点源理论上声压波前应当是一个完美的圆柱波传播速度应为水声速1500 m/s。我实际跑这个案例时用的是500×500网格Courant数取0.45。连续观察不同时间步的波场快照能看到一个圆环从源点均匀向外扩散且半径随时间线性增长。量一下不同方向的波前半径偏差在1%以内说明离散格式的数值各向异性被控制得很好。更贴近实际应用的是聚焦换能器案例在计算域的左侧布置一条120个阵元的线阵每个阵元按焦点位置计算延时激励信号使用中心频率1 MHz的Ricker子波。跑完2000个时间步后焦点区域会出现一个明显的声压增强点焦点位置和理论延时计算的位置相差不超过一个网格。同时焦点周围会自然出现旁瓣结构这正是后续评价换能器性能时最重要的输出信息。3.3 数值结果怎么验证数值程序最容易出现的问题是“看起来有图有波动但不知道对不对”。我的做法是至少做两项验证。第一项是解析衰减曲线对比。二维情况下无限介质中的点源辐射远场幅度近似按 ( 1/\sqrt{r} ) 衰减。在波前还没有到达边界的时间窗口内沿波前径向取声压峰值并绘制它相对于传播距离的曲线如果和 ( 1/\sqrt{r} ) 拟合曲线高度重合说明传播幅度是可信的。第二项是能量守恒检查。在没有PML吸收、且脉冲已经完全离开声源之后整个计算域内的声压平方和大致应该保持不变。如果总能量明显单调下降说明数值耗散过大通常和网格太粗或时间步长过大有关。如果能量异常增长则基本可以判断格式已处于不稳定状态。写程序时保留一个每100步计算一次总能量的钩子对早期调试非常有帮助。4. 调试与避坑记录4.1 数值发散的现象和真正原因最经典的现象是运行到某一步突然出现棋盘格状的交替正负大数值然后几步之内数值涨到天文数字。遇到这种情况优先检查Courant数。二维线性声波方程的显式格式稳定要求通常写作 ( c\Delta t / \Delta x \le 1/\sqrt{2} )我没有一次例外凡是突破这个上限的算例不管源信号多温和最终都发散了。建议初始调试阶段把Courant数设置在0.4以下先确保稳定再逐步调大看计算效率收益。另一个多发原因是声速赋值错误。如果某个网格点的声速因为代码逻辑问题被赋成极大值那局部的Courant数就会超限导致从该点向外发散。排查时可以在每个时间步检查声压数组的最大值如果发现某一个固定坐标点附近率先出现异常优先检查那里的材料参数。4.2 边界反射怎么判断和处理判断边界反射最直观的办法是让源信号发射完足够长时间后观察靠近计算域边缘的波场快照。如果看到圆弧形波前在边界处“弹”回内部说明边界吸收没有生效。最常见原因是PML厚度太薄或者衰减系数梯度太陡。我曾经为了省内存把PML设成6个网格结果反射波只被削弱了大概10 dB波前还是清晰可见。后来老老实实加到20层反射波幅度降到了背景噪声水平以下。PML内部还要注意一个问题介质声速不均匀时PML的衰减系数最好按当地声速来归一化否则不同方向上的吸收效果会出现不一致。这个细节在均匀介质验证时发现不了一旦换成多层介质模型就很容易露馅。4.3 网格分辨率不足的假象网格太粗不会让程序报错但结果会误导人。一个典型表现是圆形波前呈现出明显的方向性轴向和斜向传播速度不一致波前失去了光滑圆弧的形态。另一个表现是窄带脉冲经过一定距离传播后主波后面拖了一条振荡尾巴像梳子齿一样排列这就是数值色散导致的寄生振荡。遇到这类现象先把网格细化到每个波长至少20个点看现象是否消失。如果细化后结果变化很大说明之前的确是网格不足。早期做快速验证时可以容忍粗网格但正式的参数研究一定要在细化网格后重新确认结论。4.4 排查顺序速查表症状优先检查项处理建议几百步内数值爆炸Courant数、声速赋值降到0.4以内检查材料数组某固定点向外发散声源附近材料参数检查该点声速和密度赋值边缘反射明显PML厚度、衰减系数PML加到20层σ渐变波前不圆、有方向性网格分辨率细化到 λ/15~λ/20波后拖尾振荡网格色散细化网格降低时间步焦点位置偏移源延时、声速取值复查各阵元延时计算这个排查顺序是我自己反复用过的。每调一个参数只改一处重新跑一次对比别同时动网格又动声源又改边界否则问题源头根本定位不出来。5. 一点个人经验最后说一个我自己用过才知道的体会FDTD程序不是“能跑起来”就完事的。最开始我图省事把网格设到 ( \lambda/6 )算出来的焦斑位置偏移了将近半个波长当时一度以为是声源延时算错了。后来逐项排查发现单纯是网格色散导致的等效相速度偏差累积把网格细化到 ( \lambda/15 ) 之后问题自动消失。这说明网格分辨率、时间步长、边界厚度这些参数不是简单地“设一个数”就行它们共同决定了结果的可信度。如果你也想自己写一版超声声场FDTD程序我的建议是别急着直接上复杂的换能器阵列模型。先用水介质中的点源把圆柱波传播验证好再考虑加层状介质、加PML、加阵列聚焦。前面几步走稳了后面的功能扩展就是水到渠成的事。计算域规模也可以从小到大慢慢试先在低频小模型上把整个流程跑通再逐步加大算例这样定位问题会轻松很多。本文还有配套的精品资源点击获取