技术领域
本发明涉及超声波回波成像领域，尤其是一种实时超声弹性成像方法和系 统。
背景技术
超声波回波成像技术目前已经被广泛应用于军事、医疗等领域。超声波回 波成像中关于弹性成像的概念最早由Ophir等人于1991年首次提出。之后，弹 性成像技术在近二十年里得到了迅猛发展，它被称为继A型、B型、D型、M型 超声之后的E型模式。超声弹性成像是通过超声成像系统进行生物组织弹性参 数成像，超声弹性图能够提供传统B超图像所无法反映的生物组织弹性特征， 对于肿瘤检测等临床应用有非常大的帮助。由于具有易于实现、适用于实时诊 断和对组织无侵害性等优点，超声弹性成像技术受到了广泛的关注。
现有的弹性成像方法中，计算位移场的方法大多基于互相关算法和改进的 互相关算法，其使用的数据为RF射频数据，而不能直接使用B模式图像，由于 处理RF射频数据需要很大的运算，因此造成了现有的方法计算效率低下的问题； 光流法是指空间运动物体在观察成像平面上的像素运动的瞬时速度，利用图像 序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧 之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。计 算应变的方法大多基于最小二乘拟合算法和梯度法，使用这些方法计算出的应 变，在精度方面还有提高的空间。
1998年，Negahdaripour将光流(OpticalFlow)定义为动态图像的几何 变化和辐射度变化的全面表示。光流法利用图像序列中的像素强度数据的时域 变化和相关性来确定各自像素位置的运动，表达了图像灰度在时间上的变化与 景象中物体结构及其运动的关系。用于评价光流法精度的一个经典标准来自于 Middlebury，其用于光流法评价的数据集有12张图像，除此之外还有Sintel 和KITTI。
光流法的基本原理如下：
令I(x,y,t)表示在t时刻图像上像素(x,y)处的亮度(或颜色)，那么光流法 的目的就是求出在t+1时刻的图像上，该像素相对于原来(x,y)的位移量(u,v)， 用方程表示即：
I(x+u,y+v,t+1)＝I(x,y,t)(1)
其中u和v是待求解的位移量。该方程被称为BrightnessConstancyModel (亮度恒定模型)。
在经典的光流法中，一般利用一阶泰勒展开作为工具来建立图像梯度和位 移之间的关系，这一步骤通常被称为线性化。具体原理如下：
假设图像的亮度是连续的，如图2的一维示例，曲线1表示frame1(帧1)中的图像，曲线2表示frame2(帧2)中的图像，待求的位移是向右的箭头对曲线进行一阶泰勒展开，其实就是假设曲线的局部是线性的，这样可以考察如图2的向右箭头、向上箭头、粗线段组成的三角形。并非是粗线段的长度，而是其斜率。这样可以得到图中所示的关系，注意负号是因为斜率其实表示的是钝角的tan值。这样一来就建立了图像的导数和位移之间的关系，注意是图像在空间上的导数，是图像在时间上的导数。
将图2中的方程用在二维的图像上，对于每个像素，可以写出以下方程：
Ixu+Iyv+It＝0(2)
其中，Ix和Iy是图像沿空间的x和y两个方向的导数，即图像梯度的两个 分量；It是图像沿时间变化的导数，可以用两帧图像的差异来近似；u和v是该 像素沿x和y方向的位移，也就是待求的光流未知量。该方程是Brightness ConstancyModel(亮度恒定模型)线性化的结果，被称为GradientConstraint Equation(梯度约束方程)。
需要注意的是，上面这个模型是建立在两个假设基础上：第一，图像变化 是连续的；第二，位移不是很大。如果这两点假设不成立，即图像不连续且位 移很大，则无法将位移和图像梯度联系起来，泰勒展开得到的近似结果将很差。 在实际操作中，可以很容易使上述两个假设其成立。关于第一点假设，一般可 以预先对图像进行高斯平滑，使其变化较为平缓；关于第二点假设，一般对图 像降分辨率建立金字塔，通过Coarse-To-Fine(由粗到精)的方式去求解。
基于上面的方程，产生了最经典的两个光流法：Lucas-Kanade方法和 Horn-Schunck方法，他们分别从不同的角度增加了求解该方程的稳定性。 Lucas-Kanade方法是将每个像素周围的一些像素考虑进来，每个像素的未知量 单独求解，是一种局部光流方法；而Horn-Schunck方法是将上面的方程纳入到 一个规则化的框架中，将局部平滑的优先级加入光流计算，所有像素的未知位 移量之间相互依赖，需要用全局优化的方法求解。
当前如何优化现有弹性成像系统，是技术人员需要考虑的研发方向，现有 超声系统大都使用CPU进行全局的运算，本发明运用GPU进行弹性成像系统中 图像处理运算，能够提高现有超声弹性成像系统的运算速度，减少CPU的负载 从而提高了原有系统的稳定性。同时，随着日益增长的诊断需求，如何快速准 确进行辅助医生进行疾病的诊断，也是医疗领域的发展方向。
发明内容
本发明的第一目的在于提供一种实时超声弹性成像方法，运算速度快，精 度较高，同时具备较强的鲁棒性。本发明还同时提出了一种实时超声弹性成像 系统。本发明采用的技术方案是：
一种实时超声弹性成像方法，包括下述步骤：
S10.对生物组织进行缓速挤压时对生物组织目标区域进行超声波扫查并接 收回波信号；
S20.对步骤S10中接收的回波信号进行处理，形成线数据；
S30.对步骤S20得到的线数据进行处理，形成生物组织受外力缓速挤压后 产生形变的不同状态下的B模式图像；
S40.对生物组织受外力缓速挤压后产生形变的不同状态下的两帧B模式图 像选定两个在不同帧图像的同一位置区域，即选定两帧图像的ROI区域；
S50.通过第一CPU数据处理模块对该两帧图像的ROI区域使用光流法求取 位移场；并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数 据缓存模块，进行缓存；
S60.第一GPU数据缓存模块将接收到的子步骤的参数信息传输给各自的GPU 数据处理模块进行数据处理，S50中的含重复计算的各子步骤的运算分配给至少 一个GPU核工作单元，由至少一个GPU核工作单元组成GPU工作组；
S70.将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存 中，等待所有GPU工作组的参数信息映射完成，然后各GPU工作组进行步骤S50 中的各子步骤的运算处理；
S80.将GPU工作组的运算结果信息传输至第二GPU数据缓存模块；
S90.第二GPU数据缓存模块将运算结果信息传输至第二CPU数据处理模块， 第二CPU数据处理模块读出运算结果信息并带入光流法运算过程，最终求出ROI 区域的位移场即光流场；
S100.使用低通滤波器对轴向光流场求取应变，得到图像ROI区域的轴向应 变场；
S105.对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理；
S110.对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物 组织彩色弹性图像。
所述步骤S20中，通过波束合成模块对回波信号具体进行A/D转换、同相 位叠加处理；
所述步骤S30中，通过B模式图像处理模块对步骤S20得到的线数据具体 进行降噪、滤波、提高信噪比处理。
进一步地，所述步骤S50中，含重复计算的子步骤的参数信息是指：构建 图像金字塔、中心梯度计算、计算目标图像及其导数的形变、估计双变量、估 计位移这些计算的参数信息。
进一步地，所述步骤S100中，低通滤波器的核如下列公式所示：
y ( M + k ) = 25 * ( 3 M 4 + 6 M 3 - 3 M + 1 ) * k - 35 * ( 3 M 2 + 3 M - 1 ) * k 3 ( 2 M + 3 ) * ( 2 M + 1 ) * ( 2 M - 1 ) * ( M + 2 ) * ( M + 1 ) * M * ( M - 1 ) y ( M - k ) = - 25 * ( 3 M 4 + 6 M 3 - 3 M + 1 ) * k - 35 * ( 3 M 2 + 3 M - 1 ) * k 3 ( 2 M + 3 ) * ( 2 M + 1 ) * ( 2 M - 1 ) * ( M + 2 ) * ( M - 1 ) * M * ( M - 1 ) y ( M ) = 0 ]]>
该滤波核的大小为N行*1列，N为正奇数,其中，M＝(N-1)/2，k为大于0 且不大于M的整数，1≤k≤M，y代表滤波核不同位置的值。
更进一步地，滤波核的N取值不大于31。
进一步地，所述步骤S110具体为：首先将由轴向应变场构成的灰度图像经 波段合成转换成彩色图像，再将该彩色图像映射到自然颜色系统、孟塞尔颜色 系统、PANTONE色卡颜色系统、TILO管理颜色系统，或者自定义的颜色系统的 至少一种，形成彩色化的生物组织弹性图像。
进一步地，步骤S110之后还包括步骤S120，对ROI区域中的对象区域进行 弹性分级，当生物组织弹性低于预设门限时显示相应等级的报警图标。
步骤S120中的弹性分级具体包括：
首先，提取出ROI区域中的对象区域；
然后，统计对象区域所有像素值分别位于各弹性值区间的个数；各弹性值 区间依弹性值由低到高分布；
再计算出每段弹性值区间的像素个数占该对象区域总像素的比例；
将占对象区域总像素的比例最高的弹性值区间的弹性级别判为该对象区域 生物组织的弹性级别。
本发明提出的一种实时超声弹性成像系统，包括：
换能器，用于对生物组织进行缓速挤压时对生物组织目标区域进行超声波 扫查并接收回波信号；
波束合成模块，用于对接收的回波电信号进行A/D转换、同相位叠加处理， 形成线数据；
B模式图像处理模块，用于对上述线数据进行处理，形成生物组织受外力缓 速挤压后产生形变的不同状态下的B模式图像；
第一CPU数据处理模块，用于对生物组织受外力缓速挤压后产生形变的不 同状态下的两帧B模式图像的同一位置区域即ROI区域使用光流法求取位移场； 并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数据缓存模 块；
第一GPU数据缓存模块，用于缓存接收到的含重复计算的子步骤的参数信 息并传输给各自的GPU数据处理模块进行数据处理；
一个或多个GPU数据处理模块，用于将上述各子步骤的运算分配给至少一 个GPU核工作单元，将至少一个GPU核工作单元组成GPU工作组，将每个GPU 工作组计算需要的参数信息映射到每个GPU工作组的显存中，等待所有GPU工 作组的参数信息映射完成，然后各GPU工作组进行光流法中含重复计算的各子 步骤的运算处理；
第二GPU数据缓存模块，用于缓存GPU数据处理模块的运算结果信息并传 输至第二CPU数据处理模块；
第二CPU数据处理模块，用于读出运算结果信息并带入光流法运算过程， 最终求出ROI区域的位移场即光流场；使用低通滤波器对轴向光流场求取应变， 得到图像ROI区域的轴向应变场；对经过滤波得到的图像ROI区域的轴向应变 场信息进行降噪声处理；
彩色定量显示模块，用于对降噪后的ROI区域轴向应变场信息进行可视化、 彩色处理得到生物组织彩色弹性图像；
图像辅助预诊断模块，用于对ROI区域中的对象区域采用统计的方法进行 弹性分级，当生物组织弹性低于预设门限时控制显示单元显示相应等级的报警 图标；
显示单元，用于显示彩色化的生物组织弹性图像，以及弹性等级对应的报 警图标。
本发明的优点在于：
1.采用本发明的光流法求取的位移场精度高，能够显著提高图像处理的精 度。
2.采用GPU加速的方法解决了光流法效率低下的问题，因此本发明既能求 得准确的生物组织弹性，又满足了实际使用中弹性计算的实时性要求。
3.使用对生物组织的弹性进行分级的方法，相对于直接使用阈值判定法更 科学合理。
4.采用GPU加速运算大量重复子步骤信息，能够提高系统的鲁棒性，减少 系统卡顿的现象，提高系统的流畅性；同时降低了CPU的工作负荷，有利于低 配置的超声诊断设备进行超声弹性成像处理或版本的升级。
附图说明
图1为本发明的结构组成示意图。
图2为求取位移场的光流法的原理示意图。
图3为本发明的滤波核示意图。
图4为本发明的流程图。
具体实施方式
下面结合具体附图和实施例对本发明作进一步说明。
本发明所提出的实时超声弹性成像方法，由CPU与多核GPU(Graphic ProcessingUnit)联合实现，该方法步骤如下：
S10.通过超声换能器对生物组织进行缓速挤压时对生物组织目标区域进行 超声波扫查并接收回波信号；
S20.通过波束合成模块对步骤S10中的接收的回波电信号进行A/D转换、 同相位叠加等处理，形成线数据；
S30.通过B模式图像处理模块对步骤S20得到的线数据进行降噪、滤波、 提高信噪比等处理，形成生物组织受外力缓速挤压后产生形变的不同状态下的B 模式图像；
S40.通过使用ROI框(感兴趣区域)对生物组织受外力缓速挤压后产生形 变的不同状态下的两帧B模式图像选定两个在不同帧图像的同一位置区域；
ROI(RegionOfInterest)，即感兴趣区域，在图像处理中是指从待处理 图像中提取出的要处理的区域。ROI的使用一方面可以突出框内区域的特征，另 一方面又可以提高图像处理的速度。本系统预设了一个ROI框，用户可以调整 ROI框的大小和位置以满足不同的要求。
S50.通过第一CPU数据处理模块对该两帧图像的ROI区域使用光流法求取 位移场；并在计算过程中将需要加速运算的含大量重复计算的子步骤的参数信 息传输至第一GPU数据缓存模块，进行缓存；其中需加速运算的含大量重复计 算的子步骤的参数信息是指：构建图像金字塔、中心梯度计算、计算目标图像 及其导数的形变、估计双变量、估计位移等计算的参数信息。
光流法求取位移场公式如下：
E ( U , V ) = ∫ ∫ ψ ( | I ( x + u , y + v , t + 1 ) - I ( x , y , t ) | 2 ) + λψ ( | ▿ u | 2 + | ▿ v | 2 ) dxdy ]]>
其中小写的u和v分别表示每个像素沿x和y方向的位移；大写的U和V分别表示由所有像素的u和v组成的位移场；t代表时间参数；λ代表一个常量，代表权重系数，λ≤6；10-6为此函数的一个优化系数。
S60.第一GPU数据缓存模块将接收到的需加速运算的子步骤的参数信息传 输给各自的GPU数据处理模块进行数据处理，S50中需加速运算的各子步骤的运 算分配给至少一个GPU核工作单元，由至少一个GPU核工作单元组成GPU工作 组。需加速运算的含大量重复计算的子步骤的参数信息是指：构建图像金字塔、 中心梯度计算、计算目标图像及其导数的形变、估计双变量、估计位移这些计 算的参数信息。例如：
构建图像金字塔的参数信息，包括：上下文、上一层图像的宽与高、下一 层图像的宽与高、上一层图像数据、下一层图像数据、图像通道数、图像位深、 图像步长等。
中心梯度计算的参数信息，包括：图像金字塔该层的图像数据、上下文、 图像宽与高、图像步长、x方向导数、y方向导数、x方向导数步长等。
计算目标图像及导数的形变的参数信息，包括：图像金字塔该层的图像数 据、上下文、图像宽与高、图像纹理、图像步长、x与y方向位移量等。
估计双变量的参数信息，包括：初步估计的位移数据、上下文、全局线程、 本地线程、tau值、x与y方向位移量等。
估计位移的参数信息，包括：上下文、初步估计的位移数据、全局线程、 本地线程、theta值、图像宽与高、图像步长、容错数据、x与y方向位移量等。
S70.将每个GPU工作组计算需要的参数信息映射到每个GPU工作组的显存 中，等待所有GPU工作组的参数信息映射完成，然后各GPU工作组进行步骤S50 中需加速运算的各子步骤的运算处理；GPU数据处理模块单元由至少一个的GPU 工作组组成，且GPU工作组包含至少一个GPU工作单元。
S80.将GPU工作组的运算结果信息传输至第二GPU数据缓存模块；
此步骤将需要GPU加速的各子步骤经过S60、S70运算的运算结果信息传输 至第二GPU数据缓存模块；
S90.第二GPU数据缓存模块将运算结果信息传输至第二CPU数据处理模块， 第二CPU数据处理模块读出运算结果信息并带入光流法运算过程，最终求出ROI 区域的位移场即光流场；(用光流法计算出的位移场也叫光流场)。
将第二GPU数据缓存模块中的信息由第二CPU数据处理模块读出，并将处 理结果带入光流法的运算过程。待所有需要GPU加速的处理子步骤均完成并将 各自的信息带入光流法，且完成计算后，图像ROI区域的光流场即求得。
在实际处理中，上述方法由第一CPU数据处理模块、第一GPU数据缓存模 块、GPU数据处理模块、第二GPU数据缓存模块、第二CPU数据处理模块共同进 行处理。这样，将光流法中构建图像金字塔、中心梯度计算、计算目标图像及 其导数的形变、估计双变量、估计位移等巨大繁重而又重复的工作，交由GPU 工作组的至少一个的GPU核工作单元来处理，当然一般而言GPU核工作单元的 数量大于2个，通过这样的方法大大提高了运行速度，对于提高算法的实时性、 扩展算法的应用领域具有非常重要和积极的作用，减少了系统卡顿现象，从而 提高了超声仪器的整体系统的稳定性，提高了使用人员操作感观的舒适度。
S100.使用低通滤波器对轴向光流场求取应变，得到图像ROI区域的轴向应 变场；此步骤在第二CPU数据处理模块中进行。
将光流法得到的位移场分解成X方向位移场和Y方向位移场，即横向位移 场与轴向位移场。由于轴向位移场对应着探头按压、组织受力的方向，因此本 发明更关注轴向位移场和轴向应变场。使用一个低通滤波器对轴向位移场进行 卷积，计算轴向位移场的偏导，得到图像ROI区域的轴向应变场，从而反映ROI 区域的组织弹性。
本发明采用的低通滤波器的核如下列公式所示：
y ( M + k ) = 25 * ( 3 M 4 + 6 M 3 - 3 M + 1 ) * k - 35 * ( 3 M 2 + 3 M - 1 ) * k 3 ( 2 M + 3 ) * ( 2 M + 1 ) * ( 2 M - 1 ) * ( M + 2 ) * ( M + 1 ) * M * ( M - 1 ) y ( M - k ) = - 25 * ( 3 M 4 + 6 M 3 - 3 M + 1 ) * k - 35 * ( 3 M 2 + 3 M - 1 ) * k 3 ( 2 M + 3 ) * ( 2 M + 1 ) * ( 2 M - 1 ) * ( M + 2 ) * ( M - 1 ) * M * ( M - 1 ) y ( M ) = 0 ]]>
如图3所示，该滤波核的大小为N行*1列(N为正奇数),其中，M＝(N-1)/2， k为大于0且不大于M的整数(1≤k≤M)，y代表滤波核不同位置的值。
由于本发明所采用的光流法的精度高，离群值较少，因此该滤波核中的N 可以取较小的值，即滤波核可以设计得较小，这样既不会降低应变计算的效果， 又能提高效率，还能保存一些图像细节。本发明采用的滤波核的N取值不大于 31。
S105.对经过滤波得到的图像ROI区域的轴向应变场信息进行取对数法降噪 声处理；此步骤在第二CPU数据处理模块中进行。
对经过滤波得到的图像ROI区域的轴向应变场信息进行降噪声处理。首先 对上个步骤得到的图像ROI区域轴向应变场取绝对值abs(dVR)；再对绝对值使 用取对数的方法log(1+abs(dVR))去除噪声；然后将其归一化到0～255区间的 整数，得到降噪后的图像ROI区域的轴向应变场。
S110.对降噪后的ROI区域轴向应变场信息进行可视化、彩色处理得到生物 组织彩色弹性图像。此彩色化定量显示的步骤在彩色定量显示模块中完成。
首先将由轴向应变场构成的灰度图像经波段合成转换成彩色图像，再将该 彩色图像映射到自然颜色系统、孟塞尔颜色系统、PANTONE色卡颜色系统、TILO 管理颜色系统，或者自定义的颜色系统中的至少一种，形成彩色化的生物组织 弹性图像。
S120.对ROI区域中的对象区域采用统计的方法进行弹性分级，当生物组织 弹性低于预设门限时显示相应等级的报警图标。弹性分级在图像辅助预诊断模 块中进行，显示单元可显示相应的报警图标。
生物组织弹性低于预设门限时，分为4个等级，即严重、中等、轻微、警 告，测量值达到相应等级即显示相应的图标。
首先，使用阈值分割方法或其他图像分割方法提取出ROI区域中组织弹性 较小的区域，该组织弹性较小的区域为疑似病变区域，也就是关注的对象区域。 该疑似病变区域通常是ROI区域的一部分子区域，也有可能充满整个ROI区域 (这时病变部位非常大)。然后，根据系统预设值Value_verylow、Value_low、 Value_middle、Value_high，统计该对象区域所有像素弹性值分别位于[0, Value_verylow)、[Value_verylow,Value_low)、[Value_low,Value_middle)、 [Value_middle,Value_high)四段区间的个数，再计算出每段区间的像素个数 占该组织弹性较小的区域总像素的比例。若[0,Value_verylow)区间的比例最 高，则判定该ROI区域中疑似病变区域的组织弹性非常低，系统上的图标显示 严重；若[Value_verylow,Value_low)区间的比例最高，则判定该ROI区域中 疑似病变区域的组织弹性很低，系统上的图标显示中等；若[Value_low, Value_middle)区间的比例最高，则判定该ROI区域中疑似病变区域的组织弹性 较低，系统上的图标显示轻微；若[Value_middle,Value_high)区间的比例最 高，则判定该ROI区域中疑似病变区域的组织弹性有些低，系统上的图标显示 警告。使用该统计的方法对生物组织的弹性进行分级，相对于直接使用阈值判 定法更科学合理。
本发明提出的实时超声弹性成像系统如图1所示，包括：
换能器，用于对生物组织进行缓速挤压时对生物组织目标区域进行超声波 扫查并接收回波信号；
波束合成模块，用于对接收的回波电信号进行A/D转换、同相位叠加处理， 形成线数据；
B模式图像处理模块，用于对上述线数据进行处理，形成生物组织受外力缓 速挤压后产生形变的不同状态下的B模式图像；
第一CPU数据处理模块，用于对生物组织受外力缓速挤压后产生形变的不 同状态下的两帧B模式图像的同一位置区域即ROI区域使用光流法求取位移场； 并在计算过程中将含重复计算的子步骤的参数信息传输至第一GPU数据缓存模 块；
第一GPU数据缓存模块，用于缓存接收到的含重复计算的子步骤的参数信 息并传输给各自的GPU数据处理模块进行数据处理；
一个或多个GPU数据处理模块，用于将上述各子步骤的运算分配给至少一 个GPU核工作单元，将至少一个GPU核工作单元组成GPU工作组，将每个GPU 工作组计算需要的参数信息映射到每个GPU工作组的显存中，等待所有GPU工 作组的参数信息映射完成，然后各GPU工作组进行光流法中含重复计算的各子 步骤的运算处理；
第二GPU数据缓存模块，用于缓存GPU数据处理模块的运算结果信息并传 输至第二CPU数据处理模块；
第二CPU数据处理模块，用于读出运算结果信息并带入光流法运算过程， 最终求出ROI区域的位移场即光流场；使用低通滤波器对轴向光流场求取应变， 得到图像ROI区域的轴向应变场；对经过滤波得到的图像ROI区域的轴向应变 场信息进行降噪声处理；
彩色定量显示模块，用于对降噪后的ROI区域轴向应变场信息进行可视化、 彩色处理得到生物组织彩色弹性图像；
图像辅助预诊断模块，用于对ROI区域中的对象区域采用统计的方法进行 弹性分级，当生物组织弹性低于预设门限时控制显示单元显示相应等级的报警 图标；
显示单元，用于显示彩色化的生物组织弹性图像，以及弹性等级对应的报 警图标。
显示单元可以是常规的台式机显示器，也可以是触摸屏显示器或手机等终 端接收显示单元。