CN106214182A - 基于lk光流法的hifu损伤剪切波弹性特性估计方法 - Google Patents

Abstract

本发明公开一种基于LK光流法的HIFU损伤剪切波弹性特性估计方法，在HIFU作用于组织仿体的过程中用高速摄影机获得损伤形成及振动的连续图像，运用图像金字塔的Lucas‑Kanade光流法获得损伤的振动位移曲线，再运用峰值时间法,即TTP(time‑to‑peak)法计算得到损伤上在HIFU作用下产生的声辐射力剪切波的传播速度，根据介质密度、介质中剪切波波速和介质剪切弹性模量之间的关系式即可对损伤的剪切弹性模量进行定量估计，进而实现对治疗过程中靶区组织的弹性特性进行实时监控并评价治疗效果。

Description

基于LK光流法的HIFU损伤剪切波弹性特性估计方法

技术领域

本发明属于超声成像技术领域，特别涉及一种基于LK金字塔光流法跟踪剪切波估计高强度聚焦超声(HIFU)热损伤弹性特性的系统及方法。

背景技术

高强度聚焦超声(HIFU)治疗是最具有潜在应用前景的肿瘤治疗新方法之一，因其非侵入性、针对深部靶组织的精确治疗等优点，已在各种肿瘤治疗方面获得了临床运用。HIFU治疗将超声换能器发射的高能超声束聚焦于生物体目标靶区内，在超声聚焦焦点处产生60℃以上的瞬态高温,通过热效应作用几秒钟后破坏靶区组织而不损伤周围组织来达到治疗的目的。HIFU对肿瘤组织的病理改变以不可逆的凝固性坏死为主，在治疗过程中，靶区组织的硬度会随着治疗进程发生变化，其变化范围为初始硬度的4～12倍。为了达到较好的治疗效果、确定合适的治疗剂量，临床上需要对靶区组织的弹性特性进行实时监控并评价治疗效果。

超声弹性成像通过获取相关组织的弹性信息进行成像，提供了一种独特的成像观察方法。弹性成像主要分为静态/准静态法和动态法(即声辐射力剪切波弹性成像)。和静态/准静态法采用外部机械施压使组织发生形变，再通过超声回波技术或核磁共振技术检测该形变来获取组织的弹性参数不同，声辐射力剪切波弹性成像通过超声波聚焦到组织内部来对组织施加压力，声辐射力可深入组织内部进行局部施压而不受组织边界的影响，可实时定量地对目标区域的弹性特性进行评估。

现有的声辐射力成像技术是通过发射聚焦超声脉冲在被测生物组织靶区内诱导出剪切波，再通过高帧率的超声平面波对剪切波进行跟踪来获得剪切波的相关参数,如：中国发明专利申请CN201310404937.0，专利名称“声辐射力脉冲成像估算方法和系统”，采用了超声成像监测剪切波的方法；中国发明专利申请CN201310558087.X，专利名称“弹性模量测量方法和系统”，通过发射超声射频波束跟踪剪切波；中国发明专利申请CN201410091601.8，专利名称“使用自适应时间阶段的超声ARFI位移成像”，使用超声扫描跟踪在焦点区周围的位置处的组织随时间的位移；中国发明专利申请CN201410514340.6，专利名称“医学超声成像中的剪切波检测”，利用超声来检测剪切波。然而运用超声跟踪声辐射力剪切波存在一些问题：1.超声激励诱发的剪切波振动位移在微米数量级，其回波易受噪声及生物体自身运动的干扰；2.现有的超声设备信噪比较低，所测得的振动位移结果精度较低；3.超声影像设备的图像分辨率较低，难以保证弹性图像的分辨率和对比度；4.利用超声对剪切波进行跟踪的实时性较差。

发明内容

本发明的目的在于提供一种基于LK光流法的HIFU损伤剪切波弹性特性估计系统及方法，以克服现有弹性成像方法中的不足；本发明方法具有良好的精确性和实时性，可定量监控治疗过程中损伤区域弹性特性的变化。

为实现上述目的，本发明采用如下的技术方案：

基于LK光流法的HIFU损伤剪切波弹性特性估计方法，包括以下步骤：

1)采集被测介质损伤振动的N张连续高速摄影图像；

2)采用LK法对步骤1)采集的N张连续高速摄影图像进行处理，得到通过两标记点标记的感兴趣区域的各时刻位移，进而绘制位移曲线；

3)对步骤2)处理后获得的标记点位移曲线，运用TTP算法求得损伤图像中两标记点之间的平均剪切波波速；

4)通过步骤3)获得的平均剪切波波速，获得被测介质两标记点之间的区域的剪切弹性模量。

进一步的，步骤1)具体包括以下步骤：任意波形发生器输出的信号经过射频功率放大器放大后激励HIFU换能器，对LED光源照射下的水箱中的被测介质施加作用；高速摄影设备同时被任意波形发生器触发，实现高帧率的图像采集。

进一步的，采用LK法对步骤1)采集的N张连续高速摄影图像进行处理前，先对图像进行降采样，形成图像金字塔。

进一步的，采用LK法对步骤1)采集的N张连续高速摄影图像进行处理具体包括以下步骤：

步骤S1、确认高速摄影所得连续高速摄影图像的总帧数N；

步骤S2、读入第i帧图像数据；

步骤S3、读入第i+1帧图像数据；

步骤S4、通过LK法计算出i和i+1两帧图像之间的位移场P0(i)；

步骤S5、将位移场P0(i)和P0(i-1)矢量相加，取各像素点的标量位移矩阵，得到叠加位移场P1(i)；

步骤S6、将当前图像帧数i与图像总帧数N进行比较，若i<N，进行步骤S7；若i＝N，则进行步骤S8，取出感兴趣区域的各时刻位移，绘制位移曲线。

进一步的，步骤3)中TTP算法具体为：在剪切波传播的横向路径上选取两个距离已知的标记点，作出这两个标记点的剪切振动位移曲线，在位移曲线上测量出这两个标记点在剪切振动过程中依次达到位移峰值的时间间隔，由两点的距离间隔除以时间间隔即可得到剪切波在这两点之间的平均传播速率。

进一步的，被测介质的杨氏弹性模量E和该被测介质中剪切波波速c_t的关系如下：

$c_t=\sqrt{\frac{\mathrm{\&mu;}}{\mathrm{\&rho;}}}---\left(1\right)$

$c_t=\sqrt{\frac{E}{3\mathrm{\&rho;}}}---\left(2\right)$

其中：

c_t——剪切波速/m·s^-1；

ρ——介质的密度/kg·m^-3；

μ——剪切弹性模量/kPa；

E——杨氏弹性模量/kPa。

相对于现有技术，本发明具有以下有益效果：

本发明为了克服超声跟踪剪切波的不足，提出了一种新的剪切波跟踪方法，即基于光流法的HIFU热损伤剪切波弹性特性估计方法。光流法可在估计区域得到稠密的位移场，具有灵敏度高、鲁棒性好、实时性好等优点，其检测精度达到微米级，达到了对剪切波进行跟踪的精度要求。本发明采用LK金字塔光流法来处理高速摄影采集的HIFU热损伤形成过程的图像序列，对HIFU作用下热损伤中的剪切波进行实时跟踪，进而定量估计目标区域的弹性特性。

附图说明

图1是本发明采用的计算介质中剪切波波速的TTP算法示意图。

图2是本发明中LK金字塔光流算法的图像降采样流程图。

图3是本发明对高速摄影设备所得的一系列高频图像帧进行位移估计的算法流程图。

图4是本发明中任意波形发生器作用于HIFU换能器和高速摄影设备的一组激励脉冲的构成图。

图5是本发明一种基于LK光流法的HIFU损伤剪切波弹性特性估计系统的框图。

图6是本发明中高速摄影设备拍摄所得的一帧BSA仿体中HIFU热损伤图像。

图7(a)是HIFU损伤形成初期在BSA仿体上选取的两个标记点，其坐标分别为(138，37)、(138，45)；图7(b)是HIFU损伤形成中期在BSA仿体上选取的两个标记点，其坐标分别(184，30)、(184，48)；图7(c)是HIFU损伤形成末期在BSA仿体上选取的两个标记点，其坐标分别为(205，68)、(205，89)。

图8(a)是HIFU损伤形成初期的标记点位移曲线图；图8(b)是HIFU损伤形成中期的标记点位移曲线图；图8(c)是HIFU损伤形成末期的标记点位移曲线图。

具体实施方式

已知含有牛血清蛋白(BSA)的聚丙烯酰胺凝胶仿体的剪切弹性模量μ和该仿体中剪切波波速c_t的关系、BSA仿体的杨氏弹性模量E和该仿体中剪切波波速c_t的关系如下：

$c_t=\sqrt{\frac{\mathrm{\&mu;}}{\mathrm{\&rho;}}}---\left(1\right)$

$c_t=\sqrt{\frac{E}{3\mathrm{\&rho;}}}---\left(2\right)$

其中：

c_t——剪切波速/m·s^-1；

ρ——介质的密度/kg·m^-3；

μ——剪切弹性模量/kPa；

E——杨氏弹性模量/kPa。

由上述两式可见，在BSA仿体密度已知的前提下，测得仿体中的剪切波波速即可对仿体的剪切弹性进行定量估计。

为获得介质中的剪切波波速，本发明采用了峰值时间算法，即TTP(time-to-peak)算法：在剪切波传播的横向路径上选取两个距离已知(通常为0.5mm-1mm)的标记点，作出这两个标记点的剪切振动位移曲线，在位移曲线上测量出这两个标记点在剪切振动过程中依次达到位移峰值的时间间隔，由两点的距离间隔除以时间间隔即可得到剪切波在这两点之间的平均传播速率。由TTP算法原理可见，该方法只关注两个标记点相继到达位移峰值的时刻，与峰值的具体数值无关，故用此法求剪切波波速误差较小。

为获得剪切波传播横向路径上所选两标记点的位移曲线，本发明采用如下方法：高速摄影设备对HIFU作用下BSA仿体中的损伤形成及振动过程进行视频采集，即任意波形发生器输出的信号经过射频功率放大器放大后激励HIFU换能器，对LED光源照射下的水箱中的BSA仿体施加作用，高速摄影设备和HIFU换能器同时被任意波形发生器触发，对BSA仿体振动过程进行背光投射的高帧率图像采集。得到一系列连续图像帧后，本发明用光流法对图像帧进行处理，得到仿体目标区域的位移曲线。

光流法具有高精度、高灵敏度等优点，可实时跟踪剪切波，故本发明采用迭代高斯金字塔的Lucas-Kanade(LK)光流法对高速摄影设备所得图像中目标区域的振动位移进行跟踪。

LK光流法是由Bruce D.Lucas和Takeo Kanade提出的一种基于图像序列的时空梯度，通过Newton-Raphson迭代直接寻找符合最佳匹配位置的图像匹配技术。相比于其他光流法，LK法的运算速度较快，且在图像中亮度梯度较小的位置仍能保持较高的准确性。

LK法基于3个假设：

1.亮度恒定，即物体上每个点的灰度是恒定的。这是基本光流法的假定，用于得到光流法基本方程；

2.小运动，即时间的变化不会引起位置的剧烈变化,保证了图像的灰度是可微的；

3.空间一致性，即一个场景上邻近的点投影到图像上也是邻近点，且邻近点速度一致。LK法这一独有的假定使其所得的估计结果是局部最优的。

将图像平面点(x,y)处、t时刻的灰度用E(x,y,t)表示，若当时间变化δt时，该点在x、y方向的位移分别为δx、δy，因其亮度恒定，故有：

E(x,y,t)＝E(x+δx,y+δy,t+δt) (3)

在LK法中，函数F(x)和G(x)是两张连续图像中的曲线，算法期望得到一个可以将F(x+h)和G(x)之间差异降到最小的视差向量h。先假定视差向量h的初始值为0，然后基于图像各点处的灰度梯度信息不断修正h的当前值以达到更好的匹配。在匹配过程中，误差的表达式为：

E≈∑_x[F(x+h)-G(x)]² (4)

在二维的线性近似情况下，为获得最佳匹配，即使误差关于h最小，有：

$h\&ap;\frac{{\mathrm{\&Sigma;}}_x{F}^{\&prime;}\left(x\right)\&lsqb;G\left(x\right)-F\left(x\right)\&rsqb;}{{\mathrm{\&Sigma;}}_x{F}^{\&prime;}{\left(x\right)}^2}---\left(5\right)$

此时视差向量h迭代形式为：

h₀＝0

$h_{k+1}=h_k+\frac{{\mathrm{\&Sigma;}}_{x}w\left(x\right)F^{\&prime;}(x+h)\&lsqb;G\left(x\right)-F(x+h_k)\&rsqb;}{{\mathrm{\&Sigma;}}_{x}w\left(x\right)F^{\&prime;}{(x+h_k)}^2}---\left(6\right)$

计算二维的h估计值，需要计算区域R内五个量的加权和：(G-F)F’(x)、(G-F)F’(y)、(F’(x))²、(F’(y))²及F’(x)*F’(y)。

为了提高算法对光照变化、图像区域移动的跟踪能力的稳定性及精确性，本发明选用了迭代的高斯图像金字塔来实现经典的LK光流算法，即通过对图像进行降采样来达到光流计算中“小位移”的要求。常规的图像金字塔在对图像进行降采样时，在图像长度方向与宽度方向上使用相同的降采样系数，该方法适用于长宽比接近于1的图像，而HIFU热损伤呈梭形，图像长度与宽度相差较大，为提高结果的准确性，本发明采用了非对称的图像金字塔来进行损伤图像的降采样。

在得到金字塔某层图像的位移信息后，先将此位移信息的尺寸扩展为下一层图像的尺寸，在扩展的同时，位移信息也会增大相应的倍数。用扩展后的位移信息对下一层图像进行插值运算，将得到的新图像与该层另一帧图像进行运算，即用低分辨率匹配约束在高分辨率情况下检测得到的可能匹配区域。本发明采用双立方插值对图像进行插值。

金字塔跟踪算法总体流程如下：

1.计算金字塔顶端图像I^Lm的光流；

2.由最上层光流的计算结果估计次上层图像L_m-1光流的初始值，再计算次上层图像L_m-1图像的光流精确值；

3.将次上层的光流结果作为下一层图像L_m-2光流的初始值，计算其精确值后再带入下一层图像L_m-3；

4.不断迭代直至原始图像。

采用图像金字塔的LK算法对运动图像进行处理，则在计算一个较大的视差向量h时，每一层的视差向量h^L始终保持为一个微小的值，提高了算法的鲁棒性和精确性，在对仿体HIFU热损伤的运动估计中能获得更加理想的局部跟踪定位精度。

以上图像金字塔的LK法得出了相邻两帧运动图像帧之间的位移场，对本发明中HIFU作用下BSA仿体中的损伤形成及振动过程所得的一系列图像帧，需要进行位移场的叠加来获取不同时刻的位移场：对于总帧数为N的一系列高速摄影图像，第一步读入第i帧和第i+1帧图像数据，按照上述LK法计算出两帧图像之间的位移场P0(i)；第二步将位移场P0(i)和P0(i-1)矢量相加，取各像素点的标量位移矩阵，得到叠加位移场P1(i)；第三步将i值与N值进行比较，若此时i<N，则将此刻的i赋值为i+1，重返第一步进行迭代；若i＝N，则取出感兴趣区域的各时刻位移，进而绘制出HIFU作用下BSA仿体中的损伤形成及振动过程的整体位移曲线。

下面结合附图对本发明做详细描述。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

请参阅图5所示，本发明一种基于LK光流法的HIFU损伤剪切波弹性特性估计系统，包括任意波形发生器1、功率放大器2、高强度聚焦换能器(HIFU)3、光源4、水箱5、BSA防体6、告诉摄影机7和PC机8。

任意波形发生器1(AWG420,Tektronix Inc.,US)输出的信号经过射频功率放大器2(AG1017,T&C Power Conversion Inc.,US)放大后激励HIFU换能器3，对LED光源4照射下的水箱中的BSA仿体6施加作用，HIFU换能器为中心频率为1.06MHz的单阵元换能器，换能器直径95mm，焦距100mm，焦区长度18mm(Chongqing Haifu Medical Technology Co.,Ltd,CN)；高速摄影设备7(MotionPro Y3,Integrated Design Tools Inc.,US)是基于MotionPro X高速数字摄影机，和HIFU换能器同时被任意波形发生器触发，实现高帧率的图像采集；采集得到的高速摄影图像由PC 8进行保存和处理。

BSA仿体6为质量分数为7％的牛血清蛋白(BSA)聚丙烯酰胺凝胶仿体，表1为BSA仿体的配方，其中TRIS用于调节仿体的PH至8，丙烯酰胺用于增加仿体硬度，TEMED用于催化仿体凝结。

表1 100mlBSA仿体溶液配方

本发明一种基于LK光流法的HIFU损伤剪切波弹性特性估计，包括以下步骤：任意波形发生器发出的驱动信号用于形成损伤和产生脉冲声辐射力，同时发射触发信号触发高速摄影进行图像采集。驱动信号如图4所示，由10组脉冲构成，每一组由100个脉冲长度为PD＝500μs、脉冲重复频率(PRF)为100Hz的脉冲串组成。为防止温度过高，每组脉冲串之间间隔30ms来保证散热。整个驱动信号的占空比(DC)为4.85％。由AG420的Marker输出一个脉冲上升作为触发信号来同步高速摄影。功率放大器的输出功率为50W。高速摄影的帧率为3000fps，对仿体进行背光投射的拍摄，得到的成熟的损伤图像如图6所示，图像中的黑色投影为损伤，未形成损伤的正常仿体为图像中的明亮背景，图像中一个像素点的尺寸大小为20um*20um。

得到损伤振动的一系列高速摄影图像后，用LK光流法对损伤图像上选取的标记点进行运动跟踪，在不同损伤形成时期所选择的标记点如图7(a)、图7(b)和图7(c)所示，记损伤中心点为“标记点1”，损伤下边的点为“标记点2”。

LK法对一系列高速摄影图像帧的处理流程如图3所示，步骤S1确认高速摄影所得系列图像的总帧数N；步骤S2读入第i帧图像数据；步骤S3读入第i+1帧图像数据；步骤S4通过LK法计算出i和i+1两帧图像之间的位移场P0(i)；步骤S5将位移场P0(i)和P0(i-1)矢量相加，取各像素点的标量位移矩阵，得到叠加位移场P1(i)；步骤S6将当前图像帧数i与图像总帧数N进行比较，若i<N，进行步骤S7，即将此刻的i赋值为i+1，重返步骤S2进行循环；若i＝N，则进行步骤S8，取出感兴趣区域的各时刻位移，进而绘制位移曲线。

由于此处针对的是较大的位移状况，为了保证光流法的估计精度，需要先对图像进行降采样，形成图像金字塔，过程如图2所示。首先基于LK光流法对顶层图像进行处理，得到其运动信息。然后用得到的运动信息对初始图像进行卷积，并将卷积后的图像尺寸扩展为下一层图像的尺寸，与另外一幅相同尺寸的图像再次进行运算，即以尺寸较小的图像对的计算结果为初始估计值，提高其底层图像的计算精度及计算速度。循环此过程，直至计算到原始尺寸的图像。

对图像进行处理后获得标记点的位移曲线，运用TTP算法即可求得损伤图像中两标记点之间的平均剪切波波速。TTP算法如图1所示，即剪切波传播的横向路径上选取两个距离已知的探测点，测量出这两点依次达到位移峰值的时间间隔，由两点的距离间隔除以时间间隔即可得到剪切波在这两点之间的平均传播速率。图8(a)为损伤形成初期的标记点位移曲线图，可见两个标记点到达峰值的时间间隔为10帧，因采用的高速摄像的帧频为3000帧，可知两点到达峰值的时间间隔为10*(1/3000)＝3.33ms，两点的距离为0.16mm，则由TTP算法可得剪切波波速为0.16mm/3.33ms＝0.048m/s，再由介质中剪切波波速和介质剪切弹性模量的关系式可知标记点之间的区域的剪切弹性模量为2.44Pa；图8(b)为损伤形成中期的标记点位移曲线图，可见两个标记点到达位移峰值的间隔为9帧，同理可得这一时期的剪切波速为0.12m/s，两个标记点之间的区域的剪切模量为15.26Pa。图8(c)为损伤形成末期的标记点位移曲线图，可见两个标记点到达位移峰值的间隔为3帧，同理可得这一时期的剪切波波速为0.42m/s，剪切模量为186.98Pa。此处所描述的具体实施例分析的时段约在20-30帧图像，但在调节高速摄影的帧频后，可以在更长或者更短的时间段内进行分析；在调节图像中选取的标记点的坐标后，可以分析图像中任何像素位置的位移。

表2所示为不同损伤时期对仿体目标区域进行处理后得到的剪切波波速及其对应的剪切模量。由表可见，随着损伤的不断成熟，剪切波波速逐渐增大，两个标记点之间区域的剪切模量也逐渐增大。在实际HIFU的作用下，仿体热损伤的硬度会随着治疗时间的加长而变大，使得剪切波的传导速度变快，这一实际和本发明实施过程所得的结果一致。

表2 不同损伤时期所得剪切波速及剪切模量

以上实例描述了本发明的具体实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，在不脱离本发明构思的前提下做出的若干变形和改进，都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (6)

1.基于LK光流法的HIFU损伤剪切波弹性特性估计方法，其特征在于，包括以下步骤：

1)采集被测介质损伤振动的N张连续高速摄影图像；

2)采用LK法对步骤1)采集的N张连续高速摄影图像进行处理，得到通过两标记点标记的感兴趣区域的各时刻位移，进而绘制位移曲线；

3)对步骤2)处理后获得的标记点位移曲线，运用TTP算法求得损伤图像中两标记点之间的平均剪切波波速；

4)通过步骤3)获得的平均剪切波波速，获得被测介质两标记点之间的区域的剪切弹性模量。

2.根据权利要求1所述的基于LK光流法的HIFU损伤剪切波弹性特性估计方法，其特征在于，步骤1)具体包括以下步骤：任意波形发生器输出的信号经过射频功率放大器放大后激励HIFU换能器，对LED光源照射下的水箱中的被测介质施加作用；高速摄影设备同时被任意波形发生器触发，实现高帧率的图像采集。

3.根据权利要求1所述的基于LK光流法的HIFU损伤剪切波弹性特性估计方法，其特征在于，采用LK法对步骤1)采集的N张连续高速摄影图像进行处理前，先对图像进行降采样，形成图像金字塔。

4.根据权利要求1所述的基于LK光流法的HIFU损伤剪切波弹性特性估计方法，其特征在于，采用LK法对步骤1)采集的N张连续高速摄影图像进行处理具体包括以下步骤：

步骤S1、确认高速摄影所得连续高速摄影图像的总帧数N；

步骤S2、读入第i帧图像数据；

步骤S3、读入第i+1帧图像数据；

步骤S4、通过LK法计算出i和i+1两帧图像之间的位移场P0(i)；

步骤S5、将位移场P0(i)和P0(i-1)矢量相加，取各像素点的标量位移矩阵，得到叠加位移场P1(i)；

步骤S6、将当前图像帧数i与图像总帧数N进行比较，若i<N，进行步骤S7；若i＝N，则进行步骤S8，取出感兴趣区域的各时刻位移，绘制位移曲线。

5.根据权利要求1所述的基于LK光流法的HIFU损伤剪切波弹性特性估计方法，其特征在于，步骤3)中TTP算法具体为：在剪切波传播的横向路径上选取两个距离已知的标记点，作出这两个标记点的剪切振动位移曲线，在位移曲线上测量出这两个标记点在剪切振动过程中依次达到位移峰值的时间间隔，由两点的距离间隔除以时间间隔即可得到剪切波在这两点之间的平均传播速率。

6.根据权利要求1所述的基于LK光流法的HIFU损伤剪切波弹性特性估计方法，其特征在于，被测介质的杨氏弹性模量E和该被测介质中剪切波波速c_t的关系如下：

$c_t=\sqrt{\frac{\mathrm{\&mu;}}{\mathrm{\&rho;}}}---\left(1\right)$

$c_t=\sqrt{\frac{E}{3\mathrm{\&rho;}}}---\left(2\right)$

其中：

c_t——剪切波速/m·s^-1；

ρ——介质的密度/kg·m^-3；

μ——剪切弹性模量/kPa；

E——杨氏弹性模量/kPa。