超声平面波成像仿真（ultrasound plane_wave imaging simulation）

超声平面波复合成像仿真（ultrasound plane_wave compound imaging simulation）

​ 与传统超声聚焦成像相比，平面波成像通过一次发射可以获取整个感兴趣区域图像信息，帧率>1000Hz以上，但由于无聚焦，图像对比度和分辨率降低。为了克服这一限制，可以通过发射N个不同角度的平面波，通过相干叠加得到完整图像，在一定的角度范围内，平面波数量越多，图像质量越好。

​

​ 图像来自Ultrafast imaging in biomedical ultrasound

​ 上图是传统聚焦成像与平面波复合成像示意图们可以明显看出平面波的优缺点。

​ 本文通过Filed II对线阵和凸阵进行平面波复合成像仿真，通过发射延时控制平面波的形成，对于发射延时的计算可以看上节超声成像发射声场仿真（Ultrasound Emit Field Simulation）。

​ 本文仅对波束合成进行仿真，不涉及额外处理，如发射变迹、接收变迹等。

Field II仿真包含四个部分，

  1. 参数配置
  2. AD数据获取
  3. 波束合成
  4. 图像显示

在线阵中会对上述部分进行说明，凸阵设置基本一样

一 线阵平面波复合成像（plane-wave compound of linear array）

​

波束合成中的延时由上图中的两部分构成发射延时$Txdelay$和接收延时$Rvdelay$，在仿真过程中需要减去数据采样起始时间$stime$，在实际工程中$stime$可以代表波束合成参数准备时间、声束在透镜中传播时间、校正发射脉冲产生时间等，具体根据工程实现的方式去计算。对于延时时间最终转化为：
$$
delay =Txdelay+Rvdelay-stime
$$
在Field II中$stime$由下面的函数返回

1
[scat, stime] = calc_scat_multi (Th1, Th2, points, amplitudes);

$$
Txdelay = ((x_ntan(\theta)+z_n)cos(\theta) = x_nsin(\theta)+z_ncos(\theta))/c
$$

$$
Rvdelay = (sqrt((x_n-x)^2+(z_n-z)^2))/c
$$

上述即延时计算公式，$c$为声速，$(x_n,y_n)$为扫描区域位置，$(x,z)$为阵元坐标位置。

Field II仿真代码：

  1. 参数设置

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    clear all
    clc
    close all
    set(0,'defaultfigurecolor','w')
    addpath('../TransmitField');

    wavetype = 'plane_wave' ;%plane_wave \diverging_wave\fcous_wave
    %初始化
    field_init(0)
    %探头参数
    probe.type = 'linear';
    c = 1540;
    ch =128;

    f0 = 7500000;
    N_elements = 128;
    width = 0.17e-3;
    height = 5e-3;
    kerf = 0.03e-3;
    pitch = width+kerf;
    focus = 20e-3;
    fs = 100e6;                    %不要设置太大
    Th = xdc_linear_array (N_elements, width, height, kerf, 1, 10,[0,0,focus]);
    Rh = xdc_linear_array (N_elements, width, height, kerf, 1, 10,[0,0,focus]);

    %设置2个周期高斯脉冲相应、1个周期激励脉冲
    dt  = 1/fs;
    t0 = (-1/f0): dt:(1/f0);
    impulse_response = gauspuls(t0, f0);
    impulse_response = impulse_response-mean(impulse_response);
    pulse_duration = 1;
    te = 0:dt:pulse_duration/f0;
    excitation = square(2*pi*f0*te);
    %设置激励脉冲
    xdc_excitation (Th, excitation);
    %设置脉冲相应
    xdc_impulse (Th, impulse_response);
    xdc_impulse (Rh, impulse_response);

    %阵元参数

    x_ele = ([0:N_elements-1]-(N_elements-1)/2).*pitch;
    z_ele  = zeros(1,length(x_ele));
    probe.N_elements = N_elements;
    probe.pitch = pitch;


    %发射偏转角度
    steer_Num = 15;
    steer_angle = linspace(-15,15,steer_Num);
    %仿体设置
    point_position(1,:) = [0 0 10e-3];
    point_position(2,:) = [0 0 30e-3];
    point_position(3,:) = [-5e-3 0 30e-3];
    point_position(4,:) = [5e-3 0 30e-3];
    point_position(5,:) = [0 0 40e-3];
    point_amplitudes = ones(size(point_position,1),1);

  2. AD数据获取

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    % 获取AD数据
    ADInfo.data =zeros(round(2*60e-3/c/dt),N_elements,steer_Num);  
    for i = 1:steer_Num
        %实际工程中发射变迹通过发射波形去控制，以后再实施，这里不做发射变迹
        xdc_apodization(Th,0,ones(1,N_elements));
        %设置偏转发射延时
        xdc_center_focus(Th,[0 0 0]);
        emit_delay = plane_wave_tranmit_delay(probe,steer_angle(i)*pi/180,c);
        xdc_times_focus(Th,0,emit_delay);
        %接收变迹-使用矩形窗不做变迹
        xdc_apodization(Rh,0,ones(1,N_elements));
        %接收不聚焦
        xdc_center_focus(Rh,[0 0 0]);
        xdc_focus_times(Rh,0,zeros(1,N_elements))
        %获取AD数据
        [v,t]=calc_scat_multi(Th, Rh, point_position, point_amplitudes);
        ADInfo.data(1:size(v,1),:,i)=v;
        ADInfo.time(i) = t;
    end

  3. 波束合成

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    interp =0;%是否插值
    [pt,line,steer_Num] = size(ADInfo.data);
    %图像成像区域
    x = interp1(1:line,x_ele,linspace(1,line,line*2));
    z = [0:pt-1]*(c/fs/2);

    % x = interp1(1:line,x_ele,linspace(1,line,line*2));
    % z = linspace(5e-3,60e-3,256);

    [x,z] = meshgrid(x,z);
    scan_x = x(:);
    scan_z = z(:);
    %接收动态聚焦
    for i = 1:steer_Num
        %选取数据
        data = ADInfo.data(:,:,i);
        %数据插值
        steer = steer_angle(i)*pi/180;
        %发射距离，扫描区域点的位置相对于发射位置的距离
        t_dis = scan_x*sin(steer)+scan_z*cos(steer);
        %接收距离，扫描区域点的位置相对于阵元坐标距离
        r_dis = sqrt((scan_x-x_ele).^2+(scan_z(:)-z_ele).^2);
        %延时
        delay  = (t_dis+r_dis)/c-ADInfo.time(i);
        %是否插值，延时转化为点
     
        if ~interp
            pt_loc  = round(delay*fs);
            pt_loc(pt_loc<1) = 1;
            pt_loc(pt_loc>pt) = pt;
            pt_loc = pt_loc+[0:N_elements-1]*pt;
        
            for j =1:length(pt_loc)
                
                loc = pt_loc(j,:);
                RF(j,i) =  sum(data(loc));
               
            end
            
        else
            RF(:,i) =  data_interp(data,delay,fs,downsampleM,[size(x,1),size(x,2)]);
        end
        
    end
    %相干复合
    RFdata =  reshape(sum(RF,2),[size(x,1),size(x,2)]);

  4. 图像显示

1
2
3
4
5
6
7
envelope = abs(hilbert(RFdata));
envelope = db(envelope./max(envelope(:)));
figure;
imagesc(figure);
colormap gray
caxis(axis_handle,[-60 0]);

结果为：

从仿真效果上看15个角度旁瓣少，主瓣窄，提高平面波发射角度数量可以有效提升图像的分辨率和信噪比。

二 凸阵平面波复合成像（plane-wave compound of curve array）

​ 仿真流程与线阵一样，不同地方在于发射延时计算方式与线阵不一样，波束合成中延时的距离计算需要由极坐标系转化为笛卡尔坐标系，坐标系转化后应该如下图：

​

​ 按照线阵的流程配置，最终结果为：

明显可以看出15个角度的图像质量更好。

三 相控阵平面波复合成像（plane-wave compound of phase array）

相控阵发射设置与线阵完全一样，唯一不同的地方在于接收位置需要由极坐标转化为笛卡尔坐标系。坐标系转化后应该如下图：

​

按照线阵的流程配置，最终结果为：

总结

​ 通过Filed II仿真了线阵、凸阵、相控阵的平面波波束合成，从仿真结果上印证了增加复合角度数目，可以有效提升图像的分辨率和信噪比。

​ 在波束合成过程中如果假如发射、接收变迹，图像的质量可以进一步提升，变迹会在后续讨论，感兴趣小伙伴的可以尝试一下