CN103431874B - 声辐射力脉冲成像估算方法和系统 - Google Patents

Abstract

本发明涉及一种声辐射力脉冲成像估算方法和系统。所述方法包括：读取输入数据，并将所述输入数据分成预定数量数据矩阵；分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；对所述解析数据进行互相关运算，求取得到预定数量位移矩阵；对每个位移矩形累加，取累积和最大的列位移，生成新的位移矩阵；对得到的新的位移矩阵进行时间方向上的三次样条插值；对插值结果进行拉东变换，得到拉东变换矩阵；获取所述声辐射力脉冲成像的参数，根据所述参数及拉东变换矩阵求取剪切波速率。减小了互相关计算位移时产生的误差，计算结果更加准确，且对微小位移能精确估算，提高了灵敏度。

Description

声辐射力脉冲成像估算方法和系统

技术领域

本发明涉及医学影像领域，特别是涉及一种声辐射力脉冲成像估算方法和系统。

背景技术

声辐射力脉冲成像（Acoustic Radiation Force Impulse，ARFI）是通过利用医学超声功率范围内的聚焦超声波束产生的辐射力，使生物粘弹性组织局部区域产生微小变形并形成沿侧向传播的剪切波，利用超声成像技术对微小形变进行追踪并计算出剪切波的侧向传播速度，最后利用弹性重构算法由剪切波速度估算出组织的弹性分布的一种成像技术。

声辐射力脉冲成像计算是指估算组织内部剪切波的位移，后期通过对位移的判断了解组织情况。位移估算方法有多种，例如基于通用GPU（GraphicsProcess Unit，图形处理器）编程的实时超声成像监测方法，步骤包括：1）展开数据；2）加窗数据；3）互相关处理数据；4）估计位移；5）在整个时间深度上重复步骤3和4；6）中值滤波；7）估计速率。该监测方法使用了普通的基于时域互相关算法，针对微小位移的估算灵敏度不高，准确度低。

发明内容

基于此，有必要针对微小位移估算灵敏度低的问题，提供一种灵敏度高的声辐射力脉冲成像估算方法。

此外，还有必要提供一种灵敏度高的声辐射力脉冲成像估算系统。

一种声辐射力脉冲成像估算方法，包括：

数据划分步骤，读取通过声辐射力脉冲成像所采集的输入数据，并将所述输入数据分成预定数量数据矩阵；

数据滤波变换步骤，分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；

位移求取步骤，对所述解析数据进行互相关运算，求取得到预定数量位移矩阵；

位移合并步骤，对每个位移矩形进行时间累加，取累积和最大的列位移，再将预定数量位移矩阵抽取的列位移合并，生成新的位移矩阵；

插值步骤，对得到的新的位移矩阵进行时间方向上的三次样条插值；

拟合步骤，对插值结果进行拉东变换，得到拉东变换矩阵；

估算步骤，获取所述声辐射力脉冲成像的参数，根据所述参数及拉东变换矩阵求取剪切波速率。

一种声辐射力脉冲成像估算系统，包括：

数据划分模块，用于读取通过声辐射力脉冲成像所采集的输入数据，并将所述输入数据分成预定数量矩阵；

数据滤波变换模块，用于分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；

位移求取模块，用于对所述解析数据进行互相关运算，求取得到预定数量位移矩阵；

位移合并模块，用于对每个位移矩形进行时间累加，取累积和最大的列位移，再将预定数量位移矩阵抽取的列位移合并，生成新的位移矩阵；

插值模块，用于对得到的新的位移矩阵进行时间方向上的三次样条插值；

拟合模块，用于对插值结果进行拉东变换，得到拉东变换矩阵；

估算模块，用于获取所述声辐射力脉冲成像的参数，根据所述参数及拉东变换矩阵求取剪切波速率。

上述声辐射力脉冲成像估算方法，通过对输入数据进行前向滤波处理和逆向滤波处理，消除了滤波造成的相对偏移，并对滤波结果进行希尔伯特变换，减小了后续互相关计算位移时产生的误差，使得计算结果更加准确，且对微小位移能精确估算，提高了灵敏度。

另外，估算过程中采用CUDA并行计算实现，提高了效率，且在并行实现上，不增加重复计算量，不使用线程同步和共享内存，实现线程高度并行，提高了处理效率。

附图说明

图1为一个实施例中声辐射力脉冲成像估算方法的流程图；

图2为该位移求取步骤的流程图；

图3为拉动变换示意图；

图4为一个实施例中声辐射力脉冲成像估算系统的结构框图；

图5为该位移求取模块的内部结构示意图。

具体实施方式

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

如图1所示，为一个实施例中声辐射力脉冲成像估算方法的流程图。该声辐射力脉冲成像估算方法可运行于计算机上。该声辐射力脉冲成像估算方法，包括：

步骤102，数据划分步骤，读取通过声辐射力脉冲成像所采集的输入数据，并将所述输入数据分成预定数量数据矩阵。

具体的，将输入数据按侧向采集位置分为m个数据矩阵。m为侧向采集位置的条数。每个数据矩阵大小为row×col，表示包含row帧数据（时间方向），每帧有col个元素（深度方向）。

步骤104，数据滤波变换步骤，分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据。

具体的，对数据矩阵进行前向滤波处理和逆向滤波处理，可消除滤波造成的相对偏移。

本实施例中，1）该数据滤波变换步骤中前向滤波处理的公式为：

当i≤iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{i}{\mathrm{\Σ}}}\left(\mathrm{param}\[j\]*\mathrm{input}(k,i-j)\right)+\underset{j=i+1}{\overset{\mathrm{iParaLen}-2}{\mathrm{\Σ}}}\left(\mathrm{param}{\[j\]}^{*}\mathrm{input}(k,0)\right)---\left(1\right)$

当i＞iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{\mathrm{iParaLen}}{\mathrm{\Σ}}}\left(\mathrm{param}\[j\]*\mathrm{input}(k,i-j)\right)---\left(2\right)$

其中，param是滤波器，iParalen是滤波器长度，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，tOutPut是滤波结果矩阵，*表示乘法。

该前向滤波处理采用图形处理器（GPU，Graphics Process Unit）并行处理。该前向滤波处理的内核函数Bandpass_front（即图形处理器调用的内核函数）调用的线程数设计：Bandpass_front<<<blocks，threads>>>(input，aram，iParaLen，tOutPut)，此处blocks=dim3(input->rows)；threads=dim3(input->cols);其中，input->rows是输入数据矩阵的行数，input->cols是输入数据矩阵的列数。该输入数据矩阵即为待滤波矩阵。前向滤波器参数由Matlab等工具计算求得。其中，<<<、>>>是CUDA语法符号，表示GPU启动的线程，<<<表示启动block数量，>>>表示block里面线程数。input为一个结构类型，->表示取结构的一个元素。

GPU为面向计算密集型处理，适用数据并行化计算，GPU与CPU（CentralProcessing Unit，中央处理器）相比具有更多的计算单元晶体管和处理核心，可实现高度并行化计算，如Kepler GK110由71亿个晶体管组成，拥有15个流处理器核心，每个核心有196个流处理器，提供每秒超过1万亿次双精度浮点计算的吞吐量。

2）该数据滤波变换步骤中逆向滤波处理的公式为：

当i≥iParalen时：

$\mathrm{output}(k,\mathrm{Width}-1-i)=\underset{j=0}{\overset{\mathrm{iParalen}-1}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\&lsqb;j\&rsqb;*\mathrm{input}(k,\mathrm{Width}-1+j-i)\right)---\left(3\right)$

其中，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，output为逆向滤波的滤波结果矩阵，Width为滤波器宽度，*表示乘法。output的初始内容与input相同。

该逆向滤波处理采用图形处理器处理。该逆向滤波处理的内核函数Bandpass_back（即图形处理器调用的内核函数）调用的线程数设计：Bandpass_back<<<blocks，threads>>>()，其中，blocks=dim3(input->rows);threads=dim3(input->cols)，滤波器和前向滤波器相同，input->rows是输入数据矩阵的行数，input->cols是输入数据矩阵的列数。该输入数据矩阵即为待滤波矩阵。内核函数Bandpass_back将一个元素按照公式（3）处理，将获取的结果存放在output对应的位置。

3）该数据滤波变换步骤中希尔伯特变换的公式为：

r^(m,t)＝r(m,t)*(1/(πt))    （4）

式（4）中，r^(m,t)为信号r(m,t)的希尔伯特变换，m为采样点的坐标，t是时间，*是卷积运算。

在算法上采用傅里叶变换和逆变换来实现希尔伯特变换，具体为：对初始的数据矩阵进行傅里叶变换，将变换结果的实部乘-1后和虚部交换，再进行傅里叶归一化逆变换，再交换实虚部，再和初始的矩阵相加得到解析数据。其中，归一化是指变换结果除以变换元素个数。

希尔伯特变换采用GPU优化设计，其中，傅里叶变换和逆变换均使用CUDA（即CPU+GPU编程模型）的CUFFT库函数cufftExecC2C。

其中，两次实虚部交换以及矩阵相加用2个内核函数来实现，第一个内核函数kernel_1实现第一次的傅里叶变换前的实部虚部交换，第二个内核函数kernel_2实现第二次实部虚部交换，并且交换后相加，因cufftExecC2C结果没有归一化处理，在第二次实虚部交换时候处理。

内核函数调用顺序：1）cufftExecC2C（正变换）；2）kernel_1；3）cufftExecC2C（逆变换）；4）kernel_2。

两次内核函数的blocks数量和threads数量相同，均为blocks=dim3(input->rows)；threads=dim3(input->cols)，input->rows是输入数据矩阵的行数，input->cols是输入数据矩阵的列数。

在一个实施例中，该数据滤波变换步骤中，采用多个流按顺序依次对预定数量数据矩阵进行前向滤波处理、逆向滤波处理和希尔伯特变换处理。具体的，因GPU计算的一个特点是待处理数据需要从内存传送到显存，结果数据需要从显存传回内存这两个数据传输过程。由CPU调用GPU函数接口把待处理数据传输到显存，GPU通过显存读取数据，计算数据，把结果存放显存，CPU再调用数据传输函数接口将显存数据传回内存，进行下一步逻辑处理。数据在CPU和GPU之间的传输过程可与GPU计算过程并行。可使用多个流（stream）实现数据传输与GPU计算过程并行处理。流与流之间是并行运行的，因此虽然GPU一次只能处理一个内核函数，而当多个流存在时，GPU在处理内核函数的同时，显存可以传输数据，从而提高设备的利用率。

步骤106，位移求取步骤，对该解析数据进行互相关运算，求取得到预定数量位移矩阵。

如图2所示，该位移求取步骤包括：

步骤202，选取窗口，该窗口长度为预定值。

具体的，需要选取合适的窗口，窗口太小会导致计算距离不准确，太大会增加不必要的计算量。本实施例中，窗口长度为M，M为正数。

步骤204，根据该窗口长度按步长滑动的方式对该解析数据计算每个窗口的互相关系数，并形成互相关系数矩阵。

本实施例中，采用加窗互相关求取位移，加窗互相关是指一次只针对长度为窗口长度的一段求互相关，窗口位置随步长滑动，直到数据末位为止。

每个窗口的互相关系数计算公式为：

$k(q,s)=\underset{u=M/2}{\overset{M/2-1}{\mathrm{\&Sigma;}}}r^(m+u,s)\stackrel{\&OverBar;}{r^(m+u+q,s+1)}$ （5）

式（5）中，k(q,s)是互相关系数，r^(m+u,s)是解析数据，是解析数据r^(m+u+q,s+1)求复共轭，s是采样时间，q为[-M,M]。

q每一次变化称一个lag（即变化间隔参数），lag的范围[-M,M]。当估算位移比较小时，lag取一个范围[-N,N]，N＜M。

采用内核函数ComputeXcross计算所有窗（WinNum个）中一次lag的互相关系数。WinNum=(InputWidth－WindowHW)/Step+1，其中，WindowHW是窗口宽度，步长为Step，InputWidth为输入数据矩阵的列数，为了减少重复计算提高时间效率，WinNum个互相关系数计算方法如下：计算数据矩阵相邻两行Row(i)，Row(i+1)互相关系数，则先计算第0个窗口的互相关系数。第n个窗口的互相关系数计算为：第n-1个窗口的互相关系数减去第n-1个窗口前Step个元素的互相关系数，加上第n个窗口最后Step个元素的互相关系数。例如，第一个窗口的互相关系数为：第0个窗口的互相关系数减去第0个窗口前Step个元素的互相关系数，加上第一个窗口最后Step个元素的互相关系数；第二个窗口的互相关系数：第一个窗口的互相关系数减去第一个窗前Step个元素的互相关系数，加上第二个窗口最后Step个元素的互相关系数，依次计算到最后一个窗口为止。

内核函数ComputeXcross计算1个lag互相关系数，设计1个block的线程数等于lag总数，可以计算出相邻两行所有窗口的互相关系数。一个矩阵需要计算（rows－1）次加窗互相关（rows为矩阵行数），则由（rows-1）个block实现，每个block计算一个相邻两行的互相关系数。因此内核函数线程数如下：内核函数的线程数ComputeXcross<<<blocks，threads>>>()，其中，blocks=dim3(rows-1);threads=dim3(lagNum)；lagNum是lag（变化间隔参数）总数。

假设一个窗口的互相关系数（有lagNum个互相关系数）用Cij表示，通过上述互相关计算后，得到互相关系数矩阵包含（rows-1）*WinNum个Cij，i＝0，1...rows-2,j＝0,1...WinNum-1。

步骤206，获取该互相关系数矩阵中最大互相关系数对应的变化间隔参数值。

具体的，采用内核函数MatDelay对互相关系数矩阵求极值和位移。求一个Cij内最大相关系数对应的lag值lag_max。

步骤208，根据该变化间隔参数值求取位移。

具体的，由lag_max求位移delay。delay的求取方法如公式（6）：

delay＝（lag_max+delta_t）*factor    （6）

其中，delta_t采用公式（7）获取。factor＝c/(2*fs)，c是声速1540m/s，fs采样频率，*表示乘法。

当微小位移引起的超声回波信号时移δt小于采样周期时，可以通过互相关系数求出，公式如下：

$\mathrm{\&delta;t}=\frac{2\&angle;k^(n_{\max },s)}{\&angle;k^(n_{\max }+1,s)-\&angle;k^(n_{\max }-1,s)}---\left(7\right)$

式（7）中，n_max是最大互相关系数对应的坐标，∠是角度运算，即∠(a+jb)＝arctan(b/a)，得到了时移值δt。本实施例中，delta_t即为时移δt可通过式（7）计算得到。

通过上述计算得到位移矩阵。输入数据矩阵大小为row×col，则得到的位移矩阵为(row-1)×WinNum，WinNum为窗口数。

此外，微小位移d(m,s)可通过超声波速c求得:

$d(m,s)=\mathrm{\&delta;t}\frac{c}{2}---\left(8\right)$

本实施例中，内核函数线程数：blocks=dim3(rows-1)，threads=dim3(WinNum)，MatDelay<<<blocks,threads>>>(…)，其中WinNum窗个数，通过这些线程一次计算后，所有窗都计算了位移delay。

互相关系数求位移可采用CUDA的流方式求取。该CUDA的流方式可以隐藏数据传输时间，进一步缩短程序运行时间。例如在求位移的流处理示例如下：

同样，可采用该流的方式进行前向滤波处理、逆向滤波处理和希尔伯特变换。

步骤108，位移合并步骤，对每个位移矩形进行时间累加，取累积和最大的列位移，再将预定数量位移矩阵抽取的列位移合并，生成新的位移矩阵。

具体的，对得到的m个位移矩阵分别进行处理：进行时间（列方向）累加，取累积和最大的列位移，再将m个位移矩阵抽取的列位移合并，生成新的位移矩阵，即m×(row-1)。

步骤110，插值步骤，对得到的新的位移矩阵进行时间方向上的三次样条插值。

该插值步骤包括：采用中央处理器求取每个插值段之间的三次样条插值函数的系数，再采用图形处理器根据该三次样条插值函数并行计算每段的插值倍数减1个插值点的值，其中，图形处理器所调用的内核函数的blocks数量为待插值矩阵的行数，threads数量为待插值矩阵的列数。三次样条插值倍数可根据需要设定。用三次样条插值方法在相邻两点间补充p个点。本实施例中，采用的插值倍数为20，即相邻两点中间补充19个点，p=19。

具体的，插值段[x_i,x_i+1]的三次样条插值函数s_i(x)如下：

s_i(x)＝a_i(x-x_i)³+b_i(x-x_i)²+c_i(x-x_i)+d_i（9）

式（9）中，a_i、b_i、c_i、d_i为插值函数系数，可通过解线性方程获得，因线性方程是不定解方程，可以采用自然插值（Natural splines）确定方程解。

求插值系数的线性方程矩阵为三对角矩阵，直接LU分解求解。由于LU直接分解是不断迭代过程，采用CPU计算出全部插值函数系数，再用GPU并行计算出所有插值点的值。

设计内核函数spline_kern计算一个[x_i,x_i+1]插值段的所有值。对输入数据矩阵每行插值，需要线程数如下：spline_kern<<<blocks,threads>>>()，其中，blocks=dim3(inputMat->rows)；threads=dim3(inputMat->cols-1)，inputMat->rows是待插值矩阵行数，inputMat->cols是待插值矩阵的列数。

步骤112，拟合步骤，对插值结果进行拉东变换，得到拉东变换矩阵。

如图3所示，为拉动变换示意图。图3可知，依次求一点到对端一点的累加和，考虑实际应用，终点的序号不能小于起点的序号。假设矩阵宽度为Q，需要计算Q*(Q+1)/2次斜对角线的和。为了存取数据方便，通常grid和block的维数直接用矩阵的长宽相关值，这里经过插值后的输入矩阵，Q的值通常较大比如1960，那么拉东变换矩阵大小为1960*1960，而计算能力为3.5的GPU支持block内部最大线程数是1024，因此block内部线程数不能用Q定义。

定义grid（X*X），block（Y*Y）。其中，Y取16或者32，Y取满足X*Y>=Q最小值即可。索引值按如下计算：

index_start=blockIdx.x*blockDim.x+threadIdx.x；

index_end=blockIdx.y*blockDim.y+threadIdx.y;

其中，*表示乘法。

即grid的x维和block的x维组合表示起始端点的序号index_start（0到1959），grid的y维和block的y维组合表示终点端的序号index_end。

内核函数实现一个斜对角求和（从初始点index_start，到终点index_end的对角连线上的数据求和）功能，得到的结果Sum（index_start，index_end）保存到输出数据中。

步骤114，估算步骤，获取该声辐射力脉冲成像的参数，根据该参数及拉东变换矩阵求取剪切波速率。

该估算步骤中该声辐射力脉冲成像的参数包括探头中两个阵元的间距、探头个数、输入数据矩阵的行数、帧频和采样倍数。

该估算步骤中求取剪切波速率的公式为：

vs＝pitch*pitchNum*(dataRow-1)/(L*fs*spline)

其中，vs为剪切波速率，pitch是探头中两个阵元的间距，pitchNum是探头个数，dataRow是输入数据矩阵的行数，fs是帧频，spline是采样倍数，*表示乘法，L＝(Rmax).y-(Rmax).x，L是位移，Rmax是拉东变换矩阵最大值，L是最大值对应的y坐标减去最大值对应的x坐标。

上述声辐射力脉冲成像估算方法，通过对输入数据进行前向滤波处理和逆向滤波处理，消除了滤波造成的相对偏移，并对滤波结果进行希尔伯特变换，减小了后续互相关计算位移时产生的误差，使得计算结果更加准确，且对微小位移能精确估算，提高了灵敏度。

另外，估算过程中采用CUDA并行计算实现，提高了效率，且在并行实现上，不增加重复计算量，不使用线程同步和共享内存，实现线程高度并行，提高了处理效率。

如图4所示，为一个实施例中声辐射力脉冲成像估算系统的结构框图。该声辐射力脉冲成像估算系统中未详尽描述的细节，可参考上述对应的方法实施例描述。该声辐射力脉冲成像估算系统，包括数据划分模块410、数据滤波变换模块420、位移求取模块430、位移合并模块440、插值模块450、拟合模块460和估算模块470。其中：

数据划分模块410用于读取通过声辐射力脉冲成像所采集的输入数据，并将该输入数据分成预定数量矩阵。

具体的，将输入数据按侧向采集位置分为m个数据矩阵。m为侧向采集位置的条数。每个数据矩阵大小为row×col，表示包含row帧数据（时间方向），每帧有col个元素（深度方向）。

数据滤波变换模块420用于分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据。

本实施例中，1）该数据滤波变换模块420进行前向滤波处理的公式为：

当i≤iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{i}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\&lsqb;j\&rsqb;*\mathrm{input}(k,i-j)\right)+\underset{j=i+1}{\overset{\mathrm{iParalen}-2}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\&lsqb;j\&rsqb;*\mathrm{input}(k,0)\right)$

当i＞iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{\mathrm{iParalen}-2}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\&lsqb;j\&rsqb;*\mathrm{input}(k,i-j)\right)$

其中，param是滤波器，iParalen是滤波器长度，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，tOutPut是滤波结果矩阵，*表示乘法。

该前向滤波处理采用GPU并行处理。该前向滤波处理的内核函数Bandpass_front（即图形处理器调用的内核函数）调用的线程数设计：Bandpass_front<<<blocks，threads>>>(input，aram，iParaLen，tOutPut)，此处blocks=dim3(input->rows)；threads=dim3(input->cols);其中，input->rows是输入数据矩阵的行数，input->cols是输入数据矩阵的列数。该输入数据矩阵即为待滤波矩阵。

2）该数据滤波变换模块420进行逆向滤波处理的公式为：

当i≥iParalen时：

$\mathrm{output}(k,\mathrm{Width}-1-i)=\underset{j=0}{\overset{\mathrm{iParalen}-1}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\&lsqb;j\&rsqb;*\mathrm{input}(k,\mathrm{Width}-1+j-i)\right)$

其中，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，output为逆向滤波的滤波结果矩阵，Width为滤波器宽度，*表示乘法。

该逆向滤波处理采用图形处理器处理。该逆向滤波处理的内核函数Bandpass_back（即图形处理器调用的内核函数）调用的线程数设计：Bandpass_back<<<blocks，threads>>>()，其中，blocks=dim3(input->rows);threads=dim3(input->cols)，滤波器和前向滤波器相同，input->rows是输入数据矩阵的行数，input->cols是输入数据矩阵的列数。该输入数据矩阵即为待滤波矩阵。内核函数Bandpass_back将一个元素按照公式（3）处理，将获取的结果存放在output对应的位置。

3）该数据滤波变换模块420进行希尔伯特变换公式为：

r^(m,t)＝r(m,t)*(1/(πt))    （4）

式（4）中，r^(m,t)为信号r(m,t)的希尔伯特变换，m为采样点的坐标，t是时间，*是卷积运算。

在算法上采用傅里叶变换和逆变换来实现希尔伯特变换，具体为：对初始的矩阵进行傅里叶变换，将变换结果的实部乘-1后和虚部交换，再进行傅里叶归一化逆变换，再交换实虚部，再和初始的矩阵相加得到解析数据。

希尔伯特变换采用GPU优化设计，其中，傅里叶变换和逆变换均使用CUDA（即CPU+GPU编程模型）的CUFFT库函数cufftExecC2C。

其中，两次实虚部交换以及矩阵相加用2个内核函数来实现，第一个内核函数kernel_1实现第一次的傅里叶变换前的实部虚部交换，第二个内核函数kernel_2实现第二次实部虚部交换，并且交换后相加，因cufftExecC2C结果没有归一化处理，在第二次实虚部交换时候处理。

内核函数调用顺序：1）cufftExecC2C（正变换）；2）kernel_1；3）cufftExecC2C（逆变换）；4）kernel_2。

两次内核函数的blocks数量和threads数量相同，均为blocks=dim3(input->rows)；threads=dim3(input->cols)，input->rows是输入数据矩阵的行数，input->cols是输入数据矩阵的列数。

在一个实施例中，该数据滤波变换模块420还采用多个流按顺序依次对预定数量数据矩阵进行前向滤波处理、逆向滤波处理和希尔伯特变换处理。具体的，因GPU计算的一个特点是待处理数据需要从内存传送到显存，结果数据需要从显存传回内存这两个数据传输过程。由CPU调用GPU函数接口把待处理数据传输到显存，GPU通过显存读取数据，计算数据，把结果存放显存，CPU再调用数据传输函数接口将显存数据传回内存，进行下一步逻辑处理。数据在CPU和GPU之间的传输过程可与GPU计算过程并行。可使用多个流（stream）实现数据传输与GPU计算过程并行处理。流与流之间是并行运行的，因此虽然GPU一次只能处理一个内核函数，而当多个流存在时，GPU在处理内核函数的同时，显存可以传输数据，从而提高设备的利用率。

位移求取模块430用于对该解析数据进行互相关运算，求取得到预定数量位移矩阵。

如图5所示，该位移求取模块430包括选取子模块432、互相关系数计算子模块434、比较子模块436和位移计算子模块438。其中：

选取子模块432用于选取窗口，该窗口长度为预定值。本实施例中，窗口长度为M，M为正数。

互相关系数计算子模块434用于根据该窗口长度按步长滑动的方式对该解析数据计算每个窗口的互相关系数，并形成互相关系数矩阵。

本实施例中，采用加窗互相关求取位移，加窗互相关是指一次只针对长度为窗口长度的一段求互相关，窗口位置随步长滑动，直到数据末位为止。

每个窗口的互相关系数计算公式为：

$k(q,s)=\underset{u=M/2}{\overset{M/2-1}{\mathrm{\&Sigma;}}}r^(m+u,s)\stackrel{\&OverBar;}{r^(m+u+q,s+1)}$ （5）

式（5）中，k(q,s)是互相关系数，r^(m+u,s)是解析数据，是解析数据r^(m+u+q,s+1)求复共轭，s是采样时间，q为[-M,M]。

q每一次变化称一个lag（即变化间隔参数），lag的范围[-M,M]。当估算位移比较小时，lag取一个范围[-N,N]，N＜M。

采用内核函数ComputeXcross计算所有窗（WinNum个）中一次lag的互相关系数。WinNum=(InputWidth－WindowHW)/Step+1，其中，WindowHW是窗口宽度，步长为Step，InputWidth为输入数据矩阵的列数，为了减少重复计算提高时间效率，WinNum个互相关系数计算方法如下：计算数据矩阵相邻两行Row(i)，Row(i+1)互相关系数，则先计算第0个窗口的互相关系数。第n个窗口的互相关系数计算为：第n-1个窗口的互相关系数减去第n-1个窗口前Step个元素的互相关系数，加上第n个窗口最后Step个元素的互相关系数。例如，第一个窗口的互相关系数为：第0个窗口的互相关系数减去第0个窗口前Step个元素的互相关系数，加上第一个窗口最后Step个元素的互相关系数；第二个窗口的互相关系数：第一个窗口的互相关系数减去第一个窗前Step个元素的互相关系数，加上第二个窗口最后Step个元素的互相关系数，依次计算到最后一个窗口为止。

内核函数ComputeXcross计算1个lag互相关系数，设计1个block的线程数等于lag总数，可以计算出相邻两行所有窗口的互相关系数。一个矩阵需要计算（rows－1）次加窗互相关（rows为矩阵行数），则由（rows-1）个block实现，每个block计算一个相邻两行的互相关系数。因此内核函数线程数如下：内核函数的线程数ComputeXcross<<<blocks，threads>>>()，其中，blocks=dim3(rows-1);threads=dim3(lagNum)；lagNum是lag（变化间隔参数）总数。

假设一个窗口的互相关系数（有lagNum个互相关系数）用Cij表示，通过上述互相关计算后，得到互相关系数矩阵包含（rows-1）*WinNum个Cij，i＝0，1...rows-2,j＝0,1...WinNum-1。

比较子模块436用于获取该互相关系数矩阵中最大互相关系数对应的变化间隔参数值。具体的，采用内核函数MatDelay对互相关系数矩阵求极值和位移。求一个Cij内最大相关系数对应的lag值lag_max。

位移计算子模块438用于根据该变化间隔参数值求取位移。

具体的，由lag_max求位移delay。delay的求取方法如公式（6）：

delay＝（lag_max+delta_t）*factor    （6）

其中，delta_t采用公式（7）获取。factor＝c/(2*fs)，c是声速1540m/s，fs采样频率，*表示乘法。

当微小位移引起的超声回波信号时移δt小于采样周期时，可以通过互相关系数求出，公式如下：

$\mathrm{\&delta;t}=\frac{2\&angle;k^(n_{\max },s)}{\&angle;k^(n_{\max }+1,s)-\&angle;k^(n_{\max }-1,s)}---\left(7\right)$

式（7）中，n_max是最大互相关系数对应的坐标，∠是角度运算，即∠(a+jb)＝arctan(b/a)，得到了时移值δt。本实施例中，delta_t即为时移δt可通过式（7）计算得到。

通过上述计算得到位移矩阵。输入数据矩阵大小为row×col，则得到的位移矩阵为(row-1)×WinNum，WinNum为窗口数。

此外，微小位移d(m,s)可通过超声波速c求得:

$d(m,s)=\mathrm{\&delta;t}\frac{c}{2}---\left(8\right)$

本实施例中，内核函数线程数：blocks=dim3(rows-1)，threads=dim3(WinNum)，MatDelay<<<blocks,threads>>>()，其中WinNum窗个数，通过这些线程一次计算后，所有窗都计算了位移delay。

互相关系数求位移可采用CUDA的流方式求取。

位移合并模块440用于对每个位移矩形进行时间累加，取累积和最大的列位移，再将预定数量位移矩阵抽取的列位移合并，生成新的位移矩阵。具体的，对得到的m个位移矩阵分别进行处理：进行时间（列方向）累加，取累积和最大的列位移，再将m个位移矩阵抽取的列位移合并，生成新的位移矩阵，即m×(row-1)。

插值模块450用于对得到的新的位移矩阵进行时间方向上的三次样条插值。

本实施例中，该插值模块450还用于采用中央处理器求取每个插值段之间的三次样条插值函数的系数，再采用图形处理器根据该三次样条插值函数并行计算每段的插值倍数减1个插值点的值，其中，图形处理器所调用的内核函数的blocks数量为待插值矩阵的行数，threads数量为待插值矩阵的列数。三次样条插值倍数可根据需要设定。用三次样条插值方法在相邻两点间补充p个点。本实施例中，采用的插值倍数为20，即相邻两点中间补充19个点，p=19。

具体的，插值段[x_i,x_i+1]的三次样条插值函数s_i(x)如下：

s_i(x)＝a_i(x-x_i)³+b_i(x-x_i)²+c_i(x-x_i)+d_i    （9）

式（9）中，a_i、b_i、c_i、d_i为插值函数系数，可通过解线性方程获得，因线性方程是不定解方程，可以采用自然插值（Natural splines）确定方程解。

求插值系数的线性方程矩阵为三对角矩阵，直接LU分解求解。由于LU直接分解是不断迭代过程，采用CPU计算出全部插值函数系数，再用GPU并行计算出所有插值点的值。

设计内核函数spline_kern计算一个[x_i,x_i+1]插值段的所有值。对输入数据矩阵每行插值，需要线程数如下：spline_kern<<<blocks,threads>>>()，其中，blocks=dim3(inputMat->rows)；threads=dim3(inputMat->cols-1)，inputMat->rows是待插值矩阵行数，inputMat->cols是待插值矩阵的列数。

拟合模块460用于对插值结果进行拉东变换，得到拉东变换矩阵。

如图3所示，为拉动变换示意图。图3可知，依次求一点到对端一点的累加和，考虑实际应用，终点的序号不能小于起点的序号。假设矩阵宽度为Q，需要计算Q*(Q+1)/2次斜对角线的和。为了存取数据方便，通常grid和block的维数直接用矩阵的长宽相关值，这里经过插值后的输入矩阵，Q的值通常较大比如1960，那么拉东变换矩阵大小为1960*1960，而计算能力为3.5的GPU支持block内部最大线程数是1024，因此block内部线程数不能用Q定义。

定义grid（X*X），block（Y*Y）。其中，Y取16或者32，Y取满足X*Y>=Q最小值即可。索引值按如下计算：

index_start=blockIdx.x*blockDim.x+threadIdx.x；

index_end=blockIdx.y*blockDim.y+threadIdx.y;

其中，*表示乘法。

即grid的x维和block的x维组合表示起始端点的序号index_start（0到1959），grid的y维和block的y维组合表示终点端的序号index_end。

内核函数实现一个斜对角求和（从初始点index_start，到终点index_end的对角连线上的数据求和）功能，得到的结果Sum（index_start，index_end）保存到输出数据中。

估算模块470用于获取该声辐射力脉冲成像的参数，根据该参数及拉东变换矩阵求取剪切波速率。

本实施例中，该声辐射力脉冲成像的参数包括探头中两个阵元的间距、探头个数、输入数据矩阵的行数、帧频和采样倍数；

该估算模块用于求取剪切波速率的公式为：

vs＝pitch*pitchNum*(dataRow-1)/(L*fs*spline)

其中，vs为剪切波速率，pitch是探头中两个阵元的间距，pitchNum是探头个数，dataRow是输入数据矩阵的行数，fs是帧频，spline是采样倍数，*表示乘法，L＝(Rmax).y-(Rmax).x，L是位移，Rmax是拉东变换矩阵最大值，L是最大值对应的y坐标减去最大值对应的x坐标。

上述声辐射力脉冲成像估算系统，通过对输入数据进行前向滤波处理和逆向滤波处理，消除了滤波造成的相对偏移，并对滤波结果进行希尔伯特变换，减小了后续互相关计算位移时产生的误差，使得计算结果更加准确，且对微小位移能精确估算，提高了灵敏度。

另外，估算过程中采用CUDA并行计算实现，提高了效率，且在并行实现上，不增加重复计算量，不使用线程同步和共享内存，实现线程高度并行，提高了处理效率。

上述声辐射力脉冲成像估算方法和系统，模拟结果：实验中给定仿体的剪切波速度为3m/s（米/秒），计算结果为3.002m/s，误差仅为0.002m/s。

GPU优化效果：本估算方法以4组100*512的采集数据作为一组输入数据，设置窗口长度100，lag总数为81，一次计算耗时为80ms。该耗时时间包括了GPU与CPU之间数据传输时间和GPU计算时间，不包括GPU申请显存时间。由于GPU申请显存比较慢（申请足够的内存用时大约500ms），考虑到申请的显存可以重复使用，可在系统启动时一次性申请足够显存，在系统关闭时释放；申请内存时间不作为运行时间计算是合理的。对比GPU优化前，利用配置为双核主频3.3G、3.29G内存2G的个人计算机运行算法，耗时为2000ms左右。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（Read-Only Memory，ROM）或随机存储记忆体（Random Access Memory，RAM）等。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (12)

1.一种声辐射力脉冲成像估算方法，包括：

数据划分步骤，读取通过声辐射力脉冲成像所采集的输入数据，并将所述输入数据分成预定数量数据矩阵；

数据滤波变换步骤，分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；

位移求取步骤，对所述解析数据进行互相关运算，求取得到预定数量位移矩阵；

位移合并步骤，对每个位移矩形进行时间累加，取累加和最大的列位移，再将预定数量位移矩阵抽取的列位移合并，生成新的位移矩阵；

插值步骤，对得到的新的位移矩阵进行时间方向上的三次样条插值；

拟合步骤，对插值结果进行拉东变换，得到拉东变换矩阵；

估算步骤，获取所述声辐射力脉冲成像的参数，根据所述参数及拉东变换矩阵求取剪切波速率；

所述估算步骤中所述声辐射力脉冲成像的参数包括探头中两个阵元的间距、探头个数、输入数据矩阵的行数、帧频和采样倍数；

所述估算步骤中求取剪切波速率的公式为：

vs＝pitch*pitchNum*(dataRow-1)/(L*fs*spline)

其中，vs为剪切波速率，pitch是探头中两个阵元的间距，pitchNum是探头个数，dataRow是输入数据矩阵的行数，fs是帧频，spline是采样倍数，*表示乘法，L＝(Rmax).y-(Rmax).x，L是位移，Rmax是拉东变换矩阵最大值，L是最大值对应的y坐标减去最大值对应的x坐标。

2.根据权利要求1所述的声辐射力脉冲成像估算方法，其特征在于，所述数据滤波变换步骤中，采用多个流按顺序依次对预定数量数据矩阵进行前向滤波处理、逆向滤波处理和希尔伯特变换处理。

3.根据权利要求1所述的声辐射力脉冲成像估算方法，其特征在于，所述数据滤波变换步骤中前向滤波处理的公式为：

当i≤iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{i}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,i-j))+\underset{j=i+1}{\overset{\mathrm{iParalen}-2}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,0))$

当i＞iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{\mathrm{iParalen}-2}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,i-j))$

其中，param是滤波器，iParalen是滤波器长度，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，tOutPut是滤波结果矩阵；

所述数据滤波变换步骤中逆向滤波处理的公式为：

当i≥iParalen时：

$\mathrm{output}(k,\mathrm{Width}-1-i)=\underset{j=0}{\overset{\mathrm{iParalen}-1}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,\mathrm{Width}-1+j-i))$

其中，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，output为逆向滤波的滤波结果矩阵，Width为滤波器宽度；

所述前向滤波处理的内核函数和逆向滤波处理均采用图形处理器并行处理，所述图形处理器调用的内核函数的blocks数量和threads数量相同，所述blocks数量为待滤波矩阵的行数，所述threads数量为待滤波矩阵的列数。

4.根据权利要求1所述的声辐射力脉冲成像估算方法，其特征在于，所述数据滤波变换步骤中希尔伯特变换包括：

对初始的数据矩阵进行傅里叶变换，将变换结果的实部乘-1后和虚部交换，再进行傅里叶归一化逆变换，再交换实虚部，再和初始的数据矩阵相加得到解析数据。

5.根据权利要求1所述的声辐射力脉冲成像估算方法，其特征在于，所述位移求取步骤包括：

选取窗口，所述窗口长度为预定值；

根据所述窗口长度按步长滑动的方式对所述解析数据计算每个窗口的互相关系数，并形成互相关系数矩阵；

获取所述互相关系数矩阵中最大互相关系数对应的变化间隔参数值；

根据所述变化间隔参数值求取位移。

6.根据权利要求1所述的声辐射力脉冲成像估算方法，其特征在于，所述插值步骤包括：

采用中央处理器求取每个插值段之间的三次样条插值函数的系数，再采用图形处理器根据所述三次样条插值函数并行计算每段的插值倍数减1个插值点的值，其中，图形处理器所调用的内核函数的blocks数量为待插值矩阵的行数，threads数量为待插值矩阵的列数。

7.一种声辐射力脉冲成像估算系统，其特征在于，包括：

数据划分模块，用于读取通过声辐射力脉冲成像所采集的输入数据，并将所述输入数据分成预定数量矩阵；

数据滤波变换模块，用于分别对每个数据矩阵进行前向滤波处理和逆向滤波处理，并对滤波处理后的数据矩阵进行希尔伯特变换得到解析数据；

位移求取模块，用于对所述解析数据进行互相关运算，求取得到预定数量位移矩阵；

位移合并模块，用于对每个位移矩形进行时间累加，取累加和最大的列位移，再将预定数量位移矩阵抽取的列位移合并，生成新的位移矩阵；

插值模块，用于对得到的新的位移矩阵进行时间方向上的三次样条插值；

拟合模块，用于对插值结果进行拉东变换，得到拉东变换矩阵；

估算模块，用于获取所述声辐射力脉冲成像的参数，根据所述参数及拉东变换矩阵求取剪切波速率；

所述声辐射力脉冲成像的参数包括探头中两个阵元的间距、探头个数、输入数据矩阵的行数、帧频和采样倍数；

所述估算模块用于求取剪切波速率的公式为：

vs＝pitch*pitchNum*(dataRow-1)/(L*fs*spline)

其中，vs为剪切波速率，pitch是探头中两个阵元的间距，pitchNum是探头个数，dataRow是输入数据矩阵的行数，fs是帧频，spline是采样倍数，*表示乘法，L＝(Rmax).y-(Rmax).x，L是位移，Rmax是拉东变换矩阵最大值，L是最大值对应的y坐标减去最大值对应的x坐标。

8.根据权利要求7所述的声辐射力脉冲成像估算系统，其特征在于，所述数据滤波变换模块还用于采用多个流按顺序依次对预定数量数据矩阵进行前向滤波处理、逆向滤波处理和希尔伯特变换处理。

9.根据权利要求7所述的声辐射力脉冲成像估算系统，其特征在于，所述数据滤波变换模块进行前向滤波处理的公式为：

当i≤iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{i}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,i-j))+\underset{j=i+1}{\overset{\mathrm{iParalen}-2}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,0))$

当i＞iParalen时：

$\mathrm{tOutPut}(k,i)=\underset{j=0}{\overset{\mathrm{iParalen}-2}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,i-j))$

其中，param是滤波器，iParalen是滤波器长度，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，tOutPut是滤波结果矩阵，*表示乘法；

所述数据滤波变换模块进行逆向滤波处理的公式为：

当i≥iParalen时：

$\mathrm{output}(k,\mathrm{Width}-1-i)=\underset{j=0}{\overset{\mathrm{iParalen}-1}{\mathrm{\&Sigma;}}}\left(\mathrm{param}\right[j]*\mathrm{input}(k,\mathrm{Width}-1+j-i))$

其中，input为待滤波矩阵，k、i分别表示待滤波矩阵的第k行和第i列，output为逆向滤波的滤波结果矩阵，Width为滤波器宽度，*表示乘法；

所述前向滤波处理的内核函数和逆向滤波处理均采用图形处理器并行处理，所述图形处理器调用的内核函数的blocks数量和threads数量相同，所述blocks数量为待滤波矩阵的行数，所述threads数量为待滤波矩阵的列数。

10.根据权利要求7所述的声辐射力脉冲成像估算系统，其特征在于，所述数据滤波变换模块进行希尔伯特变换包括：

对初始的数据矩阵进行傅里叶变换，将变换结果的实部乘-1后和虚部交换，再进行傅里叶归一化逆变换，再交换实虚部，再和初始的数据矩阵相加得到解析数据。

11.根据权利要求7所述的声辐射力脉冲成像估算系统，其特征在于，所述位移求取模块包括：

选取子模块，用于选取窗口，所述窗口长度为预定值；

互相关系数计算子模块，用于根据所述窗口长度按步长滑动的方式对所述解析数据计算每个窗口的互相关系数，并形成互相关系数矩阵；

比较子模块，用于获取所述互相关系数矩阵中最大互相关系数对应的变化间隔参数值；

位移计算子模块，用于根据所述变化间隔参数值求取位移。

12.根据权利要求7所述的声辐射力脉冲成像估算系统，其特征在于，所述插值模块还用于采用中央处理器求取每个插值段之间的三次样条插值函数的系数，再采用图形处理器根据所述三次样条插值函数并行计算每段的插值倍数减1个插值点的值，其中，图形处理器所调用的内核函数的blocks数为待插值矩阵的行数，threads数为待插值矩阵的列数。