前言
Field_II是丹麦超声专家J. A. Jensen等利用声学原理设计的一个超声系统。它可以仿真超声探头所形成的声场和超声图像等。虽然其在Matlab上进行仿真，但是它的建模计算主要采用封装的C语言来实现，这部分代码并没有开源，matlab代码只是用来配置参数和调用函数。用户手册只是对命令流进行简要介绍，缺少必要的解释或案例，加之网上的教程很少，因此该软件对我个人而言简直是晦涩难懂。本贴结合了网上零星的案例、本人购买的教程以及自己的学习心得，一是为了记录自己的学习过程，二是为了分享交流和讨论。如果把官网的用户手册比作是心法，希望本贴能够成为你学习超声仿真的招式，祝你在超声领域有所成就。
正文
0.参数配置
首先是所有仿真模拟的第一步，对参数进行配置。
c = 1540; % S声速
f0 = 2.5e6; % 换能器中心频率 [Hz]
fs = 100e6; % 采样频率 [Hz]
lambda = c/f0; % 波长
element_height = 13/1000; % 阵元高度 [m] (垂直方向)
element_width = 18.5/1000; % 阵元宽度 [m] (方位方向)
kerf = 0; %阵元间隔
focus = [0 0 60]/1000; % Fixed emitter focal point [m] (与单阵元换能器无关)
_elements = 1; % 阵列中物理阵元的数量
N_sub_x = 1; % x方向阵元的细分，即x方向子阵元数量
N_sub_y = 1; % y方向阵元的细分，即y方向子阵元数量
1 设置发射孔径
在Field II中，有关换能器的命令通常是以xdc_开头，这里使用xdc_linear_array命令创建一个线阵换能器。其中的focus参数一直令我费解，暂时按下不表，有理解的同学欢迎不吝赐教，谢谢。
emit_aperture = xdc_linear_array (N_elements, element_width, element_height, kerf, N_sub_x, N_sub_y, focus);
图1 用户手册关于xdc_linear_array的注释
图2 16阵元矩形线阵换能器
接下来，将对不同参数下的换能器阵列进行配置和绘图显示，以期能够直观的理解发射孔径的创建。
1.1 单一阵元换能器
定义单一阵元换能器，阵元个数为1，x和y方向上的子阵个数也设置为1：N_elements = 1, N_sub_x=1; N_sub_y=1。
用show_xdc命令绘制换能器，结果如图所示。
图3 单阵元换能器，N_sub_x=1，N_sub_y=1
需要注意的是，N_sub_x和N_sub_y是对阵元的进一步细分，使计算结果更加精确，并不影响仿真计算的结果。如图4，按照N_sub_x=10，N_sub_y=10进行设置。
图4 单阵元换能器，N_sub_x=10，N_sub_y=10
1.2 创建阵列换能器
在这里，使用指令创建聚焦换能器阵列，命令的用户手册说明如图5所示，阵列换能器的部分参数可以参考图6，代码如下：
element_height = 13/1000; % Height of element [m] (elevation direction)
pitch = 0.290/1000; % Distance between element centers
kerf = 0.025/1000; % Width of fill material between the ceramic elements
element_width = pitch-kerf; % Element width [m] (azimuth direction)
Rfocus = 60/1000; % Elevation lens focus (or radius of curvature, ROC)
focus = [0 0 60]/1000; % Fixed emitter focal point [m] (irrelevant for single element transducer)
N_elements = 64; % Number of physical elements in array
N_sub_x = 1; % Element sub division in x-direction
N_sub_y = 2; % Element subdivision in y-direction
emit_aperture = xdc_focused_array (N_elements, element_width, element_height, kerf, Rfocus, N_sub_x, N_sub_y, focus);
图5 xdc_focused_array用户手册注释
图6 线阵探头示意图
其中，Rfocus参数定义了换能器的曲率半径，同样用show_xdc命令绘制创建的换能器，如图7所示。
图7 创建的仰角聚焦线阵换能器，N_elements = 64; N_sub_x = 1; N_sub_y = 2
2 设置换能器的激励和脉冲响应
在定义好发射孔径之后，接下来使用函数设置发射孔径的激励脉冲
ex_periods = 1.5;
t_ex=(0:1/fs:ex_periods/f0);
excitation=square(2*pi*f0*t_ex);
xdc_excitation (emit_aperture, excitation);之后，使用函数设置换能器发射孔径的脉冲响应
t_ir = -2/f0:1/fs:2/f0;
Bw = 0.6; %带宽设为0.6
impulse_response=gauspuls(t_ir,f0,Bw);
set_sampling(fs);
xdc_impulse (emit_aperture, impulse_response);注意，这里设置的是发射孔径的激励脉冲和发射孔径的脉冲响应，而二者的卷积是我们要计算和得到的空间脉冲响应Spatial impulse response，不知道我理解的对不对。绘制激励脉冲和脉冲响应如图8所示。
图8 激励脉冲
图9 脉冲响应
3 计算并绘制发射声场强度
Field II可以计算空间中某点、某线以及某平面上的空间脉冲响应和声场强度。使用的是calc_h和calc_hp函数，如图10和11所示。
图10 计算孔径空间脉冲响应的程序
图11 计算发射声场的程序
这里以聚焦阵列换能器为例，对其发射声场的声场强度进行绘制显示计算结果。
3.1 空间某点计算结果
绘制空间中（0，0，60mm）处的空间脉冲响应和声压，那么用到的代码如下：
points=[0,0,60/1000];
[h,start_time]=calc_h(emit_aperture,points);
figure
plot((start_time+(0:length(h)-1)/fs)*1000000,h)
xlabel('time [us]'); title('点（x=0mm,y=0mm,z=60mm）处的空间脉冲响应');
[hp,start_time]=calc_hp(emit_aperture,points);
figure
plot((start_time+(0:length(hp)-1)/fs)*1000000,hp)
xlabel('time [us]'); title('点（x=0mm,y=0mm,z=60mm）处的声压');
图12 （0，0，60）mm处的空间脉冲响应
图13 （0，0，60）mm处的声压
3.2 空间某线上的计算结果
选择(-20,0,60)mm到(20,0,60)mm处选择101个测量点形成一条线，计算聚焦换能器阵列在60mm深度处的空间脉冲响应和声场强度。代码如下：
xstart=-20/1000;xend=20/1000;
x=linspace(xstart,xend,100)';
y=zeros(length(x),1);
z=60/1000*ones(length(x),1);
points=[x,y,z];
[hp,start_time]=calc_hp(emit_aperture,points);
图14 仰角聚焦阵列换能器60mm深度位置测量示意图
图15 60mm深度位置的声场强度
比较图13和图15我们可以发现，声场强度都是时间的函数，其实这很好理解，因为声波是从孔径向外传播的，因此在空间中的某点或某线的声强，是从无到有的，当声波还未传播至此则声强为0，随着时间的推移，声波也慢慢传播过来了。细心的朋友可能会发现，无论是点还是线的声场强度hp的时间轴都不是从0开始的，而是从start_time开始的，并且记录到的声场强度的持续时间也是有限的，这其实就会带来一些问题，比如无法绘制在整个计算空间中声场强度随时间的变化。
3.3 平面上的场强
计算切面方向上的声场强度，其中的变化范围为-15mm至15mm，的变化范围为5mm至150mm。照葫芦画瓢，平面中的声场强度也比较好计算，代码如下：
measDepthStart = 5/1000; % Start depth along z-axis to place measurement points
measDepthEnd = 150/1000; % End depth along z-axis to place measurement points
xStart = -10/1000; % Start position of measurement points in x direction
xEnd = 10/1000; % End position of measurement points in x direction
Nmpx = 81;
Nmpz = 59;
mx = linspace(xStart,xEnd,Nmpx)';
my = zeros(Nmpx*Nmpz, 1);
mz = linspace(measDepthStart,measDepthEnd,Nmpz)';
[X,Z] = meshgrid(mx,mz);
measurement_points = [X(:),my,Z(:)];
[pressure_tx, startTime_tx] = calc_hp(emit_aperture, measurement_points);
其实计算平面场强并不难，当时困扰我的并不是计算，而是声场的绘制，一直得不到跟网上案例一致的结果，纠结了很久，后来发现其实我的计算没有问题，主要是绘制的方式不得其法，有两个要点：1是归一化，2是将声强用对数表示。图16得到的是归一化的计算结果，图17是每个深度点归一化的计算结果。
在这里再啰嗦几句，不同于点和线的声场强度的展示方式，平面声场强度里并没有时间轴，我个人推测对于发射声场更关注的是声场分布，而不是声场的变化。因此真实的场强，应该是随时间在空间中移动变化的。
图16 归一话的场强
图17 聚焦换能器在xz平面每个深度位置归一化的场强分布
