去年我着手做超声阵列相关的预研第一步就卡在“怎么把一套换能器方案说清楚”上。直接买探头做水听器实验当然最真实但一轮打样加实验的周期和成本都太高。后来我把目光转到Field_II仿真才算是把“探头设计-声场分析-波束合成效果验证”这条链路完整跑通。这篇就是我自己的Field_II学习笔记包括环境配置、最小示例、核心原理理解以及我在内存爆炸、参数设置上踩过的坑。如果你正准备做超声声场仿真、脉冲回波成像仿真或者只是想找一个在MATLAB里好上手的换能器声场计算工具这篇内容应该能帮你少走不少弯路。1. 超声仿真工具不少为什么我最终选Field_II不吃灰1.1 Field_II在超声仿真版图里到底处在什么位置先说结论Field_II适合解决“给定探头参数怎么算出它在空间里的声场分布怎么模拟一个脉冲回波成像过程”这一类问题。它是丹麦技术大学Jørgen Arendt Jensen团队从90年代开始持续维护的MATLAB工具箱已经被大量学术论文引用算得上是超声成像仿真的“老牌基础设施”。超声仿真这几年能用的工具其实不少我按自己了解的程度整理过一张对比表工具/方法基本方法擅长场景主要限制Field_II空间冲激响应 线性系统叠加换能器声场、脉冲回波、波束合成不擅长非线性传播、热效应、组织结构弹性k-Wave基于k空间伪谱法求解声波动方程非线性声场传播、不均匀介质、热疗/空化换能器细节建模麻烦算得慢PZFlex / Sim4Life有限元/有限差分多物理场压电振子本身、声-固-电耦合商业化license贵建模门槛高自研射线模型几何声学近似稀疏阵列快速布局、远场近场/衍射/旁瓣基本没法看我当时的目标很明确验证一个1D线阵的阵元间距、阵元宽度、焦点位置对成像的影响先不关心换能器壳体振动、不关心非线性组织谐波。这种情况下Field_II是性价比最高的选择——完全基于MATLAB接口设计很统一官方示例代码覆盖了从“单阵元声场”到“完整B模式图像仿真”的所有层次。1.2 开始前必须接受的三个限制Field_II不是万能的网上有些朋友一上来就指望它能仿真“组织剪切波”“微泡造影振荡”那肯定要失望。我给它划定的能力边界是这样的第一它是基于线性声学假设的。声波在介质中传播被当成线性系统所有叠加都用空间冲激响应实现。对常规B模式成像、10MHz以下频率、中等声压范围这个假设基本够用但如果你要研究谐波成像的多次谐波分量或者高机械指数下微泡的非线性响应Field_II会直接失真。第二它仿真的是声场和射频信号不包含“探头结构振动”和“声-电-机械耦合”。你设置的是孔径几何、激励波形和声学参数而不是压电陶瓷材料参数。想模拟晶片本身的振型、背衬层影响那要上PZFlex这类有限元工具。第三计算速度很依赖你的实现方式。官方提供的例子里有不少是“小而美”的演示但只要你把阵列规模改成128阵元、把散射体加到几万个时间复杂度立刻爆炸。这个问题我在后面专门写一节先说结论必须理解它的计算复杂度来源才能把工具用好。2. 从下载到第一张声场图环境配置与最小示例2.1 安装和路径那些事Field_II的获取方式是直接去它的官方页面下载压缩包解压后你会看到一堆源代码文件和编译好的mex文件。它在MATLAB里没有自带安装器所有安装动作其实就是“把路径加进去”。顺手提醒一句解压后不要放在有中文路径的目录下这算是我个人踩过的最小坑但也是最容易让人崩溃的坑。如果你用Windows系统还有可能遇到编译版本兼容问题——官方提供的mex文件未必能在你的MATLAB版本里直接运行。遇到报错提示找不到field_II函数时很可能是mex二进制不兼容建议确认是否需要在当前MATLAB里重新编译。安装路径设置我习惯写成这样addpath(D:\tools\field_II\); savepath;然后全局初始化field_init(-1);field_init(-1)表示每次仿真时关闭所有旧的Field_II状态打印版本信息。脚本跑完用field_end收尾这个用法几乎是所有例子的开头和结尾你也可以养成习惯。在配置这一步最需要记住的一点Field_II的“全局属性”基本都是用过一次就不要再动的。比如set_field设定声速、采样频率、衰减系数这些一进入主循环就不要反复改否则很容易出现计算结果和参数对不上的诡异问题。2.2 一个12行代码就能跑通的连续波声场示例我这里给一个非常简化的连续波声场仿真脚本目标是让一个32阵元的1D线阵聚焦到(0, 0, 40mm)处的声场网格并绘制轴向声场分布。不追求物理上的绝对精确只用来打通“Field_II 配置孔径 配置激励 计算响应”这条主线。field_init(-1); fs 100e6; % 采样频率 100 MHz set_field(fs, fs); set_field(c, 1480); % 水中声速 % 32阵元阵元宽度0.3mm间隙0.05mm高度5mm Th xdc_linear_array(32, 0.30e-3, 0.05e-3, 5e-3); % 5MHz中心频率2周期正弦激励 t0 0:1/fs:2/(5e6); excitation sin(2*pi*5e6*t0); xdc_excitation(Th, excitation); xdc_apodization(Th, ones(1,32)); % 焦点设置在z40mm xdc_focus(Th, 0, [0 0 40e-3]); % 计算坐标轴上的点 z linspace(1e-3, 80e-3, 401); points zeros(length(z), 3); points(:,3) z(:); [hp, h] calc_h(Th, points); p_axis abs(hp); plot(z*1e3, p_axis/max(p_axis)); xlabel(轴向距离/mm); ylabel(归一化幅度); field_end;这里我用calc_h计算的是“空间冲激响应对应的时间序列”没做完整的连续波积分所以幅度只是一个相对值但你已经可以通过这个脚本看到声场能量沿轴向的变化趋势。更重要的是如果这一段跑通了后续所有复杂的脉冲回波仿真都只是在这个框架上叠加。2.3 fs和c这两个参数决定了整个仿真的可信度很多人一开始不在意fs和c直接沿用别人的配置结果换了个探头参数就怎么都对不上。先说fs。它是整个时间离散化的基准。Field_II内部会把所有时间序列按采样间隔1/fs离散如果你的激励是5MHz正弦至少要满足100MHz以上的采样率才能把脉冲包络和相位关系看清楚。过低会导致波形严重失真过高会导致内存和计算时间指数级增长。我个人的习惯是先取中心频率的20到30倍做初始验证比如5MHz探头取100MHz~150MHz信任计算结果后再决定是否降采样。再说c。这个参数直接影响延迟计算、焦点位置、波长尺度。水的声速通常取1480m/s生物软组织取1540m/s。忘了改c是最隐蔽的错误它会让你在“看上去完全正常”的结果里得到系统性偏差比如焦点的实际位置和预期差了百分之几旁瓣结构也不对。还有一个容易被忽略的参数是set_field(freq_att)和衰减相关设置。如果想完全避免介质衰减的影响可以简单用set_field(att, 0)我和绝大部分快速验证场景都是这么做的如果要加衰减后面会单独提到三种设置的含义。3. 空间冲激响应Field_II唯一的“秘密武器”3.1 用“每块小面积发声”理解空间冲激响应我第一次看Field_II的官方文档时被“空间冲激响应”这个词绕晕了。后来我用一个比较朴素的模型去理解反而一直没出过错。想象有一个矩形压电阵元被电信号激励后它的表面并不是整整齐齐地像一个活塞一样整体位移而是每一小块面积无限小的点源都在独立向空间辐射球面波。空间冲激响应这个数学工具做的事情就是把换能器表面分割成无数个小点源然后把它们在同一时刻到达某个空间点的贡献做积分叠加。这就是为什么它叫“空间”冲激响应——它描述的是给定一个点上的声压对理想脉冲激励的响应并且随位置变化。在代码层面函数calc_h(Th, points)算的就是这个输入孔径定义和一系列空间坐标输出每个坐标上的脉冲响应时间序列。得到这个时间序列后再和你的实际激励波形做卷积就得到了该点处的声压波形。整个过程可以理解为空间冲激响应 换能器几何 介质声速 - 某点的时间响应 任意激励波形 空间冲激响应 与 激励波形 的卷积3.2 为什么它比“纯射线声学”更适合换能器近场如果你只关心远场方向性完全可以用几个公式算出阵元的指向性函数。但换能器近场区域恰恰是成像应用的关键区域因为焦点、孔径、旁瓣都集中在那一段范围里。传统几何声学/射线模型在处理近场时有个致命问题它忽略了衍射效应即声波遇到孔径边缘、阵元边界后会发生“绕弯”和干涉。而空间冲激响应是直接基于瑞利积分和线性声学理论推导的天然包含了这部分信息。它把每个阵元视为扩展孔径考虑到了不同部分到观测点的距离不同、时延不同、相位不同所以能正确预测焦点附近的干涉峰值、旁瓣结构和近场“自然聚焦”效应。我做个直观类比单阵元声场在近场并不是初看上去那样“越远越弱”而是先出现若干极大极小交替的振荡区间再进入单调衰减的远场。射线声学完全看不到这层结构只有衍射积分才能复现。这也是为什么我选择Field_II做验证的原因——它能非常接近真实物理行为地把近场算对这是精确计算超声成像系统的基础。3.3 一个直观验证计算自然焦点位置并与理论比较为了验证“用Field_II仿真出来的声场具有物理意义”我做了一个非常经典的自测仿真一个圆形活塞阵元的轴向声压分布并找到声压极大值出现的位置自然焦点和理论公式比较。对于半径为a的圆形活塞频率为f介质声速为c自然焦点位置大致在z_focus ≈ a^2 / λ其中λ c / f。我用一个半径5mm的圆形阵元5MHz水中声速1480m/s那么λ约0.296mm理论自然焦点位置在z ≈ (0.005^2) / 0.000296 ≈ 84.5mm在Field_II中对应实现只需要把孔径从线阵换成圆形曲面活塞然后在轴上采样幅度寻找峰值。我得到的结果和理论值误差在百分之几以内原因主要是矩形离散近似圆形的网格精度有限。这个验证给了我很大信心Field_II的计算不只是一个“看上去很美”的黑盒它的输出是可以用经典声学反演的。如果你也刚开始学建议做一次这种自己能用理论公式检验的小实验它比跑通100个官方demo都更能帮你建立对工具的信任感。4. 从声场到射频信号做一个能出图的脉冲回波小实验4.1 散射体怎么摆、幅度怎么设理解了连续波声场之后下一个关键步骤是仿真“发射-反射-接收”的完整脉冲回波过程。Field_II里这一步基本由calc_scat或calc_scat_multi完成。官方样例里最常出现的散射体模型是均匀随机分布的点散射体回波幅度按散射强度设置。比如我要做一个直径5mm的低回声囊肿区域就可以这样做rng(2024); N 2000; x (rand(N,1)-0.5) * 24e-3; % x方向范围 ±12mm y (rand(N,1)-0.5) * 6e-3; % y方向范围 ±3mm z 25e-3 rand(N,1) * 30e-3; % z方向范围 25-55mm amp ones(N,1); % 在 (0, 0, 40mm) 附近生成无散射体的球形区域 center [0 0 40e-3]; radius 2.5e-3; inside ((x-center(1)).^2 (y-center(2)).^2 (z-center(3)).^2) radius^2; amp(inside) 0;散射体数量控制在2000个左右对演示足够计算量也能接受。如果数量少于几百图像上会看到明显的颗粒状噪声多于两三万个计算时间会让你怀疑人生。想要用Field_II模拟“均匀组织”就必须明白一个现实现实中无数个微小散射体在频率下返回的信号是相干叠加的点散射体数量足够多且间距小于分辨率时才能呈现漫散射斑纹speckle。用随机点只是取了一个统计等价的近似。4.2 发射、接收、延迟求和与B模式成像脉冲回波的发射部分和我前面连续波示例类似需要给阵元设置激励脉冲。接收部分需要考虑的是每个阵元接收到回波信号后要按焦点位置做时间延迟再叠加起来这就叫延迟求和波束成形。下面是核心代码思路省去中间辅助函数% 发射孔径和接收孔径这里用同一个阵列 xdc_excitation(Th_tx, excitation); xdc_focus(Th_tx, 0, [0 0 40e-3]); % 发射焦点在40mm % 接收时动态聚焦这里简化为逐点延迟 focus_z linspace(25e-3, 55e-3, 101); rf_lines zeros(length(focus_z), 1); for line_idx 1:21 % 21条扫描线 x_line (line_idx - 11) * 1e-3; % 计算所有散射体到每个阵元的往返时延叠加回波 [rf, start_t] calc_scat_multi(Th_tx, Th_rx, points, amp); % 添加接收延迟按焦点补偿 % ... rf_lines(:, line_idx) env(rf); % 取包络 end img db(rf_lines); % 压缩成灰度 imagesc(img);这里我没有列出所有接收延迟的代码因为太占篇幅但你要抓住的重点是Field_II输出的rf就是最原始的射频线数据它不带波束成形因为波束成形本来就是超声算法层面的东西Field_II不做这件事。这也是它作为“仿真器”让人舒服的地方——你可以自由替换自己的波束成形算法公平对比不同策略的效果。实际做的时候发射焦点一般取一个固定深度接收焦点逐点变化这样才能保证全图像都有较好的侧向分辨率。否则你在一个深度上清晰远一点就开始模糊。4.3 仿真结果怎么解读当第一张B模式图从Field_II的数据里渲染出来时我会建议你至少看三样东西第一囊肿边界是否清晰。如果周围斑纹正常而囊肿内部明显黑下去说明散射体模型和波束成形基本生效了。第二点目标回波的旁瓣结构。如果你在散射体里放了一两个强反射点幅度远高于其他散射体B模式图上它们应该呈现出和系统点扩散函数一致的小亮斑。如果亮斑旁边出现明显的横向长条纹那就是旁瓣表现偏弱或阵元数不够。第三轴向与侧向分辨率的大概对比。通常轴向分辨率高于侧向尤其是阵元数少的时候。还有一种常见的陷阱是直接用0到1的线性幅度显示图像结果整幅图黑乎乎一片这是因为射频信号动态范围非常大。做法上先取包络再转成dB显示然后设置你想要的动态范围比如imagesc(20*log10(env/max(env))); caxis([-60 0]);。5. 避开那些让我重跑一整夜的坑5.1 内存爆炸从“一个数组就崩了”说起我刚开始做阵元-散射体批量计算时MATLAB直接提示内存不足。原因是Field_II里的calc_scat_multi会为每个阵元、每个散射体、每个时间点维护中间变量复杂度大约是阵元数 × 散射体数 × 时间样本数。举个例子64阵元、5000个散射体、5000个时间采样点如果内部矩阵按复数存储一个8字节复数那么一次运算需要的内存规模大概是64 × 5000 × 5000 × 16 byte ≈ 25.6 GB这个数字对大多数个人电脑来说是灾难。所以后来的做法是把大矩阵拆散小批量处理比如一次只计算500个散射体循环10次内存占用就只有原来的十分之一代价是时间略多。另一个非常实用的内存优化技巧是在三维网格里算连续波声场时别一次性申请全覆盖的meshgrid。我以前习惯直接[X,Y,Z] meshgrid(...)再算每一点的声场结果一个略精细的网格就会让内存爆掉。后来改成逐层切片计算比如固定一个深度面算完一片保存一片内存就稳定在几个GB内。5.2 fcut衰减参数三种取值对应的三种物理学场景Field_II的衰减设置是很多新人最容易忽略的细节。set_field(att, att)这个函数参数取值不同含义完全不同设置方式含义适用场景set_field(att, 0)不开启介质衰减验证声场理论模型、算法快速迭代set_field(att, 0.5)set_field(freq_att, 0.5)set_field(att_f0, 5e6)频率相关衰减衰减系数0.5 dB/cm/MHz模拟软组织中声波随深度和频率的损耗set_field(att, 1)... 配合其他方式需要更复杂衰减模型时自定义经验公式我犯过的错误是只设置了att不设置freq_att和att_f0导致衰减跟频率无关和实际物理差很远。你需要明确一点超声在组织里的衰减近似与频率成正比所以必须告诉Field_II你的中心频率基准是多少它才知道同样0.5dB/cm/MHz系数在5MHz和10MHz下分别该衰减多少。对于最开始的算法验证我强烈建议先关掉衰减。这样你看到的效果都是纯粹声场和波束造成的结果问题排查时变量更少。等仿真框架完全跑通再加上衰减模拟真实组织。5.3 频率混叠与网格间距做3D声场计算时另一个隐蔽问题是空间网格步长。简化的经验法则是网格步长要小于最小波长的一半即dx λ_min / 2 c / (2 * f_max)如果你的脉冲含有较高频率分量比如5MHz中心频率、但谐波分量到15MHz也有能量那么最小波长是1480/15e6 ≈ 0.0987mm网格步长至少要小于约0.05mm。否则声场分布会出现明显的栅瓣或波纹状伪影而且这个伪影不是物理效应而是采样不足造成的混叠。这里有个权衡网格越细需要计算的坐标点越多。我的建议是对于初始总体观察用中间波长的一半做步长如果你特别关注焦点区域的精细结构再在关键区域局部加密。Field_II的calc_h是按点计算离散点的所以局部加密非常方便不需要把整个空间一起加密。5.4 计算时间优化哪些近似值得做实际工作中计算时间往往比内存更让人烦躁。仿真一个稍微像样的孔径跑上几小时是很常见的。好在Field_II有几条便宜的优化路线一是从三维退到二维。对狭长阵元或线阵如果在某一个方向比如阵元高度方向声场变化很慢可以先把该方向尺寸设小甚至用二维近似观察主平面内的声场结构。这在前期设计验证中足够了。二是压缩时间窗口。Field_II内部会为每个空间点计算冲激响应时间采样点数量直接影响计算量。如果你只关心焦点前后一小段深度范围就通过计算起点和终点控制时间窗口而不是从0时刻一直算到无穷远。三是减少自定义阵列的子孔径单元数。Field_II把阵元表面进一步细分细分越多计算越重。对于普通线阵把一个阵元在宽度方向细分到8~16段、高度方向细分到1~4段通常精度就足够了。再往上加点数计算时间翻倍精度提升却很微弱。四是开启MATLAB并行池。Field_II的某些内部计算在多核环境下能获得直接收益。不过要注意并行池在多散射体循环里的加速并不稳定建议先在小规模数据上测一下确认收益后再上并行避免花了大半天开池子结果加速比只有1.2。6. 学习Field_II的高效路线以及下一步怎么走6.1 一个可以复制的“三阶段”路径我在实际摸索中把学习路线拆成了三个阶段现在回头看这样走比漫无目的地刷官方demo高效很多。第一阶段是“跑通闭环”。不要贪多用最简单的单阵元、连续波、几个观察点把field_init、设置参数、计算、画图、field_end这套流程跑顺。目标是能在半小时内从零生成一条轴向声压曲线。第二阶段是“拓展参数”。把阵列规模增加阵元数、焦点、频率都改一改观察声场图相对变化。这个阶段你要训练的是“直觉”什么参数会让焦点变浅、什么参数会让旁瓣抬升、什么参数会让主瓣变窄。第三阶段是“完整成像”。用一两个强反射点和一个囊肿模拟把B模式图做出来再做分辨率测量和旁瓣评估。做到这里你才算是真正把Field_II用成了超声系统仿真工具而不是“一个能画彩色图的MATLAB脚本”。6.2 做一次“自测小目标”检验掌握程度如果你也想验证自己是不是真的入门了我建议给自己定一个小目标用Field_II仿真一个64阵元1D线阵中心频率5MHz阵元间距0.4mm焦点深度30mm仿真三个直径为2/4/8mm的囊肿最后出一张B模式图并且量出它们各自的横向和纵向尺寸。如果测量值和真实尺寸偏差在合理范围内你对Field_II的配置参数、散射体建模、波束成形和结果解读就都已经过关了。做这个自测时一个容易被忽略的细节是测量分辨率时要考虑系统的点扩散函数会“胀大”目标轮廓尤其是直径小于分辨率的囊肿。小囊肿在图像上显示出来的直径会明显大于真实直径。建议先仿真一个点目标得到系统的-6dB宽度再对囊肿结果做去卷积或至少修正口径。6.3 把Field_II和k-Wave串起来一条更接近真实世界的路径Field_II在换能器阵列建模上的优势非常强但它处理不了复杂的组织非线性传播。我现在正在做的一件事是用Field_II计算换能器表面的压力分布和焦点条件再把这个压力分布作为k-Wave仿真的输入让k-Wave去模拟声波在真实组织模型中传播、产生非线性分量、被组织散射后的结果。这种“混合仿真pipeline”在学术论文里很常见实际跑下来有个明显的好处既能用Field_II准确建模探头又能享受k-Wave处理非均匀介质和非线性声学的能力。两者的接口就是在某个平面/曲面上传递声压时间序列这个传递过程需要注意网格重采样和时间对齐不然两个工具的噪声都会叠加进去。如果你以后打算做谐波成像、对比剂微泡、或者精确组织模型加热仿真非常推荐往这个方向走。你可以先读一读k-Wave官方的声源定义文档了解它能不能接收自定义输入时变压力源然后把Field_II的输出存成MATLAB数组直接喂进去。最后分享一个我自己摸索出来的习惯每次跑新仿真前先把理论公式能算的量都算一遍比如焦点深度、主瓣宽度估算、自然焦点位置再和Field_II的结果对比。做这个动作不是为了验证Field_II对不对——它大概率是对的——而是为了验证我自己的参数设置和我想表达的问题确实一致。很多时候仿真结果“不对”并不是因为工具坏了而是因为脑袋里想的问题和代码里设的参数本来就是两个东西。