简介基于MATLAB的超声声场FDTD计算程序面向HIFU治疗规划、超声成像与无损检测领域的科研和工程人员用于模拟超声波在复杂介质中的传播、衰减与聚焦。压缩包共10个文件包含4个M脚本、3个MAT数据文件以及DLL、MEXGLX和C源码分别承担核心算法、参数存储与编译加速功能整体仅331KB轻量完整。目前已有408人学习使用。借助这套程序读者能直接运行FDTD主流程观察网格初始化、时间步进、边界处理与后处理的具体实现通过调整超声频率、功率和介质参数可拓展至多种声场场景为HIFU剂量控制和探头设计提供仿真依据。代码模块划分清晰适合具备一定仿真基础的科研人员快速开展二次开发与参数研究。1. 超声声场FDTD计算程序把 HIFU 焦点声场从偏微分方程变成一张能用的图做 HIFU高强度聚焦超声换能器设计时很多人第一反应是用瑞利积分算焦点声场解析法快几秒出图。可一旦把组织分层、颅骨、反射边界放进来解析解根本推不动你只能老老实实回到数值方法。这套超声声场 FDTD 计算程序解决的就是这个痛点在 MATLAB 里把时域有限差分方法铺到二维声压网格上完整模拟换能器激励后声压场的瞬态演化再提取焦点位置、焦域尺寸、旁瓣水平这些设计指标。适合做 HIFU 单阵元聚焦评估、多阵元相位控制验证以及想搞明白声波在非均匀介质里到底怎么走的开发者。2. 技术选型与原理为什么 FDTD 能算超声解析法却在 HIFU 上失灵HIFU 声场计算的核心困难在于介质不均匀、边界形状复杂、还要看瞬态过程。解析法假设均匀介质、无限大空间一旦模型里出现分层或倾斜边界公式就变成一摊没法解的积分。FDTD 的思路是把连续偏微分方程离散到网格上一步一步往前推时间每个网格点的声压都算出来边界和介质变化只是改几个系数的事。这套程序选择二维声场作为切入点网格规模可控在普通笔记本上就能跑同时保留了 FDTD 处理复杂模型的核心能力。2.1 从波动方程到差分格式速度-压力交错网格线性声学里声压 p 和质点振速 u 满足一阶偏微分方程组∂p/∂t -ρc²(∇·u) ∂u/∂t -(1/ρ)∇pFDTD 的做法是把这两个方程在空间和时间上交错离散。声压 p 放在整数网格点上x 方向振速 ux 放在 x 方向偏移半格的位置y 方向振速 uy 放在 y 方向偏移半格的位置时间上也是错半步更新形成一个蛙跳推进格式。这样做的好处是二阶精度、数值稳定性好而且不需要解大型线性方程组每一步只做局部的差分运算。% 交错网格参数示例 dx 0.1e-3; % 网格步长 0.1 mm dy dx; c0 1500; % 水中声速 m/s rho 1000; % 密度 kg/m^3 CFL 0.5; % CFL 数必须小于 1/sqrt(2) dt CFL * dx / c0; % 时间步长约 33 ns这里的 CFL 数是稳定性条件对应二维情况下 c*dt/dx 必须小于 1/√2实际程序里取 0.5 留出余量。如果 CFL 数大于临界值迭代几步后数值就会发散声压场直接溢出成 NaN这是 FDTD 新手最容易遇到的第一道坎。2.2 PML 吸收边界与源项注入边界不反射能量进得来有限计算域的边界必须处理否则声波撞到边界会反射回来污染场内声压。这套程序在计算域四周加了 PML完美匹配层吸收边界层内引入衰减系数 σ振速更新时乘以 exp(-σdt)让波在边界附近指数衰减。激励源则作为质量源项注入在换能器表面的网格点上每步叠加一个源信号值模拟换能器向介质中辐射声波。% PML 衰减系数生成一维示例二维同理 d_pml 20; % PML 层数 m 3; % 多项式阶次 R 1e-4; % 理论反射系数 sigma_max -log(R) * c0 * (m 1) / (2 * d_pml * dx); sigma_x zeros(Nx, 1); for i 1:d_pml sigma_x(i) sigma_max * ((d_pml - i 0.5) / d_pml)^m; sigma_x(Nx - i 1) sigma_max * ((d_pml - i 0.5) / d_pml)^m; endPML 系数的核心是 σ 从内边界向外边界逐渐增大这样阻抗渐变反射最小。阶次 m 取 3 是常见做法反射系数 R 设 1e-4 已经比实际需求的精度高很多再小数值上容易出问题。2.3 程序目录结构与数据流向这套程序按功能拆成了几个脚本每个文件只干一件事fdtd_ultrasound/ ├── main.m % 程序入口按顺序调用各模块 ├── setup_grid.m % 网格尺寸、介质参数、CFL 校验 ├── setup_source.m % 换能器几何、激励信号、时间延迟 ├── setup_pml.m % PML 层数与衰减系数 ├── fdtd_step.m % 单步时间迭代核心函数 ├── compute_field.m % 焦点提取、场分布指标计算 └── plot_results.m % 声压场可视化数据流是单向的setup_grid 先生成网格坐标系和介质数组setup_source 根据换能器几何算出哪些网格点属于激励面以及每个点的激励延迟setup_pml 生成吸收系数然后 main.m 循环调用 fdtd_step 推进时间步最后 compute_field 从完整场数据里提取指标plot_results 画图。这种一文件一职责的结构后期改成三维或者加非线性项只需要替换对应的模块不用重写整套流程。3. 核心模块拆解一个声脉冲从激励面到焦点的完整链路这一章直接进代码。我按程序实际执行顺序把每个关键模块挑出来讲清楚参数配什么值、激励怎么加、迭代循环里发生了什么、最后怎么把焦点数据抠出来。3.1 setup_grid网格、介质与 CFL 稳定性校验网格参数决定了计算精度和内存开销。这个程序的默认参数按 1 MHz 水中聚焦超声设计波长约 1.5 mm网格步长 0.1 mm 意味着一个波长覆盖 15 个网格点对线性声场来说精度足够。% setup_grid.m Nx 400; Ny 400; % 网格点数对应 40 mm x 40 mm 区域 dx 0.1e-3; dy 0.1e-3; % 空间步长 0.1 mm c0 1500; rho 1000; Nt 6000; % 时间步数对应约 200 μs 物理时长 % CFL 稳定性校验 cfl c0 * dt / dx; if sqrt(cfl^2 * 2) 1 error(CFL 条件不满足请减小 dt 或增大 dx); end参数说明Nx 和 Ny 取 400 是平衡点网格太少焦点区分辨率不够太多则内存和时间成本翻倍。Nt 取 6000 对应物理时长约 200 μs足以让声波从换能器传播到焦点再稳定下来。cfl 必须满足二维稳定性条件程序里直接写死校验防止改参数时把时间步长改出问题。3.2 setup_source圆弧换能器的几何聚焦与时间延迟换能器用圆弧表面模拟几何聚焦靠的是凹面到焦点距离不同导致的时间延迟差。程序对每个激励网格点计算它到几何焦点的距离再把距离差换算成时间延迟激励信号按延迟依次发出波前在焦点处自然叠加增强。% setup_source.m src_radius 80e-3; % 凹面曲率半径 80 mm aperture 50e-3; % 换能器孔径 50 mm f0 1e6; % 中心频率 1 MHz focus [Nx/2, Ny/2]; % 几何焦点位于网格中心 src_mask zeros(Nx, Ny); % 滑动式判断网格点是否在换能器表面 for i 1:Nx for j 1:Ny x (i - focus(1)) * dx; y (j - focus(2)) * dy; if abs(sqrt(x^2 y^2) - src_radius) dx/2 sqrt(x^2 y^2) aperture/2 src_mask(i, j) 1; end end end % 计算每个激励点的时间延迟 delay zeros(Nx, Ny); for i 1:Nx for j 1:Ny if src_mask(i, j) 0 dist sqrt(((i - focus(1))*dx)^2 ((j - focus(2))*dy)^2); delay(i, j) (src_radius - dist) / c0; % 距离差转时间差 end end end时间延迟的本质是让凹面上各点发出的波在几何焦点处同相叠加。曲率半径越大聚焦越浅孔径越大焦点越尖。注意程序里用 abs(...) dx/2 判断换能器表面这个容差直接决定激励面的厚度太薄会漏网格点太厚会让换能器带宽变宽。3.3 fdtd_step主迭代循环里的 PML 与源项这是整个程序最核心的部分。每次迭代做三件事更新两个方向的振速、更新声压、注入激励源。振速更新时乘上 PML 系数声压更新用交错差分源项在声压数组的激励面上叠加当前时刻的信号值。% fdtd_step.m核心迭代 function p fdtd_step(p, p_old, ux, uy, sx_face, sy_face, params) dt params.dt; dx params.dx; dy params.dy; rho params.rho; c0 params.c0; % 1. 更新 uxx 方向振速乘 PML 衰减 ux sx_face .* (ux - dt/(rho*dx) * (p(2:end,:) - p(1:end-1,:))); % 2. 更新 uyy 方向振速乘 PML 衰减 uy sy_face .* (uy - dt/(rho*dy) * (p(:,2:end) - p(:,1:end-1))); % 3. 更新声压 p中心网格点 p(2:end-1, 2:end-1) p(2:end-1, 2:end-1) ... - rho*c0^2*dt * ( ... (ux(2:end-1, 2:end-1) - ux(1:end-2, 2:end-1)) / dx ... (uy(2:end-1, 2:end-1) - uy(2:end-1, 1:end-2)) / dy ); end这段代码里交错网格的维度关系是难点。ux 矩阵尺寸为 (Nx-1)×Nyuy 为 Nx×(Ny-1)两者都不包括边界声压中心点求导刚好用到两个振速的内部点索引对齐要小心。PML 系数 sx_face、sy_face 在 setup_pml 中生成尺寸与对应振速矩阵一致。主循环里注入激励源时时域包络要加斜坡防止突然激励产生宽带分量污染结果% main.m 中的激励注入 ramp_time 2e-6; % 2 μs 斜坡时间约 2 个周期 for n 1:Nt t (n-1) * dt; envelope min(1, t / ramp_time); % 线性斜坡 src_wave 1.0 * sin(2*pi*f0*t) * envelope; p fdtd_step(...); % 执行一步迭代 p(src_mask 0) p(src_mask 0) src_wave * dt; % 压电面加源 end斜坡时间很关键取 2 个周期左右可以平滑启动。如果直接上满幅正弦换能器表面会辐射出很宽的频谱焦点区会出现非预期的旁瓣。3.4 compute_field 与 plot_results焦点提取与场分布输出迭代完成后compute_field 从声压场数据里提取峰值声压分布和焦点指标。通常做法是取每个网格点在全部时间步上的绝对值最大值得到空间峰值声压分布图再找全局最大值的位置作为实际焦点。% compute_field.m peak_p max(abs(p_hist), [], 3); % p_hist 是所有时间步的累积场 [max_val, idx] max(peak_p(:)); [fx_idx, fy_idx] ind2sub(size(peak_p), idx); % 计算 -3 dB 焦域尺寸相对峰值下降 3 dB 的区域直径 threshold max_val / sqrt(2); [rows, cols] find(peak_p threshold); focal_width (max(cols) - min(cols)) * dx; % 焦域横向宽度 focal_depth (max(rows) - min(rows)) * dy; % 焦域纵向深度焦点偏移量是 (fx_idx-Nx/2)*dx 和 (fy_idx-Ny/2)*dy对于几何聚焦的换能器这个值应该接近零。如果偏移超过一个波长说明激励延迟或者网格映射有问题。plot_results 用 imagesc 加坐标轴画峰值声压云图焦点区在图上表现为一个亮斑旁瓣是环绕焦点的若干暗环。4. 运行与验证在 MATLAB 里跑通 HIFU 水中聚焦场原理和代码都过了一遍这一章把程序跑起来并验证结果是否合理。运行环境很常规MATLAB R2019b 以上不需要任何额外工具箱脚本按顺序执行即可。4.1 环境要求与脚本入口程序是纯 MATLAB 脚本不依赖 Signal Processing Toolbox 或 Parallel Computing Toolbox向量化写得好单核跑就够了。入口是 main.m它会依次调用 setup_grid、setup_source、setup_pml进入时间迭代循环最后调用 compute_field 和 plot_results。# 在 MATLAB 命令行中运行 cd fdtd_ultrasound main第一次跑建议把 Nt 从 6000 降到 2000先验证流程能通再逐步加长模拟时间。主循环里可以在每 500 步输出一次当前进度方便观察计算是否正常推进。4.2 一组可以直接抄的标准算例参数我常用来做基准测试的参数组合如下这个配置在 8GB 内存的机器上约 8 分钟跑完峰值内存约 2GB参数取值说明网格点数400 × 40040 mm × 40 mm 计算域网格步长0.1 mm每波长 15 个点时间步数6000物理时长 200 μs时间步长33 nsCFL 0.5孔径50 mm换能器口径曲率半径80 mm几何焦距频率1 MHz水中波长 1.5 mmPML 层数20边界吸收厚度跑完看两个数全局最大值位置和焦点峰值声压。焦点峰值声压理论上应该比换能器表面声压高出一个数量级这是几何聚焦的增益。如果峰值位置偏离网格中心超过 1 mm优先检查时间延迟计算的符号这是最容易翻车的点。4.3 结果验证焦点位置、焦域半宽和旁瓣验证算得对不对不能只看图要看数值指标。我一般对比三样东西焦点横向尺寸、焦点深度、主旁瓣比。对于 1 MHz、f/数约 1.6 的换能器焦点横向宽度应该在 1.5 mm 左右纵向深度约 5 mm主旁瓣比通常大于 10 dB。% 验证指标输出 fprintf(焦点位置偏移: %.2f mm, %.2f mm\n, fx_offset*1e3, fy_offset*1e3); fprintf(焦点横向宽度: %.2f mm\n, focal_width*1e3); fprintf(焦点纵向深度: %.2f mm\n, focal_depth*1e3); fprintf(主旁瓣比: %.1f dB\n, 20*log10(max_val / sidelobe_max));如果算出来的焦域尺寸偏差超过 30%先看网格分辨率再确认激励斜坡时间。网格步长从 0.1 mm 加到 0.15 mm 后焦点尺寸误差能明显变大这就是色散误差在作怪。焦域的数值模拟结果和解析解对不上时八成不是 FDTD 公式问题而是激励源或者边界处理的细节没到位。5. 参数调优与常见问题排查那些容易翻车的细节用这套程序跑了几百组参数后我把最常踩的坑整理在这里。每一条都是真实翻车经验按现象、原因、解决的顺序写排查时直接对着看。5.1 网格步长太粗场图出现方向性畸变现象焦点周围出现沿网格轴方向的条纹焦点形状不圆左右上下不对称。原因网格步长超过波长的 1/10 后FDTD 的数值色散开始显现不同传播方向的波速不一样导致原本圆对称的场图出现方向性畸变。1 MHz 水中波长 1.5 mm0.2 mm 步长就是 1/7.5 波长畸变已经肉眼可见。解决把 dx 和 dy 同时改成 0.1 mm 或更小。注意只改一个轴没用两个方向必须一致否则椭圆畸变更严重。内存吃紧时优先减小计算域而不是放大网格步长因为色散误差没法靠后处理补偿。5.2 PML 层太薄或衰减系数不对边界反射污染焦点现象场图边界附近出现明显的同心圆条纹焦点旁边多出虚假的旁瓣峰值位置轻微偏移。原因PML 层数少于 10 层且 σ 曲线太陡时数字反射比理论上大得多边界反射波回到焦点区域和正常波场叠加成干涉条纹。解决把 d_pml 从 10 调到 20σ 的阶次 m 保持 3。判断 PML 是否有效的方法是把计算域扩大一倍再跑一遍对比焦点区的声压曲线如果两次结果差异大于 2%说明边界反射还在场里《指南》里讲这个阈值时我还没意识到后来对比了几组才确认。5.3 激励源斜坡太陡连续波算不出稳定焦点现象焦点声压随模拟时间波动不收敛频谱图上出现除 1 MHz 外的额外峰。原因连续正弦波如果从第一拍就满幅注入等效于一个阶跃激励频谱从中心频率向两侧展宽低频分量穿透力强破坏了聚焦的相位关系。解决斜坡时间从 0 改到至少 2 个周期用 min(1, t/ramp_time) 线性斜坡就可以。想更平滑可以用余弦平方斜坡但对 HIFU 场景线性斜坡的精度已经够用。验证方法是看焦点声压包络是否在 50 μs 后进入稳定震荡如果还在波动就继续加长斜坡。5.4 计算时间和内存失控现象把网格从 400×400 改成 800×800 后内存直接爆掉一个场的历史缓存就占了几个 GB。原因如果保存场的历史数据用三维数组 Nx×Ny×Nt 累积800×800×12000 个 double 数组约 60 GB自然必崩。很多新手在这里摔过认为 FDTD 是内存怪兽其实是保存时机不对。解决不要全时间步都存在内存里。常见做法是每个网格点只保留峰值用一个 Nx×Ny 的数组边迭代边更新 max 值内存占用直接降到几 MB。另外可以用单精度 float 代替 double误差在可接受范围内内存再减一半。时间不够就从 Nt 下手先跑 2000 步确认流程再跑完整 6000 步。6. 二次开发思路把单焦点几何聚焦改造成多焦点相控阵基础的单阵元流程跑通后二次开发的第一站是做多焦点相控阵。HIFU 临床上常需要同时加热多个靶点或者避开危险区域这时候换能器从单凹面变成平面相控阵阵元每个阵元独立控制相位焦点由一个变成多个。6.1 从时间延迟到相位叠加单焦点靠凹面几何形状形成天然延迟多焦点靠阵元发射信号叠加。N 个阵元同时发波焦点位置 F 处的声压等于每个阵元在该点的贡献求和p(F, t) ΣA_i·sin(2πf₀t φ_i)。要实现两个焦点就把两套相位分布叠加通过加权控制两个焦点的能量占比。% 阵元相位计算 phase1 2*pi*f0 * (dist_to_focus1 - dist_to_center) / c0; % 焦点1 phase2 2*pi*f0 * (dist_to_focus2 - dist_to_center) / c0; % 焦点2 phase_total angle(0.7*exp(1j*phase1) 0.3*exp(1j*phase2)); % 加权叠加6.2 最小改动点只改 setup_source这套程序的多焦点改造只需动 setup_source.m 一个文件。把 src_mask 从凹面网格改成阵元阵列把 delay 换成 phase_total 对应的延迟时间其他模块全部复用。迭代循环里源项注入从单一表面改成多个阵元面每个阵元的振幅和相位可以独立控制。验证时重点看两个焦点的峰值是否达到预期比例以及两个焦点之间的旁瓣水平是否超标。改完代码后我习惯先跑小网格快速验证相位是否正确再上完整算例。从那以后每次改激励源我都强制先跑一版 200×200、1000 步的快速测试确认焦点数量和位置对了才放完整计算这动作帮我省下了大量反复调参的时间。希望帮到你。本文还有配套的精品资源点击获取