企业官网建设流程全解析

简介：基于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/ρ)∇p

FDTD 的做法是把这两个方程在空间和时间上交错离散。声压 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; end

PML 系数的核心是 σ 从内边界向外边界逐渐增大，这样阻抗渐变，反射最小。阶次 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. 更新 ux（x 方向振速），乘 PML 衰减 ux = sx_face .* (ux - dt/(rho*dx) * (p(2:end,:) - p(1:end-1,:))); % 2. 更新 uy（y 方向振速），乘 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)×Ny，uy 为 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：

跑完看两个数：全局最大值位置和焦点峰值声压。焦点峰值声压理论上应该比换能器表面声压高出一个数量级，这是几何聚焦的增益。如果峰值位置偏离网格中心超过 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 mm，0.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 步的快速测试，确认焦点数量和位置对了才放完整计算，这动作帮我省下了大量反复调参的时间。希望帮到你。

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