Field II 超声成像仿真 1--得到Bmode图像

本文单纯记录学习内容，记录跑通思路~

平面波 M-mode 成像仿真

一.使用工具：

        USTB (UltraSound ToolBox, MATLAB)：【Git 仓库克隆：把整个 USTB 仓库 clone 下来。这是推荐方式，因为可以方便更新。】

在终端运行

git clone https://github.com/ustb/ustb.git

二.STAI_L11_probe_sim.m 记录

1. 逻辑：

这段代码的目的是使用 Field II 软件模拟一个超声合成发射孔径 (STA) 数据集，然后使用 USTB (Ultrasound Toolbox) 对其进行波束成形，最终生成一幅超声图像

      2.代码：

---

第一部分：初始化与设置

c0=1540;     % 声速 [m/s]
fs=100e6;    % 采样频率 [Hz]
dt=1/fs;     % 采样时间间隔 [s]

• Field II 初始化：

field_init(0);  %初始化
set_field('c',c0); %设置声速
set_field('fs',fs); %设置采样率
set_field('use_rectangles',1); % 元件用矩形模型（而不是椭圆）

---

第二部分：定义超声探头

• 定义探头：定义一个 128 阵元的线阵探头。

这里是在打造一把虚拟的超声探头。你精确地定义了它的物理规格：它有 128 个小小的听筒（阵元），每个听筒多宽 (element_width)，多高 (element_height)，听筒之间隔多远 (pitch)，以及它们之间的缝隙有多窄 (kerf)。中心频率 f0 决定了这把探头发射超声波的音调高低

probe = uff.linear_array();

f0 = 5.1333e+06;         % 中心频率 [Hz]
lambda = c0/f0;          %  计算波长 ≈ 0.3 mm

probe.element_height = 5e-3;  %阵元高
probe.pitch = 0.300e-3;  %阵元间距 pitch
kerf = 0.03e-03;  % 阵元之间的间隙 0.03 mm

probe.element_width = probe.pitch-kerf;  %阵元宽度 = pitch - kerf
lens_el = 20e-3;
probe.N = 128; 
pulse_duration = 2.5;    % 脉冲长度 [cycles]

---

第三部分：定义发射的超声脉冲

• 定义发射脉冲

pulse = uff.pulse(); %定义一个 USTB 的 pulse 对象
pulse.center_frequency = f0;
pulse.fractional_bandwidth = 0.65; %带宽 = 0.65 

• 模拟脉冲的产生

%然后生成冲激响应和激励信号

t0 = (-1/pulse.fractional_bandwidth/f0): dt : (1/pulse.fractional_bandwidth/f0);

impulse_response = gauspuls(t0, f0, pulse.fractional_bandwidth);  %gauspuls函数生成高斯包络的载波信号

impulse_response = impulse_response-mean(impulse_response); %直流分量去除

再生成方波激励，并和冲激响应卷积

te = (-pulse_duration/2/f0): dt : (pulse_duration/2/f0); %% 发射激励时间轴
excitation = square(2*pi*f0*te+pi/2); % 生成一个方波激励电信号

one_way_ir = conv(impulse_response,excitation);  %​​单向脉冲响应​
two_way_ir = conv(one_way_ir,impulse_response); %双向脉冲响应

lag = length(two_way_ir)/2 + 1; % 时间轴中心点

单向脉冲响应： 模拟单个换能器发射时的电-声转换过程

双向脉冲响应：超声波遇到目标散射回来，又被同一个探头接收。这个过程相当于脉冲经历了“发射-接收”两个过程，所以用发射脉冲和接收脉冲再做一次卷积来模拟。这被称为“双向脉冲响应”

% We display the pulse to check that the lag estimation is on place

% (and that the pulse is symmetric)

figure;

plot((0:(length(two_way_ir)-1))*dt -lag*dt,two_way_ir); hold on; grid on; axis tight

plot((0:(length(two_way_ir)-1))*dt -lag*dt,abs(hilbert(two_way_ir)),'r')

plot([0 0],[min(two_way_ir) max(two_way_ir)],'g');

legend('2-ways pulse','Envelope','Estimated lag');

title('2-ways impulse response Field II');

第四部分：设置 Field II 的孔径模型

• 创建发射和接收孔径

noSubAz = round(probe.element_width/(lambda/8)); % 将每个阵元在 azimuth 方向细分
noSubEl = round(probe.element_height/(lambda/8));% 将每个阵元在 elevation 方向细分

Th = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]);

Rh = xdc_linear_array (probe.N, probe.element_width, probe.element_height, kerf, noSubAz, noSubEl, [0 0 Inf]);

为了计算更精确，Field II 会把探头的一个物理阵元细分成很多个更小的点（子元素） 来进行模拟。这里根据波长（λ）来决定细分的程度（通常每个波长划分 8 个点以上）。Th (Transmit handle) 和 Rh (Receive handle) 就像是两个虚拟的探头，一个专门负责发射，一个专门负责接收。

xdc_linear_array → Field II 函数，定义线阵探头几何和子阵元划分

• 为孔径设置属性

xdc_excitation(Th, excitation); % 设置发射孔径的激励信号
xdc_impulse(Th, impulse_response); % 设置发射孔径的脉冲响应
xdc_impulse(Rh, impulse_response); % 设置接收孔径的脉冲响应
...

将我们之前定义好的激励电信号和脉冲响应“加载”到虚拟探头 Th 和 Rh 上。

% We also set the excitation, impulse response and baffle as below:

xdc_excitation (Th, excitation);

xdc_impulse (Th, impulse_response);

xdc_baffle(Th, 0); %0 = 无限大刚性挡板，防止背面辐射（只向前发声）

xdc_center_focus(Th,[0 0 0]); %设置发射阵元的“几何中心”对准 (0,0,0)。

xdc_impulse (Rh, impulse_response);  %接收也需要一样设置

xdc_baffle(Rh, 0);

xdc_center_focus(Rh,[0 0 0]);

第五部分：定义仿体并计算 STA 数据

•  定义仿体 (Phantom)

sca = [0 0 20e-3]; % 散射点坐标 [x, y, z] = [0, 0, 20mm]
amp = 1;            % 散射点强度

这里创造了一个非常简单的仿体——在探头正前方 20mm 深度处有一个孤零零的、强度为 1 的点状散射物（可以想象成一个微小的玻璃珠）。

• 计算 STA 数据集 (核心循环)

for n=1:probe.N % 遍历每一个阵元
% 发射孔径：只激活第 n 个阵元
xdc_apodization(Th, 0, [zeros(1,n-1) 1 zeros(1,probe.N-n)]);  %第n个为1

% 接收孔径：所有 128 个阵元都接收
xdc_apodization(Rh, 0, ones(1,probe.N));  %第 n 个发射 -> 所有阵元接收 的数据

% 计算散射回声：计算第 n 个阵元发射，所有阵元接收到的信号
[v,t] = calc_scat_multi(Th, Rh, sca, amp);

% 存储数据：将数据存入一个大矩阵 STA 中
STA(1:size(v,1),:,n) = v;

% 记录发射序列信息：记录第 n 次发射的物理参数（如哪个位置发的、延迟等）
seq(n) = uff.wave();
...
end

详细代码

%% Output data
% 
% We define the variables to store our output data
cropat=round(1.1*2*sqrt((max(sca(:,1))-min(probe.x))^2+max(sca(:,3))^2)/c0/dt);   % maximum time sample, samples after this will be dumped
STA=zeros(cropat,probe.N,probe.N);    % impulse response channel data%% Compute STA signals
% 
% Now, we finally reach the stage where we generate a STA (Synthetic
% Transmit Aperture) dataset with the help of Field II.disp('Field II: Computing STA dataset');
wb = waitbar(0, 'Field II: Computing STA dataset');
for n=1:probe.Nwaitbar(n/probe.N, wb);% transmit aperturexdc_apodization(Th, 0, [zeros(1,n-1) 1 zeros(1,probe.N-n)]);xdc_focus_times(Th, 0, zeros(1,probe.N));% receive aperture    xdc_apodization(Rh, 0, ones(1,probe.N));xdc_focus_times(Rh, 0, zeros(1,probe.N));% do calculation[v,t]=calc_scat_multi(Th, Rh, sca, amp);% build the datasetSTA(1:size(v,1),:,n)=v;% Sequence generation  seq(n)=uff.wave();seq(n).probe=probe;seq(n).source.xyz=[probe.x(n) probe.y(n) probe.z(n)];seq(n).origin.x = probe.x(n);seq(n).sound_speed=c0;seq(n).delay = probe.r(n)/c0-lag*dt+t; % t0 and center of pulse compensation
end
close(wb);

第六部分：组织数据并波束成形

• 创建 USTB 信道数据对象

% In this part of the code, we creat a uff data structure to specifically

% store the captured ultrasound channel data.

channel_data = uff.channel_data();

channel_data.sampling_frequency = fs;

channel_data.sound_speed = c0;

channel_data.initial_time = 0;

channel_data.pulse = pulse;

channel_data.probe = probe;

channel_data.sequence = seq;

channel_data.data = STA./max(STA(:));

将原始的通道数据 channel_data 以 UFF (Ultrasound File Format) 的标准格式保存下来。这样以后就可以直接加载这个文件进行分析或使用不同的参数进行波束成形，而无需重新运行耗时的 Field II 仿真

channel_data 就像是一个 超声“原始通道数据”的容器，既保存数据，又保存所有元信息（meta-data）

• 定义扫描区域

%% SCAN

sca=uff.linear_scan('x_axis',linspace(-4e-3,4e-3,256).', 'z_axis', linspace(16e-3,24e-3,256).');

定义你最终想生成图像的区域和分辨率。这里定义了一个从深度 16mm 到 24mm，宽度从 -4mm 到 4mm 的矩形区域，并将其划分为 256x256 个像素点

• 设置波束成形管道 (Pipeline)

pipe = pipeline(); % 创建处理管道
pipe.channel_data = channel_data; % 放入原始数据
pipe.scan = sca;                  % 放入扫描区域定义

% 设置参数：使用 Boxcar（矩形窗）函数进行变迹，F数为 1.7
pipe.receive_apodization.window = uff.window.boxcar;
pipe.receive_apodization.f_number = 1.7;
...

• 设置接收和发射孔径加权 (Apodization)

pipe.receive_apodization.window=uff.window.boxcar; %矩形窗，相当于均匀加权
pipe.receive_apodization.f_number=1.7;
pipe.transmit_apodization.window=uff.window.boxcar;
pipe.transmit_apodization.f_number=1.7;

变迹 (Apodization)：给不同阵元的信号赋予不同的权重（通常中间权重大，两边权重小），类似于加一个窗函数，可以抑制旁瓣，改善图像质量 ，

F数：类似于相机的光圈，控制接收波束的聚焦深度。F数越小，波束越窄，分辨率越高，但噪声也可能更多

• 执行波束成形并显示

b_data=pipe.go({midprocess.das() postprocess.coherent_compounding()});

pipe.go(): 启动处理管道，开始执行波束成形。

midprocess.das(): 使用延迟叠加 (Delay-And-Sum) 算法，这是最基础的波束成形算法

postprocess.coherent_compounding(): 相干复合。因为我们有 128 次发射的数据，这个步骤会将所有 128 幅低质量的单发射图像融合成一幅高质量的最终图像

• 显示成像结果

b_data.plot()

第七部分：保存数据

• 保存 UFF 文件

channel_data.write([ustb_path(),'/data/FieldII_PSF_simulation_v2.uff'],'channel_data');

---

仿真结果如下图：

总结

这段代码完成了一个完整的超声成像仿真流水线：

1. “造器械”：定义虚拟探头、脉冲等硬件参数。

2. “造目标”：定义一个简单的点状仿体。

3. “数据采集”：使用 Field II 模拟最原始的 STA 数据采集过程（循环发射-接收）。

4. “数据打包”：将原始数据整理成 USTB 的标准格式。

5. “图像处理”：使用 USTB 的波束成形算法（DAS + 复合）将原始数据重建为图像。

6. “出片”：显示并保存最终图像和原始数据。

最终，你应该会看到一幅图像，中心有一个明亮的点，周围有一些固有的旁瓣（sidelobe）噪声。这个点就是仿体中那个散射点的像，图像的质量（点的胖瘦、旁瓣的高低）体现了你所定义的探头、脉冲和波束成形参数的性能。