简介一份基于MATLAB的超声声场有限差分时域计算程序面向超声医学、高强度聚焦超声治疗和无损检测领域的研究人员与工程师用于模拟声波在介质中的传播、反射与衰减并辅助探头设计及治疗剂量优化。压缩包内含十个文件整体大小约三百三十一千字节其中四个脚本文件负责主程序与可视化三个数据文件保存介质与声场参数另有C语言源码及动态链接库文件便于跨平台运行与加速计算。目前已有四百零八人学习下载。读者可获得完整的模拟代码与配套数据掌握网格初始化、时间步进、边界处理、结果可视化等核心流程并能针对不同频率和功率的超声源调节介质参数直观分析高强度聚焦超声等应用场景下的声场分布为实验研究或临床规划提供可靠参考。1. 超声声场 FDTD 计算为什么 HIFU 模拟离不开时域有限差分做超声换能器设计或者 HIFU 治疗规划的人迟早会撞上一个问题声场到底怎么算解析法在简单几何下还能用一旦碰到多层组织、非均匀介质、不规则边界解析解就彻底失灵了。FDTD时域有限差分把波动方程在空间和时间上直接离散一个一个时间步往前推虽然算得慢但胜在能处理任意非均匀介质而且能直接拿到每个网格点的时域波形。这套 MATLAB 程序把 FDTD 的核心流程——网格生成、激励源设置、迭代计算、声压提取——全部封装好了拿到手改几个参数就能跑出自己的声场分布。对做 HIFU 焦点控制、换能器阵元布局验证、或者想搞懂 FDTD 原理的人来说这份代码值得花一下午拆一遍。2. FDTD 声场模拟的原理与网格设定从波动方程到离散迭代2.1 波动方程的 FDTD 离散形式速度与压力的交错网格FDTD 求解声场的基本思路不复杂就是把一阶压力-速度方程组在时间和空间上交错采样。常见的做法是 Yee 式交错网格压力 p 定义在网格中心速度分量 ux、uz 定义在网格边界上。这样做的原因是空间差分可以做到二阶精度时间上也是中心差分整个格式的数值色散能控制在可接受范围。MATLAB 里实现这个交错网格通常用数组索引错位的方式不需要真正建立两套网格坐标。代码里体现出来就是速度场的更新要用到压力场的相邻点% 更新 x 方向速度分量 ux(2:end-1, 2:end-1) ux(2:end-1, 2:end-1) - ... dt/rho(2:end-1, 2:end-1) .* (p(3:end, 2:end-1) - p(1:end-2, 2:end-1)) / dx; % 更新 z 方向速度分量 uz(2:end-1, 2:end-1) uz(2:end-1, 2:end-1) - ... dt/rho(2:end-1, 2:end-1) .* (p(2:end-1, 3:end) - p(2:end-1, 1:end-2)) / dz;这里rho是介质密度矩阵dx和dz是空间步长。注意数组索引的偏移p(3:end) - p(1:end-2)本质上是中心差分但起点偏移了半个网格这就是交错网格的体现。压力场的更新稍微不同需要用到体积模量K或者直接用声速c和密度rho换算% 更新压力场 p(2:end-1, 2:end-1) p(2:end-1, 2:end-1) - ... K(2:end-1, 2:end-1) * dt .* ... ( (ux(2:end-1, 2:end-1) - ux(1:end-2, 2:end-1)) / dx ... (uz(2:end-1, 2:end-1) - uz(2:end-1, 1:end-2)) / dz );参数说明dt是时间步长K rho * c^2。这套迭代每走一步就是完成一次速度场和压力场的交替更新。最外层再套一个时间循环循环次数乘以dt就是总的模拟时长。实际写代码时建议把rho、K、dx、dt都放在一个配置结构体里管理这样后面做参数扫描时只需要改结构体字段不用动核心迭代逻辑。2.2 稳定性条件与网格尺寸选择CFL 数不是摆设FDTD 的稳定性条件用脚趾头想都知道信息在一个时间步内不能跨过一个网格。具体到声学方程要求dt 1 / (c_max * sqrt(1/dx^2 1/dz^2))在均匀网格下简化成c_max * dt / dx 1/sqrt(2)。这个比值就是 CFL 数工程经验上取 0.3 到 0.5 比较稳妥取太接近上限虽然算得快但数值色散会变大波形会变形。网格尺寸的选择跟频率直接绑定。一个波长内至少要有 10 个网格这是底线。HIFU 常用的频率在 0.5 MHz 到 5 MHz 之间水中的声速约 1500 m/s对应波长 0.3 mm 到 3 mm。如果你想分辨谐波分量网格还得更细因为谐波的波长更短。代码里可以在初始化阶段自动计算网格数c0 1500; % 介质声速 freq 1e6; % 激励频率 1 MHz lambda c0 / freq; % 波长 dx lambda / 12; % 每个波长 12 个网格 Nz round(60e-3 / dx); % 模拟区域 60 mm Nx round(40e-3 / dx); % 模拟区域 40 mm dt 0.4 * dx / c0; % CFL 0.4dx的选择要同时考虑分辨率和计算量。网格细化一倍二维计算量是原来的四倍三维则是八倍这是个指数爆炸的关系所以一开始别贪心先用粗网格把流程跑通再用细网格做精算。2.3 吸收边界PML 层的设置与参数调试吸收了边界反射FDTD 声场模拟才算真正能看。最简单的吸收边界是一阶 Mur实现容易但吸收效果一般尤其是在斜入射时反射率偏高。做 HIFU 这种需要长时间模拟的场景边界反射会污染焦点的声压值所以更推荐用 PML完美匹配层。PML 在 MATLAB 里的实现不复杂本质上是给边界区域额外引入一个衰减因子让波在里面指数衰减。关键的参数是衰减系数sigma的分布形式通常用多项式渐变% PML 参数设置 pml_thickness 20; % PML 层厚度单位网格数 sigma_max 2.0; % 最大衰减系数 x_pml linspace(0, 1, pml_thickness1); sigma_x sigma_max * x_pml(1:end-1).^2; % 二次方渐变sigma_max取值太大会产生数值反射太小吸收不干净。经验上先按sigma_max (0.5~2) / dx的范围试看边界处最大反射幅值是否小于主声场的 1%。判断 PML 效果的办法很简单做一个没有散射体的均匀介质模拟看计算域边缘的残余声压幅值对比中心区域低多少低两三个量级就算及格。3. 换能器建模与声场计算把激励源“放进”计算域3.1 激励源的建模连续波、脉冲波与阵元相位控制FDTD 里的激励源本质上是给某个网格点或多个网格点施加一个随时间变化的压力扰动。激励函数的选择直接决定模拟结果的性质——想知道稳态声场分布就用连续正弦波想知道时间响应就用来几个周期的脉冲。最常见的激励是高斯包络正弦波频带窄、旁瓣低既能保证主频能量集中又不会因为突变引起数值震荡t (0:nt-1) * dt; f0 1e6; % 中心频率 1 MHz ncyc 5; % 脉冲周期数 t0 ncyc / f0 / 2; % 高斯包络中心时间 sigma_t ncyc / f0 / 4; % 包络宽度 source exp(-((t - t0).^2) / (2 * sigma_t^2)) .* sin(2 * pi * f0 * t);参数说明ncyc决定脉冲的频带宽度周期数越少频带越宽在高衰减介质中传播后波形畸变也越明显。做 HIFU 连续波模拟时源就直接用sin(2*pi*f0*t)但要等波场充分建立模拟时间要够长否则焦点处的声压没达到稳态。对于换能器阵列每个阵元相当于一个独立源用不同相位激励就能实现波束聚焦。相位的计算基于阵元位置到焦点距离的差for i 1:N_elem dist sqrt((xelem(i) - x_focus).^2 (zelem(i) - z_focus).^2); phase(i) mod(dist / c0 * 2 * pi * f0, 2*pi); source_signal(i, :) sin(2 * pi * f0 * t - phase(i)); end这里的phase就是每个阵元的激励延迟单位是弧度。实际代码里要做的是把延迟换算成时间步判断是否需要整数步插值否则相位误差会直接影响焦点质量。3.2 计算域初始化与参数设置密度、声速、吸收系数的赋值计算域里的每个网格点都要有对应的介质参数。均匀介质最简单——所有点赋同一个rho和c。但做 HIFU 仿真肯定要分层或者放个组织模型进去这时候参数矩阵就得按几何形状填值。代码里典型做法是先初始化全零矩阵然后用逻辑索引或者循环赋值rho_map 1000 * ones(Nz, Nx); % 默认水 c_map 1500 * ones(Nz, Nx); % 在某个区域内设置组织参数 mask (z 20e-3) (z 30e-3); rho_map(mask) 1050; % 组织密度 c_map(mask) 1580; % 组织声速 % 体积模量矩阵 K_map rho_map .* c_map.^2;注意mask的生成要用二维坐标网格[zz, xx] meshgrid(z, x)转置后和mask维度对齐。这里最容易出的错是行列顺序搞反Nz对应行数z 方向、Nx对应列数x 方向循环里一旦写错介质模型就变成横着躺的了。参数矩阵准备好后要单独写一段验证代码把rho_map和c_map用imagesc画出来看一眼确认形状和位置正确再开始迭代。3.3 声场后处理从原始网格数据到声压分布图算完的原始数据是每个时间步的压力场直接看的话是一堆动态图并不直观。常用做法是做两类后处理一类是取某个时刻的全场幅值快照另一类是取某个点的时域波形做频谱分析。对 HIFU 来说最有价值的展示是稳态声场的峰值包络分布也就是把每个网格点在所有时间步里的最大正压值记录下来peak_pressure max(p_history, [], 3); % p_history 是 Nz x Nx x Nt 的三维数组 imagesc(x*1000, z*1000, peak_pressure); colorbar; xlabel(x (mm)); ylabel(z (mm));这里p_history如果全存下来内存压力很大。一个1000x1000x2000的单精度数组就要 8 GB。常见做法是每隔若干步记录一次压力场或者只记录每个点的峰值和均方根值不保留完整时间序列。代码里要加一个判断if mod(n, record_interval) 0 p_history(:, :, n/record_interval) p; endrecord_interval的取值取决于你需要的波形精度做频谱分析至少要满足奈奎斯特采样率所以record_interval * dt不能超过最小关注周期的一半。4. HIFU 聚焦声场仿真换能器阵列的参数标定与焦点控制4.1 相控阵列相位延迟计算球面波聚焦的离散化实现HIFU 换能器阵列的聚焦原理说白了就是相位补偿——让每个阵元发出的波前到达焦点时同相叠加声压在焦点处增强焦点之外相互抵消。圆形活塞换能器本身就自带几何聚焦但阵列换能器可以通过电子相位控制灵活移动焦点这才是 HIFU 治疗的核心技术。相位延迟的精确计算要考虑到阵元不是点源而是一个面。近场效应会让焦点位置略偏移几何中心。代码里的简化做法是取阵元几何中心到焦点的距离差% 阵列参数 elem_radius 2e-3; % 单个阵元半径 pitch 3e-3; % 阵元间距 N_elem 32; % 阵元数量 focus [0, 40e-3]; % 焦点坐标 (x, z) % 计算每个阵元的相位延迟 x_pos ((1:N_elem) - (N_elem1)/2) * pitch; for ii 1:N_elem dist sqrt(x_pos(ii)^2 focus(2)^2); delay_time(ii) (dist - focus(2)) / c0; phase_delay(ii) 2 * pi * f0 * delay_time(ii); endphase_delay算出来后每个阵元的激励信号写成sin(2*pi*f0*t - phase_delay(ii))。有一个细节如果直接把相位取模到[0, 2*pi)会造成时间上的不连续体现在声压波形上就是初始段有瞬态畸变。更好的做法是保留绝对延迟时间让所有阵元在启动时有一段时间的静默期等最远的阵元开始工作后整个阵列才开始稳定辐射。4.2 焦点验证半波宽、焦点偏移与声场增益的量化评估跑完仿真别急着出图量化验证焦点质量至少要算三个指标焦域尺寸-6 dB 半波宽、焦点峰值压力相对阵元表面的增益、焦点位置与预设位置的偏差。用 MATLAB 提取这些指标很简单从峰值包络分布图上找最大值点再沿 x 和 z 方向做剖面线% 找焦点位置 [max_val, idx] max(peak_pressure(:)); [row, col] ind2sub(size(peak_pressure), idx); focus_x_actual x(col) * 1000; % 转成 mm focus_z_actual z(row) * 1000; % 沿 x 方向提取 -6 dB 宽度 x_profile peak_pressure(row, :); half_max max_val / 2; above_idx find(x_profile half_max); beam_width (x(above_idx(end)) - x(above_idx(1))) * 1000;beam_width就是焦域的横向宽度HIFU 的典型值是几个毫米。如果算出来的焦点比几何预估偏了超过 5%大概率是 PML 层反射或者网格色散造成的相位累积误差要回头看边界设置和网格分辨率。焦点增益聚焦增益的定义是焦点声压与阵元表面平均声压之比理论上阵元越多增益越高但受阵元间距产生的栅瓣限制。仿真里如果栅瓣出现说明阵元间距大于半波长这是阵列设计阶段就能预判的问题。4.3 参数扫描频率、孔径与聚焦深度对焦域的影响参数扫描是 FDTD 仿真最有价值的用途之一。改一个参数、重跑一次、对比指标迭代出最优设计。但盲扫的效率太低得有策略地扫。先固定焦点位置和阵列几何扫频率。频率升高焦域变小变窄但组织中的衰减也增大焦点实际声压可能反而下降。衰减系数一般按频率的 1.1 次方线性增加所以高频并不是越高越好有个最佳平衡点freq_list [0.5e6, 0.75e6, 1e6, 1.5e6, 2e6, 3e6]; for kk 1:length(freq_list) params.f0 freq_list(kk); % 重新生成激励源 lambda c0 / params.f0; params.dx lambda / 12; % 更新网格注意整个网格模型要跟着变 params rebuild_grid(params); % 跑仿真 result(kk) run_fdtd_sim(params); end这里有个坑只改频率不改网格尺寸等效于波长内网格数变了结果对比就没有意义。正确做法是每次调整网格保证网格分辨率一致。孔径和聚焦深度的扫描逻辑类似但每次改阵列几何后要重新确认阵元在网格上的位置映射是整数索引否则会出现亚网格偏差导致相位计算错误。参数扫描跑完之后把结果整理成表格——频率、孔径、焦点尺寸、峰值声压、旁瓣水平——一眼就能看出哪个参数组合最合适做后续实验验证。5. 常见问题与排查FDTD 仿真中那些“看起来对其实错”的坑5.1 问题一声压值随时间步数线性增大不收敛现象模拟跑了上千步焦点处的声压幅值还在持续涨完全没有饱和趋势。原因最常见的两种情况一是激励源是连续正弦波但模拟区域边界没有吸收或者吸收不够波在计算域里来回反射累加二是介质参数设置出错导致局部声速异常能量被限制在一个小区域里堆积。解决先用脉冲激励代替连续波如果脉冲声场在边界处有明显反射回波先优化 PML 参数再用均匀介质模拟做一个基准测试确认整洁的介质模型下声压能在几十个周期内稳定最后检查rho和K矩阵是否有零值或负值MATLAB 里除零会静默产出 Inf然后被传播到整个波场。5.2 问题二边界处出现明显反射回波但 PML 已经设了现象压力分布图上靠近边界的位置有一条明亮的弧形条纹从外往内传播污染了内部场。原因PML 层的衰减系数sigma_max取值不合适太大在层与主区域界面处阻抗失配反而产生反射太小则波还没被衰减就穿透 PML 到达截断边界。另外检查 PML 层是否覆盖了所有外边界四个角区域容易漏掉。解决做一个无源场的测试从中心放一个脉冲源看边界反射幅值。把sigma_max按(0.5/dx)、(1/dx)、(2/dx)各跑一遍选反射最低的值同时确认 PML 区域厚度不少于 10 个网格推荐 20 个网格以上太薄了吸收效果就是玄学。5.3 问题三焦点位置整体偏移但网格和源都检查过现象预设焦点在 (0, 40 mm)实际算出来的焦点在 (-2 mm, 38 mm)偏移量还不小。原因阵元相位延迟计算里忽略了阵元的几何尺寸把大面积阵元当点源处理。阵元表面是一个曲面或平面不同位置的相位其实不同。还有一个隐蔽原因是阵元所在的介质声速和计算用声速不一致如果阵元表面覆盖了匹配层或透镜声速差异引起的相位差不可忽略。解决把阵列模型细化每个阵元在网格上对应多个源点每个源点独立计算延迟或者至少用等效声程法把阵元表面上距焦点最近和最远的点到焦点的距离差平均掉。更稳妥的办法是做一个没有相位控制的测试先验证单个阵元的声场对称性再叠加相位控制看焦点偏移。5.4 问题四算得实在慢一个网格数 500x500 的二维模拟要跑一晚上现象时间循环内每一步的计算量巨大MATLAB 的循环效率奇低。原因代码里大量使用for嵌套循环更新每个网格点而 MATLAB 的性能优势在于向量化运算。二维 FDTD 的核心更新用矩阵切片完全可以向量化。解决把内部的空间循环全部改写成矩阵运算for n 1:nt % 速度场更新向量化 ux(2:end-1, 2:end-1) ux(2:end-1, 2:end-1) - ... dt ./ rho_core .* (p(3:end, 2:end-1) - p(1:end-2, 2:end-1)) / dx; % 压力场更新向量化 p(2:end-1, 2:end-1) p(2:end-1, 2:end-1) - ... K_core * dt .* ( (ux(2:end-1, 2:end-1) - ux(1:end-2, 2:end-1))/dx ... (uz(2:end-1, 2:end-1) - uz(2:end-1, 1:end-2))/dz ); end这样改写之后速度能提升一两个数量级。如果还嫌慢下一步是开parfor做参数扫描或者上 GPU 计算MATLAB 的gpuArray对 FDTD 这种规则网格计算天然友好。但 GPU 版调试起来麻烦建议先把 CPU 向量化版本调通了再考虑。6. 进阶技巧交叉验证与自检——当系统做过全局对比验证用一套自研 FDTD 跑出来的结果直接用于判定设计总觉得不踏实。一个简单但极其有效的习惯是架构性地准备一次交叉验证。找一个声誉极佳的独立求解器作为参照——比如 K-Wave 工具箱同样的几何、同样的激励参数算两遍放在一起对比声压分布图和时间波形。对比不是肉眼看形状像就收工要看具体的量化指标焦点峰值声压的相对误差是否在 5% 以内、焦域半波宽的绝对误差是否小于一个网格分辨率、焦点位置偏差是否小于两个网格步长。当一套自研代码过了一致性检验之后后续的改动就敢放手做了。从那以后我每次大改控制参数都会强制走一遍同样的对比验证流程花不了多少时间但能避免很多“看着没问题实际上边界条件都错位了”的假象希望帮到你。本文还有配套的精品资源点击获取