基于频域广义旁瓣相消和维纳后置滤波的超声成像方法

1.本发明属于超声成像技术领域，涉及一种基于频域广义旁瓣相消和维纳后置滤波的超声成像方法。


背景技术：

2.超声成像中应用最为广泛的，也是最简单的波束形成技术即延时叠加算法(delay and sum，das)，它是根据阵元通道几何位置关系对所接收的回波信号进行延时量的计算，然后对延时后的数据对齐叠加。传统das算法复杂度低，成像速度快，但由于其采用固定窗函数加权导致主瓣宽度增加，分辨率较低。
3.近年来，为了提高波束形成算法的对比度和分辨率，自适应算法被应用于超声成像领域。其中，最小方差(minimum variance,mv)波束形成算法是目前使用最为广泛的自适应算法。该方法在保持期望方向增益不变的情况下，使阵列输出能量达到最小，通过动态计算聚焦延时信号的加权矢量，对信号进行加权求和，提高了图像的分辨率，但该算法的缺点是稳健性远不如传统的延时叠加算法，而且容易使有用信号相消，这在信噪比较低的情况下对图像质量有较大影响，并且在对比度上存在明显不足。广义旁瓣相消算法(generalized sidelobe canceler,gsc)作为最小方差算法的一种鲁棒结构，具有比最小方差算法更好的背景成像质量。目前，有学者提出将广义旁瓣相消器引入频域(generalized sidelobe canceler based on frequency domain segmentation,fdsgsc)，从而改善算法适用条件，但是算法分辨率和背景质量仍然需进一步提高。
4.因此，亟需发明一种能够进一步提高超声成像广义旁瓣相消波束形成算法分辨率，对比度及稳健性的算法，促进超声成像质量进一步改善。


技术实现要素：

5.有鉴于此，本发明的目的在于提供一种基于频域广义旁瓣相消和维纳后置滤波的超声成像方法，通过在频域广义旁瓣相消算法输出的基础上，设计一个频域维纳后置滤波器对输出信号进行加权，进而显著提高超声图像的成像分辨率和对比度，能有效克服广义旁瓣相消算法对比度低的问题，改善图像整体成像质量。
6.为达到上述目的，本发明提供如下技术方案：
7.一种基于频域广义旁瓣相消和维纳后置滤波的超声成像方法，具体包括以下步骤：
8.s1：对阵元接收的采样信号进行延时处理，以获得超声回波信号；
9.s2：选取适合超声回波信号的频域分段长度，同时利用短时傅里叶变换将超声回波信号转换为频域的窄带子信号，即子频带的超声回波数据；
10.s3：对子频带的超声回波数据依次进行子阵空间平滑、对角加载等处理，得到子频带下的频域样本协方差矩阵；
11.s4：根据广义旁瓣相消原理计算出各个子频带下的波束形成权值；
12.s5：根据子频带下的频域广义旁瓣相消器的波束形成权值设计相应的维纳后置滤波器，求取维纳滤波器系数；
13.s6：将子频带下的频域广义旁瓣相消器的波束形成权值与其后置滤波器系数相乘，得到最终的复合权矢量；
14.s7：利用复合权矢量对频域下的回波数据进行加权求和，得到波束形成器的频域输出；
15.s8：利用快速傅里叶逆变换将波束形成器的频域输出值转换为时域输出值，得出最终时域输出进行成像。
16.进一步，步骤s2中，获取子频带的超声回波数据，具体包括以下步骤：
17.s21：利用短时傅里叶变换对超声回波信号x(k)进行频域分割，过程如下式所示：
[0018][0019]
其中，z(k)表示窗长和短时傅里叶变换点数均为64且无信号重叠的hanning窗，x(k)表示超声回波时域信号，s(m，ω)表示信号x(k)在时频域中的对应数值，m表示分段序号，ω表示在分段序号m下的频域序号，k表示采样时刻，k表示采样时刻总数，i表示虚数单位，e表示指数；
[0020]
s22：基于上一步骤，将每个传感器阵元的超声回波信号依次转换为若干个独立的等间隔窄带子信号，第n个孔径上的第m段窄带的子信号sn(m，ω)表达式为：
[0021]
sn(m，ω)＝[sn(m，1)，...，sn(m，w-1)，sn(m，w)]
[0022]
其中，m＝1，2，...，m，m表示第n个孔径信号一共被划分的段数；ω＝1，2，...，w，w表示各个窄子带信号的长度；m
×
w表示孔径信号在时域中的总长度，sn(m，ω)表示第n个孔径在时频域中的点(m，ω)处的信号幅值。
[0023]
进一步，步骤s3中，得到子频带下的样本协方差矩阵，具体包括以下步骤：
[0024]
s31：提取每个时频点的阵列信息，时频点(m，ω)的阵列信号x(m，ω)表示为：
[0025]
x(m，ω)＝[s1(m，ω)，s2(m，ω)...，sn(m，ω)]
[0026]
s32：对各个时频点阵列信号进行子阵划分，并基于子阵划分之后的阵列信号进行构造频域样本协方差矩阵r(m，ω)：
[0027][0028]
其中，x
l
(m，ω)＝[s
l
(m，ω)，s
l 1
(m，ω)，...，s
l l-1
(m，ω)]表示第l个子阵的频域前向平滑向量，并且l＝(1，2，...，n-l 1)，n表示超声阵列总长度，l表示子阵列长度；x
l
(m，ω)h为x
l
(m，ω)的共轭转置；
[0029]
s33：对频域样本协方差矩阵r(m，ω)进行对角加载处理，以增强算法的鲁棒性，得到对角加载后的协方差矩阵表达式：
[0030][0031]
其中，β＝trace(r(m，ω))
·
δ，trace(r(m，ω))为信号的等效功率，trace(
·
)是求矩阵迹的函数，δ为空间噪声与信号功率之比，r0为单位矩阵。
[0032]
进一步，步骤s4具体包括：根据广义旁瓣相消原理计算出各个时频点(m，ω)下的
波束形成权值，表达式为：
[0033]wfds-gsc
(m，ω)＝w
sq
(m，ω)-bw
sa
(m，ω)
[0034]
其中，w
fds-gsc
(m，ω)表示子频带广义旁瓣相消器在时频点(m，ω)的波束形成权值，w
sq
(m，ω)为子频带广义旁瓣相消器在时频点(m，ω)的主支路的非自适应权矢量，w
sa
(m，ω)为子频带广义旁瓣相消器在时频点(m，ω)辅助支路的自适应权矢量；其中，w
sq
与w
sa
的关系为：
[0035][0036]
其中，bh为阻塞矩阵b的共轭转置，为对角加载后的频域估计样本协方差矩阵，(
·
)-1
表示矩阵求逆运算。
[0037]
进一步，步骤s5具体包括：将子频带下的广义旁瓣相消器的输出作为时频点(m，ω)期望信号的估计值，设计相应的维纳后置滤波器，得到基于子频带广义旁瓣相消器的频域维纳后置滤波器系数h
fdsgsc-wiener
，表达式为：
[0038][0039]
其中，w
fds-gsc
为子频带广义旁瓣相消器的自适应波束形成权值，为w
fds-gsc
的共轭转置，sn(m，ω)为时频点(m，ω)下第n个孔径的超声频域信号，为sn(m，ω)的共轭转置。
[0040]
进一步，步骤s6具体包括：将子频带广义旁瓣相消器的波束形成权值与其后置滤波器权值相乘，得到最终的复合权值w
fdsgsc-wiener
，表达式为：
[0041]wfdsgsc-wiener
＝h
fdsgsc-wienerwfds-gsc
[0042]
其中，w
fdsgsc-wiener
表示频域广义旁瓣相消器结合维纳后置滤波器的复合权值；h
fdsgsc-wiener
表示基于频域的广义旁瓣相消器的频域维纳后置滤波器系数，w
fds-gsc
表示子频带广义旁瓣相消器的自适应波束形成权值。
[0043]
进一步，步骤s7中，得到波束形成器的频域输出的表达式为：
[0044][0045]
其中，y
fdsgsc-wiener
(m，ω)表示融合维纳后置滤波器的子频带超声广义旁瓣相消波束形成器在时频域点(m，ω)的频域输出值，w
fdsgsc-wiener
表示频域广义旁瓣相消器结合维纳后置滤波器的复合权值，表示w
fdsgsc-wiener
的共轭转置；x
l
(m，ω)表示第l个子阵在时频域点(m，ω)上的频域数据；n表示超声阵列总长度，l表示子阵列长度。
[0046]
进一步，步骤s8中，得到最终时域输出的表达式为：
[0047][0048]
其中，y
fdsgsc-wiene
r(k)表示融合维纳后置滤波器的子频带超声广义旁瓣相消波束形成器在采样时刻k下的时域输出值；istft(
·
)表示短时傅里叶逆变换，y
fdsgsc-wiener
(m，ω)表示融合维纳后置滤波器的子频带超声广义旁瓣相消波束形成器在时频域点(m，ω)的
频域输出值。
[0049]
本发明的有益效果在于：本发明首先利用短时傅里叶变换，将时域超声信号划分为短窄带信号引入时频域，并在时频域做广义旁瓣相消波束形成算法，之后根据频域下波束形成器的输出值计算得到频域下的维纳后置滤波器系数，对频域输出进一步加权优化，进一步提高算法分辨率和旁瓣抑制能力。最后利用短时傅里叶逆变换方法，将时频域输出转换为时域进行成像。因此，本发明能够进一步提高传统超声广义旁瓣相消算法的分辨率和对比度性能，对旁瓣抑制具有更高的效果。
[0050]
本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。
附图说明
[0051]
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：
[0052]
图1为本发明基于频域广义旁瓣相消和维纳后置滤波的超声成像方法的流程图；
[0053]
图2为4种算法点目标成像效果图；
[0054]
图3为4种算法50mm深度下点目标横向分辨率曲线图；
[0055]
图4为4种算法多斑成像效果图；
[0056]
图5为4种算法多斑成像50mm深度下横向分辨率曲线图。
具体实施方式
[0057]
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。
[0058]
请参阅图1～图5，图1为本发明提供的基于频域广义旁瓣相消器的维纳后置滤波器超声成像算法，包括以下步骤：
[0059]
步骤s1：对超声阵元接收的回波信号进行放大、ad转换和延时聚焦处理，以获得超声回波数据；得到延时聚焦处理之后的信号x(k)，x(k)表示为x(k)＝[x1(k)，x2(k)，...，xn(k)]，其中n表示超声阵列的阵元个数，k表示为对应采样深度的采样时刻；
[0060]
步骤s2：选取适合超声回波信号的频域分段长度，同时利用短时傅里叶变换将超声回波信号转换为频域的窄带子信号。具体包括以下步骤：
[0061]
s21：利用短时傅里叶变换对超声回波信号x(k)进行频域分割，过程如下式所示：
[0062][0063]
其中，z(k)表示窗长和短时傅里叶变换点数均为64且无信号重叠的hanning窗，x
(k)表示超声回波时域信号；s(m，ω)表示x(k)信号在时频域中的对应数值；m表示分段序号，ω则表示在分段序号m下的频域序号；k表示采样时刻，k表示采样时刻总数；i表示虚数单位，e表示指数。
[0064]
s22：基于上一步骤，将每个传感器阵元的超声回波信号依次转换为若干个独立的等间隔窄带子信号，第n个孔径上的第m段窄带的子信号sn(m，ω)表达式为：
[0065]
sn(m，ω)＝[sn(m，1)，...，sn(m，w-1)，sn(m，w)]
[0066]
其中，m＝1，2，...，m，m表示该孔径信号一共被划分的段数；ω则表示在分段序号m下的频域序号；ω＝1，2，...，w，w是各个窄子带信号的长度；m
×
w表示孔径信号在时域中的总长度；sn(m，ω)表示第n个孔径在时频域中的点(m，ω)处的信号幅值。
[0067]
步骤s3：提取每个时频点的阵列信息，对子频带的超声回波数据依次进行子阵空间平滑、对角加载等处理，得到子频带下的样本协方差矩阵。具体包括以下步骤：
[0068]
s31：提取每个时频点的阵列信息，时频点(m，ω)的阵列信号x(m，ω)可以表示为：
[0069]
x(m，ω)＝[s1(m，ω)，s2(m，ω)...，sn(m，ω)]
[0070]
s32：对各个时频点阵列信息进行子阵划分，并基于子阵划分之后的阵列信息进行构造频域样本协方差矩阵r(m，ω)：
[0071][0072]
其中，x
l
(m，ω)＝[s
l
(m，ω)，s
l 1
(m，ω)，...，s
l l-1
(m，ω)]表示第l个子阵的频域前向平滑向量，并且l＝1，2，...，n-l 1，n表示超声阵列总长度，l表示子阵列长度，x
l
(m，ω)h为x
l
(m，ω)的共轭转置。
[0073]
s33：对频域样本协方差矩阵r(m，ω)进行对角加载处理，以增强算法的鲁棒性，得到对角加载后的协方差矩阵如下：
[0074][0075]
其中，β＝trace(r(m，ω))
·
δ，trace(r(m，ω))为信号的等效功率，trace(
·
)是求矩阵迹的函数，δ为空间噪声与信号功率之比，r0为单位矩阵。
[0076]
步骤s4：依据广义旁瓣相消器原理计算出各个时频点(m，ω)下的波束形成器输出权值，如下式所示：
[0077]wfds-gsc
(m，ω)＝w
sq
(m，ω)-bw
sa
(m，ω)
[0078]
式中，w
fds-gsc
(m，ω)表示子频带广义旁瓣相消器在时频点(m，ω)的加权矢量，w
sq
(m，ω)为子频带广义旁瓣相消器在时频点(m，ω)的主支路的非自适应权矢量，w
sa
(m，ω)为子频带广义旁瓣相消器在时频点(m，ω)辅助支路的自适应权矢量。其中，w
sq
与w
sa
的关系为：
[0079][0080]
其中，bh为阻塞矩阵b的共轭转置，为对角加载后的频域估计样本协方差矩阵，(
·
)-1
表示矩阵求逆运算。
[0081]
步骤s5：将子频带下的广义旁瓣相消器的输出作为时频点(m，ω)期望信号的估计值，设计相应的维纳后置滤波器，得到基于子频带广义旁瓣相消器的频域维纳后置滤波器
系数h
fdsgsc-wiener
如下：
[0082][0083]
其中，w
fds-gsc
为子频带广义旁瓣相消器的自适应波束形成权值，为w
fds-gsc
的共轭转置，sn(m，ω)为时频点(m，ω)下，第n个孔径的超声频域信号，为sn(m，ω)的共轭转置。
[0084]
步骤s6：将子频带广义旁瓣相消器的权值与其后置滤波器权值相乘，得到最终的复合权矢量w
fdsgsc-wiener
：
[0085]wfdsgsc-wiener
＝h
fdsgsc-wienerwfds-gsc
[0086]
其中，w
fdsgsc-wiener
表示融合维纳后置滤波器的子频带超声广义旁瓣相消器算法最终求得的复合权矢量。h
fdsgsc-wiener
表示基于频域的广义旁瓣相消器算法的频域维纳后置滤波器系数，w
fds-gsc
表示原频域广义旁瓣相消波束形成器的加权系数。
[0087]
步骤s7：利用复合权矢量对频域下的回波数据进行加权求和，得到波束形成器的频域输出，如下式所示：
[0088][0089]
其中，y
fdsgsc-wiener
(m，ω)表示融合维纳后置滤波器的子频带超声广义旁瓣相消波束形成器在时频域点(m，ω)的频域输出值，w
fdsgsc-wiener
表示频域广义旁瓣相消器结合维纳后置滤波器的复合权值，且表示w
fdsgsc-wiener
的共轭转置。x
l
(m，ω)表示时频域点(m，ω)上，第l个子阵的频域数据；n表示超声阵列总长度，l表示子阵列长度。
[0090]
步骤s8：利用快速傅里叶逆变换将该波束形成器的频域输出值转换为时域输出值，得出最终时域输出进行成像：
[0091][0092]
其中，y
fdsgsc-wiener
(k)表示融合维纳后置滤波器的子频带超声广义旁瓣相消波束形成器在采样时刻k下的时域输出值，istft(
·
)表示短时傅里叶逆变换，y
fdsgsc-wiener
(m，ω)表示融合维纳后置滤波器的子频带超声广义旁瓣相消波束形成器在时频域点(m，ω)的频域输出值。
[0093]
验证实验：
[0094]
field ii是丹麦理工大学基于声学原理开发的一款超声实验仿真平台，其在理论研究上获得了广泛的认可和使用。为验证所提算法的有效性，利用field ii对超声成像中常用的点散射目标和无回声囊肿目标进行成像并利用实际实验数据进行成像对比实验。在点目标仿真实验中，设置单列纵向间隔为2.5mm的12个点目标，深度分布在40mm～70mm之间，采用发射定点聚焦和接收动态聚焦方式，发射焦点固定在50mm处，设置图像的成像动态范围为60db。
[0095]
在多斑成像仿真实验中，设置一列横向位置位于中心0mm处的，纵向位于深度32.5mm，50mm，67.5mm的3个散射点目标，并且在纵向50mm处，横向
±
5mm处另设置两个散射
点目标，用于观察各个算法的横向分辨率，采用合成孔径聚焦方式，并设置图像的成像动态范围为60db。同时，在斑点介质中，设置两个半径为3mm的无回声囊肿，圆心分别位于(-5mm，40mm)、(5mm，50mm)，以及两个半径为3mm的强散斑，圆心位于(5mm，40mm)和(-5mm，50mm)。散射点在块状囊肿和背景之间的振幅比为10倍，在无回声囊肿和背景之间的振幅比为40倍。另外，实验所采用的阵元中心频率为7mhz，阵元数目为64个，间距为0.24mm，采样频率为100mhz，声速为1540m/s，设成像动态范围为60db。
[0096]
对上述二个实验目标采用延时叠加算法(das)，广义旁瓣相消算法(gsc)，子频带广义旁瓣相消算法(fdsgsc)，融合维纳后置滤波器的子频带广义旁瓣相消算法(fdsgsc-wiener)进行对比成像实验。
[0097]
图2给出了4种算法点目标成像结果，从图2中可以看出das算法成像质量最差，旁瓣伪像严重，相比于其他算法具有更高的横向伪影和更宽的主瓣宽度，难以区分目标点。gsc算法在主瓣宽度上明显缩窄，横向分辨率明显提高，但是在旁瓣伪像的抑制上仍然表现不佳。
[0098]
fdsgsc算法在整个深度范围内分辨率进一步提高，但是伪影仍然可见。fdsgsc-wiener算法进一步提高了成像质量，主瓣宽度和旁瓣等级为4种算法最优成像质量最佳。并且其远场分辨率得到了明显的改善。
[0099]
图3给出了50mm焦点处4种算法横向分辨率曲线图，表1给出了4种成像算法在点目标情况下，40mm，50mm，60mm深度时的在-6db下的主瓣宽度数值，可以看出，fdsgsc-wiener的主瓣宽度最优，代表其具有最好的点目标分辨率。fdsgsc第二，das的效果差。相比于das算法，fdsgsc-wiener在分辨率主瓣宽度的改善最大提高了80.7％；相比于gsc算法，fdsgsc-wiener在分辨率主瓣宽度的改善最大提高了51.6％，分辨率改善明显。
[0100]
表1点目标实验中不同深度下4种算法-6db的fwhm对比
[0101][0102]
图4给出4种算法多斑成像效果图，表2给出了多斑成像实验下，4种算法不同深度对点目标横向分辨率曲线在-6db的主瓣宽度数据。结合表2和图4可以看出，das算法的分辨率最差，点成像最宽，分辨率低，并且暗斑内部具有较为明显的伪像。gsc算法相比于das算法分辨率有所提高，fdsgsc-wiener则相比gsc和fdsgsc算法在分辨率性能上进一步地改善。相比于das，和gsc算法，fdsgsc-wiener的分辨率提高了74.2％和50％，可以看到fdsgsc-wiener算法的远场分辨率改善程度非常明显。为了更直观地观察分辨率情况，图5给出了4中算法多斑成像50mm深度下横向分辨率曲线图，可以看出fdsgsc-wiener算法具有最优的主瓣宽度，和最低的旁瓣等级，这表明其具有最高的点目标成像性能。
[0103]
表2多斑成像实验中不同深度下4种算法-6db的fwhm对比
[0104][0105]
表3给出了多斑成像仿真不同成像算法成像性能指标的比较，可以看出，相比于das算法，gsc算法在对比度cr上有所损失，并且信噪比下降明显。fdsgsc算法在cr方面的改善也仅仅与das持平，而fdsgsc-wiener的cr值却达到了33.18，相比于gsc算法提高了5.75，这表明fdsgsc-wiener对于斑成像质量的提高是非常明显的，并且他的cnr和sd数值相比于gsc算法都有明显的改善。
[0106]
表3多斑成像仿真不同成像算法成像性能指标的比较
[0107][0108]
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。