一种基于单空化源分离及定位的超分辨超声被动空化成像方法及系统

技术领域

本发明属于超声检测与超声成像技术领域，涉及一种通过对单空化源进行分离和定位以获得小区域内空化源空间分布的精细解析的超分辨超声被动空化成像方法及系统。

背景技术

声空化是指在空化激励超声换能器发射的超声波的激励下介质中空化泡形成、膨胀收缩和溃灭的物理现象，其增强了超声波与生物体之间的相互作用，是肿瘤消融、血栓溶解、药物释放等治疗超声应用的核心物理基础。空化的检测有助于监控超声治疗过程并进一步澄清其背后的物理机制，最常用的手段有光学检测和声学检测方法。高速摄影、声致发光等光学检测方法可以直观地观察空化变化过程且具有高的空间分辨率，但其受到介质透明度的限制，不能用于检测生物组织中的空化。声学检测方法则可以对生物组织中空化的声学特性进行定量分析，根据检测换能器的工作模式分为主动检测和被动检测，其中后者具有在超声波激励过程中同步地对空化进行检测的优势。被动声学检测方法根据检测换能器的种类进一步分为单阵元检测和阵列检测(即超声被动空化成像)，其中后者弥补了前者无法提供空化空间位置和分布信息的缺陷，近年来得到广泛关注。

超声被动空化成像一般是对阵列超声换能器被动接收到的空化信号进行基于相对渡越时间信息的波束合成处理并根据所得信号计算声能量，其空间分辨性能取决于阵列超声换能器的衍射模式(包括孔径和接收带宽等因素)。目前实施超声被动空化成像的阵列超声换能器主要有常用于临床诊断的换能器(如线阵、相控阵等)和专用于经颅成像的半球阵换能器，其中前者适用于人体不同部位，具有更为广泛的应用前景。但诊断用换能器有限的孔径会限制空化图像的空间分辨率(尤其是轴向分辨率)，同时组织异质性、源相互作用等因素的影响使得空化图像中往往伴随着高水平的干扰伪迹。目前研究人员主要通过改进和优化成像方法来提高图像质量。一类典型的成像方法采用了自适应波束合成技术，包括基于信号二阶或高阶统计特性的波束合成和基于协方差矩阵特征分解的波束合成等，其对施加到各阵元的权重系数进行优化并据此进行波束合成，从而使得干扰伪迹得到明显消除，增强了成像分辨率。除此之外，研究人员也通过对波束合成信号施加幅度/相位相干权因子，对信号中心频率正负两侧增量相等的两个频率处的幅度进行耦合相乘，对延时信号进行符号开方、耦合相乘及叠加，以及对阵列信号进行互补变迹处理并以互相关因子对波束合成信号加权等方式对成像方法进行改进，研究表明这些方法均有助于成像分辨率的提高。

尽管目前发展的多种成像方法表现出良好的性能，但其并不能从根本上克服成像换能器衍射限制的问题，以至于成像分辨率难以得到突破性的提高。而在声空化及相关的治疗超声应用中，局部小区域空间范围内空化的检测对于研究声空化物理现象本身以及解释和澄清相关的生物物理机制至关重要。鉴于此，亟待提出一种能够对小区域内空化源空间分布进行精细解析的超分辨超声被动空化成像方法及系统，而这也是超声被动空化成像领域一直以来亟待攻克的技术难题。

发明内容

本发明的目的在于提供一种基于单空化源分离及定位的超分辨超声被动空化成像方法及系统。

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

一种基于单空化源分离及定位的超分辨超声被动空化成像方法，包括以下步骤：

1)建立成像坐标系并规划时间门控能量图的像素网格，根据该像素网格中每个像素坐标到空化激励超声换能器表面中心坐标处切线或切面的距离和每个像素坐标与阵列超声换能器各阵元坐标之间的距离计算双程延时，以双程延时为时间门控区间的初始时刻并设置时间门控区间的时间长度，得到各阵元的时间门控区间，根据各阵元的时间门控区间从阵列超声换能器被动接收所得的每帧原始射频信号中提取各阵元的时间门控信号，对各阵元的时间门控信号沿阵元方向叠加，得到时间门控合成信号，将时间门控合成信号的平方沿采样点方向叠加，得到每帧原始射频信号对应的时间门控能量图；

2)对每帧时间门控能量图进行阈值化处理后根据所有非零值像素建立全局像素集合，根据设置的闭邻域半径和闭邻域像素数目阈值，通过重复初始化像素簇及扩展像素簇过程，对全局像素集合中的像素进行聚簇，得到每帧时间门控能量图的全局像素簇；

3)针对每帧时间门控能量图的每个全局像素簇划定局部聚簇窗，利用局部聚簇窗对时间门控能量图进行框选，得到局部能量图，对局部能量图进行阈值化处理后根据所有非零值像素建立局部像素集合，根据设置的闭邻域半径和闭邻域像素数目阈值，通过重复初始化像素簇及扩展像素簇过程，对局部像素集合中的像素进行聚簇，得到局部像素簇，从局部像素簇中寻找局部像素主簇并将局部能量图中不包含在局部像素主簇中的所有像素的像素值置为零，得到单空化源能量图，从而实现单空化源的分离；

4)提取所有单空化源能量图的像素峰值，根据设置的像素峰值阈值对其中每个单空化源能量图的像素峰值进行判断，由像素峰值小于像素峰值阈值的单空化源能量图建立单空化源能量图集合；

5)在成像坐标系中建立高斯分布函数，对单空化源能量图集合中的每个单空化源能量图进行高斯拟合以定位单个空化源，得到高斯分布函数中参数的估计值，其中高斯分布峰值位置的估计值为空化源坐标的估计值，将参数的估计值代入到高斯分布函数中计算高斯拟合结果，计算给定置信水平下的空化源坐标估计不确定度，根据单空化源能量图及单空化源能量图的高斯拟合结果计算拟合优度；

6)根据所有空化源的空化源坐标估计不确定度和拟合优度建立空化源坐标估计不确定度集合和拟合优度集合，对各集合中元素分别进行升序排列并计算升序排列后的各集合中元素的四分位数和四分位距，根据四分位数和四分位距设置空化源坐标估计不确定度的阈值和拟合优度的阈值，利用设置的空化源坐标估计不确定度的阈值和拟合优度的阈值检测出正常定位的空化源；

7)规划超分辨成像的像素网格并针对每个正常定位的空化源建立定位分布函数，计算该像素网格中每个像素坐标的定位分布函数值，得到每个正常定位的空化源的定位分布图，对所有正常定位的空化源的定位分布图进行叠加，然后进行归一化及对数化处理，得到超分辨成像结果。

优选的，所述步骤1中，双程延时的计算公式表示为：

其中，

当空化激励超声换能器的中轴线在阵列超声换能器的成像平面内时，

当空化激励超声换能器的中轴线与阵列超声换能器的成像平面相交且与阵列超声换能器的阵元方向垂直时，

其中，α为空化激励超声换能器中轴线与阵列超声换能器中轴线的夹角，|·|表示取绝对值。

优选的，所述步骤1中，时间门控区间的时间长度为空化激励超声换能器发射脉冲的脉冲长度。

优选的，所述步骤1中，各阵元的时间门控信号表示为：

其中，i＝1,2,...,N

优选的，所述步骤2或步骤3中，闭邻域半径以像素数目为单位且设置为

优选的，所述步骤2或步骤3中，初始化像素簇的条件为：从全局像素集合或局部像素集合中任选一未被判断的像素，若所选像素的闭邻域中像素的数目小于闭邻域像素数目阈值，则重新选择，直至所选像素的闭邻域中像素的数目大于等于闭邻域像素数目阈值，此时所选像素为核心像素，并初始化一个像素簇。

优选的，所述步骤2或步骤3中，扩展像素簇具体包括以下步骤：

2.1)将从全局像素集合或局部像素集合中选择并确定为核心像素的像素聚入到所述像素簇中；

2.2)将所述核心像素的去心闭邻域中的所有像素归入到一待判断像素集合中；

2.3)在待判断像素集合中任选一未被判断的像素，若所选像素是核心像素，则将该像素去心闭邻域中未被判断的像素归入到所述待判断像素集合中，并转至步骤2.4；若所选像素不是核心像素，则直接转至步骤2.4；

2.4)若步骤2.3所选像素未被聚簇，则将该像素聚入到所述像素簇中；

2.5)重复步骤2.3和2.4，直至遍历完所述待判断像素集合中的所有像素，从而生成一个全局像素簇或局部像素簇。

优选的，所述步骤3中，局部聚簇窗的中心为全局像素簇的质心对应的像素，局部聚簇窗的尺寸根据点扩散函数确定。

优选的，所述步骤3中的阈值化处理所用的阈值系数小于步骤2中的阈值化处理所用的阈值系数。

优选的，所述步骤3中，从局部像素簇中寻找局部像素主簇具体包括以下步骤：设置标签像素为全局像素簇的质心对应的像素，则包含标签像素的某一局部像素簇为局部像素主簇。

优选的，所述步骤4中，像素峰值阈值为所有单空化源能量图的像素峰值的最大值的0.4～0.6倍。

优选的，所述步骤5中，高斯分布函数允许角度偏转，其表达式为：

其中，A为高斯分布峰值，μ

优选的，所述步骤5中，高斯拟合需要对高斯分布函数中的参数进行初始化：将高斯分布峰值初始化为单空化源能量图中的最大像素值，将高斯分布峰值位置初始化为单空化源能量图中最大像素值所在坐标，根据点扩散函数初始化高斯分布标准差。

优选的，所述步骤5中，高斯拟合采用基于Levenberg-Marquardt算法的非线性最小二乘法。

优选的，所述步骤6中，空化源坐标估计不确定度和拟合优度的阈值分别为：

其中，

优选的，所述步骤6中，检测出正常定位的空化源具体包括以下步骤：若x轴方向和z轴方向的空化源坐标估计不确定度δ

优选的，所述步骤7中，超分辨成像的像素网格的像素数目为时间门控能量图的像素网格的像素数目的λ

优选的，所述步骤7中，定位分布函数为允许角度偏转的高斯分布函数，该高斯分布函数的峰值位置、标准差、峰值以及偏转角分别为对单空化源能量图进行高斯拟合所得的x轴方向和z轴方向的空化源坐标的估计值、x轴方向和z轴方向的空化源坐标估计不确定度、高斯分布峰值的估计值以及高斯分布偏转角的估计值。

一种基于单空化源分离及定位的超分辨超声被动空化成像系统，该系统包括时间门控能量图计算模块、单空化源分离模块、单空化源能量图像素峰值判断模块、高斯拟合空化源定位模块、空化源定位检测模块和空化源定位分布图计算及处理模块；

所述时间门控能量图计算模块用于执行上述步骤1，主要是用于建立成像坐标系并规划时间门控能量图的像素网格、根据该像素网格中每个像素坐标到空化激励超声换能器表面中心坐标处切线或切面的距离和每个像素坐标与阵列超声换能器各阵元坐标之间的距离计算双程延时、以双程延时为时间门控区间的初始时刻并设置时间门控区间的时间长度、根据各阵元的时间门控区间从阵列超声换能器被动接收所得的每帧原始射频信号中提取各阵元的时间门控信号、对各阵元的时间门控信号沿阵元方向叠加，以及将叠加所得时间门控合成信号的平方沿采样点方向叠加，从而得到每帧原始射频信号对应的时间门控能量图；

所述单空化源分离模块用于执行上述步骤2和步骤3，具体包含像素全局聚簇子模块和像素局部聚簇及寻找局部像素主簇子模块；

所述像素全局聚簇子模块用于执行上述步骤2，主要是用于对时间门控能量图计算模块得到的每帧时间门控能量图进行阈值化处理后根据所有非零值像素建立全局像素集合，以及根据设置的闭邻域半径和闭邻域像素数目阈值，通过重复初始化像素簇及扩展像素簇过程，对全局像素集合中的像素进行聚簇，从而得到每帧时间门控能量图的全局像素簇；

所述像素局部聚簇及寻找局部像素主簇子模块用于执行上述步骤3，主要是用于针对像素全局聚簇子模块得到的每帧时间门控能量图的每个全局像素簇划定局部聚簇窗、利用局部聚簇窗对时间门控能量图进行框选、对框选所得局部能量图进行阈值化处理后根据所有非零值像素建立局部像素集合、根据设置的闭邻域半径和闭邻域像素数目阈值，通过重复初始化像素簇及扩展像素簇过程，对局部像素集合中的像素进行聚簇，以及从聚簇所得局部像素簇中寻找局部像素主簇并将局部能量图中不包含在局部像素主簇中的所有像素的像素值置为零，从而得到单空化源能量图；

所述单空化源能量图像素峰值判断模块用于执行上述步骤4，主要是用于提取像素局部聚簇及寻找局部像素主簇子模块得到的所有单空化源能量图的像素峰值、根据设置的像素峰值阈值对其中每个单空化源能量图的像素峰值进行判断，以及由像素峰值小于像素峰值阈值的单空化源能量图建立单空化源能量图集合；

所述高斯拟合空化源定位模块用于执行上述步骤5，主要是用于在成像坐标系中建立高斯分布函数、对单空化源能量图像素峰值判断模块得到的单空化源能量图集合中的每个单空化源能量图进行高斯拟合以定位单个空化源、将拟合所得高斯分布函数中参数的估计值代入到高斯分布函数中计算高斯拟合结果、计算给定置信水平下的空化源坐标估计不确定度，以及根据单空化源能量图及单空化源能量图的高斯拟合结果计算拟合优度；参数的估计值中，高斯分布峰值位置的估计值为空化源坐标的估计值；

所述空化源定位检测模块用于执行上述步骤6，主要是用于根据高斯拟合空化源定位模块得到的所有空化源的空化源坐标估计不确定度和拟合优度建立空化源坐标估计不确定度集合和拟合优度集合、对各集合中元素分别进行升序排列并计算升序排列后的各集合中元素的四分位数和四分位距、根据四分位数和四分位距设置空化源坐标估计不确定度的阈值和拟合优度的阈值，以及利用设置的空化源坐标估计不确定度的阈值和拟合优度的阈值检测出正常定位的空化源；

所述空化源定位分布图计算及处理模块用于执行上述步骤7，主要是用于针对空化源定位检测模块得到的每个正常定位的空化源建立定位分布函数、计算规划的超分辨成像的像素网格中每个像素坐标的定位分布函数值，以及对计算所得所有正常定位的空化源的定位分布图叠加后进行归一化及对数化处理，从而得到超分辨成像结果。

本发明的有益效果体现在：

本发明提出的超分辨超声被动空化成像方法，通过时间门控能量图计算、基于像素全局聚簇和像素局部聚簇及寻找局部像素主簇的单空化源分离、单空化源能量图像素峰值判断、高斯拟合空化源定位、空化源定位检测和空化源定位分布图计算及处理，解决了超声被动空化成像分辨率受阵列超声换能器衍射模式限制的问题，使得超声被动空化成像的分辨率得到大幅度提高。本发明可以对小区域内空化源的空间分布进行微米级尺度的精细解析，为空化物理现象本身以及空化介导的多种治疗超声应用(如组织毁损、药物递送等)的研究提供了重要的技术手段，并为实现超声精准诊疗奠定了基础。

进一步地，本发明通过设置一初始时刻为双程延时、时间长度为空化激励超声换能器发射脉冲的脉冲长度的时间门控区间，从原始射频信号中提取出时间门控信号，然后对时间门控信号进行处理得到时间门控能量图，有效抑制阵列超声换能器轴向方向的干扰伪迹，为后续的像素聚簇空化源分离和高斯拟合空化源定位处理提供了高质量的原始图像，从而提高空化源定位的精度。

进一步地，本发明通过对时间门控能量图进行阈值化处理(所用阈值系数较大)并对非零值像素进行全局聚簇，初步地识别出时间门控能量图中包含的空化源；再通过划定局部聚簇窗对时间门控能量图进行框选、对局部能量图进行阈值化处理(所用阈值系数较小)并对非零值像素进行局部聚簇、以及寻找局部像素主簇，准确地分离出单个空化源的区域；这使得直接对单个空化源进行定位成为可能，从而提高空化源定位的精度。

进一步地，本发明通过设置像素峰值阈值对单空化源能量图的像素峰值进行判断，并将像素峰值大于等于像素峰值阈值的单空化源能量图剔除，使得空化源粘连导致的空化源定位不准确的问题得到解决，从而提高超分辨超声被动空化成像的准确度。

进一步地，本发明通过对单空化源能量图进行高斯拟合得到空化源坐标的估计值和空化源坐标估计不确定度，可实现对单个空化源的直接定位；其中，进行高斯拟合时采用的是允许角度偏转的高斯分布函数，该函数对于不同情境下任意偏转的能量波束具有普适性，有助于空化源定位精度的提高。

进一步地，本发明采用空化源坐标估计不确定度和拟合优度作为高斯拟合空化源定位的评价指标，根据四分位数和四分位距分别设置空化源坐标估计不确定度的阈值和拟合优度的阈值，并联合空化源坐标估计不确定度的阈值和拟合优度的阈值来检测空化源的定位是否正常，使得异常定位的空化源得以剔除，从而提高超分辨超声被动空化成像的准确度。

进一步地，本发明针对正常定位的空化源建立定位分布函数并以高密度的像素网格计算定位分布图，在此基础上通过叠加得到超分辨成像结果，从而在微米级尺度上对空化源空间分布进行精细的解析。

附图说明

图1为本发明实施例中时间门控能量图计算的流程图。

图2为本发明实施例中像素全局聚簇的流程图。

图3为本发明实施例中像素局部聚簇及寻找局部像素主簇的流程图。

图4为本发明实施例中单空化源能量图像素峰值判断及单空化源能量图集合建立的流程图。

图5为本发明实施例中高斯拟合空化源定位的流程图。

图6为本发明实施例中空化源定位检测的流程图。

图7为本发明实施例中空化源定位分布图计算及处理的流程图。

图8为本发明实施例中随机分布(a、d)、直线分布(b、e)、曲线分布(c、f)三种空化源分布下的超声被动空化成像结果；其中，a、b和c使用传统基于延时叠加积分的成像方法，d、e和f使用本发明提出的基于单空化源分离及定位的超分辨成像方法。

具体实施方式

下面结合附图和实施例对本发明作进一步的详细说明。所述实施例仅用于解释本发明，而非对本发明保护范围的限制。

本发明提出的基于单空化源分离及定位的超分辨超声被动空化成像方法，具体包括以下步骤：

参见图1，通过计算双程延时并设置时间门控区间的时间长度得到时间门控区间，从原始射频信号中提取各阵元的时间门控信号，沿阵元方向叠加该信号并沿采样点方向叠加所得时间门控合成信号的平方，从而得到时间门控能量图。具体流程见下述步骤(1.1)～(1.11)；

(1.1)以阵列超声换能器(例如，线阵换能器)的中心为原点，以阵列超声换能器的阵元方向为x轴、阵列超声换能器的中轴线方向为z轴，建立成像坐标系；在成像坐标系的基础上，以与阵列超声换能器的成像平面垂直的方向为y轴，建立三维坐标系；

(1.2)在步骤(1.1)所述成像坐标系中规划时间门控能量图的像素网格，该网格由x轴方向和z轴方向的成像范围(例如，-4～4mm和36～44mm)以及x轴方向和z轴方向的像素间隔(例如，0.1mm和0.1mm)确定；

(1.3)对于步骤(1.2)中像素网格中的任意一像素坐标(x,z)，计算该像素坐标到空化激励超声换能器表面中心坐标处切线或切面的距离

当空化激励超声换能器的中轴线在阵列超声换能器的成像平面内时，像素坐标(x,z)到空化激励超声换能器表面中心坐标处切线的距离

当空化激励超声换能器的中轴线与阵列超声换能器的成像平面相交且与阵列超声换能器的阵元方向垂直时，像素坐标(x,z)到空化激励超声换能器表面中心坐标处切面的距离

其中，x

(1.4)计算步骤(1.3)中所述像素坐标(x,z)与阵列超声换能器第i个阵元坐标(x

其中，i＝1,2,…,N

(1.5)根据步骤(1.3)和(1.4)所得距离

其中，c为声传播速度；

(1.6)将阵列超声换能器第i个阵元的时间门控区间的初始时刻设置为步骤(1.5)所得双程延时τ

(1.7)根据步骤(1.6)所得第i个阵元的时间门控区间从阵列超声换能器被动接收所得的一帧原始射频信号中提取出第i个阵元的时间门控信号

其中，k′＝1,2,…,TGS

(1.8)将步骤(1.7)所得所有阵元的时间门控信号沿阵元方向叠加，得到时间门控合成信号

(1.9)将步骤(1.8)所得时间门控合成信号

(1.10)重复步骤(1.3)～(1.9)，直至计算完所有像素坐标的时间门控能量值，由所有像素坐标的时间门控能量值构成时间门控能量图(即时间门控能量图各像素坐标处的像素值为对应的时间门控能量值)；

(1.11)对空化激励超声换能器发射N

参见图2，对每帧时间门控能量图进行阈值化处理并建立全局像素集合，设置闭邻域半径和闭邻域像素数目阈值，通过重复初始化像素簇及扩展像素簇过程得到每帧时间门控能量图的全局像素簇。具体流程见下述步骤(2.1)～(2.7)；

(2.1)设置阈值系数γ

(2.2)根据步骤(2.1)所得所有非零值像素建立全局像素集合S

(2.3)设置闭邻域半径ε(以像素数目为单位，例如，

(2.4)初始化全局像素簇：

(2.4.1)从步骤(2.2)所得全局像素集合S

(2.4.2)若所述像素(u

(2.5)扩展获得全局像素簇：

(2.5.1)将步骤(2.4.2)所得核心像素

(2.5.2)将步骤(2.4.2)所得核心像素

(2.5.3)在待判断像素集合S

(2.5.4)若步骤(2.5.3)中的像素(u

(2.5.5)重复步骤(2.5.3)和(2.5.4)，直至遍历完待判断像素集合S

(2.6)重复步骤(2.4)和(2.5)，直至遍历完步骤(2.2)所得全局像素集合S

(2.7)重复步骤(2.1)～(2.6)，得到N

参见图3，针对每帧时间门控能量图的每个全局像素簇划定局部聚簇窗，对该窗框选得到的局部能量图进行阈值化处理并建立局部像素集合，对局部像素集合中的像素进行聚簇得到局部像素簇，从局部像素簇中寻找局部像素主簇，从而得到单空化源能量图。具体流程见下述步骤(3.1)～(3.9)；

(3.1)针对第p帧时间门控能量图的第q个全局像素簇(p＝1,2,…,N

(3.2)利用步骤(3.1)所得局部聚簇窗对第p帧时间门控能量图进行框选，得到包含所述全局像素簇的质心对应的像素且尺寸确定(例如，5×5)的局部能量图；

(3.3)设置阈值系数γ

(3.4)根据步骤(3.3)所得所有非零值像素建立局部像素集合S

(3.5)根据步骤(2.4)～(2.6)中的方法对步骤(3.4)所得局部像素集合S

(3.6)设置标签像素为步骤(3.1)所述全局像素簇的质心对应的像素，并根据标签像素从步骤(3.5)所得若干个局部像素簇中寻找局部像素主簇；即判断标签像素是否在某一局部像素簇中，若在，则该局部像素簇为局部像素主簇；

(3.7)将步骤(3.2)所得局部能量图中不包含在步骤(3.6)所得局部像素主簇中的所有像素的像素值置为零，得到单空化源能量图，从而实现单空化源的分离；

(3.8)重复步骤(3.1)～(3.7)，得到第p帧时间门控能量图的

(3.9)重复步骤(3.1)～(3.8)，得到N

参见图4，提取所有单空化源能量图的像素峰值，设置像素峰值阈值并对每个单空化源能量图的像素峰值进行判断，由像素峰值小于像素峰值阈值的单空化源能量图建立单空化源能量图集合。具体流程见下述步骤(4.1)～(4.5)；

(4.1)提取所有单空化源能量图的像素峰值，即分别提取步骤(3.9)所得N

(4.2)设置像素峰值阈值，例如设置该阈值为步骤(4.1)所得所有单空化源能量图的像素峰值的最大值的0.5倍；

(4.3)初始化单空化源能量图集合为一空集；

(4.4)对步骤(3.9)所得N

(4.5)重复步骤(4.4)，直至对N

参见图5，在成像坐标系中建立允许角度偏转的高斯分布函数，对单空化源能量图集合中的各单空化源能量图进行高斯拟合，计算空化源坐标估计不确定度和拟合优度。具体流程见下述步骤(5.1)～(5.7)；

(5.1)在步骤(1.1)所述成像坐标系中建立允许角度偏转的高斯分布函数：

其中，A为高斯分布峰值，μ

(5.2)针对步骤(4.5)所得单空化源能量图集合中的第n个单空化源能量图SCE

(5.3)根据步骤(5.1)所述高斯分布函数和步骤(5.2)所得各参数的初始化值，对第n个单空化源能量图SCE

(5.4)将步骤(5.3)所得各参数的估计值代入到步骤(5.1)所述高斯分布函数中，计算得到单空化源能量图SCE

(5.5)设定置信水平(例如，0.68)，根据该置信水平下x轴方向和z轴方向的空化源坐标估计的置信区间

(5.6)根据步骤(5.2)所述单空化源能量图SCE

其中，SCE

(5.7)重复步骤(5.2)～(5.6)，直至对单空化源能量图集合中的

参见图6，建立空化源坐标估计不确定度集合和拟合优度集合，根据集合中元素的四分位数和四分位距设置空化源坐标估计不确定度和拟合优度的阈值，利用所得阈值检测每帧时间门控能量图中空化源的定位。具体流程见下述步骤(6.1)～(6.7)；

(6.1)以步骤(5.7)所得

(6.2)对步骤(6.1)所得集合Δ

(6.3)计算步骤(6.2)所得升序排列后的集合Δ

其中，M表示元素数目为

令以上两式中的M为Δ

(6.4)根据步骤(6.3)所得结果分别计算集合Δ

/>

(6.5)根据步骤(6.3)所得四分位数和步骤(6.4)所得四分位距，分别设置空化源坐标估计不确定度δ

(6.6)对第n个空化源的定位(n＝1,2,...,

(6.7)重复步骤(6.6)，直至检测完

参见图7，针对正常定位的空化源建立定位分布函数并计算定位分布图，叠加所有正常定位的空化源的定位分布图后进行归一化及对数化处理，得到超分辨成像结果。具体流程见下述步骤(7.1)～(7.4)；

(7.1)在步骤(1.1)所述成像坐标系中规划超分辨成像的像素网格，该像素网格的像素数目为步骤(1.2)所述像素网格的像素数目的λ

(7.2)针对步骤(6.7)所得

该函数为允许角度偏转的高斯分布函数，其峰值位置、峰值、以及偏转角分别为对相应单空化源能量图进行高斯拟合所得的x方向和z方向的空化源坐标的估计值

(7.3)根据步骤(7.2)所得定位分布函数，计算步骤(7.1)所得像素网格中每个像素坐标的定位分布函数值，得到第s个正常定位的空化源的定位分布图；

(7.4)重复步骤(7.2)～(7.3)，得到

以下采用数值仿真方法检验以上提出的超分辨超声被动空化成像方法的性能，仿真的空化源分布有三种：随机分布、直线分布和曲线分布。图8a、图8b和图8c分别为随机分布、直线分布和曲线分布下使用传统基于延时叠加积分的成像方法得到的超声被动空化成像结果(以线性尺度显示)，图8d、图8e和图8f分别为随机分布、直线分布和曲线分布下使用本发明提出的基于单空化源分离及定位的超分辨成像方法得到的超声被动空化成像结果(以对数尺度显示)。图8a、图8b和图8c所示成像结果中出现若干个尺寸很大的亮斑，成像分辨率很差，无法有效地反映空化源空间分布信息；相比而言，图8d、图8e和图8f所示成像结果中空化源的空间分布在一个非常小的尺度下(可达微米级)被精细地解析出来，成像分辨率大幅度提高，有效地反映了空化源空间分布信息。

本发明具有以下优点：

(1)超声被动空化成像存在空间分辨率差的问题，尽管目前已发展的多种改进成像方法能在一定程度上提高图像质量，但并不能从根本上克服成像换能器衍射模式对成像分辨率的限制。而本发明通过时间门控能量图计算、基于像素全局聚簇和像素局部聚簇及寻找局部像素主簇的单空化源分离、单空化源能量图像素峰值判断、高斯拟合空化源定位、空化源定位检测和空化源定位分布图计算及处理，实现了对空化源空间分布的超分辨成像，从根本上解决了成像分辨率受换能器衍射模式限制的问题，使得超声被动空化成像的分辨率得到大幅度提高。

(2)传统超声被动空化成像方法生成的图像中会出现沿阵列超声换能器轴向方向的高水平干扰伪迹，利用该图像对空化源进行基于高斯拟合的定位会造成较大的定位误差。而本发明提取出初始时刻为双程延时、时间长度为空化激励超声换能器发射脉冲的脉冲长度的时间门控区间内的信号，据此计算的时间门控能量图中干扰伪迹得到有效抑制，为后续的像素聚簇空化源分离和高斯拟合空化源定位处理提供了高质量的原始图像，从而提高空化源定位的精度。

(3)本发明首先通过高阈值化处理时间门控能量图和像素全局聚簇对时间门控能量图中包含的空化源进行初步的识别，然后通过划定局部聚簇窗框选时间门控能量图、低阈值化处理局部能量图、像素局部聚簇和寻找局部像素主簇对单个空化源所在的区域进行准确的分离，从而允许直接对单个空化源进行定位，提高了空化源定位的精度；本发明进一步通过对单空化源能量图的像素峰值进行判断，使得空化源粘连导致的空化源定位不准确的问题得到解决，从而提高超分辨超声被动空化成像的准确度。

(4)本发明对单空化源能量图进行高斯拟合时采用的是允许角度偏转的高斯分布函数，该函数对于不同情境(例如，空化激励超声换能器与阵列超声换能器以不同的方式放置等)下任意偏转的能量波束具有普适性，从而有助于空化源定位精度的提高。

(5)本发明采用空化源坐标估计不确定度和拟合优度来评价高斯拟合空化源定位并基于四分位法设置这两种指标的阈值，将异常定位的空化源检测出来并在叠加空化源定位分布图时将异常定位的空化源剔除(即只叠加正常定位的空化源的定位分布图)，从而使得超分辨超声被动空化成像的准确度得到提高。

(6)本发明对正常定位的空化源建立定位分布函数并在高密度的像素网格中计算定位分布图，进一步通过叠加生成超分辨成像结果，从而使得空化源空间分布得到微米级尺度的精细解析。

(7)本发明为空化激励超声换能器激励下产生的空化源的空间分布的解析，特别是为小区域内空化源的空间分布的微米级尺度的精细解析提供强有力的技术手段，可以对组织毁损、药物递送等多种高强度或低强度治疗超声应用中的空化活动进行精细的监控成像，并进一步深入地澄清这些应用背后的生物物理机制，有利于促进超声精准诊疗及诊疗一体化技术的发展。