鲍喜荣, 沈晓燕, 张石, 苏婷
东北大学 计算机科学与工程学院, 辽宁 沈阳 110169
收稿日期：2018-06-21
基金项目：中央高校基本科研业务费专项资金资助项目(N171604011)。
作者简介：鲍喜荣(1978-), 男, 湖北当阳人, 东北大学讲师, 博士;
张石(1963-), 男, 辽宁抚顺人, 东北大学教授，博士生导师。

摘要：针对基于DMAS波束合成的平面波超声成像算法的图像对比信噪比偏低及算法复杂度大, 无法实现实时成像的问题, 提出了一种具有较高成像质量和较小计算复杂度的平面波超声成像算法DSBMGCF(delay sum before multiply and generalized coherence factor), 借助FieldII仿真工具在Matlab上进行点目标和囊肿目标的仿真实验, 验证了该算法的成像质量.同时, 对所提出的新算法进行了并行化研究和改进, 得到了一种适合在GPU上并行实现的平面波超声成像并行算法PDMASGCF(parallel delay multiply and sum generalized coherence factor), 并在实验室戴尔T7810普通工作站上进行了平面波超声成像时间和质量的验证实验, 获得了较高的成像帧频, 并保证了较好的成像性能.
关键词：平面波DMAS对比信噪比GPU并行
Plane Wave Ultrasonic Imaging Algorithm Based on Improved DMAS and Its GPU Implementation
BAO Xi-rong, SHEN Xiao-yan, ZHANG Shi, SU Ting
School of Computer Science & Engineering, Northeastern University, Shenyang 110169, China
Corresponding author: SHEN Xiao-yan, E-mail: 15566037485@163.com
Abstract: The contrast signal-to-noise ratio(SNR) of the image synthesized by the plane wave ultrasonic imaging algorithm based on DMAS beamforming is low and the real-time imaging problem cannot be realized due to the complexity of the algorithm.A new plane wave ultrasonic imaging algorithm named DSBMGCF(delay sum before multiply and generalized coherence factor) with higher imaging quality and less complexity of calculation was proposed. Then, the simulation experiment of point target and cyst target were done on Matlab by FieldII simulation tool. The imaging quality of the algorithm was verified. At the same time, the proposed algorithm was parallelized and improved, and a parallel algorithm for plane wave ultrasound imaging, PDMASGCF(parallel delay multiply and sum generalized coherence factor), which is suitable for parallel implementation on GPU was obtained. Finally, on the ordinary workstation Dell T7810 in the laboratory, the verification experiments on time and quality of imaging were carried out. Higher imaging frame rate was obtained, while better imaging performance was guaranteed.
Key words: plane waveDMAS(delay multiply and sum)contrast signal-to-noise ratiographic processing unit(GPU)parallel
传统B超系统需要上百次的波束发射和合成才能形成一幅超声图像, 帧频只有35 Hz左右.而为了满足更多器官的实时成像需求, 需要提高超声成像帧率, 但成像质量会受影响.因此, 如何在成像质量和高帧频成像之间作好权衡成为医学超声成像领域内的研究热点之一.目前, 已有学者提出了很多方法.并行接收是在发射端发射一个较宽的波束, 可以实现一次性接收到多个波束, 可以大幅度提高帧频^[1].平面波超声成像方法通过对感兴趣区域进行一次扫描即获得一帧图像, 但因没有进行波束聚焦, 成像不佳.为了解决这个问题, 学者们将发射的不同角度的平面波进行空间复合(coherent, CO), 复合角度个数越多, 成像质量越好, 但帧频会随之降低^[2-3].于是, 一些学者将最小方差波束合成法、特征空间法等应用到平面波发射的快速成像中用于改善成像质量, 但是算法复杂度太高, 成像帧频依然很低^[4].随后, 有学者将延时乘累加波束形成法(delay multiply and sum, DMAS)应用在乳腺癌的检测研究^[5], DMAS便陆续被应用在B超系统的研究中且效果显著, 但在对比噪声比和工程运用方面表现欠佳^[6-11].
为加快超声算法的工程化实现, GPU被越来越多地应用在医学超声领域, 针对声辐射力脉冲弹性成像提出相应的GPU并行算法, 实现了实时处理^[12].2016年, 尉明望利用CUDA实现了平面波相干叠加算法、傅里叶成像算法和F-K迁移算法, 但效果不是太理想^[13].
为改善平面波成像质量, 本文提出了一种基于改进DMAS的平面波超声成像算法; 为提高成像帧频, 加快算法的工程化实现, 本文利用并行化硬件平台GPU对该新型算法作并行化改进和硬件加速, 最后通过执行时间测试和成像性能测试验证了算法的有效性.
1 方法1.1 DMAS平面波复合成像法如图 1所示, DMAS平面波合成算法是对M个激励阵元施加一定的延时时间线, 以发射一定角度的平面波, 通过采集回波信号, 获得不同通道所对应的射频回波信号s_m(t), 将其两两组合, 对组合信号分别进行符号运算和取绝对值后求取算术平方根的操作, 再将两个结果相乘, 最后再将各组合的结果叠加求和, 得到对应位置的像素点.
图 1(Fig. 1)

图 1 DMAS算法框图Fig.1 DMAS algorithm diagram

由此, 可得到基于DMAS的平面波超声成像法的每条波束合成的表达式为

(1)

DMAS算法中的组合相乘再叠加的操作与空间自相关函数的计算步骤等价, 所以, DMAS波束合成算法充分考虑了回波信号之间的相干性, 而噪声之间的这种相干性特别弱, 故可以在很大程度上去除噪声, 获得高质量的超声图像, 但其至少进行C_M²次乘法运算, 成像时间相对较长.
1.2 DSBM平面波复合成像法从工程实现角度来看, 过多的乘法运算次数会增加硬件电路的复杂程度.而从Matlab编程实现层面来讲, 对于高精度的数据运算, 过多的乘法运算次数会在很大程度上延长计算耗时, 从而拉低成像帧频.因此, 为了减少乘法的运算次数, 对DMAS算法公式作以下变形, 即
其中, x_m(t)为延迟时间后的第m通道回波数据.因变形后该算法要先进行累加和, 再进行乘法运算, 称之为延时累加乘波束合成算法(delay sum before multiply, DSBM).从变形公式可以看到, DSBM方案将M-1个公共项提取出来, 分别与对应的不同长度的和因式相乘, 仅需要M-1次乘法, 因此, 从仿真实现角度上来讲, 完成一个像素点的重建工作, DSBM算法具有较少的乘法运算量, 成像时间平均可缩短87%.同时, 从理论上讲, DSBM也在一定程度上保留了延时乘累加波束合成算法利用回波信号之间较强空间相关性去除噪声的性质, 可保证良好的成像质量.
1.3 DSBMGCF平面波超声成像算法基于DSBM波束合成的平面波成像方法通过减少乘法次数缩短了计算时间, 牺牲了一定的抗噪能力.为此, 提出了一种基于DSBM与广义相干系数(generalized coherence factor, GCF)相融合的平面波超声成像方法(DSBMGCF).通过引入GCF, 突出强调成像过程的高相干区域, 消除因声速不均性所导致的相位畸变, 有效改善图像的对比度和分辨率.
现简要介绍本文所采用的GCF的方法:对回波数据组成的二维矩阵X(n)作二维傅里叶变换, 获取其二维空间频谱^[8]p(f).其中上标和下标分别表示复合角度个数和阵列维度对应的频率分量.

(2)

与一维广义相干系数的定义相类似, 定义GCF为低频能量占总能量之比:

(3)

其中, M₁, M₂分别代表合适的低频分量的取值.M₁, M₂都为0, 指的是不同角度的平面波发射的情况下, 各阵元所接收到的所有数据中共同的成分.对于点扩散函数, 令M₁, M₂都为0, 可获得最佳分辨率^[9].
图 2为DSBMGCF平面波超声成像算法的实现流程图.此外, 为了进一步验证引入GCF法可以明显改善图像的分辨率以及对比DSBMGCF平面波成像方法的优越性, 类似地, 将GCF引入DMAS平面波成像方法, 记为DMASGCF, 并进行散射点和囊肿目标点的仿真实验.
图 2(Fig. 2)

图 2 DSBMGCF平面波超声成像算法的实现流程图Fig.2 Flow chart of the new plane wave imaging method DSBMGCF

2 实验及分析2.1 Matlab仿真借助FieldII在Matlab上进行以上几种算法成像下的点目标和囊肿目标的仿真实验.点目标设置10对点, 分布深度为20~60 mm, 横向距离为4 mm, ,纵向距离为10 mm.复合角度间隔为0.25°, 阵元个数为128.首先, 将复合角度个数设为1, 不同成像方法的散射点仿真结果(取深度25~45 mm)如图 3所示.散射点越小越亮, 即散射点的聚焦效果越好, 表明该波束合成算法的散射点成像质量越理想.明显可见, 引入GCF之后, DSBMGCF和DMASGCF平面波成像方法具有最好的散射点聚焦效果, 它可以极好地消除CO算法在水平方向以及DMAS和DSBM算法在纵向方向的伪影现象.
图 3(Fig. 3)

图 3 不同成像方法的散射点仿真结果Fig.3 Results of different imaging methods on scattered point (a)—CO；(b)—DMAS；(c)—DMASGCF；(d)—GCF；(e)—DSBM；(f)—DSBMGCF.

图 4展示了散射点目标分别在30 mm和60 mm深度处的横向剖面图.从旁瓣上看，DSBMGCF算法不论是在近场还是在远场, 旁瓣明显低于其他算法, 且主瓣宽度也明显窄于其他算法.
图 4(Fig. 4)

图 4 不同深度处横向剖面图Fig.4 Transverse map at different depths (a)—30 mm；(b)—60 mm.

在囊肿仿真实验中, 设置500 000个强度服从高斯分布的点, 成像深度为30~80 mm, 设置一个直径为10 mm的圆形囊肿, 内部点强度均为0.为了清晰看到成像细节, 设置数据动态压缩范围为60 dB.不同算法的囊肿仿真结果如图 5所示.
图 5(Fig. 5)

图 5 不同算法的囊肿仿真结果Fig.5 Simulation result of different methods on cyst (a)—CO；(b)—DMAS；(c)—DMASGCF；(d)—GCF；(e)—DSBM；(f)—DSBMGCF.

由图 5可见, GCF算法的囊肿内外对比强烈, 但边缘毛刺较多.DMAS和DSBM算法的囊肿横向分辨率很好, 但从纵向分辨率上看, 毛刺较多, 使得囊肿内外界限不清晰.而DSBMGCF算法的囊肿具有较清晰的囊肿边缘.现取55 mm深度处, 囊肿内部和外部各一块(边长为7 mm的矩形, 接近囊肿内接矩形, 如图中白色方框), 对比度CR和对比信噪比CNR的计算公式为
其中：μ_b, μ_c分别表示背景区域和囊肿内部区域的平均强度, dB; σ_b, σ_c分别代表背景区域和囊肿区域的方差. CR绝对值越大, 图像对比度越好.而CNR代表囊肿内部和外部噪声变化的比值, CNR数值越大意味着囊肿更容易被可视化以及亮度分布更均匀.表 1给出了复合角度个数为1时, 不同算法下囊肿内部CR和CNR的统计结果.
表 1(Table 1)

表 1 不同算法下, 囊肿内部CR和CNR的统计结果Table 1 CR and CNR of cysts for different methods/dB 性能指标 CO GCF DMAS DSBM DMASGCF DSBMGCF μc -45 -75.5 -35.4 -41.8 -68 -80.2 μb -15.5 -25.9 -27.3 -25.5 -52.4 -47.1 CR 29.6 49.5 18.1 17.3 15.5 33.2 CNR 3.5 3.5 0.8 2.15 0.8 2.2

表 1 不同算法下, 囊肿内部CR和CNR的统计结果 Table 1 CR and CNR of cysts for different methods/dB

由表 1中可以看到:DSBMGCF算法μ_c最小, 说明囊肿内部噪声最小, 同时对比度CR大于CO算法, 远大于DMAS算法, 仅次于GCF算法.由于囊肿仿体的采样点数目有限, 与真实人体组织有一定的差距, 所以在囊肿的外部区域会有一定的黑色斑点填充, 导致背景区域平均强度较小, 使DSBMGCF算法的CNR值略受影响, 不过, 相对DMAS要改进很多.
2.2 GPU实现2.2.1 算法的并行化GPU并行计算是一种单指令多数据的计算模式, 所以在算法并行化设计上的主要思路是降低或者去除波束合成算法中数据的依赖性.DSBM算法较DMAS减少了乘法次数, 缩减了计算时间, 但仍存在计算冗余.为此, 本文提出了一种更简单的计算方法, 推导如下.因为
类比多项式求和公式为

(4)

(5)

有

(6)

所以,
可以看到合成一个点需要一次乘法运算(忽略加权过程所消耗的N次乘法运算量).并行改进后的DMAS复合算法(parallel DMAS, PDMAS)具有较低的算法复杂度, 且从并行角度分析, 它可对输入的各通道回波数据单独处理, 具有较强的数据并行性.在单个角度平面波发射模式下, 采用PDMAS波束合成法合成像素点, 且所有像素点共用一个延时计算模型, 像素点之间具有相互独立性, 不同深度、不同通道的像素点的合成工作具备充分的并行条件.因此, 令一个工作项负责一个像素点的合成工作, 所有工作项并发执行, 完成一帧图像的重建工作.PDMAS平面波波束合成并行方案如图 6所示.
图 6(Fig. 6)

图 6 PDMAS平面波波束合成并行方案Fig.6 Parallel implementation scheme of PDMAS plane wave beamforming

而对于广义相干系数部分的并行化实现, 重点是将回波数据的FFT运算移植到GPU.通过设置INIT_kernel, 结合本实验平台为AMD GCN架构以及补零到2的整数次幂之后, 回波数据大小为1 048 576, 令各个工作项访问处理16个元素的FFT, 因此共需要256个工作组, 65 536个工作项, GCF并行方案如图 7所示.
图 7(Fig. 7)

图 7 GCF并行任务分配示意图Fig.7 Assignment diagram of GCF parallel task

首先, 工作项将输入地址按位反转, 接着按照新的地址从全局内存中加载4个float2型数据(每个float2包含两个float型, 对应实数和虚数部分), 并对其作4点FFT.然后, 将4元素归并为8元素的FFT, 循环迭代, 直到工作项所能处理的元素数量达到上限.再设置Stage_kernel, 将16个元素的FFT归并32个元素的FFT序列.同样经过在更大FFT序列中分配一个位置, 完成相应元素的FFT, 这个过程循环往复直到整个FFT计算结束.综上, 将并行化改进后的算法称作PDMASGCF算法, 其成像步骤主要包括3个阶段:阶段一, 进行PDMAS波束合成; 阶段二, 按照空间二维广义相干系数法, 通过并行FFT计算得到回波数据的频域变换, 进而计算得到广义相干系数GCF; 阶段三, 用(1+GCF)作用于阶段一合成的图像, 得到最终的超声图像.
2.2.2 GPU实现在实验室戴尔T7810工作站进行了算法的GPU实现.所采用的硬件平台见表 2.
表 2(Table 2)

表 2 GPU硬件资源Table 2 GPU hardware resources 显卡型号 显存 带宽 CU 单精度浮点TFLOPs OpenCL GB GB·s-1 个 w7100 8 160 1 792 3.3 2.0

表 2 GPU硬件资源 Table 2 GPU hardware resources

图 8为每个工作项的数据计算流程图.可以看到, GPU实现要依赖CPU, GPU和Matlab 3个平台完成.CPU+GPU异构平台主要是利用OpenCL编程, 在CPU端完成波束合成以及回波数据的FFT等计算任务.Matlab端则是进行平面波超声数据的产生、回波数据的采集和广义相干系数的计算以及后续的成像工作.
图 8(Fig. 8)

图 8 每个工作项的数据计算流程图Fig.8 Data calculation flowchart for each work item

因Matlab 2016a在底层实现中对矩阵运算作了优化, 且采用了GPU加速, 所以将GPU实现结果与Matlab的执行结果对比具有一定的参考意义.分别在Matlab和GPU上进行了10次关于执行时间的测试, 结果如表 3所示, GPU和Matlab的平均执行时间分别为16 ms和18.3 s, 相对Matlab来说, 本设计的GPU并行方案具有高达1 129倍的加速比, 且在有限的硬件资源条件下, 成像帧频可达62帧/s, 远超过实时成像的帧频要求.
表 3(Table 3)

表 3 Matlab和GPU的计算时间Table 3 Calculation time on Matlab & GPU 次数 1 2 3 4 5 6 7 8 9 10 16 16 16 16 16 16 16 16 16 16 19 18 19 19 18 18 19 18 17 18

表 3 Matlab和GPU的计算时间 Table 3 Calculation time on Matlab & GPU

接下来借助Matlab进行了成像性能测试.从图 9的单点成像和图 10的囊肿成像的测试结果中可直观看到, 引入GCF方法之后, PDMASGCF算法的Matlab成像质量明显变好, 仅次于DSBMGCF算法.PDMASGCF算法分别在Matlab和GPU上执行后, 发现经GPU加速后, 成像效果基本相同, 且都明显好于DMAS算法.
图 9(Fig. 9)

图 9 单点成像Fig.9 Single point imaging (a)—Matlab上实现DMAS；(b)—Matlab上实现DSBM；(c)—Matlab上实现DSBMGCF；(d)—Matlab上实现PDMAS；(e)—Matlab上实现PDMASGCF；(f)—CPU上实现PDMASGCF.

图 10(Fig. 10)

图 10 囊肿成像Fig.10 Cyst imaging (a)—Matlab上实现DMAS；(b)—Matlab上实现DSBM；(c)—Matlab上实现DSBMGCF；(d)—Matlab上实现PDMAS；(e)—Matlab上实现PDMASGCF；(f)—CPU上实现PDMASGCF.

3 结论本文提出的平面波超声成像算法在对比噪声比和复杂度方面都显著优于DMAS平面波超声成像算法, 在一定程度上解决了平面波成像算法的质量与帧频不能兼得的问题.借助GPU完成了本文算法的并行化改进与硬件加速, 通过执行时间和成像性能的测试, 验证了该算法具有优良的成像性能.同时, 可在普通GPU工作站上实现实时成像, 具有很好的工程应用前景.
参考文献

[1]	Shattuck D P, Weinshenker M D, Smith S W, et al. Explososcan:a parallel processing technique for high speed ultrasound imaging with linear phased arrays[J].The Journal of the Acoustical Society of America, 1984, 75(4): 1273–1282.DOI:10.1121/1.390734
[2]	Tanter M, Fink M. Ultrafast imaging in biomedical ultrasound[J].IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2014, 61(1): 102–119.DOI:10.1109/TUFFC.2014.2882
[3]	Montaldo G, Tanter M, Bercoff J, et al. Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[J].IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2009, 56(3): 489–506.DOI:10.1109/TUFFC.2009.1067
[4]	Deylami A M, Jensen J A, Asl B M.An improved minimum variance beamforming applied to plane-wave imaging in medical ultrasound[C]//2016 IEEE International Ultrasonics Symposium(IUS).Tours, 2016: 1-4.https://www.researchgate.net/publication/309772902_An_improved_minimum_variance_beamforming_applied_to_plane-wave_imaging_in_medical_ultrasound
[5]	Lim H B, Nhung N T T, Li E P, et al. Confocal microwave imaging for breast cancer detection:delay-multiply-and-sum image reconstruction algorithm[J].IEEE Transactions on Biomedical Engineering, 2008, 55(6): 1697–1704.DOI:10.1109/TBME.2008.919716
[6]	赵金鑫, 汪源源, 郭威, 等. 基于广义相干系数的超声平面波复合成像[J].仪器仪表学报, 2014, 35(sup2): 186–190. ( Zhao Jin-xin, Wang Yuan-yuan, Guo Wei, et al. Ultrasonic plane wave composite imaging based on generalized coherence coefficient[J].Chinese Journal of Scientific Instrument, 2014, 35(sup2): 186–190.)
[7]	Matrone G, Savoia A S, Caliano G, et al. The delay multiply and sum beamforming algorithm in ultrasound B-mode medical imaging[J].IEEE Transactions on Medical Imaging, 2015, 34(4): 940–949.DOI:10.1109/TMI.2014.2371235
[8]	Matrone G, Savoia A S, Caliano G, et al.Ultrasound plane-wave imaging with delay multiply and sum beamforming and coherent compounding[C]//The 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society(EMBC).Orlando, 2016: 3223-3226.https://www.ncbi.nlm.nih.gov/pubmed/28268994
[9]	Li P C, Li M L. Adaptive imaging using the generalized coherence factor[J].IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2003, 50(2): 128–141.DOI:10.1109/TUFFC.2003.1182117
[10]	Denarie B, Tangen T A, Ekroll I K, et al. Coherent plane wave compounding for very high frame rate ultrasonography of rapidly moving targets[J].IEEE Transactions on Medical Imaging, 2015, 32(7): 1265–1276.
[11]	Lu J Y. 2D and 3D high name rate imaging with limited diffraction beams[J].IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 1997, 44(4): 839–856.DOI:10.1109/58.655200
[12]	Rosenzweig S, Palmeri M, Nightingale K. GPU-based real-time small displacement estimation with ultrasound[J].IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2011, 58(2): 399–405.DOI:10.1109/TUFFC.2011.1817
[13]	尉明望.超快速超声成像方法研究及其CUDA实现[D].哈尔滨: 哈尔滨工业大学, 2016. ( Yu Ming-wang.Ultra-fast ultrasonic imaging and implementation on CUDA[D].Harbin: Harbin Institute of Technology, 2016.http://cdmd.cnki.com.cn/Article/CDMD-10213-1016774465.htm)