一种基于GPU的声辐射力成像组织位移估计算法

谭勇 王丛知

【摘要】位移估计是弹性成像的第一步，它决定着整个成像的时间以及质量，目前位移估计算法有两大类：传统的归一化互相关和相位检测算法。相位检测法由于计算速度快被广泛运用，通过对相位检测算法的深入研究提出一种新的基于原始RF信号的Doppler位移估计算法，该方法只需一次积和运算以及不用构建IQ信号，大大缩短了计算时间。通过生物仿体研究证明，该方法比传统相位检测法和Doppler速度估计法更快，从性能上看，RF-Doppler法得到应变图像质量与相位检测法相似，但是速度却快了接近100倍。

【关键词】声辐射力成像;位移估计;GPU;并行算法

1.引言

声辐射力弹性成像（ARFI）是一种非侵害性的影像技术，将局部组织的弹性模量数值映射成彩色编码的图像信息[1]。它的原理是，首先利用聚焦超声波束产生的声辐射力在聚焦位置造成组织偏移，并形成向侧向传播的剪切波，然后由高帧频的探测超声波束追踪该剪切波的传播，从射频（RF）数据中计算出剪切波引起的组织位移变化，进而计算出剪切波传播速度，最后根据剪切波速度计算出组织的弹性模量值。它可以提供一些以生物软组织（如乳房、肝脏、肾、肌肉、眼睛、胰腺等）的弹性病变为基础的诊断信息，具有重要的临床价值和广阔的应用前景。

目前超声弹性成像所使用的组织位移估计算法主要包括两大类：传统的归一化互相关算法和基于IQ信号的Doppler相位检测算法。传统的归一化互相关计算需要在一定范围内确定互相关系数的极大值位置，但直接的搜索结果只能是采样周期的整数倍，因此，通常需要对互相关函数极大值点的前后两点之间进行抛物线插值，然后在插值范围内重新寻找更为精确的极大值点位置。插值倍数越大，搜索极大值点位置的精度就越高，但耗时也越长。基于IQ信号的Doppler相位检测算法的原理是，对原始的超声射频（RF）信号进行IQ分解后（其中一个步骤是让信号通过低通滤波器去除载波成分），将所得到的正交IQ信号进行互相关计算，则其互相关系数极大值点的相位，再除以相邻两帧数据间的时间间隔，对应于组织的相对平均移动速度。据此便可进一步计算出精确的组织相对位移距离。该方法不需要进行插值，减少了计算时间并提高了计算精度，已被广泛采用。

要解决计算量和计算速度之间的矛盾，可以通过使用GPU进行并行计算加速实现。NVIDIA公司于2006年推出的“统一编程架构”使得拥有NVIDIA GPU的用户能轻松使用GPU进行编程[2]。GPU具有比CPU更多的运算单元，更高的带宽以及浮点运算能力，因此GPU特别适合高密集型的运算。目前，GPU已经被广泛用在计算机断层扫描、核磁共振等成像上并取得了显著效果，但在超声成像，特别是声辐射力超声弹性成像上的应用还很少见于报道。

本文提出一种应用于ARFI中的特别适合于利用GPU进行加速计算的快速组织位移估计算法，并通过实验证明该算法有较快的计算速度以及不错的应变图像信噪比。

2.组织位移估计算法介绍

目前在商用设备和科学研究中被广泛采用的组织位移估计算法主要包括基于RF信号的改进互相关算法和基于IQ信号的Doppler相移检测算法，本章将分别对现有算法和我们提出的新算法进行介绍。

2.1 基于RF信号的改进互相关算法

传统的互相关算法需要进行大量的错位积和计算和插值计算以提高定位精度。计算量大且组织位移的结果不够准确。为解决这些问题，Cabot在1981年提出了一种通过构造互相关系数的解析信号并对其相位进行分析来提高定位精度的新算法[3]。首先对原始的RF数据分段并对相邻两帧中的对应数据段进行互相关计算，再对所得的互相关系数进行希尔伯特变换，构造其解析信号。该解析信号在互相关系数极大值附近的相位过零点位置，被证明精确地对应于互相关系数的真实极大值点的位置。该算法具体步骤：

1）先对回波信号分段，本文每段长为100个数据点。

2）对每段进行时域互相关：

（1）

式中：m窗口索引;s为帧序列;L为窗长;n是延迟间隔。

3）对互相关函数构造解析信号。

（2）

（3）

h（n）是希尔伯特变换函数[4]，是的希尔伯特变换，为互相关函数的解析信号。

4）寻找相位过零点：

（4）

使用牛顿迭代法[5]寻找到τ使得，即该延时精确对应时域互相关函数最大值的延时。

该方法的弊端除了前面提到的多次错位运算外，运用牛顿迭代并不适合在GPU上并行运算，因为后一次的运算需要前一次运算得出的结果。

2.2 基于IQ信号的Doppler相移检测算法

此算法由Kasai于1985年提出，最初用于Doppler血流速度估计[6]，后来又被用于组织位移估计。该算法首先对原始的RF信号进行IQ分解，在对得到的IQ信号进行分段，并计算相邻两帧中对应数据段之间的点积和，并计算其相对于中心频率的平均相位移动，从而计算出组织的相对位移，其公式如（5）。c是波速，f是RF信号中心频率。

由于ARFI数据中两帧相邻信号间的位移幅度很小，所以运用该方法时基本不会出现相位卷绕的问题。但由于IQ分解所引起的系统复杂度上升和信号信噪比下降的问题依然存在，并最终影响到组织位移估计结果的精度。

（5）

2.3 特别适合于GPU加速的快速组织位移估计新算法

以上两种方法虽然比传统归一化互相关算法高效和精确，但是第一种算法需要多次错位运算和迭代运算，计算量仍旧很大，第二种构建IQ信号需要专门的电路，若用软件实现，计算量大而且繁琐。通过分析发现，构建解析信号与IQ信号虽然从构建方法上看完全不同，但利用其相位进行位移或者速度估计的方法从本质上说是相似的。用软件实现希尔伯特变换来构建解析信号比用软件构建IQ信号简单的多，而且可以利用傅立叶正、反变换在频域上很方便地实现，不再需要完成IQ分解的硬件电路，简化了系统设计，也不再需要低通滤波，保留了完整的超声回波信息，具有很大优势。所以我们提出，综合上述两种方法的优势，直接利用原始的RF数据进行计算，直接对整帧的RF数据进行解析信号的构建，然后再进行分段等后续计算，这一步骤可以大大提高利用GPU进行并行计算加速的效果。算法的大致步骤是：首先对原始RF信号构建希尔伯特变换，分段后，对应数据段内的对应数据点进行一次求积和，公式如下：

（6）

（7）

（8）

是相邻两帧解析信号对应窗口的积和值，m是窗口索引，L是窗长，*是共轭标志，f是RF信号中心频率，c代表波速。

3.基于GPU的快速组织位移估计算法实

3.1 基于RF信号的改进互先关算法

一条A-Line信号总共有K帧每帧有N个采样点，为了达到并行运算，我们分配K-1个线程块，每个线程块有M个线程（M是延迟数量），每个线程块处理S帧和（S+1）帧信号的互相关运算。

由于做互相关时窗口覆盖率达到98%，若使用传统的互相关运算，即S和（S-1）的各个对应窗口进行相关运算，这样会导致大量的重复运算，因此我们使用滑窗的办法。首先每个线程计算第一个窗口的所有数据，计算完后，进行下一个窗口的计算时只需计算新增的点积和，然后再减去被移出的点的点积和，这样至少可以节约90%的运算量。有线程存取共享内存比全局内存有更高的访问速度，所以在做互相关之前先将信号存入共享内存，经验证能提高15%的运算速度。

3.2 基于IQ信号的Doppler相移检测算法

为了要将原始信号的每一帧构造IQ信号，需要将原始信号拷入两个矩阵I和Q。在这两个矩阵里分别对其进行I转换和Q转换。在做互相关时，需要K个线程块，每个线程块有M个线程，K代表需要做互相关的信号，M代表窗口数量。K线程块的序号0、1、2…，分别对应I、Q矩阵的第0、1、2…信号帧，即线程块0处理I信号和Q信号的第0帧信号，线程块1处理I信号和Q信号的第1帧信号…。线程序号对应每个窗口序号，相当于公式（8）中的m，每个线程负责计算该窗口内所有采样点的点积和。为了能提高效率，同样先将要进行互相关运算的I、Q信号存入共享内存，当所有点积和运算完后，使用库函数里的atan2f对该和求出相位值，并将其写回GPU全局内存。

3.3 特别适合于GPU加速的快速组织位移估计新算

该方法在GPU上的实现相当简单。首先我们利用CUDA自带的CUFFT库中的傅立叶变换函数对整帧信号做希尔伯特变换。如果依然采用传统的先分段再构造解析信号的方式，希尔伯特变换会反复多次调用cufftExecC2C函数，且很多数据点会被反复多次计算（因为步长仅为2%，相邻窗口之间有98%的数据点重叠），效率低、耗时长。根据Nvidia的CUFFT库函数使用指南，由于其采用了高效的FFTW快速傅立叶变换算法，若需要进行变换的数据长度固定，可以用cufftPlan1d函数进行预优化处理，并且信号长度越长，计算速度提升的幅度就越大。因此，采用对整帧信号进行解析信号构造的方式，会极大地提高计算效率。同时，经过公式推导和实验验证，我们确认，对于ARFI中所使用的超声RF回波信号来说，这种调换计算顺序的方式对最终的组织位移估计结果的影响极小，可以忽略不计。

4.实验与讨论

4.1 实验

本文使用NIVIDIA GT630M，它有96个CUDA核和1GB显存，编程平台为Visual Studio2010。实验数据采集自我们自制的弹性仿体，三个算法分别对两个区域进行检测。区域1是一个较软与较硬区域位置的交界面的位置，区域（2）是一个硬度分布均匀且较软的区域。以上两个区域大小为10*10cm，区域内每个推动位置的大小为1*1mm，即每一排每一列的推动位置数量为10。为了比较各个算法成像的差异，选择声辐射力弹性成像后的图像质量进行评价。评价弹性图质量的方式主要有信噪比（SNRe）[7]和对比度噪声比（CNRe）：

S是所测区域弹性系数的平均值，δe是所测区域弹性系数的标准差。st、sb分别是硬度区和背景区的均值，δt、δb分别为硬度区和背景区的标准差。信噪比的计算方式是选取硬度均匀的区域，如图2-2，对比度噪声比的计算需选择包含不同硬度值的区域，如图2-1。信噪比值越大说明算法抑制噪声的能力越强，对比度噪声比值越大说明算法弹性检测能力越强。

在编程平台上实现算法时我们采用了在CPU上单独编程和基于CPU和GPU的混合编程，并且对基于RF信号的改进互相关算法采取了两种实现手段：

（1）先构造解析信号，再分段互相关。

（2）先分段互相关，再由互相关系数构造解析信号。

4.2 结果与讨论

表1 各算法的时间对比

算法 时间

NCC_1

NCC_2

New

Doppler CPU/ms

1701

2920

141

5290 GPU/ms

266

1619

2.4

1630

表1给出了各算法在不同窗长下的运行时间。其中NCC_1代表先构造解析信号，再分段互相关的实现方式，NCC_2代表先分段互相关，最后用互相关系数构造解析信号， Doppler代表基于IQ信号的位移计算时间，New代表基于RF信号的我们所提出的新算法的计算时间。从上表可以发现算法在GPU上的实现方式对计算速度有很大影响，先对原始RF信号进行希尔伯特变换再分段互相关比先分段再互相关最后进行希尔伯特变化速度有明显提升，大概为6倍。Doppler无论在CPU上还是在GPU上，速度都比其它算法慢很多。本文提出的新方法，相对CPU速度提升了近58倍，并且无论在CPU还是在GPU上速度均是最快的。

图2-1、2-2、2-3分别是NCC_1、New、Doppler算法对图2中区域1成像（硬度分布不均匀）的比较;图2-1、2-2、2-3是对区域2（弹性分布稳定且较软）的成像图。图3给出了各算法在图1中区域各个推动位置测得的具体弹性模量值。从分布来看NCC_1与New比较吻合，而Doppler差别较明显，这主要是由于我们程序中使用了低通滤波器，但是很难达到完全理想的滤波效果，这就造成信号波形和相位的偏移和扭曲，进而影响估算Doppler频率和组织位移的准确程度。NCC_1、New、Doppler算法对区域2所得弹性图像的SNRe分别为：8.67、7.5、4.26。三种算法的对比度噪声比分别为：0.991、0.953、0.684。

从信噪比的结果来看，NCC和New算法所得应变图像都有较高的图像信噪比，而Doppler算法所得应变图像信噪比略低。并且前两种算法均比Doppler算法有更高的对比度噪声比值，即NCC_1和New有更好的弹性检测能力。

图3 各算法对区域1各推动位置的弹性值

5.结论

本文通过对传统超声弹性成像位移算法的研究以及结合声辐射力位移的特点，提出了一种特别适合于GPU加速的快速组织位移估计新算法。并介绍了基于解析信号的互相关相位改进算法、Kasais的Doppler运动估计算法和本文提出的算法在GPU上实现的一些技巧。最后通过实验表明本文提出的新算法能取得与常用算法相似的应变图像信噪比以及对比度噪声比，但是速度却提高了近100倍。除此以外本文提出的新算法并不局限于被应用在基于声辐射力的二维超声弹性成像系统，也可应用于任何需要利用超声回波射频信号进行组织位移估计的设备或方法中，比如准静态超声弹性成像（quasi-static elastography）、瞬态超声弹性成像（transient elastography）等其他超声弹性成像方法，或者超声温度成像（temperature imaging）、超声热应变成像（thermal strain imaging）等方法。

参考文献

[1]彭博，谌勇，流动权.基于GPU的超声弹性成像并行研究[J]，光电工程，2013，40（5）：97-105.

[2]CUDA Tookit Document.NVIDIA.

[3]RICHARD C.CABOT.A Note on the Application of the Hilbert Transform to time delay estimation[J].IEEE TRANSICTION ACOUSTICS，SPEECH，AND SIGNAL PROCESSINF，1981，29（3）：607-609.

[4]SCHWARTZ M，BENNRTT W R，STEIN S.Communicaion Systems and Techniques[M].New-York：McGraw-Hill，1965.

[5]PESAVENTO A，PERREY C，KRUEGER M，et at.A Time-Efficient and Accurate Strain Estimation Concept for Ultrasonic Elastography Using Iterative Phase Zero Estimation[J].IEEE TRANSAXTIONS ON ULTRASONICS，FERROELECTICS，AND FREQUENCY CONTROL，1999，46（5）：1057-1067.

[6]DE C，NAMEKAWA K，KOYANO K，et al.Real-Time two-dimensionalblood flow imaging using autocorrelation technique[J].IEEETRANSACTION SONICS ULTRASON，1985，42（4）：458-463.

[7]万明习等.生物医学超声学[M].北京：科学出版社，2010：

483-486.

作者简介：

谭勇（1985—），男，成都理工大学硕士研究生在读，研究方向：嵌入式系统研究。

通讯作者：王丛知（1977—），男，博士，中国科学院深圳先进技术研究院助理研究员，研究方向：生物医学信号处理。