CN104535645B - 微秒分辨空化时空分布的三维空化定量成像方法 - Google Patents

Abstract

本发明提供一种微秒分辨空化时空分布的三维空化定量成像方法，使用宽波束检测空化后，阵列换能器移动一个单元位置，待空化核分布恢复，再在同样的空化能量激励下，再次宽波束检测空化，得到不同单元位置的空间系列二维空化原始射频数据，通过改变空化源能量作用时间、能量源装置激励与阵列换能器发射宽波束之间的时间延迟以及脉动泵与能量源装置激励之间的时间延迟，可得到时间序列二维空化原始射频数据，然后结合宽波束最小方差自适应波束合成、Nakagami参量成像、三维重建算法得到微秒分辨空化三维时空分布图像和空化微泡密度的定量图像，其时间分辨率可达到微秒，该方法有潜力发展成为类似于声场测量的三维空化成像标准方法。

Description

微秒分辨空化时空分布的三维空化定量成像方法

技术领域

本发明涉及空化物理与应用及超声成像技术领域，该方法结合阵列Plane-by-plane宽波束空化检测，宽波束最小方差自适应波束合成以及Nakagami参量成像算法，实现稳态自由场和脉动流条件下微秒分辨空化三维时空分布成像和空化微泡密度的定量成像。

背景技术

空化是指液体中的空化核在外加能量(热/力)的作用下被激活，出现微小泡核的振荡、生长、收缩乃至崩溃等一系列动力学过程，是生物医学领域中药物释放、基因转染、体外碎石、溶栓、止血、热疗以及肿瘤热消融等方面的主要机制。空化的过程包括以下几个阶段：空化成核、空化泡线性和非线性振动，空化泡生长、空化急速收缩至坍塌破裂以及空化泡消散，可分为以非惯性空化为特点的稳态空化和惯性空化为特点的瞬态空化。在液体介质中，产生空化的最小能量值称为空化起始阈值，其大小取决于媒介液体静态压、初始温度、液体本身的结构状态以及液体中外加的多样性空化核，因此液体媒介中空化的产生具有一定的随机性，但相同环境以及空化能量作用时，其空化泡群形状及分布具有可重复性。目前，为更好地研究不同媒介空化的产生机制以便更好地控制和利用空化，需要研究有效的空化检测与成像方法。

现有的空化检测与成像主要有光学和声学方法。光学检测成像主要包括通过高速/超高速摄影、声致发光以及声致化学发光等，可拍摄观察空化泡的行为以及时空动态特性，具有直观、同步性好、时间分辨率高的优点，缺点是一方面对媒介透光性要求很高且不适用于原位研究，另一方面所得图像是沿光穿透方向信息的重叠。声学检测方法是基于空化过程中或空化微泡产生的声信息，包括谐波、次谐波、超谐波和宽带噪声等，其中得到最广泛应用的就是被动空化检测(PassiveCavitationDetection,PCD)和主动空化检测(ActiveCavitationDetection,ACD)。PCD利用换能器被动接收由空化微泡所产生的声散射信号，而ACD采用低压脉冲回波探测可能发生空化的区域，但PCD和ACD由于一般采用单阵元换能器，受限于有限的空间检测区域，无法提供空化微泡的空间分布。

在PCD和ACD基础上，使用二维阵列换能器作为空化检测换能器，发展出被动空化成像(PassiveCavitationImaging,PCI)和主动空化成像(ActiveCavitationImaging,ACI)。由于超声空化具有瞬态特性，空化微泡的振动、坍塌破裂以及消散的时间都是微秒级，因此空化成像方法的时间分辨率需要达到微秒。同时，针对空化的瞬态特性，有必要得到空化微泡的时空分布，包括不同空化能量源作用时间和空化消散随时间的序列时空空化分布。PCI通过阵列换能器被动接收和通道信号源重建得到空化泡的二维空间分布，重建算法复杂且空间分辨率不高。ACI包括常规的B超成像和超快速主动空化成像方法。由于B超图像是通过逐线扫描得到的，同一帧图像不同扫描线之间存在时间差，且时间分辨率无法达到微秒级。而超快速主动空化成像由于发射的是平面波，其灵敏度以及横向分辨率有待改进，且其时间分辨率为几百个微秒，无法满足研究空化瞬态分布的要求。

在空化成像的基础上需要对空化泡进行定征，包括空化量化、空化尺寸及密度分布等。当前的空化量化方法主要有惯性空化剂量和非惯性空化剂量，是通过计算特定频段内宽带噪声或次谐波幅度的均方根值作为空化强度的一种相对度量，可分别衡量瞬态空化和稳态空化的相对大小，但一般针对PCD所采集到的一维射频数据，这种量化方法无法反应空化强度分布。现有的空化密度检测方法有激光相位多普勒法，该方法主要针对空化泡在不同尺寸上的分布，而关于空化泡在不同空间位置的密度分布还没有研究，无法提供空间信息。

目前已有的空化检测与成像方法局限于一维和二维，实际中空化泡的分布区域遍布整个焦域甚至更大，而且在临床应用中如聚焦超声治疗时，其声波传播路径上可能存在其他组织介质，使得声场分布发生变化而出现不对称性，因此有必要发展一种微秒分辨空化三维时空分布成像和空化微泡密度的定量成像方法。此外，对于流动条件尤其是脉动流条件下的空化研究相对较少，而人体的血流是一种脉动流，因此有必要研究其条件下的三维时空空化分布尤其是脉动流周期内不同时间点的三维空化分布。

发明内容

针对上述现有技术的缺陷以及微秒分辨空化时空分布的三维定量成像的必要性，本发明的目的在于提供一种稳态自由场和脉动流条件下具有微秒分辨空化时空分布的三维空化定量成像方法。

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

一种稳态自由场下微秒分辨空化时空分布的三维空化成像和空化密度定量成像方法，其原理是：

在稳态自由场条件下，采用Plane-by-plane宽波束检测空化，克服同帧空化图像空间不同步的缺点，每次宽波束检测空化后，阵列换能器移动一个单元位置，并等待足够长的时间使得介质空化核分布恢复到空化能量源作用之前的初始状态，再在同样的空化能量激励下，再次宽波束检测空化，以此逐步得到不同单元位置的系列二维空化原始射频数据，再结合宽波束最小方差自适应波束合成(Synnevag,J.F.,A.Austeng,etal.(2007)."Adaptivebeamformingappliedtomedicalultrasoundimaging."IeeeTransactionsonUltrasonicsFerroelectricsandFrequencyControl54(8):1606-1613)、Nakagami参量成像和三维重建算法，得到稳态自由场下微秒分辨空化时空分布的三维空化图像和和空化微泡密度的定量图像。

步骤一、采用阵列Plane-by-plane宽波束检测空化，在源能量温度或压力连续可调的情况下激励空化的产生，采集空化信号：空化的产生装置包括产生能量场的能量源装置和控制时序的同步信号发生器；空化信号的检测装置包括可编程发射宽波束的阵列换能器和并行通道数据采集及存储单元；同步信号发生器产生同步信号分别控制能量源装置和阵列换能器，能量源装置产生连续可变能量激励空化的产生，阵列换能器发射宽波束对空化进行检测，得到的空化回波信号由并行通道数据采集及存储单元采集存储；等待一个足够长的时间使得媒介空化核分布恢复到初始状态，通过三维机械扫描装置控制阵列换能器沿垂直于阵列放置方向移动一个单元位置，相同参数的空化能量源再次激励产生空化，同步阵列换能器发射宽波束采集空化射频数据；重复以上过程，可得到相同参数的空化能量源条件下，沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像原始射频数据。

步骤二、采用宽波束最小方差自适应波束合成对采集到的二维空化分布图像原始射频数据进行处理，得到沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像：首先选定某一单元位置二维空化成像区域中的某个目标点，根据目标点的位置计算有效孔径及延时得到对目标点进行接收聚焦后的通道信号；然后对通道信号进行最小方差自适应波束合成，实现目标点通道信号的最优幅度变迹并得到最优输出；遍历该二维空化成像区域所有目标点得到该单元位置二维空化成像波束合成后的射频数据；以上过程遍历所有单元位置得到一系列二维空化分布的射频数据；最后通过射频成像算法对各二维空化射频数据进行成像，得到一系列高分辨、高信噪比的二维空化分布图像。

步骤三、对步骤二中得到的沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像一方面进行三维重建，得到三维空化分布图像，另一方面在步骤二中得到的波束合成后的射频数据取包络后进行Nakagami参量提取得到空化密度分布：利用三维重建算法，将不同单元位置的序列二维空化分布图像进行三维显示；对波束合成后的射频数据先进行Nakagami参量计算，得到不同单元位置的序列二维空化密度定量图像，再经过三维重建算法，得到空化密度定量图像三维显示。

步骤四、改变产生空化的能量源装置的参数以及媒介，重复步骤一到三，得到不同条件下的空化以及空化密度定量三维图像：改变空化能量源作用时间，其时间分辨可达到几个微秒，重复步骤一到三，得到随时间演化的空化以及空化密度定量三维序列图像；改变空化能量源能量大小，重复步骤一到三，得到随能量变化的空化以及空化密度定量三维序列图像；改变能量源装置激励与阵列换能器发射宽波束之间的时间延迟，重复步骤一到三，得到随时间消散的空化以及空化密度定量三维序列图像；改变空化能量源作用的媒介，重复步骤一到三，得到不同媒介下空化以及空化密度定量三维序列图像。

步骤三中所述Nakagami参量提取，具体方法为：

1)将射频数据rf进行相关去噪处理，具体为：

a)选取一定尺寸大小的背景信号区域，计算该区域平均能量P；

b)对rf分别叠加平均能量为P的随机高斯白噪声n₁和n₂，得到S₁和S₂；

c)计算S₁和S₂的相关系数，并给定阈值Th，对相关系数进行阈值化处理后与rf加权，得到rf_denoise如下，其中corrcoef(S₁,S₂)表示信号相关系数：

$r{f}_{\mathrm{denoise}}=\left\{\begin{array}{cc}\mathrm{rf},& \mathrm{corrcoef}(S_1,S_2)>\mathrm{Th}\\ 0,& \mathrm{corrcoef}(S_1,S_2)\≤\mathrm{Th}\end{array}\right.---\left(1\right)$

2)将去噪后的空化射频数据rf_denoise进行Hilbert解调，得到包络信号，记为R；

3)计算Nakagami参量如下：

$m=\frac{{\left[E\left(R^2\right)\right]}^2}{E{[R^2-E\left(R^2\right)]}^2}---\left(3\right)$

其中E(·)代表统计均值。

一种脉动流条件下的微秒分辨空化时空分布的三维空化成像和空化密度定量成像方法，其原理是在脉动流条件下，由于脉动周期时间尺度远大于空化瞬态改变时间尺度，通过脉动周期倍数的等待时间可使得空化核分布恢复到初始状态，采用Plane-by-plane宽波束检测空化，改变脉动泵和空化能量源装置的时间延迟，可得到脉动周期内不同时刻，而改变阵列换能器位置，可得到不同单元位置的系列二维空化分布图像原始射频数据，然后再结合宽波束最小方差自适应波束合成、Nakagami参量成像和三维重建算法，得到脉动流下微秒分辨空化时空分布的三维空化图像和空化微泡密度的定量图像，包括以下步骤：

步骤一、利用仿血管组织体模模拟血管及其周围的组织，管道两边连接橡皮管，并与脉动泵相连，连接仿血管组织体模与三维机械扫描装置，通过三维移动将其放置于空化源作用区域，并使空化产生在仿血管组织体模的管道内。

步骤二、测试不同脉动泵参数设置下的流动规律：脉动液体为生理盐水造影剂微泡溶液，设置脉动泵参数，包括每分钟脉动次数以及流量，观察多普勒谱并记录脉动周期和规律，最后管道内流动去离子水将管道内的生理盐水造影剂微泡溶液冲洗干净。

步骤三、采用阵列Plane-by-plane宽波束检测空化，在源能量温度或压力连续可调的情况下激励空化的产生，采集空化信号：设置脉动周期和流量，并以此编辑时序，控制脉动流装置、空化源装置和空化信号检测装置，脉动流装置包括产生脉动流的脉动泵、仿血管组织体模和控制其工作的同步信号发生器；空化源装置包括产生能量场的能量源装置和控制时序的同步信号发生器，空化信号的检测装置包括可编程发射宽波束的阵列换能器和并行通道数据采集及存储单元；同步信号发生器产生同步信号分别控制脉动泵、能量源装置和阵列换能器，脉动泵使液体按照一定脉动周期和流量流入流出仿血管组织体模，能量源装置产生连续可变能量激励空化的产生，阵列换能器发射宽波束对空化进行检测，得到的空化回波信号由并行通道数据采集及存储单元采集存储。

步骤四、阵列换能器位置不变，改变脉动泵触发与空化能量源装置激励之间的时间延迟，重复步骤三，得到脉动周期内各时刻的序列二维空化分布图像原始射频数据：固定阵列换能器，依据脉动流多普勒谱记录的脉动周期和规律，改变脉动泵与能量源装置的时间延迟，等待一个足够长的时间使得媒介空化核分布恢复到初始状态，重复步骤三，得到脉动周期内不同时间点的二维空化分布图像原始射频数据。

步骤五、改变阵列换能器位置，重复步骤三、四，得到不同空间位置及其脉动周期内各时刻的序列二维空化分布图像原始射频数据：通过三维机械扫描装置控制阵列换能器沿垂直于阵列放置方向移动一个单元位置，等待一个足够长的时间使得媒介空化核分布恢复到初始状态，相同参数的空化能量源再次激励产生空化，阵列换能器发射宽波束采集空化射频数据，重复以上过程，可得到相同参数的空化能量源条件下，沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像原始射频数据。不同单元位置时，重复步骤四，可得到该位置下脉动周期内不同时间点的二维空化分布图像原始射频数据。

步骤六、采用宽波束最小方差自适应波束合成对采集到的二维空化分布图像原始射频数据进行处理，得到沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像：首先选定某一单元位置二维空化成像区域中的某个目标点，根据目标点的位置计算有效孔径及延时得到对目标点进行接收聚焦后的通道信号；然后对通道信号进行最小方差自适应波束合成，实现目标点通道信号的最优幅度变迹并得到最优输出；遍历该二维空化成像区域所有目标点得到该单元位置二维空化成像波束合成后的射频数据；以上过程遍历所有单元位置得到一系列二维空化分布的射频数据；最后通过射频成像算法对各二维空化射频数据进行成像，得到一系列高分辨、高信噪比的二维空化分布图像。

步骤七、对步骤六中得到的沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像一方面进行三维重建，得到三维空化分布图像，另一方面在步骤六中得到的波束合成后的射频数据取包络后进行Nakagami参量提取得到空化密度分布：利用三维重建算法，将不同单元位置的序列二维空化分布图像进行三维显示；对波束合成后的射频数据先进行Nakagami参量计算，得到不同单元位置的序列二维空化密度定量图像，再经过三维重建算法，得到空化密度定量图像三维显示。

步骤八、改变产生空化的能量源装置的参数以及媒介，重复步骤三到七，得到脉动流条件下不同参数时的空化以及空化密度定量三维图像：改变空化能量源作用时间，其时间分辨可达到几个微秒，重复步骤三到七，得到随空化源作用时间演化的脉动流条件下空化以及空化密度定量三维序列图像；改变空化能量源能量大小，重复步骤三到七，得到随能量变化的脉动流条件下空化以及空化密度定量三维序列图像；改变能量源装置激励与阵列换能器发射宽波束之间的时间延迟，重复步骤三到七，得到随时间消散的脉动流条件下空化以及空化密度定量三维序列图像；改变空化能量源作用的媒介，重复步骤三到七，得到脉动流条件下不同媒介时空化以及空化密度定量三维序列图像。

针对稳态自由场和脉动流条件下的空化时空分布研究，本发明首先使用阵列Plane-by-plane宽波束检测空化，采集得到沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像原始射频数据，再结合宽波束最小方差自适应波束合成得到不同单元位置的二维空化分布序列图像，然后使用Nakagami参量算法提取得到反映空化密度的定量图像，最后使用三维重建算法得到微秒分辨空化的三维定量图像，改变空化源能量作用时间以及其与宽波束检测之间的时间延迟，可得到时间序列空化图像，而改变产生空化的能量源装置的参数以及媒介，可得到空化以及空化密度定量三维序列图像，该方法有潜力发展成为类似于声场测量的三维空化成像标准方法。

本发明与现有技术相比，具有下列优点：

与现有技术相比，本发明首先采用Plane-by-plane宽波束检测空化，具有微秒时间分辨率，可观察空化泡群瞬态分布；其次结合宽波束最小方差自适应波束合成得到的空化图像具有较高的空间分辨率和信噪比；然后通过Nakagami参量提取可得到空化微泡密度分布的定量图像；再次，使用三维重建算法将不同单元位置的二维空化分布序列图像转化为空化以及空化密度定量三维图像；最后，改变产生空化的能量源装置作用时间，可得到随时间演化的空化以及空化密度定量三维时空分布图像，其时间分辨可达到几个微秒；改变能量源装置激励与阵列换能器发射宽波束之间的时间延迟，可得到随时间消散的空化以及空化密度定量三维时空分布图像；改变脉动与能量源装置激励之间的时间延迟，可得到脉动周期内不同时间点的空化及空化密度定量时空分布图像；改变空化能量源能量大小，可得到随能量变化的空化以及空化密度定量三维序列图像；改变空化能量源作用的媒介，得到不同媒介时空化以及空化密度定量三维序列图像。

附图说明

图1是本发明阵列Plane-by-plane宽波束空化检测装置示意图；

图2是本发明的微秒分辨空化的三维定量成像流程图；

图3是本发明的Nakagami参量成像流程图；

图4是本发明脉动流微秒空化成像装置示意图；

图5是本发明的脉动流微秒空化三维定量成像流程图；

图6、图7、图8是本发明以高强度聚焦超声换能器作为能量源激励自来水中空化以及空化密度定量三维图像结果，高强度聚焦超声脉冲长度10us-200ms，电功率100W，箭头所指为超声波作用方向。图6(a)为二维时空分布空化图像，图6(b)为二维时空空化微泡密度定量图像，图7是脉冲长度为20ms时阵列换能器位置从-5mm到6mm的空间序列空化成像，图8(a)是脉动长度为20ms时三维空化图像，图8(b)是脉动长度为20ms时三维空化微泡密度图像；

图9是本发明以高强度聚焦超声换能器作为能量源激励脉动流时空化以及空化密度定量三维图像结果，高强度聚焦超声脉冲长度20us，电功率100W，箭头所指为超声波作用方向。

图中所示的序号为：空化能量源装置1，功率放大器2，波形发生器3，全数字化超声设备4，第一三维机械扫描装置5，第二三维机械扫描装置5'，阵列换能器6，水槽7，吸声材料8，脉动泵9，仿血管组织体模10，储存池11，废液池12，乳胶管13。

具体实施方式

以下结合附图和实施例对本发明作详细说明。

参见图1，一种阵列Plane-by-plane宽波束空化检测装置，包括空化产生装置和空化信号检测装置，前者包括产生空化的能量源装置1、功率放大器2和控制时序的同步波形发生器3，后者包括全数字化超声设备4的可编程发射宽波束阵列换能器6和并行通道数据采集及存储单元。波形发生器3的第一通道驱动功率放大器2，激励空化能量源装置产生能量场，一方面通过波形发生器3的波形编辑可控制能量源的参数，包括能量作用时间以及工作模式，其时间分辨率可达到微秒，另一方面通过功率放大器2的面板可控制能量源的大小。波形发生器3的第二通道驱动全数字化超声设备4发射宽波束和采集原始空化射频数据，其发射接收参数可通过对全数字化超声设备4编程实现。波形发生器3通过第一通道和第二通道实现空化产生装置和空化信号检测装置之间的时间同步，可根据需求设置二者之间的时间延迟。当能量源产生空化和空化信号采集之后，等待一个足够长的时间使得媒介空化核分布恢复到初始状态，通过三维机械扫描装置5沿垂直于阵列换能器6的方向移动一个单元位置，单元位置可通过三维机械扫描装置5进行设置，相同参数的空化能量源装置1再次激励产生空化，同步阵列换能器6发射宽波束采集空化射频数据，重复以上过程，得到相同参数的空化能量源条件下，沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像原始射频数据。以上操作在水槽7中进行，并在水槽7的底部和侧壁放置吸声材料8。

参见图2，微秒分辨空化的三维定量成像流程，其步骤如下：

①采用阵列Plane-by-plane宽波束检测空化，可编程全数字化超声设备4发射宽波束检测空化，得到空化信号原始射频数据，接着等待一个足够长的时间使得媒介空化核分布恢复到初始状态，通过三维机械扫描装置5沿垂直于阵列换能器6的方向移动一个单元位置，单元位置大小可通过三维机械扫描装置5进行设置，相同参数的空化能量源装置1再次激励产生空化，同步阵列换能器6发射宽波束采集空化射频数据，重复以上过程，得到相同参数的空化能量源条件下，沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像原始射频数据；

②采用最小方差自适应波束合成算法对所述二维空化原始射频数据进行波束合成，得到沿垂直于阵列放置方向不同单元位置的一系列二维空化射频数据，随后包络检波、对数压缩以及坐标变换得到一系列二维空化分布图像；

③采用Nakagam参量成像算法对一系列二维空化射频数据进行定量成像，提取空化密度参量，得到沿垂直于阵列放置方向不同单元位置的一系列二维空化密度图像；

④采用三维重建算法，对沿垂直于阵列放置方向不同单元位置的序列二维空化分布图像和二维空化密度图像进行三维重建，得到空化能量源作用下的三维空化分布图像和三维空化密度图像；

⑤改变产生空化的能量源装置的参数以及媒介，重复步骤①-④，得到空化以及空化密度定量三维序列图像。

参见图3，Nakagami参量成像算法流程，其步骤如下：

①取宽波束最小方差自适应波束合成后的射频数据记为rf，计算背景区域平均能量，记为P；

②对rf叠加平均能量均为P的高斯白噪声信号得到加噪信号S₁和S₂；

③给定阈值Th，对S₁和S₂进行相关性分析并对相关系数阈值化处理得到二值化相关系数矩阵Coef_corr；

④将rf与③中得到的二值化矩阵Coef_corr进行加权得到相关去噪后的射频数据rf_denoise：

$r{f}_{\mathrm{denoise}}=\left\{\begin{array}{cc}\mathrm{rf},& \mathrm{corrcoef}(S_1,S_2)>\mathrm{Th}\\ 0,& \mathrm{corrcoef}(S_1,S_2)\≤\mathrm{Th}\end{array}\right.$

⑤对④中得到的rf_denoise进行Hilbert变化得到包络信号R；

由Nakagami统计模型得到R的概率密度函数为：

$f\left(R\right)=\frac{2m^m{R}^{2m-1}}{\mathrm{\Γ}\left(m\right){\mathrm{\Ω}}^m}\exp (-\frac{m}{\mathrm{\Ω}}{R}^2)U\left(R\right)$

其中Γ(·)和U(·)分别代表gamma函数和单位阶跃函数，m和Ω分别为Nakagami统计模型的Nakagami参量和尺度参量；

⑥利用公式计算Nakagami参量。

参见图4，一种脉动流的微秒分辨空化三维定量成像装置，包括脉动流装置、空化产生装置和空化信号检测装置。脉动流装置包括脉动泵9、仿血管组织体模10和触发其工作的波形发生器3。空化产生装置包括产生空化的能量源装置1、功率放大器2和控制时序的同步波形发生器3，而空化信号检测装置包括全数字化超声设备4的可编程发射宽波束阵列换能器6和并行通道数据采集及存储单元。波形发生器3的第一通道驱动功率放大器2，激励空化能量源装置产生能量场，一方面通过波形发生器3的波形编辑可控制能量源的参数，包括能量作用时间以及工作模式，其时间分辨率可达到微秒，另一方面通过功率放大器2的面板可控制能量源的大小。波形发生器3的第二通道驱动全数字化超声设备4的发射宽波束和采集原始空化射频数据，其发射接收参数可通过对全数字化超声设备4编程实现。波形发生器3的Marker通道触发脉动泵9。波形发生器3通过第一通道、第二通道和Marker通道实现脉动流装置、空化产生装置和空化信号检测装置之间的时间同步，可根据需求设置时间延迟。将仿血管组织体模10的管道两端连接乳胶管13，并与脉动泵9相连，将液体从储存池11经管道流入废液池12，以此模拟脉动血流。连接仿血管组织体模10与第二三维机械扫描装置5'，通过三维移动将其放置于空化源作用区域，并使空化产生在仿血管组织体模10的管道内。全数字化超声设备4采用多普勒模式，测试不同脉动泵9参数下的多普勒谱，记录脉动周期和规律，以此编辑波形发生器3的同步波形时序。全数字化超声设备4的可编程发射宽波束阵列换能器6与第一三维机械扫描装置5相连，当能量源产生空化和空化信号采集之后，等待一个足够长的时间使得媒介空化核分布恢复到初始状态，通过第一三维机械扫描装置5沿垂直于阵列换能器6的方向移动一个单元位置，单元位置大小可通过第一三维机械扫描装置5进行设置，相同参数的空化能量源装置1再次激励产生空化，同步阵列换能器6发射宽波束采集空化射频数据，重复以上过程，得到相同参数的空化能量源条件下，沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像原始射频数据。以上操作在水槽7中进行，并在水槽7的底部和侧壁放置吸声材料8。

参见图5，脉动流微秒分辨空化三维定量成像流程，其步骤如下：

①制备仿血管组织体模10，将凝胶溶液倒入模具，待冷却凝固后慢慢取出管子，在凝固凝胶内部形成管道。

②定位仿血管组织体模10，先连接仿血管组织体模10和第二三维机械扫描装置5'，通过三维移动将其放置于空化源作用区域，并使空化产生在仿血管组织体模10的管道内。

③搭建脉动流装置，仿血管组织体模10的管道两边连接乳胶管13，并与脉动泵9相连，从储存池11中将液体经管道流入废液池12，模拟脉动血流。

④测量脉动流流动，储存池11中液体为生理盐水造影微泡溶液，将全数字化超声设备4改为脉冲多普勒方式，记录不同脉动泵参数下的脉动周期和规律。

⑤根据记录的脉动周期，编辑波形发生器3的通道1、通道2和Marker波形，分别输入给功率放大器2激励空化源能量装置1产生空化、全数字化超声设备4控制阵列换能器6发射接收平面波声信号、以及脉动泵9控制脉动流周期内各时刻与空化能量源装置1的时间延迟。

⑥采用阵列Plane-by-plane宽波束检测空化，可编程全数字化超声设备4发射宽波束检测空化，得到空化信号原始射频数据，接着等待一个足够长的时间使得媒介空化核分布恢复到初始状态，通过三维机械扫描装置5沿垂直于阵列换能器6的方向移动一个单元位置，单位位置大小可通过三维机械扫描装置5进行设置，相同参数的空化能量源装置1再次激励产生空化，同步阵列换能器6发射宽波束采集空化射频数据，重复以上过程，得到相同参数的空化能量源条件下，沿垂直于阵列放置方向不同单元位置的一系列二维空化分布图像原始射频数据；

⑦采用最小方差自适应波束合成算法对序列二维空化原始射频数据进行波束合成，得到沿垂直于阵列放置方向不同单元位置的序列二维空化射频数据，随后包络检波、对数压缩以及坐标变换得到序列二维空化分布图像；

⑧采用Nakagami参量成像算法对序列二维空化射频数据进行定量成像，提取空化密度参量，得到沿垂直于阵列放置方向不同单元位置的序列二维空化密度图像；

⑨采用三维重建算法，对沿垂直于阵列放置方向不同单元位置的序列二维空化分布图像和二维空化密度图像进行三维重建，得到空化能量源作用下的三维空化分布图像和三维空化密度图像；

⑩改变产生空化的能量源装置的参数以及媒介，重复步骤⑤-⑨，得到空化以及空化密度定量三维序列图像。

实施例1.以介质自来水为例,高强度聚焦超声脉冲长度10us-200ms，电功率100W，对该实验条件下产生的空化泡群进行空化以及空化密度定量三维成像：

(1)采用图1的实验装置，水槽中注入自来水，高强度聚焦超声换能器固定于水槽侧壁，采用Sonix-Touch全数字化超声设备发射宽波束探测空化微泡，阵列换能器沿着高强度聚焦超声传播方向放置，且使得其中心位置对应于高强度聚焦超声换能器的焦域，通过双通道任意波形发生器AWG420实现两者之间的同步；

(2)通过波形发生器设置高强度聚焦超声(HIFU)脉冲长度为10us-200ms，功率放大器设置电功率为100W，波形编辑设置触发功率放大器的第一通道和触发Sonix-Touch全数字化超声设备的第二通道的时间延迟为1ms；

(3)采用阵列Plane-by-plane宽波束检测空化，通过三维机械扫描装置Multiscan5800沿垂直于阵列放置方向移动，移动单元设置为0.5mm，得到不同单元位置的序列二维空化原始射频数据；

(4)采用最小方差自适应波束合成算法对序列二维空化原始射频数据进行波束合成，得到不同单元位置的序列二维空化射频数据；

(5)针对序列二维空化射频数据，一方面采用包络检波、对数压缩以及坐标变换得到序列二维空化分布图像，一方面采用Nakagami参量成像算法得到序列二维空化密度分布图像；

(6)采用三维重建算法，分别对HIFU脉冲长度为20ms的序列二维空化分布图像和空化密度分布图像进行三维重建，得到自来水中三维空化分布图像和三维空化密度图像。

当阵列换能器位置为0mm时，HIFU脉冲长度为10us-200ms的空化时空和空化密度时空分布结果参见图6，根据分析，空化微泡出现在焦域，这是因为焦域声压最高，但泡群中心位置并不是严格位于0mm焦点位置，而是偏向焦后区域，这是因为声场产生的辐射力会对空化微泡产生作用力，使其沿声传播的方向发生移动，当脉冲长度为5ms时空化泡群逐渐在焦后形成枝状结构，并逐渐增大，当脉冲长度为200ms时，声辐射力太大破坏空化泡群枝状结构。图7给出脉冲长度为20ms时，阵列换能器位置从-5到6mm的各个空化泡群断面分布。图7给出脉冲长度为20ms时，阵列换能器位置从-5到6mm的空间序列空化分布。图8给出了脉冲长度为20ms时，三维空化图像(a)和三维空化微泡密度图像(b)。

实施例2.脉动流条件下，高强度聚焦超声脉冲长度20us，电功率100W，脉动频率80次/秒，流速为5cm/s，管道直径为4mm，对该实验条件下产生的空化泡群进行空化以及空化密度定量三维成像：

(1)制备琼脂明胶仿血管组织体模，将琼脂(3％)、明胶(12％)和除气水(85％)混合加热搅拌至完全溶解，放入除气箱除气30分钟，倒入模具待冷却凝固后慢慢抽出管子，形成管道，管道两端与乳胶管相连，并连接到脉动泵，脉动泵将液体从储存池经管道流入废液池，在管道内形成脉动流。

(2)采用图4的实验装置，水槽中注入自来水，高强度聚焦超声换能器固定于水槽侧壁，采用Sonix-Touch全数字化超声设备发射宽波束探测空化微泡，阵列换能器垂直于高强度聚焦超声传播方向放置，琼脂明胶仿血管组织体模与三维移动装置相连，通过三维移动将高强度聚焦超声聚焦在管道内部，在管道内产生空化。

(3)脉动频率设置为80次/秒，流速为5cm/s，储存液为生理盐水造影微泡溶液，将Sonix-Touch全数字化超声设备设置为脉冲多普勒模式，记录该脉动参数下的多普勒谱，分析其脉动周期和规律，以此编辑波形发生器AWG420的波形，其CH1、CH2和Marker通道分别输入给功率放大器、Sonix-Touch和脉动泵，实现三者之间的同步。

(4)固定阵列换能器位置，采用阵列Plane-by-plane宽波束检测空化，改变脉动泵与功率放大器之间的时间延迟，分别设置为0-12ms，单元时间为1ms，每次停止2s使得液体空化核分布恢复初始状态，得脉动周期内各时刻的序列二维空化分布图像原始射频数据。

(5)通过三维机械扫描装置Multiscan5800沿垂直于阵列放置方向移动，移动单元距离设置为1mm，得到沿着得到不同单元位置的序列二维空化原始射频数。

(6)采用最小方差自适应波束合成算法对脉动周期内各时刻的序列二维空化分布图像原始射频数据和序列二维空化原始射频数据进行波束合成，得到脉动周期内各时刻和不同单元位置的序列二维空化射频数据。

(7)针对序列二维空化射频数据，一方面采用包络检波、对数压缩以及坐标变换得到脉动周期内各时刻和不同单元位置的序列二维空化分布图像，一方面采用Nakagami参量成像算法得到脉动周期内各时刻和不同单元位置的序列二维空化密度分布图像；

(8)采用三维重建算法，分别对序列二维空化分布图像和空化密度分布图像进行三维重建，得到高强度聚焦超声脉冲长度20us，电功率100W时脉动流下管道内三维空化分布图像和三维空化密度图像。

图9给实施例2条件下三维空化图像(a)和三维空化微泡密度图像(b)，可见空化微泡充满了仿血管组织体模管道，并沿着管道两端延伸。

使用宽波束检测空化后，阵列换能器移动一个单元位置，待空化核分布恢复，再在同样的空化能量激励下，再次宽波束检测空化，得到不同单元位置的空间系列二维空化原始射频数据，通过改变空化源能量作用时间、能量源装置激励与阵列换能器发射宽波束之间的时间延迟以及脉动泵与能量源装置激励之间的时间延迟，可得到时间序列二维空化原始射频数据，然后结合宽波束最小方差自适应波束合成、Nakagami参量成像、三维重建算法得到微秒分辨空化三维时空分布图像和空化微泡密度的定量图像，该方法结合阵列Plane-by-plane宽波束检测空化，宽波束最小方差自适应波束合成以及Nakagami参量算法成像，不仅具有空间高分辨、高信噪比的成像特点，而且可实现稳态自由场和脉动流条件下微秒分辨空化三维时空分布成像和空化微泡密度的定量成像，有潜力发展成为类似于声场测量的空化成像标准方法。

Claims (5)

1.微秒分辨空化时空分布的三维空化定量成像方法，其特征在于：包括以下步骤：采用宽波束检测空化得到二维空化原始射频数据，每次宽波束检测空化后，将用于进行宽波束检测的阵列换能器沿垂直于该阵列换能器放置方向移动一个单元位置，待空化核分布恢复到初始状态，再在同样的空化能量激励下，再次采用宽波束检测对应时间的空化，以此逐步得到对应该阵列换能器不同放置位置的一系列二维空化原始射频数据；

然后采用宽波束最小方差自适应波束合成算法对所述二维空化原始射频数据进行波束合成，得到沿垂直于所述阵列换能器放置方向不同放置位置的一系列二维空化射频数据，根据所述二维空化射频数据获得一系列二维空化分布图像，利用三维重建算法将所述二维空化分布图像进行三维显示，得到三维空化图像；

或者，

然后采用宽波束最小方差自适应波束合成算法对所述二维空化原始射频数据进行波束合成，得到沿垂直于所述阵列换能器放置方向不同放置位置的一系列二维空化射频数据，对所述二维空化射频数据进行Nakagami参量提取后得到一系列二维空化密度定量图像，利用三维重建算法将所述二维空化密度定量图像进行三维显示，得到空化微泡密度的定量三维图像。

2.根据权利要求1所述微秒分辨空化时空分布的三维空化定量成像方法，其特征在于：空化的产生装置包括能量源装置和同步信号发生器；空化的检测装置包括可编程的发射宽波束的阵列换能器和并行通道数据采集及存储单元；同步信号发生器产生同步信号分别控制所述能量源装置和阵列换能器，能量源装置产生能量激励空化的产生，且能量源装置产生的能量连续可调，所述阵列换能器发射宽波束对空化进行检测，得到的空化回波信号由并行通道数据采集及存储单元采集存储。

3.根据权利要求1所述微秒分辨空化时空分布的三维空化定量成像方法，其特征在于：所述Nakagami参量提取包括以下步骤：

1)将二维空化射频数据rf进行去噪处理：

a)选取一定尺寸大小的背景信号区域，计算该区域平均能量P；

b)对rf分别叠加平均能量为P的随机高斯白噪声n₁和n₂，得到S₁和S₂；

c)计算S₁和S₂的相关系数，并给定阈值Th，然后对相关系数进行阈值化处理后与rf加权，得到rf_denoise：

${\mathrm{rf}}_{denoise}=\left\{\begin{array}{cc}rf,& corrcoef(S_1,S_2)>Th\\ 0,& corrcoef(S_1,S_2)\≤Th\end{array}\right.---\left(1\right)$

其中corrcoef(S₁,S₂)为S₁和S₂的相关系数；

2)将rf_denoise进行Hilbert解调，得到包络信号R；

3)计算Nakagami参量：

$m=\frac{{\[E\left(R^2\right)\]}^2}{E{\[R^2-E\left(R^2\right)\]}^2}---\left(3\right)$

其中E(·)代表统计均值。

4.根据权利要求1所述微秒分辨空化时空分布的三维空化定量成像方法，其特征在于：在所述阵列换能器的同一个放置位置，根据脉动流多普勒谱记录的脉动周期和脉动规律，通过改变空化能量激励与脉动起始之间的时间延迟，得到脉动周期内不同时刻对应该放置位置的二维空化原始射频数据，在所述阵列换能器移动一个单元位置后，得到脉动周期内不同时刻对应新的放置位置的二维空化原始射频数据，最终得到脉动周期内任一时刻对应所述陈列换能器不同放置位置的一系列二维空化原始射频数据。

5.根据权利要求1所述微秒分辨空化时空分布的三维空化定量成像方法，其特征在于：改变空化能量激励时间，得到随时间演化的三维空化图像和空化微泡密度的定量三维图像；改变空化能量大小，得到随能量变化的三维空化图像和空化微泡密度的定量三维图像；改变空化能量激励与所述阵列换能器发射宽波束之间的时间延迟，得到随时间消散的三维空化图像和空化微泡密度的定量三维图像；改变空化能量作用的媒介，得到不同媒介下三维空化图像和空化微泡密度的定量三维图像。