引言
超声成像是生物医学影像领域中的主要方法之一。与其它高分辨率的影像技术相比,超声成像存在图像质量差的缺点,有必要研究高对比度和空间分辨率的超声图像重建方法。因此,本文基于超声延迟叠加算法进行波束幅度的变迹加权提升成像质量,并在算法方面进行速度优化提升成像速度。
延迟叠加算法是超声图像重建中最常见的波束形成算法。由于其原理简单,可以很容易地应用到实时超声成像应用中。但使用延迟叠加算法重建的图像中通常会出现高旁瓣和强伪影,并且成像的速度不高,这是由于其波束成形原理造成的。
本文在算法方面使用逆向运算对程序运行速度进行了优化。通过实验分析验证,发现使用反向延迟叠加算法可以在保证成像质量不变的前提下,至少有1.25倍的速度提升。
同时,为了提高延迟叠加聚焦超声波束成像质量,引入了波达方向窗函数作为优化因子,并对比研究了使用传统固定加权窗函数的主瓣宽度和旁瓣等级。通过实验,验证发现使用波达方向窗函数可以比固定窗函数更有效地抑制旁瓣并减小主瓣宽度。
关键词:超声成像;延迟叠加;窗函数;反向延迟叠加
一、延迟叠加算法
延迟叠加算法(Delay and sum,DAS),核心在于延迟的计算。延迟叠加算法的原理是通过计算不同通道回波数据的延迟,然后对每个通道数据进行延迟补偿,实现数据时域对齐,最后将对齐后的数据叠加,使得波束具有指向性,在目标点处的能量达到最大。
基本的叠加公式表示如下:
对于换能器接收到的回波信号s,可以用一个三维数组s[j][e][k]表示,其中的字母分别表示第e次发射、第k个阵元接收到的第j个信号。算法的输入是三维回波信号,输出是叠加后的射频信号RF。定义换能器的发射阵元数量为N,算法中换能器的发射阵元数量与接收阵元数量一致。定义射频信号的大小为RF[W][H],其中W表示信号列数即回波信号波束数量,大小与接收阵元数量大小一致(W=N),且一一对应,又被称为扫描线;H表示与采样点数或者采样深度,代表了扫描区域内不同深度位置的回波幅度。
传统延迟叠加算法主要由下面两个公式实现:
算法流程:
二、反向延迟叠加算法
在上面介绍的延迟叠加算法是通过计算声束传播距离,进而计算得到数组下标的方式完成叠加,其叠加过程是根据射频信号中某一点的坐标计算该点在回波信号中的位置实现的。反向延迟叠加的思路是通过回波信号中某一点的坐标逆向计算该点在射频信号中的位置来实现波束叠加(如下图所示)。
在第一节中,叠加算法是先计算距离,再计算出射频信号某一点在回波信号中对应的位置,然后将回波信号的对应位置叠加。本节中将这个过程逆向,将距离计算公式中的一个射频信号采样点下标j作为代求未知量,计算公式如下:
式中,由于fs、c、w、h均为常数,所以cm/(2fsh)、fw^2/(2ch)均为常数。计算中仅包含两次乘法、一次除法,还去除了复杂的开平方操作,这就极大提升了运行速度。
三、结果展示
四、算法关键源码
4.1 延迟叠加算法
%% DAS波束合成ticfor i=1:W for j=1:H for k=1:N d=j*h+sqrt(((k-i)*w)^2+(j*h)^2); % 计算发射到接收传感器的距离,基于两点距离公式,探头与图像中间对齐(为保证中间较强信号保留)% d=sqrt(((i-W/2)*w-(k-N/2)*pitch)^2+(j*h)^2)+j*h; %计算延迟 time=round(d*fs/c); if time==0 time=time+1; end if(time<=size(sensor_data,1)) RF(i,j) = RF(i,j) + sensor_data(time,k,i);%可以在这加窗.*window_function% RF(i,j) = RF(i,j) + sensor_data(time,k,i).*Lateral_window(k); end end endendRF=RF'; % 图像旋转% imshow(RF)tocdisp(['运行时间: ',num2str(toc)]);
4.2 反向延迟叠加算法
%% RDAS波束合成tic%计时for i=1:W for j=1:N for k=1:H index=(c/(fs*2*h))*k-(fs*w^2/(c*2*h))*(i-j)^2/k; index=round(index); if index<=0 index=1; elseif index>H index=H; end if(index<=size(sensor_data,1)) %DOA RF(i,index) = RF(i,index) + sensor_data(k,j,i);%可在此处加窗% RF(i,index) = RF(i,index) + sensor_data(k,j,i).*Lateral_window(j); end end endendRF=RF'; % 图像旋转tocdisp(['运行时间: ',num2str(toc)]);
工程源码详见:https://github.com/Lxp181417/DAS_beamforming.git
可私信论文原文(仅具参考解释意义)!!!