简介这是一份围绕Field II软件开展的PSF点散射成像仿真实例包面向光学系统仿真、图像处理及医学超声成像领域的研究人员和工程师既可用于教学演示也适合成像算法的前期验证与参数探索。资源以20个散射点的成像仿真为主线完整覆盖参数设定、点源图像生成、结果分析等环节能够帮助使用者从实际数据中理解点扩展函数PSF如何影响成像分辨率以及散射现象对成像质量的作用。压缩包共16个文件全部为MATLAB脚本.m整体大小仅4KB轻量且模块清晰脚本分工明确涵盖主仿真流程、数值孔径与波前计算、矩阵运算辅助、图像合成等功能模块便于按需调用和二次开发。目前已有282人学习下载。通过运行这些脚本可直观对比不同点间距与位置下的成像特征快速掌握Field II建模的关键步骤同时为光学设计优化、散射分析以及PSF测量提供可复用的代码基础。1. 这个 tar.gz 里的 PSF是超声成像质量的答案拿到psf_example.tar.gz的时候大多数人其实不是想研究压缩包而是想在一套超声仿真流程里确认一个基本问题这套参数下系统到底能把一个点散射体拍成什么样。Field_II 是超声仿真领域被引用最多的工具箱核心是空间冲激响应模型PSFPoint Spread Function正是这个模型下点散射目标的成像结果。一个理想点经过发射、传播、散射、接收、波束合成之后变成带有主瓣和旁瓣的二维斑。斑越窄侧向和轴向分辨率越好旁瓣越低动态范围越干净。点散射是这一切的输入散射成像则是把成千上万个点散射体叠加后的宏观表现。这篇文章从解压 tar.gz 开始走到 Field_II 初始化、单点散射体仿真、PSF 成像最后落到量化分辨率特别适合刚接触 Field_II 又想快速跑通 PSF 流程的人。2. Linux 解压 tar.gz 与 Field_II 初始化从压缩包到能跑仿真2.1 用 tar 命令解压 psf_example.tar.gz 的最小步骤很多现场错误跟仿真本身没关系卡在第一步解压。最常见的是用户把 tar.gz 下载到了~/Downloads却在 MATLAB 的工作目录里执行解压命令然后得到tar: 没有那个文件或目录。先确认文件位置、再解压是最不容易出错的顺序。cd ~/Downloads ls -lh psf_example.tar.gz file psf_example.tar.gz mkdir -p ~/ultrasound tar -xzf psf_example.tar.gz -C ~/ultrasound第一行切换目录第二行用ls -lh确认文件存在且大小不为 0第三行用file检查文件类型这一步能直接分辨出“真 tar.gz”和“被浏览器改名成 tar.gz 的 HTML 文件”。第四行用-C指定解压目标目录避免把内容散在Downloads里。tar -xzf的含义是x解压、z走 gzip 解压、f指定文件名。如果file命令输出显示gzip compressed data说明压缩包本身没问题。VSCode 用户常在集成终端里操作这里容易踩一个坑如果用了 Remote-SSHtar.gz 可能在本机而终端在远端tar当然找不到文件。先pwd看看终端所在的机器和目录再决定是上传还是本地操作。解压之后目录里一般会有psf_example.m和一个Field_II子目录或若干数据文件以下都以这个目录为工作区。2.2 Field_II 在 MATLAB 与 Octave 中的安装动作Field_II 不是 pip install 就能装好的东西它是一组 .m 文件外加需编译的 mex 内核。在 MATLAB 中关键是让field_init能被找到。addpath(/home/user/ultrasound/Field_II); addpath(/home/user/ultrasound/psf_example); savepath; cd(/home/user/ultrasound/psf_example); field_init;这段命令把 Field_II 的主目录和示例目录加入 MATLAB 搜索路径savepath让路径配置在下次启动时仍然有效。cd之后调用field_init它负责加载 mex 内核、设置默认参数结束时应当用field_end释放资源。如果不执行field_init直接调用calc_scat会得到类似Undefined function calc_scat的错误这不代表代码写错只是初始化被跳过了。Octave 用户需要额外注意一件事如果下载的是 Windows 版预编译 mex在 Linux 上直接field_init会报.mexw64无法加载。这时候要在 Field_II 的目录里重新执行mex编译常见做法是看README里给定的 MEX 源文件名逐个编译。不要指望官方一次提供所有平台的二进制这是 Field_II 使用中非常普遍的认知差。2.3 解压没有那个文件或目录时的三种查法遇到tar: 没有那个文件或目录按顺序排查三个点第一是当前路径第二是文件名拼写第三是压缩包内容结构。pwd ls -la *.tar.gz tar -tvzf psf_example.tar.gz | head -20第一条看当前目录第二条用通配符列出所有 tar.gz 文件第三条不解压直接查看包内目录结构。如果第三条输出里所有路径都带一层psf_example/前缀解压时用tar -xzf自然会生成这个目录不需要手动 mkdir。如果输出里直接是散落的.m文件解压前先建目录再进目录解压更干净。一个很隐蔽的路径问题是 MATLAB 的cd和系统 shell 的cd不共享。在 VSCode 终端里解压成功了回到 MATLAB 执行psf_example依然报文件不存在多半是 MATLAB 当前目录还在别处。用cd(/home/user/ultrasound/psf_example)显式切换比pwd之后猜路径更直接。3. 点散射体 Field_II 仿真最小流程从空间冲激响应到 PSF 成像3.1 双程空间冲激响应为什么 PSF 是超声成像的固有属性Field_II 的理论基础是 Tupholme 与 Stepanishen 的空间冲激响应方法。换能器面上每个微小面元对空间中某个点的声场贡献可以写成该点到面元的距离相关的随时间变化的响应整个阵元则是面元响应的积分。发射时这个积分形成入射声场接收时散射点产生的回波再次经过接收阵元的空间冲激响应形成一个双程系统。假设空间中只有一个位于r0的点散射体其散射强度用幅度amp表示那么接收端电压信号的频域表达是V(ω) H_t(ω, r0) · H_r(ω, r0) · amp(ω)其中H_t和H_r分别是发射和接收的空间冲激响应传递函数。这个式子说明点散射体回波里同时编码了发射孔径和接收孔径的衍射特征而 PSF 就是这个双程响应的空间分布。Field_II 的calc_scat函数做的正是这件事对给定孔径、焦点、散射体位置计算时域回波。理解这一点对后面调参数很关键——孔径决定 PSF 旁瓣方向焦点决定主瓣汇聚位置点散射体只是探测这些系统特性的探针。3.2 定义线阵、聚焦与单点散射体的可运行代码下面这段代码是 PSF 仿真的最小骨架直接在解压出的工作区里新建psf_min.m运行即可。% psf_min.m -- 用 Field_II 计算单点散射体的二维 PSF field_init; set_field(c, 1540); % 声速 1540 m/s对应软组织典型值 set_field(fs, 100e6); % 采样率 100 MHz远高于奈奎斯特要求 set_field(use_rectangles, 1); % 用矩形近似阵元计算量显著下降 set_field(threads, 4); % 4 线程并行减少多散射体计算时间 f0 3e6; % 中心频率 3 MHz lambda 1540 / f0; % 波长约 0.513 mm width lambda; % 阵元宽度取一个波长 height 5e-3; % 阵元高度 5 mm kerf 0.05e-3; % 阵元间切割缝 0.05 mm N 64; % 阵元数 64 focus [0; 0; 60e-3]; % 焦点位置 (0, 0, 60 mm) emit xdc_linear_array(N, width, height, kerf); receive xdc_linear_array(N, width, height, kerf); emit xdc_focus(emit, 0, focus); receive xdc_focus(receive, 0, focus); apod hamming(N); emit xdc_apodization(emit, 0, apod); receive xdc_apodization(receive, 0, apod); % 以焦点为中心横向 ±3 mm、轴向 58~62 mm 的成像区域 Nz 81; Nx 41; z_axis linspace(58e-3, 62e-3, Nz); x_axis linspace(-3e-3, 3e-3, Nx); psf zeros(Nx, Nz); for ix 1:Nx for iz 1:Nz scat [x_axis(ix); 0; z_axis(iz)]; [rf, ~] calc_scat(emit, receive, scat, 1); env abs(hilbert(rf)); % 解析信号取包络 psf(ix, iz) max(env); % 取该位置回波的包络峰值 end end field_end;这段代码里最需要解释的是双循环的意义。calc_scat每次计算一个点散射体的完整回波理论上把散射体放在焦点位置时双程系统对焦点的响应最大放在偏离焦点的位置时回波峰值逐渐变小。把散射体在空间网格上逐个移动并记录峰值就得到了这个系统对整个空间分布的响应也就是二维 PSF。每一次calc_scat都包含发射、传播、散射、接收全过程所以得到的是双程 PSF而不是单纯发射声场。xdc_linear_array的三个几何参数width、height、kerf单位都是米。width取一个波长会让阵元方向性图比较宽适合作为基线kerf过大会产生栅瓣通常控制在 0.1 个波长以内。hamming窗抑制旁瓣代价是主瓣略微展宽这是 PSF 仿真里绕不开的权衡点。3.3 呈现 PSF 成像结果包络、动态范围与坐标轴计算完成之后psf矩阵里保存的是 41 行乘 81 列的空间响应幅度。直接imagesc会有两个问题幅度动态范围大导致低旁瓣不可见坐标轴单位不直观。常用做法是对峰值归一化后取 dB再限制显示动态范围。figure; psf_db 20 * log10(psf / max(psf(:)) eps); imagesc(z_axis * 1e3, x_axis * 1e3, psf_db, [-40 0]); colormap gray; axis image; xlabel(轴向距离 (mm)); ylabel(侧向距离 (mm)); title(Field_II 双程 PSF);psf / max(psf(:))将峰值归一化到 0 dBeps是为了防止零幅度位置取对数得到负无穷。imagesc的显示范围设为[-40 0]意思是只显示主瓣以下 40 dB 以内的响应低于 -40 dB 的旁瓣统一显示为黑色。这样能在同一张图里同时观察主瓣宽度、近场旁瓣和离轴响应。从图中能直观看到三个特征焦点中心处亮度最高沿轴向的主瓣宽度明显小于沿侧向的主瓣宽度这是典型的超声成像现象轴向分辨率由发射脉冲带宽决定侧向分辨率由孔径的衍射极限决定主瓣两侧存在周期性旁瓣这是有限孔径截断产生的衍射环如果设置里kerf过大还会在远离主瓣的位置出现与栅瓣对应的亮带。4. 影响 PSF 成像形态的 4 个参数频率、孔径、F 数与采样率4.1 中心频率与波长分辨率刻度的基准在傅里叶成像的意义上超声 PSF 的分辨率刻度和频率严格绑定。中心频率f0决定波长lambda c / f0而波长是孔径尺寸、阵元宽度、聚焦深度的公共标尺。理论上一个圆形聚焦孔径的双程侧向分辨率约等于lambda * F_number轴向分辨率则由带宽反比决定。把f0从 3 MHz 提高到 6 MHz波长减半同样孔径下 PSF 的侧向主瓣宽度近似减半。但频率不是白给的。3 MHz 在软组织中的衰减约 0.5 dB/cm/MHz穿透深度和频率成反比。这个矛盾是 PSF 仿真里最不该忽略的仿真里改一个频率参数非常容易但真实系统的信噪比和成像深度边界也随之改变。Field_II 本身不模拟非线性衰减默认介质是均匀无损的所以仿真 PSF 更接近“系统响应”而非“真实图像”。4.2 孔径与 F 数决定焦点深度处的声场形态F 数定义为焦距与孔径宽度的比值。上一节的例子中64 个阵元、阵元宽度一个波长孔径宽度约 32.8 mm焦距 60 mmF 数约 1.8这属于较强的聚焦。F 数越小焦点处声束越窄但焦点附近声场变化越剧烈焦深越短。把阵元数从 64 减到 32孔径减半F 数变为 3.6PSF 主瓣明显变宽旁瓣结构也随之改变。4.3 采样率与带宽时间域量化对旁瓣的影响fs 100 MHz对 3 MHz 中心频率意味着每周期约 33 个采样点这远高于奈奎斯特率主要目的是保证延时叠加时的延迟精度。Field_II 计算得到的rf信号在时间上被离散化如果fs太低波形的时间量化误差会直接转化为 PSF 的旁瓣抬升和主瓣位置偏移。参数变量对 PSF 的影响常用起点中心频率f0决定波长主瓣宽度随频率升高而收窄2.5 ~ 7.5 MHz阵元数N决定孔径宽度阵元越多侧向主瓣越窄64 ~ 128F 数z / (N * width)F#越小焦点处声束越窄、焦深越短1.5 ~ 4采样率fs时间量化误差过低会抬高旁瓣8 * f0以上幅度加权apod抑制旁瓣代价是主瓣略微展宽Hamming / Hann4.4 参数对比实验的设计改一个、锁住其他调参实验最容易犯的错误是一次改多个参数导致无法归因。常见做法是把上一章的代码包成一个函数输入阵元数和加权窗输出 PSF 矩阵。function psf compute_psf(N, apod_type) % 复用第 3 章初始化代码只把 N 和 apod 变成入参 field_init; set_field(c, 1540); set_field(fs, 100e6); set_field(use_rectangles, 1); f0 3e6; lambda 1540 / f0; focus [0; 0; 60e-3]; emit xdc_linear_array(N, lambda, 5e-3, 0.05e-3); receive xdc_linear_array(N, lambda, 5e-3, 0.05e-3); emit xdc_focus(emit, 0, focus); receive xdc_focus(receive, 0, focus); switch apod_type case hamming apod hamming(N); case hanning apod hanning(N); case boxcar apod ones(1, N); end emit xdc_apodization(emit, 0, apod); receive xdc_apodization(receive, 0, apod); % 只计算焦点附近横向 41 点、轴向 81 点 [psf, x_axis, z_axis] scan_psf(emit, receive, focus); field_end; end这里的核心是switch分支保证每次只改加权类型其他几何和焦点参数全部固定。实际跑对比时分别在 64 阵元下用 boxcar 和 hamming 各算一次观察 -6 dB 主瓣宽度和第一旁瓣高度。boxcar 的第一旁瓣大约在 -13 dB 附近hamming 能把旁瓣压到 -40 dB 以下但主瓣宽度增加约 30%。这个数字不是理论空谈而是能在仿真结果上直接量出来的。5. 从单点散射到散射成像RF 数据、散斑与相位恢复的区别5.1 用随机点散射体生成散射成像仿真数据单点散射体只能给出系统响应真实的散射成像对象是大量随机分布的点散射体。把上一章的发射、接收孔径定义保持不变把单个散射体替换成数百个随机位置点就可以模拟组织散射。% 生成 500 个随机点散射体分布在 10 mm x 10 mm 区域内 N_scat 500; x_pos (rand(1, N_scat) - 0.5) * 10e-3; z_pos 55e-3 (rand(1, N_scat) - 0.5) * 10e-3; scat_pos [x_pos; zeros(1, N_scat); z_pos]; scat_amp randn(1, N_scat); % 幅度随机模拟不同散射强度 [rf, t0] calc_scat_multi(emit, receive, scat_pos, scat_amp); rf_env abs(hilbert(rf)); figure; plot(t0 * 1e6, 20 * log10(rf_env / max(rf_env) eps)); xlabel(时间 (μs)); ylabel(幅度 (dB));calc_scat_multi和循环调用calc_scat的区别在于前者把多个散射体的回波在时间上叠加后一次性输出避开了大量 MATLAB 循环开销这也是处理散射成像生成 RF 数据时的首选方式。randn幅度符合高斯分布对应 Rayleigh 散射假设这是组织散射仿真最常用的统计模型。生成的 RF 信号不再是单一的 PSF 响应而是 PSF 与多个散射体的卷积叠加包络信号呈现颗粒状的散斑纹理。5.2 散射成像读图时的概念边界PSF 仿真不涉及相位恢复看到散射成像和相位恢复这两个词一起出现时容易把问题想偏。Field_II 仿真的是相干超声成像里的正向过程已知散射体分布和系统 PSF计算接收 RF 信号。相位恢复通常出现在无透镜成像、相干衍射成像的场景描述的是从强度测量中恢复相位信息的问题属于逆问题的范畴。对这个标题里的 PSF 仿真来说接收孔径处的采集本身就保留了相位信息不需要做相位恢复。散斑统计特性反而更值得关注。当散射体数量足够大且独立分布时包络幅度服从 Rayleigh 分布均值与标准差之比约 1.91。这个比值可以当作检查仿真是否正常的指标如果计算出的比值明显偏离说明散射体数量不足或分布不再均匀。实际操作中把 500 个散射体增加到 2000 个散斑纹理更稳定但计算时间按散射体数线性增长这时候set_field(threads, 4)和calc_scat_multi的加速作用就很明显。5.3 从示例回归数据文件校验 psf_example 的输出Field_II 示例psf_example在跑完后通常会保存 RF 数据或 PSF 矩阵供后续处理复用。解压出来的psf_example.tar.gz里如果带有.mat数据文件可以对照自己的计算结果检查偏差。load(psf_example_data.mat); size(psf_ref) max(psf_ref(:))这里的检查逻辑是尺寸是否与自己计算的 41 x 81 网格一致峰值是否落在焦点位置。如果尺寸对不上大概率是网格范围或采样间距设置不同如果峰值位置偏移检查焦点定义和坐标轴方向。Field_II 的坐标约定是 z 轴沿声束传播方向x 轴沿阵元排列方向这一点在对比数据时最容易绕晕。6. 收尾技巧从 2D PSF 图量化分辨率与旁瓣6.1 画 -6 dB 等值线观察主瓣形状成像系统的分辨率通常以 -6 dB 主瓣宽度为定义对应幅度降为一半处的宽度。直接用imagesc观察很难精确读出边界等值线图是更客观的方式。psf_norm psf / max(psf(:)); db_psf 20 * log10(psf_norm eps); figure; contour(x_axis * 1e3, z_axis * 1e3, db_psf, [-6 -10 -20], LineWidth, 1.5); xlabel(轴向位置 (mm)); ylabel(侧向位置 (mm)); grid on;contour的第三个参数是等值线层级[-6 -10 -20]分别对应主瓣半功率边界、-10 dB 边界和旁瓣包络。db_psf的转置是为了匹配contour对矩阵维度的要求第一个向量对应矩阵的行坐标第二个向量对应列坐标。画出来后-6 dB 等值线的横向宽度就是侧向分辨率纵向宽度就是轴向分辨率。6.2 数值提取 FWHM做参数横向对比等值线适合观察数值量化需要直接对剖面线操作。取焦点所在位置的行和列分别提取侧向和轴向剖面计算超过 -6 dB 的宽度。[~, iz] min(abs(z_axis - 60e-3)); lateral_profile db_psf(:, iz); above lateral_profile -6; idx find(above); fwhm_lateral (numel(idx) - 1) * (x_axis(2) - x_axis(1)); [~, ix] min(abs(x_axis - 0)); axial_profile db_psf(ix, :); above_axial axial_profile -6; idx_axial find(above_axial); fwhm_axial (numel(idx_axial) - 1) * (z_axis(2) - z_axis(1)); fprintf(侧向 FWHM %.3f mm轴向 FWHM %.3f mm\n, ... fwhm_lateral * 1e3, fwhm_axial * 1e3);find(above)返回所有超过 -6 dB 的索引numel(idx) - 1是这段连续区域的采样间隔数乘上x_axis或z_axis的均匀间距得到实际物理宽度。这个方法对单峰 PSF 可靠如果旁瓣也超过了 -6 dB会导致宽度偏大这种情况下先缩小成像区域范围再用find找最大连续段。侧向 FWHM 在 3 MHz、F 数 1.8 的条件下通常在 0.8~1.2 mm 范围约两倍波长轴向 FWHM 则由脉冲持续时间决定一般在 0.5 mm 附近这个数值可以直接作为判断仿真参数是否合理的快速基准。本文还有配套的精品资源点击获取