CN111415408B - 一种超声空化的微秒级多尺度时空成像及特征图谱计算方法与系统 - Google Patents

Abstract

本发明公开了一种超声空化的微秒级多尺度时空成像及特征图谱计算方法与系统：沿X和Y方向规划感兴趣成像位置；对二维超声面阵换能器检测的时空三维空化信号进行延时处理，沿阵元方向叠加后得到空化背向散射信号，对空化背向散射信号进行小波包分解得到空化背向散射有效尺度信号，对其进行希尔伯特变换后计算空化瞬时强度，遍历所有感兴趣成像位置得到X和Y方向的微秒级空化时空成像结果；在此基础上计算空化特征指标沿X和Y方向的分布，遍历所有有效尺度得到X和Y方向的空化特征图谱。本发明可实现微秒时间分辨下超声空化的多尺度时空成像及特征图谱的计算，适用于超声空化的瞬态物理机制研究及精细调控。

Description

一种超声空化的微秒级多尺度时空成像及特征图谱计算方法 与系统

技术领域

本发明属于超声检测与超声成像技术领域，具体涉及一种超声空化的微秒级多尺度时空成像及特征图谱计算方法与系统。

背景技术

聚焦超声波作用于生物组织中会激发出空化效应，即从微泡成核、生长到最后坍塌的物理过程。低强度下空化微泡的振动可使得血脑屏障可逆地开放，有助于药物传输；高强度下空化的强烈机械作用可直接粉碎细胞，从而实现目标组织的均匀毁损，有助于肿瘤组织的精准切除。然而空化效应本身的随机性可能会造成意外的组织损伤，这就需要对空化活动进行准确的实时影像监控。在聚焦超声尤其是脉冲聚焦超声作用下，单个脉冲内空化的动力学过程以及相邻脉冲之间的交互作用机制是脉冲聚焦超声治疗方案优化的重要物理基础；而由于聚焦超声频率本身为兆赫兹级，因此微秒级时间分辨率下空化随时间和空间变化的瞬态物理过程的监控至关重要。

目前，最为常见的空化声学监控手段是一种使超声换能器工作在不发射只接收模式的被动检测方法，该方法由于不受到聚焦超声声场的干扰，可实现超声空化的实时监控。根据所使用的超声换能器类型又可分为单阵元被动检测和多阵元被动成像，其中单阵元被动检测可对空化强度的时间变化进行量化，但其不能提供空化的空间分布；在单阵元被动检测的基础上发展而来的多阵元被动成像通过对各阵元同时接收到的空化信号进行波束合成处理，从而可对空化的空间位置及空间分布进行定征。然而，多阵元被动成像主要通过对聚焦超声辐照一段时间内的空化能量进行累计，从而得到这段时间内空化活动的总体空间位置，但无法对该段时间内空化的时空瞬态物理过程进行成像和定征。

中国发明专利ZL201410834392.1提供了一种微秒分辨空化时空分布的三维空化定量成像方法，该方法中超声成像换能器单次发射宽波束，待介质中空化核分布恢复之后移动超声成像换能器再次发射宽波束，从而得到空化微泡(空化成核后演变而来)的空间分布。但该方法在实时性、检测灵敏度、可操作性以及介质适用性方面存在缺陷：(1)为了避开聚焦超声信号的干扰，超声成像换能器必须在聚焦超声作用停止一段时间(例如，1ms)之后发射宽波束，此时检测到的是已经形成的空化微泡，而并非聚焦超声作用过程中产生的空化活动，即不能实现空化的实时监控；(2)超声成像换能器发射宽波束后，在检测到空化微泡的同时也会检测到生物组织的信号，当空化微泡信号比较微弱时，该方法的检测灵敏度会降低；(3)该专利中提到利用Plane-by-plane宽波束检测空化以实现微秒时间分辨率；然而，该专利中超声成像换能器每次发射宽波束后需要等待一段足够长的时间(例如，2秒)以恢复介质中的空化核分布，而实际聚焦超声作用过程中介质中的空化核分布瞬息万变且可能难以恢复到初始状态，因此该方法并不能实现真正意义上的瞬态监控；(4)该方法通过机械扫描装置带动一维超声成像换能器来实现空化微泡的空间定征，可操作性较差，在实际应用中具有局限性；(5)该方法要求介质中的空化核分布可通过一定方式恢复到初始状态，适用于液体组织或充满液体的腔室但不适用于其他组织。

此外，空化信号中的次谐波、谐波、超谐波以及宽带噪声分别表征了不同的空化类型(稳态空化或惯性空化)，不同频率成分(即不同尺度)下的空化活动所发挥的作用不同，因此，多尺度下空化的时空成像也非常重要。

发明内容

本发明的目的在于提供一种超声空化的微秒级多尺度时空成像及特征图谱计算方法与系统。

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

一种超声空化的微秒级多尺度时空成像及特征图谱计算方法，包括以下步骤：

1)根据空化源位置，沿X方向规划若干个感兴趣成像位置，该成像位置的Y坐标和Z坐标分别与空化源位置的Y坐标和Z坐标相同；沿Y方向规划若干个感兴趣成像位置，该成像位置的X坐标和Z坐标分别与空化源位置的X坐标和Z坐标相同；

2)根据步骤1)规划的某一感兴趣成像位置到二维超声面阵换能器阵元位置的超声波传播时间，对二维超声面阵换能器检测的时空三维空化信号进行延时处理，将经过延时的空化信号沿阵元方向叠加，得到空化背向散射信号；选择小波包分解的小波基、小波包分解层数并确定小波包分解的有效尺度，对所得空化背向散射信号进行小波包分解，得到不同有效尺度下的空化背向散射有效尺度信号，根据所得任意有效尺度下的空化背向散射有效尺度信号及该有效尺度信号的希尔伯特变换结果，计算得到所述感兴趣成像位置在所述有效尺度下的空化瞬时强度；

3)针对步骤1)中沿X方向和Y方向规划的各感兴趣成像位置，通过重复步骤2)分别计算沿X方向和Y方向规划的各感兴趣成像位置在所述有效尺度下的空化瞬时强度，得到所述有效尺度下X方向和Y方向的微秒级空化时空成像结果；

4)根据步骤3)计算不同有效尺度下的X方向和Y方向的微秒级空化时空成像结果；根据X方向和Y方向的不同有效尺度下的微秒级空化时空成像结果，计算X方向和Y方向的对应空化特征图谱。

优选的，所述时空三维空化信号是通过二维超声面阵换能器被动接收(即二维超声面阵换能器不主动地向外部发射探测脉冲)的包含两个空间维度(X和Y)及一个时间维度的空化采样信号。

优选的，所述步骤2)中，空化背向散射信号表示为：

其中，CB(x,y,z,t)为空化背向散射信号，Nx和Ny分别为二维超声面阵换能器X方向和Y方向的阵元数目，ch_i,j(t)为所述时空三维空化信号，t为二维超声面阵换能器阵元接收时空三维空化信号的时刻，τ_i,j(x,y,z)为某一感兴趣成像位置(x,y,z)到二维超声面阵换能器阵元位置(xe_i,ye_j,⁰)的超声波传播时间。

优选的，所述步骤2)中，小波包分解的小波基选自小波消失矩阶数L较大的dbL小波或symL小波，L为6～10。

优选的，所述步骤2)中，根据二维超声面阵换能器的接收带宽f₁～f₂，确定小波包分解的有效尺度：

其中，f₁和f₂分别为二维超声面阵换能器的接收带宽的下限和上限，f_Nyquist＝f_Sample/2为Nyquist频率，f_Sample为采样频率，ceil(·)表示向上取整，floor(·)表示向下取整，p为小波包分解层数。

优选的，所述步骤2)中，小波包分解层数大于3。

优选的，所述步骤2)中，空化背向散射信号的小波包分解表示为：

其中，CB(x,y,z,t)为空化背向散射信号，CBE_k(x,y,z,t)和CBI_m(x,y,z,t)分别为空化背向散射有效尺度信号和空化背向散射无效尺度信号，k＝1,2,...,K，k从小到大依次表征二维超声面阵换能器的接收带宽内从低频到高频的不同的有效尺度；m＝1,2,...,2^p-K，表征无效尺度。

优选的，所述步骤2)中，空化瞬时强度按照以下公式计算：

其中，CIA_k(x,y,z,t)为某一感兴趣成像位置(x,y,z)在第k个有效尺度下的空化瞬时强度，CBE_k(x,y,z,t)为第k个空化背向散射有效尺度信号，CBEH_k(x,y,z,t)为第k个空化背向散射有效尺度信号所对应的希尔伯特变换结果。

优选的，所述步骤4)中，X方向和Y方向的空化特征图谱的计算具体包括以下步骤：

4.1)根据步骤3)所得第k个有效尺度下X方向的微秒级空化时空成像结果分别计算对应有效尺度下空化时间指标、空化能量指标和空化峰态指标沿X方向的分布，k＝1,2,...,K；根据K个有效尺度下对应的空化时间指标、空化能量指标和空化峰态指标沿X方向的分布，得到X方向的空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱；

4.2)根据步骤3)所得第k个有效尺度下Y方向的微秒级空化时空成像结果分别计算对应有效尺度下空化时间指标、空化能量指标和空化峰态指标沿Y方向的分布；根据K个有效尺度下对应的空化时间指标、空化能量指标和空化峰态指标沿Y方向的分布，得到Y方向的空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱。

优选的，所述步骤4.1)中，第k个有效尺度下空化时间指标沿X方向的分布为沿时间轴计算的步骤3)所得对应有效尺度下X方向的微秒级空化时空成像结果经二值化处理后所得的结果的和：

所述步骤4.2)中，第k个有效尺度下空化时间指标沿Y方向的分布为沿时间轴计算的步骤3)所得对应有效尺度下Y方向的微秒级空化时空成像结果经二值化处理后所得的结果的和：

其中，CTX_k(x,y_c,z_c)和CTY_k(x_c,y,z_c)分别为第k个有效尺度下空化时间指标沿X方向和Y方向的分布，

和

分别为第k个有效尺度下X方向和Y方向的微秒级空化时空成像结果按照对应阈值ΔX_k和ΔY_k进行二值化处理后所得的结果，TSX_k(x,y_c,z_c,t)和TSY_k(x_c,y,z_c,t)分别为第k个有效尺度下X方向和Y方向的微秒级空化时空成像结果。

优选的，所述阈值ΔX_k和ΔY_k分别设置为第k个有效尺度下X方向和Y方向的微秒级空化时空成像结果的最大值的一半。

优选的，所述步骤4.1)中，第k个有效尺度下空化能量指标沿X方向的分布为沿时间轴计算的步骤3)所得对应有效尺度下X方向的微秒级空化时空成像结果的均方值：

所述步骤4.2)中，第k个有效尺度下空化能量指标沿Y方向的分布为沿时间轴计算的步骤3)所得对应有效尺度下Y方向的微秒级空化时空成像结果的均方值：

其中，CEX_k(x,y_c,z_c)和CEY_k(x_c,y,z_c)分别为第k个有效尺度下空化能量指标沿X方向和Y方向的分布，T为所述时空三维空化信号的信号长度。

优选的，所述步骤4.1)中，第k个有效尺度下空化峰态指标沿X方向的分布为沿时间轴计算的步骤3)所得对应有效尺度下X方向的微秒级空化时空成像结果的最大值与均方根值的比值的平方：

所述步骤4.2)中，第k个有效尺度下空化峰态指标沿Y方向的分布为沿时间轴计算的步骤3)所得对应有效尺度下Y方向的微秒级空化时空成像结果的最大值与均方根值的比值的平方：

其中，CPX_k(x,y_c,z_c)和CPY_k(x_c,y,z_c)分别为第k个有效尺度下空化峰态指标沿X方向和Y方向的分布。

一种超声空化的微秒级多尺度时空成像及特征图谱计算系统，该系统包括二维超声面阵换能器、感兴趣成像位置规划模块、微秒级多尺度空化时空成像模块以及空化特征图谱计算模块；

所述感兴趣成像位置规划模块用于根据空化源位置规划微秒级多尺度空化时空成像所需的沿X方向和Y方向的感兴趣成像位置；

所述微秒级多尺度空化时空成像模块用于实现不同尺度下的微秒级空化时空成像，该模块包括延时叠加处理子模块、小波包分解子模块、空化瞬时强度计算子模块以及微秒级空化时空成像子模块；所述延时叠加处理子模块用于根据某一感兴趣成像位置到二维超声面阵换能器阵元位置的超声波传播时间，对二维超声面阵换能器检测的时空三维空化信号进行延时处理，及将经过延时的空化信号沿阵元方向叠加，从而得到空化背向散射信号；所述小波包分解子模块用于根据小波包分解的小波基、小波包分解层数以及小波包分解的有效尺度，对所述空化背向散射信号进行小波包分解，以生成不同有效尺度下的空化背向散射有效尺度信号；所述空化瞬时强度计算子模块用于根据任意有效尺度下的空化背向散射有效尺度信号及该有效尺度信号的希尔伯特变换结果计算所述感兴趣成像位置在所述有效尺度下的空化瞬时强度；所述微秒级空化时空成像子模块用于参照所述感兴趣成像位置规划模块规划的X方向和Y方向的感兴趣成像位置，并利用所述空化瞬时强度计算子模块，获得所述有效尺度下X方向和Y方向的微秒级空化时空成像结果；

所述空化特征图谱计算模块用于根据X方向和Y方向的微秒级空化时空成像结果计算不同指标(空化时间指标、空化能量指标、空化峰态指标)下的空化特征图谱，该模块包括空化时间特征图谱计算子模块、空化能量特征图谱计算子模块以及空化峰态特征图谱计算子模块；所述空化时间特征图谱计算子模块用于计算不同有效尺度下空化时间指标沿X方向和Y方向的分布，从而生成X方向和Y方向的空化时间特征图谱，其中某一有效尺度下空化时间指标沿X方向和Y方向的分布为沿时间轴计算的对应有效尺度下X方向和Y方向的微秒级空化时空成像结果经二值化处理后所得的结果的和；所述空化能量特征图谱计算子模块用于计算不同有效尺度下空化能量指标沿X方向和Y方向的分布，从而生成X方向和Y方向的空化能量特征图谱，其中某一有效尺度下空化能量指标沿X方向和Y方向的分布为沿时间轴计算的对应有效尺度下X方向和Y方向的微秒级空化时空成像结果的均方值；所述空化峰态特征图谱计算子模块用于计算不同有效尺度下空化峰态指标沿X方向和Y方向的分布，从而生成X方向和Y方向的空化峰态特征图谱，其中某一有效尺度下空化峰态指标沿X方向和Y方向的分布为沿时间轴计算的对应有效尺度下X方向和Y方向的微秒级空化时空成像结果的最大值与均方根值的比值的平方。

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

本发明提出的微秒级多尺度空化时空成像及特征图谱计算方法，将小波包分解、希尔伯特变换引入到对时间分辨率与时空三维空化信号采样频率(一般为10～60MHz)相对应的空化背向散射信号的处理中，从而可以获得多尺度下微秒级时间分辨的空化时空成像结果；根据该空化时空成像结果可分别得到空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱。本发明可以对聚焦超声辐照过程中产生的空化活动的瞬态物理过程进行有效定征，并为超声空化的物理机制研究及精细调控提供了有力手段。

进一步地，本发明中的二维超声面阵换能器可一次性得到包含两个空间维度及一个时间维度的时空三维空化信号；且本发明中时空三维空化信号是通过使二维超声面阵换能器被动接收得到的，一方面可以实现空化的实时检测，另一方面可以提高对微弱空化信号的检测灵敏度。

进一步地，本发明中选择与时空三维空化信号较为相似的小波基对空化背向散射信号进行小波包分解，有助于获得多个尺度下的微秒级空化时空成像结果。

进一步地，本发明中以二维超声面阵换能器的接收带宽为依据来确定小波包分解的有效尺度，可从空化背向散射信号的小波包分解结果中遴选出与接收带宽相对应的空化背向散射有效尺度信号，并去除掉超出接收带宽的空化背向散射无效尺度信号，从而可以降低后续空化瞬时强度计算以及空化特征图谱计算所需的计算量。

进一步地，本发明根据某一有效尺度下X方向和Y方向的微秒级空化时空成像结果分别计算得到某一有效尺度下空化时间指标、空化能量指标和空化峰态指标沿X方向和Y方向的分布；然后遍历所有有效尺度，从而得到X方向和Y方向的空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱，可直观地反映在不同空间位置(X方向和Y方向)和不同尺度下空化活动的发生时间、时间累积能量以及冲击特征。

附图说明

图1为本发明实施例中沿X方向和Y方向的感兴趣成像位置的规划示意图。

图2为本发明实施例中X方向和Y方向的微秒级多尺度空化时空成像的流程图。

图3为本发明实施例中微秒级多尺度空化时空成像中小波包分解所用到的小波基的小波函数图。

图4为本发明实施例中微秒级多尺度空化时空成像中小波包分解的示意图。

图5为本发明实施例中X方向(a)和Y方向(b)的微秒级多尺度空化时空成像结果。

图6为本发明实施例中X方向和Y方向的空化时间特征图谱(a和d)、空化能量特征图谱(b和e)和空化峰态特征图谱(c和f)的计算流程图。

图7为本发明实施例中计算得到的X方向和Y方向的空化时间特征图谱(a和d)、空化能量特征图谱(b和e)和空化峰态特征图谱(c和f)。

具体实施方式

下面结合附图和实施例对本发明作进一步的详细说明。

本发明针对聚焦超声空化瞬态物理过程缺乏微秒级多尺度空化时空成像监控的问题，首先沿二维超声面阵换能器的X方向和Y方向规划感兴趣成像位置；然后计算某一感兴趣成像位置到二维超声面阵换能器阵元位置的超声波传播时间，在此基础上对二维超声面阵换能器检测的时空三维空化信号进行延时处理并沿阵元方向叠加后得到空化背向散射信号；接着根据合适的小波包分解的小波基、小波包分解层数以及小波包分解的有效尺度，对空化背向散射信号进行小波包分解得到空化背向散射有效尺度信号；之后根据空化背向散射有效尺度信号及其希尔伯特变换结果计算得到不同有效尺度下的空化瞬时强度；最后遍历X方向和Y方向规划的所有感兴趣成像位置，分别得到不同有效尺度下X方向和Y方向的微秒级空化时空成像结果，从而实现微秒级时间分辨率下多尺度的空化时空成像。在此基础上计算各有效尺度下空化特征指标(时间、能量、峰态)沿X和Y方向的分布，得到X和Y方向的空化特征图谱(时间、能量、峰态)。本发明的具体步骤和结果举例说明如下。

(1)利用二维超声面阵换能器(例如，X方向的阵元数目Nx为32，Y方向的阵元数目Ny为32)被动地接收(即不主动地向外部发射探测脉冲)时空三维空化信号(三个维度分别为X方向的阵元、Y方向的阵元和时间)，接收到的信号记为ch_i,j(t)，其中i＝1,2,...,Nx，j＝1,2,...,Ny，t为时空三维空化信号接收时刻，信号长度为T；

(2)参见图1，根据空化源所在位置(x_c,y_c,z_c)，沿X方向规划NPX(例如，NPX＝100)个感兴趣成像位置，其X坐标为x₁,...,x_c,...,x_NPX(例如，间隔设置为0.048mm)，Y坐标固定为y_c，Z坐标固定为z_c；同样地，沿Y方向规划NPY(例如，NPY＝100)个感兴趣成像位置，其X坐标固定为x_c，Y坐标为y₁,...,y_c,...,y_NPY(例如，间隔设置为0.048mm)，Z坐标固定为z_c；

参见图2，步骤(3.1)～(3.10)为微秒级多尺度空化时空成像的具体流程；

(3.1)计算从步骤(2)规划的某一感兴趣成像位置(x,y,z)到二维超声面阵换能器阵元位置(xe_i,ye_j,0)的超声波传播时间τ_i,j(x,y,z)：

其中，c为超声波的传播速度；

(3.2)根据步骤(3.1)所得超声波传播时间τ_i,j(x,y,z)对步骤(1)所得时空三维空化信号ch_i,j(t)进行延时处理，然后沿阵元方向叠加后得到空化背向散射信号CB(x,y,z,t)：

(3.3)选择小波包分解的小波基(例如选择db10小波，小波函数参见图3)；

(3.4)选择小波包分解层数p(例如，p＝5)，在第p层共有2^p个尺度，其中包含有K个有效尺度和2^p-K个无效尺度；

(3.5)根据步骤(1)所述二维超声面阵换能器的接收带宽f₁～f₂，确定小波包分解的有效尺度(数目为K)为：

其中，f₁和f₂分别为二维超声面阵换能器的接收带宽的下限和上限(例如，4MHz和11MHz)，f_Nyquist＝f_Sample/2，为Nyquist频率，f_Sample为采样频率(例如，50MHz)，ceil(·)表示向上取整，floor(·)表示向下取整；

(3.6)参见图4所示的小波包分解示意图(实线框表示有效尺度，虚线框表示无效尺度)，根据步骤(3.3)选择的小波包分解的小波基和步骤(3.4)选择的小波包分解层数对步骤(3.2)所述空化背向散射信号CB(x,y,z,t)进行小波包分解；根据步骤(3.5)所得K个有效尺度分别得到K个空化背向散射有效尺度信号CBE_k(x,y,z,t)和2^p-K个空化背向散射无效尺度信号CBI_m(x,y,z,t)：

其中，k＝1,2,...,K，k从小到大依次表征二维超声面阵换能器的接收带宽内从低频到高频的不同的有效尺度；m＝1,2,...,2^p-K，表征无效尺度，空化背向散射无效尺度信号CBI_m(x,y,z,t)不进行后续步骤的计算；

(3.7)对步骤(3.6)得到的第k个空化背向散射有效尺度信号CBE_k(x,y,z,t)进行希尔伯特变换：

其中，

表示信号的卷积运算；

(3.8)根据步骤(3.6)及(3.7)所得的第k个空化背向散射有效尺度信号CBE_k(x,y,z,t)及其对应的希尔伯特变换结果CBEH_k(x,y,z,t)计算感兴趣成像位置(x,y,z)在第k个有效尺度下的空化瞬时强度CIA_k(x,y,z,t)：

(3.9)针对步骤(2)中X方向规划的NPX个感兴趣成像位置重复步骤(3.1)～(3.8)，得到第k个有效尺度下X方向的微秒级空化时空成像结果TSX_k(x,y_c,z_c,t)(该结果为二维矩阵)：

TSX_k(x,y_c,z_c,t)＝[CIA_k(x₁,y_c,z_c,t)；CIA_k(x₂,y_c,z_c,t)；...；CIA_k(x_NPX,y_c,z_c,t)]

(3.10)针对步骤(2)中Y方向规划的NPY个感兴趣成像位置重复步骤(3.1)～(3.8)，得到第k个有效尺度下Y方向的微秒级空化时空成像结果TSY_k(x_c,y,z_c,t)(该结果为二维矩阵)：

TSY_k(x_c,y,z_c,t)＝[CIA_k(x_c,y₁,z_c,t)；CIA_k(x_c,y₂,z_c,t)；...；CIA_k(x_c,y_NPY,z_c,t)]

参见图5，图5(a)从上向下依次显示了K个有效尺度(K＝9)下X方向所得的微秒级空化时空成像结果，图5(b)从上向下依次显示了K个有效尺度(K＝9)下Y方向所得的微秒级空化时空成像结果，表明成像结果时间分辨率可达微秒级；每个尺度下的结果反映了该尺度下空化出现的时间点和空间位置以及该时间点和该空间位置下空化的强度(图像中像素值越大，空化的强度越大)；不同尺度下的结果反映了不同频率成分下(9个有效尺度对应的频段分别为[3.91MHz,4.69MHz]、[4.69MHz,5.47MHz]、[5.47MHz,6.25MHz]、[6.25MHz,7.03MHz]、[7.03MHz,7.81MHz]、[7.81MHz,8.59MHz]、[8.59MHz,9.38MHz]、[9.38MHz,10.16MHz]和[10.16MHz,10.94MHz])的空化活动的强度的时空演变过程；

参见图6，图6(a)、图6(b)、图6(c)分别显示了X方向的空化时间特征图谱、空化能量特征图谱、空化峰态特征图谱的计算流程，具体计算过程见下述步骤(4.1)～(4.7)；

(4.1)对步骤(3.9)所得第k个有效尺度下X方向的微秒级空化时空成像结果TSX_k(x,y_c,z_c,t)按照阈值ΔX_k进行二值化处理：

其中，阈值ΔX_k为TSX_k(x,y_c,z_c,t)的最大值的一半；

(4.2)沿时间轴计算步骤(4.1)所得第k个有效尺度下X方向的微秒级空化时空成像结果TSX_k(x,y_c,z_c,t)的二值化结果

的和(即空化时间指标)，得到第k个有效尺度下空化时间指标沿X方向的分布CTX_k(x,y_c,z_c)：

(4.3)重复步骤(4.2)，直至获得K个有效尺度下空化时间指标沿X方向的分布，从而得到X方向的空化时间特征图谱CTX(x,y_c,z_c,k)：

CTX(x,y_c,z_c,k)＝[CTX₁(x,y_c,z_c)；CTX₂(x,y_c,z_c)；...；CTX_K(x,y_c,z_c)]

(4.4)沿时间轴计算步骤(3.9)所得第k个有效尺度下X方向的微秒级空化时空成像结果TSX_k(x,y_c,z_c,t)的均方值(即空化能量指标)，得到第k个有效尺度下空化能量指标沿X方向的分布CEX_k(x,y_c,z_c)：

(4.5)重复步骤(4.4)，直至获得K个有效尺度下空化能量指标沿X方向的分布，从而得到X方向的空化能量特征图谱CEX(x,y_c,z_c,k)：

CEX(x,y_c,z_c,k)＝[CEX₁(x,y_c,z_c)；CEX₂(x,y_c,z_c)；...；CEX_K(x,y_c,z_c)]

(4.6)沿时间轴计算步骤(3.9)所得第k个有效尺度下X方向的微秒级空化时空成像结果TSX_k(x,y_c,z_c,t)的最大值与均方根值的比值的平方(即空化峰态指标)，得到第k个有效尺度下空化峰态指标沿X方向的分布CPX_k(x,y_c,z_c)：

(4.7)重复步骤(4.6)，直至获得K个有效尺度下空化峰态指标沿X方向的分布，从而得到X方向的空化峰态特征图谱CPX(x,y_c,z_c,k)：

CPX(x,y_c,z_c,k)＝[CPX₁(x,y_c,z_c)；CPX₂(x,y_c,z_c)；...；CPX_K(x,y_c,z_c)]

参见图6，图6(d)、图6(e)、图6(f)分别显示了Y方向的空化时间特征图谱、空化能量特征图谱、空化峰态特征图谱的计算流程，具体计算过程见下述步骤(5.1)～(5.7)；

(5.1)对步骤(3.10)所得第k个有效尺度下Y方向的微秒级空化时空成像结果TSY_k(x_c,y,z_c,t)按照阈值ΔY_k进行二值化处理：

其中，阈值ΔY_k为TSY_k(x_c,y,z_c,t)的最大值的一半；

(5.2)沿时间轴计算步骤(5.1)所得第k个有效尺度下Y方向的微秒级空化时空成像结果TSY_k(x_c,y,z_c,t)的二值化结果

的和(即空化时间指标)，得到第k个有效尺度下空化时间指标沿Y方向的分布CTY_k(x_c,y,z_c)：

(5.3)重复步骤(5.2)，直至获得K个有效尺度下空化时间指标沿Y方向的分布，从而得到Y方向的空化时间特征图谱CTY(x_c,y,z_c,k)：

CTY(x_c,y,z_c,k)＝[CTY₁(x_c,y,z_c)；CTY₂(x_c,y,z_c)；...；CTY_K(x_c,y,z_c)]

(5.4)沿时间轴计算步骤(3.10)所得第k个有效尺度下Y方向的微秒级空化时空成像结果TSY_k(x_c,y,z_c,t)的均方值(即空化能量指标)，得到第k个有效尺度下空化能量指标沿Y方向的分布CEY_k(x_c,y,z_c)：

(5.5)重复步骤(5.4)，直至获得K个有效尺度下空化能量指标沿Y方向的分布，从而得到Y方向的空化能量特征图谱CEY(x_c,y,z_c,k)：

CEY(x_c,y,z_c,k)＝[CEY₁(x_c,y,z_c)；CEY₂(x_c,y,z_c)；...；CEY_K(x_c,y,z_c)]

(5.6)沿时间轴计算步骤(3.10)所得第k个有效尺度下Y方向的微秒级空化时空成像结果TSY_k(x_c,y,z_c,t)的最大值与均方根值的比值的平方(即空化峰态指标)，得到第k个有效尺度下空化峰态指标沿Y方向的分布CPY_k(x_c,y,z_c)：

(5.7)重复步骤(5.6)，直至获得K个有效尺度下空化峰态指标沿Y方向的分布，从而得到Y方向的空化峰态特征图谱CPY(x_c,y,z_c,k)：

CPY(x_c,y,z_c,k)＝[CPY₁(x_c,y,z_c)；CPY₂(x_c,y,z_c)；...；CPY_K(x_c,y,z_c)]

参见图7，图7(a)显示了X方向的空化时间特征图谱，图7(d)显示了Y方向的空化时间特征图谱，其反映了在不同空间位置(X方向和Y方向)处不同尺度下的空化活动的发生时间，从图中可直观地观察到X方向0.14mm处第三个有效尺度(对应频段为[5.47MHz,6.25MHz])下的空化活动以及Y方向-0.14mm处第一个有效尺度(对应频段为[3.91MHz,4.69MHz])下的空化活动在发生时间上占主导地位；图7(b)显示了X方向的空化能量特征图谱，图7(e)显示了Y方向的空化能量特征图谱，其反映了不同尺度下的空化活动的时间累积能量在X方向和Y方向的分布情况，从图中可直观地观察到X方向0.19mm处第四个有效尺度(对应频段为[6.25MHz,7.03MHz])下的空化活动以及Y方向0.10mm处第四个有效尺度(对应频段为[6.25MHz,7.03MHz])下的空化活动在声能量上占主导地位；图7(c)显示了X方向的空化峰态特征图谱，图7(f)显示了Y方向的空化峰态特征图谱，其反映了在不同空间位置(X方向和Y方向)处不同尺度下的空化活动的冲击特征，从图中可直观地观察到X方向1.63mm处第八个有效尺度(对应频段为[9.38MHz,10.16MHz])下的空化活动以及Y方向-1.30mm处第八个有效尺度(对应频段为[9.38MHz,10.16MHz])下的空化活动的变化更加剧烈。

本发明具有以下优点：

(1)传统的多阵元被动成像只能获得聚焦超声辐照一段时间内空化活动的总体空间位置，无法获得该段时间内的微秒级空化时空成像。而本发明提出的超声空化的微秒级多尺度时空成像方法通过对由三维空化信号延时叠加而来的空化背向散射信号进行小波包分解并对所得空化背向散射有效尺度信号进行希尔伯特变换，从而得到沿X方向和Y方向规划的感兴趣成像位置处的空化瞬时强度，可在微秒级时间分辨率下同时观察多个尺度下空化活动随时间和空间的演变规律。

(2)本发明提出的空化特征图谱计算方法根据某一有效尺度下微秒级空化时空成像结果计算得到该有效尺度下空化时间指标、空化能量指标和空化峰态指标沿X方向和Y方向的分布，进而分别得到X方向和Y方向的空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱，可用以观察空化发生时间、空化时间累积能量及空化冲击特征在不同空间位置和不同尺度下的分布。

(3)相比现有的时空分布的三维空化定量成像方法(例如，ZL201410834392.1)，本发明中采用被动地接收三维空化信号，可实现空化的实时监控并提高微弱空化信号的检测灵敏度；本发明中空化时空成像的时间分辨率与时空三维空化信号的采样频率相对应，可实现真正意义上的微秒级时间分辨的瞬态监控；本发明中采用二维超声面阵换能器一次性获得三维空化信号，避免了额外扫描装置所带来的可操作性差的不足；本发明所述方法在介质种类上不存在局限性，在实际应用中具有广泛前景。

(4)本发明既可用于连续波聚焦超声辐照过程中的微秒级多尺度空化时空成像，也可用于任意脉冲参数下脉冲聚焦超声辐照过程中的微秒级多尺度空化时空成像，这为聚焦超声治疗过程的影像监控以及超声空化的物理机制研究与精细调控提供了有效手段，同时也为实现精准和高效聚焦超声治疗奠定了基础。

Claims (9)

1.一种超声空化的微秒级多尺度时空成像方法，其特征在于：包括以下步骤：

1)根据空化源位置，沿X方向、Y方向分别规划若干个感兴趣成像位置；

2)根据步骤1)规划的某一感兴趣成像位置到二维超声面阵换能器阵元位置的超声波传播时间，对二维超声面阵换能器检测的时空三维空化信号进行延时处理，将经过延时的空化信号沿阵元方向叠加，得到空化背向散射信号；对所得空化背向散射信号进行小波包分解，得到不同有效尺度下的空化背向散射有效尺度信号，根据所得任意有效尺度下的空化背向散射有效尺度信号及该有效尺度信号的希尔伯特变换结果，计算得到所述感兴趣成像位置在所述有效尺度下的空化瞬时强度；

所述空化瞬时强度按照以下公式计算：

其中，CIA_k(x,y,z,t)为某一感兴趣成像位置(x,y,z)在第k个有效尺度下的空化瞬时强度，CBE_k(x,y,z,t)为第k个空化背向散射有效尺度信号，CBEH_k(x,y,z,t)为第k个空化背向散射有效尺度信号所对应的希尔伯特变换结果；

3)通过重复步骤2)分别计算步骤1)中沿X方向和Y方向规划的各感兴趣成像位置在所述有效尺度下的空化瞬时强度，得到所述有效尺度下X方向和Y方向的微秒级空化时空成像结果。

2.一种基于微秒级多尺度时空成像的超声空化特征图谱计算方法，其特征在于：包括以下步骤：

1)根据空化源位置，沿X方向、Y方向分别规划若干个感兴趣成像位置；

2)根据步骤1)规划的某一感兴趣成像位置到二维超声面阵换能器阵元位置的超声波传播时间，对二维超声面阵换能器检测的时空三维空化信号进行延时处理，将经过延时的空化信号沿阵元方向叠加，得到空化背向散射信号；对所得空化背向散射信号进行小波包分解，得到不同有效尺度下的空化背向散射有效尺度信号，根据所得任意有效尺度下的空化背向散射有效尺度信号及该有效尺度信号的希尔伯特变换结果，计算得到所述感兴趣成像位置在所述有效尺度下的空化瞬时强度；

所述空化瞬时强度按照以下公式计算：

其中，CIA_k(x,y,z,t)为某一感兴趣成像位置(x,y,z)在第k个有效尺度下的空化瞬时强度，CBE_k(x,y,z,t)为第k个空化背向散射有效尺度信号，CBEH_k(x,y,z,t)为第k个空化背向散射有效尺度信号所对应的希尔伯特变换结果；

3)通过重复步骤2)分别计算步骤1)中沿X方向和Y方向规划的各感兴趣成像位置在不同有效尺度下的空化瞬时强度，得到X方向和Y方向的不同有效尺度下的微秒级空化时空成像结果；

4)根据X方向和Y方向的不同有效尺度下的微秒级空化时空成像结果，计算X方向和Y方向的对应特征图谱。

3.根据权利要求1或2所述的方法，其特征在于：所述时空三维空化信号是通过二维超声面阵换能器被动接收的包含两个空间维度及一个时间维度的空化采样信号。

4.根据权利要求1或2所述的方法，其特征在于：所述步骤2)中，空化背向散射信号表示为：

其中，CB(x,y,z,t)为空化背向散射信号，Nx和Ny分别为二维超声面阵换能器X方向和Y方向的阵元数目，ch_i,j(t)为所述时空三维空化信号，t为二维超声面阵换能器阵元接收时空三维空化信号的时刻，τ_i,j(x,y,z)为某一感兴趣成像位置(x,y,z)到二维超声面阵换能器阵元位置(xe_i,ye_j,0)的超声波传播时间。

5.根据权利要求1或2所述的方法，其特征在于：所述步骤2)中，小波包分解的小波基选自小波消失矩阶数L较大的dbL小波或symL小波，L为6～10；小波包分解的有效尺度根据二维超声面阵换能器的接收带宽来确定：

其中，f₁和f₂分别为二维超声面阵换能器的接收带宽的下限和上限，f_Nyquist为Nyquist频率，ceil(·)表示向上取整，floor(·)表示向下取整，p为小波包分解层数。

6.根据权利要求2所述的方法，其特征在于：所述步骤4)具体包括以下步骤：

4.1)根据步骤3)所得第k个有效尺度下X方向的微秒级空化时空成像结果分别计算对应有效尺度下空化时间指标、空化能量指标和空化峰态指标沿X方向的分布，k＝1,2,...,K；根据K个有效尺度下对应的空化时间指标、空化能量指标和空化峰态指标沿X方向的分布，得到X方向的空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱；

4.2)根据步骤3)所得第k个有效尺度下Y方向的微秒级空化时空成像结果分别计算对应有效尺度下空化时间指标、空化能量指标和空化峰态指标沿Y方向的分布，k＝1,2,...,K；根据K个有效尺度下对应的空化时间指标、空化能量指标和空化峰态指标沿Y方向的分布，得到Y方向的空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱。

7.根据权利要求6所述的方法，其特征在于：所述第k个有效尺度下空化时间指标沿X方向的分布为沿时间轴计算的对应有效尺度下X方向的微秒级空化时空成像结果经二值化处理后所得的结果的和；所述第k个有效尺度下空化时间指标沿Y方向的分布为沿时间轴计算的对应有效尺度下Y方向的微秒级空化时空成像结果经二值化处理后所得的结果的和；

所述第k个有效尺度下空化能量指标沿X方向的分布为沿时间轴计算的对应有效尺度下X方向的微秒级空化时空成像结果的均方值；所述第k个有效尺度下空化能量指标沿Y方向的分布为沿时间轴计算的对应有效尺度下Y方向的微秒级空化时空成像结果的均方值；

所述第k个有效尺度下空化峰态指标沿X方向的分布为沿时间轴计算的对应有效尺度下X方向的微秒级空化时空成像结果的最大值与均方根值的比值的平方；所述第k个有效尺度下空化峰态指标沿Y方向的分布为沿时间轴计算的对应有效尺度下Y方向的微秒级空化时空成像结果的最大值与均方根值的比值的平方。

8.根据权利要求7所述的方法，其特征在于：所述二值化处理的阈值分别设置为第k个有效尺度下X方向和Y方向的微秒级空化时空成像结果的最大值的一半。

9.一种超声空化的微秒级多尺度时空成像及特征图谱计算系统，其特征在于：该系统包括二维超声面阵换能器、感兴趣成像位置规划模块、微秒级多尺度空化时空成像模块以及空化特征图谱计算模块；

所述感兴趣成像位置规划模块用于根据空化源位置规划微秒级多尺度空化时空成像所需的X方向和Y方向的感兴趣成像位置；

所述微秒级多尺度空化时空成像模块包括延时叠加处理子模块、小波包分解子模块、空化瞬时强度计算子模块以及微秒级空化时空成像子模块；所述延时叠加处理子模块用于根据某一感兴趣成像位置到二维超声面阵换能器阵元位置的超声波传播时间，对二维超声面阵换能器检测的时空三维空化信号进行延时处理，及将经过延时的空化信号沿阵元方向叠加，从而得到空化背向散射信号；所述小波包分解子模块用于对所述空化背向散射信号进行小波包分解，以生成不同有效尺度下的空化背向散射有效尺度信号；所述空化瞬时强度计算子模块用于根据任意有效尺度下的空化背向散射有效尺度信号及该有效尺度信号的希尔伯特变换结果计算所述感兴趣成像位置在所述有效尺度下的空化瞬时强度；所述微秒级空化时空成像子模块用于参照所述感兴趣成像位置规划模块规划的X方向和Y方向的感兴趣成像位置，并利用所述空化瞬时强度计算子模块，获得沿X方向和Y方向规划的各感兴趣成像位置在所述有效尺度下的空化瞬时强度，从而得到所述有效尺度下X方向和Y方向的微秒级空化时空成像结果；

所述空化特征图谱计算模块用于根据X方向和Y方向的不同有效尺度下的微秒级空化时空成像结果计算X方向和Y方向的空化时间特征图谱、空化能量特征图谱和空化峰态特征图谱；

所述空化瞬时强度按照以下公式计算：

其中，CIA_k(x,y,z,t)为某一感兴趣成像位置(x,y,z)在第k个有效尺度下的空化瞬时强度，CBE_k(x,y,z,t)为第k个空化背向散射有效尺度信号，CBEH_k(x,y,z,t)为第k个空化背向散射有效尺度信号所对应的希尔伯特变换结果。

Images