发明名称 一种血管内超声弹性成像二维多级混合位移估计方法
摘要 本发明公开了一种血管内超声弹性成像二维混合位移估计方法,本发明提出的二维混合位移估计方法,结合了数据点的多级搜索和单级跟踪的方法,体现了合理平衡计算量和计算精度的优点。在算法初级阶段,利用二维多级搜索的方法来获取轴向和横向偏移量,而避免在高级阶段对横向方向进行搜索。在高级阶段,为最大程度减小计算量,则利用质量指导和零点相位差跟踪方法来估计轴向位移。这种算法可以很好地适用于血管内超声成像的实时成像要求,以及成像精度。
申请公布号 CN103735287A 申请公布日期 2014.04.23
申请号 CN201310645210.1 申请日期 2013.12.05
申请人 中国科学院苏州生物医学工程技术研究所 发明人 钱馨然;崔崤峣;顾天明;焦阳;吕铁军
分类号 A61B8/08(2006.01)I;A61B8/12(2006.01)I 主分类号 A61B8/08(2006.01)I
代理机构 南京经纬专利商标代理有限公司 32200 代理人 曹毅
主权项 1.一种血管内超声弹性成像二维混合位移估计方法,其特征在于,包括以下步骤:步骤1)利用商用超声换能器得到待测生物组织的压缩前的一帧二维射频信号;步骤2)手持所述超声换能器沿着探头纵向对所述生物组织施加一个微小的挤压,得到压缩后的一帧二维射频信号;步骤3)分别对从步骤1和步骤2中得到的组织压缩前、后的射频信号进行计算,得到基带解析数据,具体步骤为先对射频数据进行下采样减少数据点,为避免下采样时组织位移信息丢失,压缩前、后的信号采样频率应尽可能高,而下采样时,在不违背麦奎斯特采样定律的前提下,下采样频率也应该尽可能高,再进行希尔伯特变换,乘以exp<sup>(-jw</sup><sub>0</sub><sup>t)</sup>,其中w<sub>0</sub>为超声换能器的中心频率,则得到基带解析信号pre(x,y)和post(x,y),其中x表示沿探头轴向方向采样点,y表示沿探头横向方向采样点;步骤4)对步骤3中得到的压缩前、后基带解析信号进行信号处理,分别得到压缩前、后信号的幅值a<sub>1</sub>(x,y)和a<sub>2</sub>(x,y);步骤5)确定第一级窗口的数量及大小,可取9个,每行每列各3个;压缩前的目标区域第一级窗口大小为轴向或横向采样点的五分之一;(2)压缩后的搜索区域第一级窗口大小略大于轴向或横向采样点的五分之一,可取四分之一;步骤6)利用步骤5中确定的参数,对步骤4中得到的压缩前、后信号的幅值a<sub>1</sub>(x,y)和a<sub>2</sub>(x,y)进行第一级窗口搜索,搜索时对压缩前、后帧进行相关系数最大值匹配,相关系数公式如式(1)所示;(1)在最低分辨率下,压缩前和压缩后帧中的每一个窗口中,每四个采样点作为一个对象点,图2中的上一图中方形点为压缩前窗口的中心,对齐到压缩后的搜索区域进行相关系数计算,找到最佳匹配的矩阵形区域并做记录;(2)将分辨率提高一级,每两个采样点作为处理对象点,上一步中找到的外圆方形点(图2中的下一图)周围的x方向三个点作为这一步的匹配候选点,计算相关系数,找到最佳匹配点,并做记录;(3)分辨率最高的一级,不间隔点进行搜索,将上一步方形点周围的三个点作为候选点,且计算相关系数时不再间隔任何采样点,找到最佳匹配点,完成第一级搜索,只有在最低分辨率时,候选匹配位置尺寸与搜索区域的大小相关;在下一步分辨率提高的情况下,对相关系数的要求高于计算量,候选匹配区域定为三个,得到9个独立的轴向-横向位移估计值;<img file="2013106452101100001DEST_PATH_IMAGE001.GIF" wi="66" he="62" /><img file="309914DEST_PATH_IMAGE002.GIF" wi="371" he="62" /><img file="976519DEST_PATH_IMAGE001.GIF" wi="66" he="62" />(1)其中<img file="2013106452101100001DEST_PATH_IMAGE003.GIF" wi="50" he="42" /><img file="55333DEST_PATH_IMAGE004.GIF" wi="41" he="21" /><img file="740261DEST_PATH_IMAGE005.GIF" wi="41" he="21" />分别是横向和轴向位移,<img file="511908DEST_PATH_IMAGE006.GIF" wi="48" he="42" /><img file="2013106452101100001DEST_PATH_IMAGE007.GIF" wi="48" he="42" />分别是压缩前和压缩后的图像幅值,<img file="400230DEST_PATH_IMAGE006.GIF" wi="48" he="42" /><img file="282735DEST_PATH_IMAGE008.GIF" wi="48" he="42" />分别是<img file="822170DEST_PATH_IMAGE006.GIF" wi="48" he="42" /><img file="764718DEST_PATH_IMAGE009.GIF" wi="48" he="42" />的内插窗口平均值,T是窗口大小;步骤7)确定第二级的窗口数量及大小,其数量与算法参数无关,其中有九个为第一级9个窗口的子集,因此第一级搜索得到的轴向和横向位移估计就可以传递到第二级,第二级其余的窗口则均匀分布整个帧,第二级运动的方法和第一级相同,但间隔因子初始即取每两个采样点作为一个对象点,压缩前的目标区域窗口大小为轴向或横向采样点的十五分之一,而压缩后的搜索区域窗口大小则略大于目标区域,可取轴向或横向采样点的十分之一;步骤8)重复步骤6中的(1)、(2)、(3)步,但第二级的匹配过程由动态跟随来完成,首先第二级的9个窗口由第一级的输出确定一个初始点p1(x0,y0),其余的窗口初始值则由第一级的结果大致定向,定位于每一个窗口中心,将p1存储于点集S中,以及一个队列L集中,第二步,计算p1的四个邻点,并将其加入到S集中,<i>S</i>2=<i>{p1,p</i>2<i>,p</i>3<i>,p</i>4<i>,p</i>5<i>}</i>,<i>p</i>2=(<i>x</i>0+1<i>,y</i>0),<i>p</i>3=(<i>x</i>0<i></i>1<i>,y</i>0),<i>p</i>4=(<i>x</i>0<i>,y</i>0+1),<i>p</i>5=(<i>x</i>0<i>,y</i>0<i></i>1);此时S集中有四个点待处理,根据式(1),选出下一个处理点P2:<i>P</i>2=max(<i>C</i>(<i>p</i>1)<i>,C</i>(<i>p</i>2)<i>,C</i>(<i>p</i>3)<i>,C</i>(<i>p</i>4)<i>,C</i>(<i>p</i>5)),其中C即为互相关系数最大,此时用P2将队列L集更新,紧接着搜索当前点P2的邻点,如此经过连续不断的递进搜索,L集队列中永远是相关系数最大值,直到判断条件S集不再有更新,则该区域搜索完毕,得到更普遍分布的第二级轴向-横向位移估计值,具体流程可参考图3,步骤9)步骤8完成后,利用Savitzky-Golay滤波方法,对得到的横向估计进行平滑;步骤10)第三级搜索之前,对第二级得到的横向位移估计值利用双线性插值法,得到第三级的横向位移估计,因为横向位移对噪声敏感,窗口越小,这种影响越大,因此,第三级不再对横向位移估计做更进一步的精确;步骤11)确定第三级窗口数量及大小,第三级需要完成整幅帧的搜索,因此数量由窗口大小确定,窗口大小在轴向方向可取四十分之一,第二级获得的轴向位移估计作为第三级窗口的子集,进行结果传递;步骤12)第三级因为窗口较小,此时默认每个窗口的弹性模量为常数,位移可以近似看做深度的线性函数,窗口的位移可以假定参考窗口中心点的位移,但由于散斑和背向散射可能会引起小窗口中心位移估计局部产生巨大变差,因此利用对数压缩方法,对基带信号进行处理,以此来减小幅值的变化,如式2和式3所示;<img file="202653DEST_PATH_IMAGE010.GIF" wi="363" he="21" />(2)<img file="560953DEST_PATH_IMAGE011.GIF" wi="381" he="21" />(3)其中,c表示压缩因子,步骤13)因为第三级窗口的数量较多,不再利用互相关系数作为匹配准则,而是改用压缩前、后信号相关性最大时,其相位差趋近于零点的特点,进行牛顿迭代快速收敛的方法,如式4所示,其中arg(pre<sub>log</sub>(x,y))和arg(post<sub>log</sub>(x,y))为步骤12中得到的压缩前、后对数压缩后的相值,W(x,y)为采样点对位移估计的贡献值,即权函数,<img file="33523DEST_PATH_IMAGE012.GIF" wi="24" he="21" /><img file="146972DEST_PATH_IMAGE013.GIF" wi="24" he="21" /><img file="259153DEST_PATH_IMAGE014.GIF" wi="39" he="21" /><img file="2013106452101100001DEST_PATH_IMAGE015.GIF" wi="15" he="21" />表示前一步得到的位移值,<img file="483461DEST_PATH_IMAGE014.GIF" wi="39" he="21" /><img file="482641DEST_PATH_IMAGE016.GIF" wi="39" he="21" /><img file="2013106452101100001DEST_PATH_IMAGE017.GIF" wi="27" he="21" />为当前所求的位移值,w<sub>0</sub>是换能器的中心频率,而且因为信号的采样点为离散点,且经过下采样处理,因此实际中很难完全趋近于零点,可以通过限制迭代次数来达到精确位移估计的目的,同时不会进入死循环,迭代次数可以由具体的信号采样点数量及所要求的精度最终确定;<img file="953943DEST_PATH_IMAGE018.GIF" wi="153" he="62" /><img file="366470DEST_PATH_IMAGE019.GIF" wi="150" he="42" /><img file="332152DEST_PATH_IMAGE020.GIF" wi="66" he="42" /><img file="513734DEST_PATH_IMAGE020.GIF" wi="66" he="42" /><img file="703407DEST_PATH_IMAGE021.GIF" wi="357" he="42" />(4)步骤14)步骤13中的W(x,y),可取1,也可取步骤4得到的压缩前、后信号的幅值a<sub>1</sub>(x,y)和a<sub>2</sub>(x,y)对应窗口的和,如式5所示;<img file="790181DEST_PATH_IMAGE020.GIF" wi="66" he="42" /><img file="356291DEST_PATH_IMAGE022.GIF" wi="293" he="42" />(5)步骤15)利用步骤8中的跟踪方法结合步骤13中提到的相位差方法,将轴向位移估计传递到剩下的每一个第三级窗口,完成算法的最后一步搜索,即可得到整幅帧图的位移估计值。
地址 215000 江苏省苏州市高新区科灵路88号