CN108670304B - 一种基于改进dmas算法的超声平面波成像方法 - Google Patents

Abstract

本发明属于超声平面波成像领域，具体涉及一种基于改进DMAS算法的超声平面波成像方法，包括以下步骤：1)通过FieldII仿真的B超设备发射一定复合角度的平面波超声信号；2)计算所述平面波超声信号的回波数据平方根

及其累加和项；3)对延时乘累加波束合成算法，即DMAS，进行改进，得出延时累加乘波束合成算法，即DSBM算法：4)重复步骤2至步骤3，得到每一帧图像数据；5)结合广义相干系数，得出DSBMGCF算法，并根据该算法得到修正后的成像结果。本发明融合了延时乘累加波束合成算法与广义相干系数的优点所在，在保证较高的成像帧频的前提下，解决了平面波空间复合成像图像质量与成像帧频不能兼得的问题，并节约了内存资源的利用率。

Description

一种基于改进DMAS算法的超声平面波成像方法

所属技术领域

本发明属于超声平面波成像领域，具体涉及一种基于改进DMAS算法的超声平面波成像方法。

技术背景

传统超声的波束形成环节通常采用单波束的延时叠加方法，对各个阵元接收的回波信号进行相应的延时后相加从而得到一条扫描线的信号，要想形成一幅超声图像就需要上百次的波束发射和处理，故，这种方法虽然简单易行，但成像帧频较低，无法满足心脏等运动型器官的成像需求。

针对帧频低的问题，在医学超声领域提出了超快速超声成像的概念，它一般是以平面波发射为基础，通过一次发射即可得到整个成像区域的数据，大大降低了成像所需的发射次数，从而提高了成像帧频，但是由于平面波超声成像发射的波束没有进行聚焦，波束较宽，每次接收的回波数据噪声过多，信噪比低，从而导致图像分辨率和对比度等都大大降低。

故，为提高成像质量，平面波空间复合成像方法被提出。它是通过改变超声换能器发射角度或者频率，从多个角度或频率获得多幅超声成像图，并将多幅图像叠加平均得到复合图像。对图像进行叠加和平均的过程可以有效平滑噪声，另外，不同的发射角度也可以有效解决成像目标边沿缺失的问题从而提高图像质量。因此，复合角度个数越多，成像质量便越好，但同时会使得成像时间变长导致一定的帧频损失，并且由于不同成像区域的复合次数不同，成像质量分布不均匀而不能有效改进图像的分辨率。

因此，改进平面波空间复合成像方法，在保持较高成像帧频的基础上提高图像质量成为快速超声成像技术的研究热点。本发明正是基于此才得以研究。

发明内容

针对上述存在的技术问题，本发明提供一种基于改进DMAS算法的超声平面波成像方法，包括以下步骤：

步骤1，通过FieldII仿真的B超设备发射一定复合角度的平面波超声信号；

步骤2，计算所述平面波超声信号的回波数据平方根

及其累加和项；

步骤3，对延时乘累加波束合成算法，即DMAS，进行改进，得出延时累加乘波束合成算法，即DSBM算法：

对所述累加和项进行合成计算，得到一条扫描线上一个点的重建结果，继而得到一条扫描线的波束合成结果y_DSBM(t)，该计算方法为：

步骤4，重复步骤2至步骤3，共得到M条扫描线的波束合成结果y_DSBM(t)，即一帧图像数据；设共有N个所述复合角度，据此方法，得到其中每个所述复合角度对应的每一帧图像数据；

步骤5，结合广义相干系数，得出DSBMGCF算法，并根据该算法得到修正后的成像结果imagesum。

所述步骤1，具体包括以下步骤：

步骤1.1，通过设置延时时间线timedelay，由FieldII仿真的B超设备的换能器发射一定复合角度的平面波超声信号；

timedelay＝(1:M-1)×(w+k)×sinθ/c，

其中w为所述换能器宽度，k为所述换能器之间间隙宽度，M为阵元个数，θ为复合角度，c为所述平面波超声信号的速度；

步骤1.2，换能器接收到射频回波信号传至后端处理。

所述步骤2，具体包括以下步骤：

步骤2.1，对射频回波信号做延时累加乘波束合成，得到回波数据x_m(t)；

步骤2.2，对x_m(t)做符号运算并求取平方根，即回波数据平方根

求前M-1项通道对应的回波数据平方根的累加和项，即：

所述步骤2.1，具体包括以下步骤：

步骤2.1.1，求每一通道上的延时时间及对应的延迟时间点数pointdelay，即

pointdelay＝round(timedelay×f_s-timesnum+1)，

式中，timesnum指盲点数；

步骤2.1.2，求得进行不同延迟时间点数后的各通道对应的回波数据x_m(t)，即

x_m(t)＝env(pointdelay,m)

式中，m表示第m条通道，M为所述通道的数量，1≤m≤M，env指的是采集到的原始回波数据矩阵。

所述步骤5，具体包括以下步骤：

步骤5.1，对步骤4得出的每一帧图像数据进行希尔伯特解调以及对数压缩得到数据矩阵env_dB，并重置为行向量，并据此获得由所有所述复合角度图像数据组成的图像矩阵P；

步骤5.2，通过低频能量占总能量的比值GCF，将所述图像矩阵P进行修正，得到成像结果imagesum。

所述步骤5.2，具体包括以下步骤：

步骤5.2.1，对所述图像矩阵P求取二维傅里叶变换，获得空间频谱矩阵p(f)，即

式中，p的上标表示所述复合角度的个数，下标表示阵列维度对应的频率分量；

步骤5.2.2，在所述空间频谱矩阵p(f)中，频率为0时的直流成分为低频能量，计算低频能量占总能量的比值GCF，该计算方法为：

式中，|P(K₁,K₂)|代表所述空间频谱矩阵p(f)中第K₁行，第K₂列的频谱幅值；M₁和M₂代表低频频率的最大值；

步骤5.2.3，对N个复合角度下得到的数据矩阵env_dB进行叠加和平均，得到平均结果imagesum(x)；

步骤5.2.4，用(1+GCF)乘以步骤5.2.3的平均结果，得到修正后的成像结果imagesum，即：

imagesum＝(1+GCF)*imagesum(x)。

本发明的有益效果：

本发明提出一种基于改进DMAS算法的超声平面波成像方法，融合了延时乘累加波束合成算法与广义相干系数的优点所在，在保证较高的成像帧频的前提下，通过提高图像对比信噪比和对比度和分辨率，解决了平面波空间复合成像图像质量与成像帧频不能兼得的问题，并通过减小算法复杂度、节约了内存资源的利用率；不论在散射点和囊肿目标仿真上，基于本发明所述方法的超声成像效果都是最好的。

本发明设计合理，易于实现，具有很好的实用价值。

附图说明

图1为本发明具体实施方式中所述不同成像算法的散射点成像结果图；

图2为本发明具体实施方式中所述散射点目标在30mm深度处的成像结果横向剖面图；

图3为本发明具体实施方式中所述散射点目标在60mm深度处的成像结果横向剖面图；

图4为本发明具体实施方式中所述不同角度个数下CO算法与DSBMGCF算法的散射点仿真比较图；

图5为本发明具体实施方式中所述不同角度个数下CO算法与DSBMGCF算法的散射点横向剖面图近场比较图；

图6为本发明具体实施方式中所述不同角度个数下CO算法与DSBMGCF算法的散射点横向剖面图远场比较图；

图7为本发明具体实施方式中所述不同成像算法的囊肿仿真效果对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施实例，对本发明做出进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明应用于FieldII仿真的B超设备，具体提出一种基于改进DMAS算法的超声平面波成像方法，包括以下步骤：

步骤1，通过FieldII仿真的B超设备发射一定复合角度的平面波超声信号，具体包括以下步骤：

步骤1.1，通过设置延时时间线timedelay，由FieldII仿真的B超设备的换能器发射一定复合角度的平面波超声信号；

timedelay＝(1:M-1)×(w+k)×sinθ/c，

其中w为所述换能器宽度，k为所述换能器之间间隙宽度，M为阵元个数，θ为复合角度，c为所述平面波超声信号的速度；

步骤1.2，换能器接收到射频回波信号传至后端处理；

步骤2，计算所述平面波超声信号的回波数据平方根

及其累加和项，具体包括以下步骤：

步骤2.1，对射频回波信号做延时累加乘波束合成，得到回波数据x_m(t)，包括以下步骤：

步骤2.1.1，求每一通道上的延时时间及对应的延迟时间点数pointdelay，即

pointdelay＝round(timedelay×f_s-timesnum+1)，

式中，timesnum指盲点数；

步骤2.1.2，求得进行不同延迟时间点数后的各通道对应的回波数据x_m(t)，即

x_m(t)＝env(pointdelay,m)

式中，m表示第m条通道，M为所述通道的数量，1≤m≤M，env指的是采集到的原始回波数据矩阵；

步骤2.2，对x_m(t)做符号运算并求取平方根，即回波数据平方根

求前M-1项通道对应的回波数据平方根的累加和项，即：

步骤3，对延时乘累加波束合成算法，即DMAS，进行改进，得出延时累加乘波束合成算法，即DSBM算法：

对所述累加和项进行合成计算，得到一条扫描线上一个点的重建结果，继而得到一条扫描线的波束合成结果y_DSBM(t)，该计算方法为：

步骤4，重复步骤2至步骤3，共得到M条扫描线的波束合成结果y_DSBM(t)，即一帧图像数据；设共有N个所述复合角度，据此方法，得到其中每个所述复合角度对应的每一帧图像数据；

步骤5，结合广义相干系数，得出DSBMGCF算法，并根据该算法得到修正后的成像结果imagesum，具体包括以下步骤：

步骤5.1，对步骤4得出的每一帧图像数据进行希尔伯特解调以及对数压缩得到数据矩阵env_dB，并重置为行向量，并据此获得由所有所述复合角度图像数据组成的图像矩阵P；

步骤5.2，通过低频能量占总能量的比值GCF，将所述图像矩阵P进行修正，得到成像结果imagesum，具体包括以下步骤：

步骤5.2.1，对所述图像矩阵P求取二维傅里叶变换，获得空间频谱矩阵p(f)，即

式中，p的上标表示所述复合角度的个数，下标表示阵列维度对应的频率分量；

步骤5.2.2，在所述空间频谱矩阵p(f)中，频率为0时的直流成分为低频能量，计算低频能量占总能量的比值GCF，该计算方法为：

式中，|P(K₁,K₂)|代表所述空间频谱矩阵p(f)中第K₁行，第K₂列的频谱幅值；M₁和M₂代表低频频率的最大值，本实施例中，M₁＝M₂＝1；

步骤5.2.3，对N个复合角度下得到的数据矩阵env_dB进行叠加和平均，得到平均结果imagesum(x)；

步骤5.2.4，用(1+GCF)乘以步骤5.2.3的平均结果，得到修正后的成像结果imagesum，即：

imagesum＝(1+GCF)*imagesum(x)。

本发明在保证较高的成像帧频的前提下，通过提高图像对比信噪比和对比度和分辨率，解决了平面波空间复合成像图像质量与成像帧频不能兼得的问题，并通过减小算法复杂度、节约了内存资源的利用率。

本发明技术融合了延时乘累加波束合成算法与广义相干系数的优点所在。延时乘累加波束合成算法，即DMAS，依据回波信号之间具有较强的空间相关性，噪声之间不具有这种高相关性，故，其通过对成对组合的回波信号乘积求平方根，可以在很大程度上去除噪声。而本发明技术则是在DMAS算法公式进行改进，得到步骤3所述延时累加乘波束合成算法，即DSBM算法；

从仿真实现角度上来讲，为完成一个像素点的重建工作，所述DSBM算法仅需要M-1次乘法，而现有DMAS算法需要至少

次乘法，因此所述DSBM算法具有较少的乘法运算量，而在高数据精度的计算过程中，过多的乘法运算次数会显著拉长计算时间，经仿真实验发现，DSBM较DMAS算法相比，成像时间缩短了87％；

此外，本发明提出DSBMGCF算法，即在所述DSBM算法的基础上依据广义相干系数能够能够降低由超声波速误差带来的聚焦误差而能够明显抑制旁瓣，突出高相关性区域，故通过步骤5，将N个角度下的平面波发射模式下的波束合成结果y_DSBM(t)按照数据压缩方法进行压缩之后，重置为N个行向量，组成图像矩阵P，计算图像矩阵P的空间二维广义相干系数，并在平面波成像过程中乘以此系数，可以在很大程度上提高图像的对比度以及分辨率；

经MATLAB仿真验证，在囊肿仿真成像时相对DMAS算法，所述DSBMGCF算法的对比度提高了3.1倍，对比信噪比提高了1.59倍，成像时间缩短了5倍，这说明本发明通过引入广义相干系数，在一定程度上弥补了DMAS算法中因回波信号自身较弱的相关性带来的能量损失，在一定程度上提高了超声图像质量，并大幅度缩短了成像时间，保证了成像帧频。

针对上述超声平面波成像方法，借助FieldII仿真工具在MATLAB上进行点目标和囊肿目标的仿真实验，具体实验方法为：

点目标设置10对点，分布深度为0～100mm，横向距离为4mm,纵向距离为10mm；复合角度间隔为0.25°,阵元个数M为128。当复合角度个数为1个的时候，不同成像算法的成像结果如图1所示，图1中distance表示所述分布深度，Depth表示所述横向距离和纵向距离；所述不同成像算法分别为：CO算法、GCF算法、DMBS算法、DMAS算法和DSBMGCF算法；

图1中可见，本发明的成像方法可以极好的消除图像中的伪影现象，具有良好的散射点聚焦效果。图2和图3为散射点目标分别在30mm和60mm深度处的成像结果横向剖面图。从旁瓣上来看，本发明不论是在近场还是在远场，旁瓣明显低于其他算法。从主瓣宽度上来看，本发明提出的DSBMGCF算法明显比其他算法窄。通过统计各算法在50mm深度处的散射点的纵向分辨力和横向分辨力得到不同成像算法横向分辨力与纵向分辨力统计表，即表1：

分辨力算法	CO	GCF	DMAS	DSBM	DSBMGCF
						横向分辨力(mm)	1.9	0.93	0.41	0.41	0.31
纵向分辨力(mm)	0.39	0.16	1	1	0.18

表1

由表1可以看到，本发明提出的DSBMGCF算法具有最小的横向分辨力和纵向分辨力，且略等于CO算法在40个复合角度的情况。由此说明，本发明提出的DSBMGCF算法在散射点聚焦成像时，可以在复合角度个数极少的情况下，表现出较好的图像分辨率，即本发明提出的DSBMGCF算法可以在实现较高成像帧频的前提下，保证较好的成像质量。

保持其他成像参数不变，通过设置平面波复合成像算法，即CO算法的复合角度个数分别为1、7、15、20和40来进行多角度的平面波空间复合成像，成像结果如图4所示，图4中distance表示所述分布深度，Depth表示所述横向距离和纵向距离；由图4可以看到，随着复合角度个数的增加，CO算法的散射点聚焦效果越来越好，当复合角度个数为40个的时候，成像效果不如DSBMGCF算法在复合角度个数为1时的效果。图5和图6为散射点目标分别在30mm和60mm深度处的横向剖面图，显而易见，本发明不论在近场还是在远场都拥有最低的旁瓣和最窄的主瓣。同样通过统计各算法下在50mm深度处的散射点的纵向分辨力和横向分辨力得到不同角度个数下CO算法与DSBMGCF算法的横向分辨力与纵向分辨力统计表，即表2：

表2

在囊肿仿真实验中，设置500000个强度服从高斯分布的点，成像深度为30～80mm，同时在55mm深度处，设置一个圆形囊肿，直径为10mm，囊肿内部点强度均为0。为了清晰看到成像细节，设置数据动态压缩范围为60dB。不同成像算法的成像结果如图7所示。从图7中可见，GCF算法的囊肿内外对比强烈，但是囊肿边缘毛刺现象较为严重；DMAS和DSBM算法下的囊肿横向分辨率很好，但是从纵向分辨率上看，毛刺较多，使得囊肿内外界限不明朗；而DSBMGCF算法的囊肿具有较清晰的囊肿边缘。现取55mm深度处，囊肿内部和外部各一块边长为7mm的矩形，接近囊肿内接矩形，如图中白色矩形所示，利用以下公式分别计算不同算法的对比度CR和对比信噪比CNR；

CR＝μ_b-μ_c

式中，μ_b和μ_c分别表示背景区域和囊肿内部区域的平均强度，δ_c和δ_b分别代表背景区域和囊肿区域的方差。CR绝对值越大，图像对比度越好。而CNR代表囊肿内部和外部噪声变化的一个比值，CNR数值越大意味着囊肿更容易被可视化以及亮度分布更均匀。表3为复合角度个数为1时，不同成像算法下囊肿成像的CR和CNR值的统计结果；

	CO	GCF	DMAS	DSBM	DMASGCF	DSBMGCF
							μ<sub>c</sub>	-45.08	-75.51	-35.42	-41.76	-67.96	-80.18
μ<sub>b</sub>	-15.51	-25.97	-27.32	-25.49	-52.4	-47.02
							CR	29.58	49.54	8.1	17.27	15.54	33.16
CNR	3.49	3.49	0.83	2.15	0.83	2.15

表3

由表3可以看到，DSBMGCF算法μ_c最小，说明囊肿内部强度最趋于零，即内部噪声最小，同时对比度CR大于CO算法，远大于DMAS算法，仅次于GCF算法。由于囊肿仿体的采样点数目有限，与真实人体组织有一定的差距，所以在囊肿的外部区域会有一定的黑色斑点填充，导致囊肿背景区域平均强度较小，DSBMGCF算法的CNR值略受影响，不过，相对DMAS算法来说要高很多。

综上，不论在散射点和囊肿目标仿真上，基于本发明所述方法的超声成像效果都是最好的。

Claims (6)

1.一种基于改进DMAS算法的超声平面波成像方法，其特征在于，包括以下步骤：

步骤1，通过FieldII仿真的B超设备发射一定复合角度的平面波超声信号；

步骤2，计算所述平面波超声信号的回波数据平方根

及其累加和项；

步骤3，对延时乘累加波束合成算法，即DMAS，进行改进，得出延时累加乘波束合成算法，即DSBM算法：

对所述累加和项进行合成计算，得到一条扫描线上一个点的重建结果，继而得到一条扫描线的波束合成结果y_DSBM(t)，y_DSBM(t)的计算方法为：

其中，M为通道的数量；

步骤4，重复步骤2至步骤3，共得到M条扫描线的波束合成结果y_DSBM(t)，即一帧图像数据；设共有N个所述复合角度，据此方法，得到其中每个所述复合角度对应的每一帧图像数据；

步骤5，结合广义相干系数，得出DSBMGCF算法，并根据DSBMGCF算法得到修正后的成像结果imagesum。

2.根据权利要求1所述的基于改进DMAS算法的超声平面波成像方法，其特征在于，所述步骤1，具体包括以下步骤：

步骤1.1，通过设置延时时间线timedelay，由FieldII仿真的B超设备的换能器发射一定复合角度的平面波超声信号；

timedelay＝(1:M-1)×(w+k)×sinθ/c，

其中w为所述换能器宽度，k为所述换能器之间间隙宽度，M为阵元个数，θ为复合角度，c为所述平面波超声信号的速度；

步骤1.2，换能器接收到射频回波信号传至后端处理。

3.根据权利要求2所述的基于改进DMAS算法的超声平面波成像方法，其特征在于，所述步骤2，具体包括以下步骤：

步骤2.1，对射频回波信号做延时累加乘波束合成，得到回波数据x_m(t)；

步骤2.2，对x_m(t)做符号运算并求取平方根，即回波数据平方根

求前M-1项通道对应的回波数据平方根的累加和项，即：

4.根据权利要求3所述的基于改进DMAS算法的超声平面波成像方法，其特征在于，所述步骤2.1，具体包括以下步骤：

步骤2.1.1，求每一通道上的延时时间及对应的延迟时间点数pointdelay，即

pointdelay＝round(timedelay×f_s-timesnum+1)，

式中，timesnum指盲点数；round()为通过四舍五入取整的函数，f_s为超声探头的工作频率；

步骤2.1.2，求得进行不同延迟时间点数后的各通道对应的回波数据x_m(t)，即

x_m(t)＝env(pointdelay,m)

式中，m表示第m条通道，M为所述通道的数量，1≤m≤M，env指的是采集到的原始回波数据矩阵。

5.根据权利要求1所述的基于改进DMAS算法的超声平面波成像方法，其特征在于，所述步骤5，具体包括以下步骤：

步骤5.1，对步骤4得出的每一帧图像数据进行希尔伯特解调以及对数压缩得到数据矩阵env_dB，并重置为行向量，并据此获得由所有所述复合角度图像数据组成的图像矩阵P；

步骤5.2，通过低频能量占总能量的比值GCF，将所述图像矩阵P进行修正，得到成像结果imagesum。

6.根据权利要求5所述的基于改进DMAS算法的超声平面波成像方法，其特征在于，所述步骤5.2，具体包括以下步骤：

步骤5.2.1，对所述图像矩阵P求取二维傅里叶变换，获得空间频谱矩阵p(f)，即

式中，p的上标表示所述复合角度的个数，下标表示阵列维度对应的频率分量；

步骤5.2.2，在所述空间频谱矩阵p(f)中，频率为0时的直流成分为低频能量，计算低频能量占总能量的比值GCF，该计算方法为：

式中，|P(K₁,K₂)|代表所述空间频谱矩阵p(f)中第K₁行，第K₂列的频谱幅值；M₁和M₂代表低频频率的最大值；

步骤5.2.3，对N个复合角度下得到的数据矩阵env_dB进行叠加和平均，得到平均结果imagesum(x)；

步骤5.2.4，用(1+GCF)乘以步骤5.2.3的平均结果，得到修正后的成像结果imagesum，即：

imagesum＝(1+GCF)*imagesum(x)。

Images