一种声辐射力双向剪切波复合成像方法及装置与流程

本发明涉及超声波成像领域，特别是一种声辐射力双向剪切波复合成像方法及装置。

背景技术：

剪切波弹性成像技术可实现生物组织实时硬度定量检测，为临床判断组织病变情况提供依据。它的基本原理如下：声辐射力聚焦冲击能在组织内部产生剪切波，由于剪切波在不同硬度的组织传播速度有差异，通过检测剪切波在不同位置的传播速度可以间接反映该位置的软硬情况。根据剪切波波速大小进行伪彩色映射即为剪切波弹性成像。因此，剪切波波速估计的准确性对组织弹性情况的判断至关重要。

剪切波在组织传播过程中会引起组织的位移形变，通过位移形变重建出观察位置的剪切波震动曲线。最后利用剪切波波速估计方法，对观察位置不同位置所重建的剪切波震动曲线进行波峰匹配，找出波峰相差时间间隔，此间隔可以理解为剪切波从该两个检测位置之间传播所用时间。而两检测位置之间的距离为已知，通过距离除以时间得出剪切波经过两个检测位置之间的平均速度。根据检测区域剪切波速度值进行伪彩色映射叠加在超声灰阶图像上，形成剪切波弹性图。

通过声辐射力聚焦冲击能够在人体组织内部产生剪切波，常见的声辐射力冲击方法有单焦点冲击或多焦点冲击(马赫锥效应)等方法。单焦点声辐射力冲击产生以焦点为源，近似球面状传播的剪切波，此方法产生的剪切波较弱，传播范围较小；多焦点声辐射力冲击利用马赫锥效应产生接近平面状形态的剪切波，使得剪切波效果增强，传播区域更大。但剪切波在人体组织中传播时，衰减非常快。因此靠近声辐射力冲击源位置的剪切波信噪比明显要比远离位置的剪切波信噪比高。导致一个弹性测量区域中，远离剪切波源的位置的剪切波波速估计准确性更低。有方法提出通过对同一观察区域，在不同位置进行声辐射力冲击，最终形成多幅剪切波弹性图，最后在复合叠加形成该观察区域的最终剪切波弹性图。该方法提高了最终剪切波弹性图的质量，但形成多幅剪切波弹性图耗时过长，而且需要多次对组织不同位置进行声辐射力冲击。实际人体组织状态是运动变化的，多幅剪切波弹性图形成过程中背景其实已经发生变化，此时最终复合的剪切波弹性图并没有准确反映组织的真实弹性情况。

技术实现要素：

针对上述问题，本发明主要提出一种声辐射力双向剪切波复合成像方法及装置，该方法先通过对人体组织选定位置进行一组声辐射力冲击，产生两种不同方向传播的剪切波。然后通过超声采集观察区域指定位置的回波信号。对各个检测位置的回波信号进行形变估计运算，并利用方向滤波技术选出不同传播方向的剪切波形变估计结果。然后得出以不同方向传播的剪切波波速估计结果，将两个结果加权复合得到最终的剪切波弹性结果。

本发明为实现其技术目的所采用的技术方案是：一种声辐射力双向剪切波复合成像方法，包括，该方法中，在选择的弹性观察区域两侧，按照预先设定的声辐射力冲击焦点位置和顺序进行声辐射力聚焦冲击；采集各横向检测位置线上各检测位置的超声回波数据；对采集的超声回波数据进行处理包括以下步骤：

步骤1、对超声回波数据进行剪切波形变估计运算的步骤；该步骤中，对各横向检测位置线上采集的超声回波数据形变估计运算，形成“形变-时间”结果矩阵；

步骤2、对形变估计结果进行方向滤波的步骤；该步骤中对各检测位置线的“形变-时间”结果矩阵进行方向滤波；

步骤3、画出经过方向滤波后的各检测位置线的“形变-时间”曲线；

步骤4、拟合“时间-距离”直线的步骤；该步骤中，根据各“形变-时间”曲线峰值对应时间值和检测位置物理距离值，形成拟合“时间-距离”直线；

步骤5、计算最终剪切波速度值的步骤；该步骤中：求取“时间-距离”直线的斜率的倒数，即可得到剪切波速度值；

步骤6、剪切波速度复合的步骤；该步骤中，将两个方向的剪切波速度值进行复合,计算如下：

剪切波速度复合值＝a*v1+b*v2

式中：a，b为正加权系数，且0＜a+b≤1；v1，v2分别为两个方向的剪切波速度值；

步骤7、生成弹性图像的步骤；重复上述步骤1至步骤6得到整个弹性成像区域的剪切波复合速度值，按照下式进行灰度映射运算：

式中：g表示结果灰度值，v表示剪切波复合速度，vmax、vmin分别表示剪切波复合速度最大值和最小值。

本发明方法通过对选定位置进行声辐射力冲击，产生两个不同源，沿不同方向传播的剪切波。合理覆盖了需要观察的组织区域。然后通过方向滤波技术选取不同方向的剪切波形变估计结果进行剪切波波速估计，得到两组不同方向，经过同一观察区域的剪切波速度结果。此方法弥补了剪切波传播时衰减导致的不同位置的剪切波波速估计信噪比差异，提升了观察区域整体的剪切波波速估计信噪比，使观察区域的弹性结果更加准确。而且本发明只需要一组冲击，即可产生两帧弹性图像进行复合优化。在保持弹性图像帧率不变的情况下提升了弹性图像的质量。

进一步的，上述的声辐射力双向剪切波复合成像方法中：冲击焦点的位置在选择观察区域左右两边界位置或者上下边界位置；冲击焦点的间隔3-10mm；声辐射力冲击脉冲长度为30-100us；声辐射力冲击脉冲重复频率为0.5-5khz；声辐射力焦点冲击顺序可以是从上到下或者从左到右。

进一步的，上述的声辐射力双向剪切波复合成像方法中：所述的步骤1中，形变估计运算对同一检测位置线，相邻次采集的数据进行下式归一化互相关运算：

ut＝cτ/2

式中：sr为当前超声射频回波信号，sd下一次超声射频回波信号，r为互相关系数，τ为对应位移值的超声回波时间延迟，c为超声传播速度，t为信号段周期长度。

进一步的，上述的声辐射力双向剪切波复合成像方法中：所述的步骤2中，方向滤波的过程如下：

步骤2.1、对各检测位置线的“形变-时间”结果矩阵进行二维傅里叶变换。

步骤2.2、与两组不同方向的方向滤波算子矩阵相乘。

步骤2.3、进行二维傅里叶反变换；选取结果实数部分作为最终滤波的结果。

进一步的，上述的声辐射力双向剪切波复合成像方法中：

所述的步骤2.2中两组不同方向的方向滤波算子矩阵为：

本发明还提供一种实现上述方法的声辐射力双向剪切波复合成像装置，包括从选择的弹性观察区域两侧，按照预先设定的声辐射力冲击焦点位置和顺序进行声辐射力聚焦冲击的声辐射力聚焦冲击的装置，采集各横向检测位置线上各检测位置的超声回波数据的超声回波数据采集装置，对超声回波数据采集装置采集的数据进行处理的数据处理装置；所述的数据处理装置包括：对超声回波数据进行形变估计运算、形成“形变-时间”结果矩阵的模块；对各检测位置线的“形变-时间”结果矩阵进行方向滤波的模块；画出经过方向滤波后的各检测位置线的“形变-时间”曲线的模块；拟合“时间-距离”直线的模块；计算最终剪切波速度值的模块；剪切波速度复合的模块；生成弹性图像的模块。

进一步的，上述的声辐射力双向剪切波复合成像装置中：所述的声辐射力聚焦冲击的装置的声辐射力冲击脉冲长度为30-100us；声辐射力冲击脉冲重复频率为0.5-5khz。

进一步的，上述的声辐射力双向剪切波复合成像装置中：所述的对各检测位置线的“形变-时间”结果矩阵进行方向滤波的模块包括：对各检测位置线的“形变-时间”结果矩阵进行二维傅里叶变换的模块；与两组不同方向的方向滤波算子矩阵相乘模块；进行二维傅里叶反变换模块。

下面结合附图和具体实施方式对本发明进行进一步的说明。

附图说明

附图1为本发明实施例1流程图。

附图2为本发明实施例1声辐射力冲击焦点轨迹示意图。

附图3为本发明实施例1剪切波传播示意图。

附图4为本发明实施例1沿检测位置线采集超声回波信号。

附图5为本发明实施例1方向滤波后的各个检测位置形变估计结果。

附图6为本发明实施例1拟合“时间-距离”直线。

具体实施方式

实施例1，本实施例是一种声辐射力双向剪切波复合成像装置，包括从选择的弹性观察区域两侧，按照预先设定的声辐射力冲击焦点位置和顺序进行声辐射力聚焦冲击的声辐射力聚焦冲击的装置，采集各横向检测位置线上各检测位置的超声回波数据的超声回波数据采集装置，对超声回波数据采集装置采集的数据进行处理的数据处理装置；所述的数据处理装置包括：对超声回波数据进行形变估计运算、形成“形变-时间”结果矩阵的模块；对各检测位置线的“形变-时间”结果矩阵进行方向滤波的模块；画出经过方向滤波后的各检测位置的“形变-时间”曲线的模块；拟合“时间-距离”直线的模块；计算最终剪切波速度值的模块；获得两个方向的剪切波速度复合值的模块；灰度映射的模块。

本实施例的声辐射力双向剪切波复合成像装置使用在生物组织弹性测量过程中，选取弹性观察区域，然后进行如下步骤：如图1所示。

步骤1：根据选择的弹性观察区域，按照预先设定的声辐射力冲击焦点位置和顺序采用声辐射力聚焦冲击的装置进行声辐射力聚焦冲击，如附图2所示。有六个声辐射力冲击焦点，分别在弹性观察区域a的两侧，从上到下依次为f1、f2、f3、f4、f5和f6，其中冲击焦点横向x向位置可选例如：观察区域a左右两边界位置；冲击焦点纵向y向位置可选例如观察区域上下边界位置，焦点纵向间隔大小，也就是f1、f3、f5或者f2、f4、f6之间的间隔可选例如：5mm，一般在3-10mm之间可行；声辐射力冲击脉冲长度可选例如：50us，一般在30-100us范围内选择。声辐射力冲击脉冲重复频率可选例如：1khz，一般在500hz-10khz之间选择；声辐射力焦点冲击顺序可选如附图2中焦点序号顺序，从上到下依次冲击。声辐射力冲击焦点数量可选例如：6个焦点；通过此步骤所产生的剪切波传播路径大致如附图3中s1和s2所示。

步骤2：完成声辐射力聚焦冲击步骤后，选定各横向检测位置线位置，如图4所示，沿检测位置线有例如：8条，实践中可选至少2条检测位置线。采集超声回波信号，并采集检测线上超声回波数据，各检测位置采集次数可选例如：50次，实践中可在30-100次中选择。采集重复频率可选例如：1khz；实践中可以在500hz到5khz之间选择。这里的检测位置是分布在检测线上的，一条线上采集数据点可以为500-2048个点,实践中可选1024个点。

步骤3：根据各检测位置回波数据，进行形变估计运算。形变估计算法可选择例如：归一化互相关运算，如下(1)式所示：对同一检测位置线，相邻次采集的数据进行归一化互相关运算，其中sr为当前超声射频回波信号，sd下一次超声射频回波信号，r为互相关系数，τ为对应位移值的超声回波时间延迟，c为超声传播速度，t为信号段周期长度；求取r的最大值所对应结果τ即为相邻回波时延，根据公式2计算得到形变值uτ；

uτ＝cτ/2(2)

步骤4：对步骤3所得的各检测位置线的“形变-时间”结果矩阵进行二维傅里叶变换，然后与两组不同方向的方向滤波算子矩阵相乘，最后进行二维傅里叶反变换，选取结果实数部分作为最终结果，得出两组不同方向的滤波结果。方向滤波算子可选如下：

上述方向选通算子矩阵中，dir_mask_l表示选取从左到右方向，dir_mask_r表示从右到左方向。其中dir_mask_l算子内部结构为:将算子矩阵四等分，第1和3象限范围部分的值为1，其余部分为0。dir_mask_l为一中心对称算子。dir_mask_r算子构造方法一样，数据值相反，这里不再重复。

步骤5：作出经过方向滤波后的各检测位置的“形变-时间”曲线，如附图5所示，作出了其中一个深度的左方向方向滤波后各检测位置的“形变-时间”曲线。搜索图中每条“形变-时间”曲线的峰值，并记录峰值所对应的横向坐标值,即时间值。

步骤6：根据步骤5中峰值对应时间值和检测位置物理距离值，作出拟合“时间-距离”直线；如图6所示。求取“时间-距离”直线的斜率的倒数，即可得到剪切波速度值。同理可得到另一个方向的方向滤波结果后最终的剪切波速度值。

步骤7：将两个方向的剪切波速度值进行复合，计算如式5所示：其中a，b为正加权系数，且0＜a+b≤1；v1，v2分别为不同方向的剪切波速度值。a，b例如可选：0.5。

剪切波速度复合值＝a*v1+b*v2(5)

步骤8：重复步骤3-步骤7即可得到整个弹性成像区域的剪切波复合速度值，根据式6的线性映射到灰度0-255之间：其中g表示结果灰度值，v表示剪切波复合速度，vmax、vmin分别表示剪切波复合速度最大值和最小值。