CN109431536B - 一种聚焦超声空化的实时高分辨时空分布成像方法与系统 - Google Patents
一种聚焦超声空化的实时高分辨时空分布成像方法与系统 Download PDFInfo
- Publication number
- CN109431536B CN109431536B CN201811083655.4A CN201811083655A CN109431536B CN 109431536 B CN109431536 B CN 109431536B CN 201811083655 A CN201811083655 A CN 201811083655A CN 109431536 B CN109431536 B CN 109431536B
- Authority
- CN
- China
- Prior art keywords
- cavitation
- image
- steady
- inertial
- frame
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000009826 distribution Methods 0.000 title claims abstract description 138
- 238000003384 imaging method Methods 0.000 title claims abstract description 118
- 230000002123 temporal effect Effects 0.000 title abstract 5
- 238000000034 method Methods 0.000 claims abstract description 44
- 238000002604 ultrasonography Methods 0.000 claims abstract description 42
- 230000008569 process Effects 0.000 claims abstract description 31
- 238000000513 principal component analysis Methods 0.000 claims abstract description 11
- 230000001427 coherent effect Effects 0.000 claims abstract description 3
- 238000002474 experimental method Methods 0.000 claims description 31
- 239000011159 matrix material Substances 0.000 claims description 30
- 238000012216 screening Methods 0.000 claims description 21
- 230000015572 biosynthetic process Effects 0.000 claims description 17
- 238000003786 synthesis reaction Methods 0.000 claims description 17
- 238000012545 processing Methods 0.000 claims description 16
- 238000000605 extraction Methods 0.000 claims description 15
- 238000001914 filtration Methods 0.000 claims description 11
- 239000013598 vector Substances 0.000 claims description 9
- 230000009466 transformation Effects 0.000 claims description 7
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 230000035945 sensitivity Effects 0.000 claims description 3
- 230000001678 irradiating effect Effects 0.000 claims description 2
- 238000013480 data collection Methods 0.000 claims 1
- 238000013500 data storage Methods 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 18
- 238000012544 monitoring process Methods 0.000 abstract description 6
- 230000001052 transient effect Effects 0.000 abstract description 5
- 238000011156 evaluation Methods 0.000 abstract description 4
- 238000004364 calculation method Methods 0.000 abstract description 3
- 238000011160 research Methods 0.000 abstract description 3
- 230000001225 therapeutic effect Effects 0.000 abstract 1
- 238000002560 therapeutic procedure Methods 0.000 abstract 1
- 230000009471 action Effects 0.000 description 18
- 238000005516 engineering process Methods 0.000 description 9
- 238000001514 detection method Methods 0.000 description 8
- 230000008859 change Effects 0.000 description 5
- 238000009210 therapy by ultrasound Methods 0.000 description 5
- 238000011282 treatment Methods 0.000 description 5
- 238000012285 ultrasound imaging Methods 0.000 description 5
- 238000002679 ablation Methods 0.000 description 3
- 230000006399 behavior Effects 0.000 description 3
- 230000002537 thrombolytic effect Effects 0.000 description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 3
- 239000011358 absorbing material Substances 0.000 description 2
- 230000008499 blood brain barrier function Effects 0.000 description 2
- 210000001218 blood-brain barrier Anatomy 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000003111 delayed effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 229940079593 drug Drugs 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 239000007788 liquid Substances 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000000451 tissue damage Effects 0.000 description 2
- 231100000827 tissue damage Toxicity 0.000 description 2
- 206010028980 Neoplasm Diseases 0.000 description 1
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 230000006378 damage Effects 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 238000012637 gene transfection Methods 0.000 description 1
- 230000012010 growth Effects 0.000 description 1
- 230000023597 hemostasis Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
- 230000001960 triggered effect Effects 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/48—Diagnostic techniques
- A61B8/481—Diagnostic techniques involving the use of contrast agents, e.g. microbubbles introduced into the bloodstream
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Biophysics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Engineering & Computer Science (AREA)
- Biomedical Technology (AREA)
- Hematology (AREA)
- Physics & Mathematics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
Abstract
本发明提供了一种聚焦超声空化的实时高分辨时空分布成像方法与系统。对空化声散射信号进行滤波，得到稳态、惯性空化信号；基于高分辨相干系数对现有被动声学成像算法进行修正，从而得到高分辨的稳态、惯性空化图像；对空化图像进行筛选并去除附加空化能量，然后通过主成分分析得到稳态、惯性空化特征图像，在此基础上得到了稳态和惯性空化的时空分布。本发明可以克服主动超声成像只能得到空化微泡分布而不能得到空化活动时空分布，以及被动声学成像分辨率差且缺乏有效的空化时空分布计算方法的问题，为空化瞬态物理过程研究和聚焦超声治疗评价提供了一种有效手段，更为推进超声引导及监控的聚焦超声治疗系统在临床中的应用奠定了基础。
Description
技术领域
本发明属于超声空化物理与应用及超声成像技术领域，具体涉及一种聚焦超声空化的实时高分辨时空分布成像方法与系统。
背景技术
空化是指液体中的空化核在外加物理场(如超声场、激光场等)等作用下所呈现出的振荡、生长、收缩以致溃灭的一系列动力学过程，在诸如基因转染、体外碎石、药物释放、止血、溶栓及肿瘤热消融等领域中均有应用。根据气泡的不同动力学行为可以将空化分为稳态空化和惯性空化两种，其中稳态空化主要是指微泡在较低声压作用下的非线性振动行为，而惯性空化主要指微泡在较高声压下瞬间破裂，并形成高温、高压和冲击波等极端物理现象。空化的检测对于聚焦超声治疗的监控和治疗效果的评价都具有重要的意义。
最常用的空化检测方法是采用声学方法来检测聚焦超声作用过程中空化的声散射信号，例如次谐波、谐波、超谐波和宽带噪声等。早期主要利用单阵元聚焦超声换能器来发射脉冲并接收回波以实现主动检测，或不发射脉冲只被动接收声散射信号以实现被动检测，但这两种方法只能提供一维的时间信息而无法提供二维的空间分布。之后基于诊断超声换能器(如线阵和相控阵)发展了基于阵列的成像技术，最常用的则是传统的B模式超声成像；但该技术是通过逐线扫描得到，同一帧图像中不同扫描线之间存在一定的时间差，导致成像帧率较低，不能很好地描述空化的瞬态行为，且聚焦波会对空化微泡造成破坏。为了实现高帧率的成像，又发展了基于平面波发射的超快速超声成像技术，该技术可大大提高成像帧率，但由于平面波本身不聚焦且声压较低，导致成像分辨率低且信噪比低。近年来，有学者提出微秒时间分辨的超声空化微泡成像技术，该技术在传统B模式超声成像的基础上通过控制同步触发每次仅发射一次声脉冲，对空化微泡的辐射力作用几乎可以忽略，同时所得到的图像扫描线之间不存在时间差，时间分辨率可达到微秒，但是该技术对介质的可重复空化微泡分布要求严格，也就是说只能应用于液体或空腔等特殊环境中而不能用于生物软组织中，且扫描过程较为复杂。值得注意的是，以上方法都是在发射/接收的基础上实现的，都属于主动超声成像，而为了避免聚焦超声信号的干扰，只能选择在聚焦超声作用完成或者聚焦超声作用间歇来进行成像，因而这些方法成像的目标是聚焦超声作用后产生的空化微泡，而并非聚焦超声作用过程中的空化活动。
事实上，聚焦超声作用过程中的空化活动直接反映了当聚焦超声施加在介质上时介质所产生的响应，对此类空化活动的检测更有必要也更有意义。近年来，有学者提出了基于阵列换能器的被动声学成像技术，该技术是被动空化检测技术的拓展，通过关闭阵列换能器的脉冲发射而只被动接收聚焦超声焦域处的空化信号，再通过图像重建算法即可实现聚焦超声作用过程中空化活动的实时监控，目前已被广泛应用于实时监控热消融、碎石、组织损毁、超声溶栓和血脑屏障开放等方面。另外，聚焦超声作用过程中的空化活动随着时间和空间的联合变化，即空化的时空分布的研究也有着重要意义，一方面空化时空分布是研究治疗中空化瞬态物理过程的有效手段，另一方面空化时空分布也可用来实时定量观测治疗区域的变化。然而，现有的基于时间曝光声学的被动声学成像算法由于微泡间相互作用、阵列换能器缺陷、有限的阵列孔径和带宽等因素，导致其空间分辨性能极差；同时目前也缺乏有效的空化活动时空分布的计算方法，无法满足聚焦超声治疗实时监控系统的临床需求。
鉴于以上内容，亟待提出一种基于被动声学成像的聚焦超声空化的实时高分辨时空分布成像方法与系统。
发明内容
本发明的目的在于提供一种聚焦超声空化的实时高分辨时空分布成像方法与系统。
为了实现上述目的，本发明采用了以下技术方案：
一种聚焦超声空化的实时高分辨时空分布成像方法，包括以下步骤：
步骤一：利用线阵换能器被动接收聚焦超声辐照介质F次过程中对应产生的F帧空化声散射信号，对每一帧空化声散射信号通过滤波提取得到稳态空化信号和/或惯性空化信号；
步骤二：对由每一帧空化声散射信号中提取得到的稳态空化信号和/或惯性空化信号进行延时及补偿，得到稳态空化延时补偿信号和/或惯性空化延时补偿信号，通过Hilbert变换分别计算稳态空化延时补偿信号和/或惯性空化延时补偿信号的瞬时相位并根据对应瞬时相位的标准差计算相位相干系数，根据利用相位相干系数加权后所得的波束合成输出计算得到高分辨相干系数，对利用高分辨相干系数加权后所得的波束合成输出的平方进行积分，得到稳态空化和/或惯性空化的高分辨被动声学成像结果，由F帧空化声散射信号对应的稳态空化和/或惯性空化的高分辨被动声学成像结果构成稳态空化图像时间序列和/或惯性空化图像时间序列。
还包括以下步骤：
步骤三：针对所述稳态空化图像时间序列和/或惯性空化图像时间序列，分别根据其中每帧空化图像的空化能量最大值所对应成像位置是否在焦域内，对对应空化图像时间序列进行筛选，然后去除由前一次辐照对后一次辐照的影响而形成的附加空化能量，得到F帧纯净稳态空化图像和/或F帧纯净惯性空化图像；
步骤四：按照步骤一至步骤三进行多次重复实验，得到多次重复实验下的纯净稳态空化图像和/或纯净惯性空化图像；对多次重复实验下的同一帧纯净稳态空化图像和/或纯净惯性空化图像分别进行主成分分析，并根据主成分分量的方差对主成分分量进行加权，得到F帧稳态空化特征图像和/或F帧惯性空化特征图像；
步骤五：对F帧稳态空化特征图像和/或F帧惯性空化特征图像，分别根据对应F帧空化特征图像中各帧图像的空化能量最大值所对应成像位置得到轴向和横向最大能量分布曲线，并根据能量分布曲线半高宽确定的两轴向坐标和横向坐标计算得到横向和轴向平均能量分布；将F帧稳态空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，得到稳态空化的横向时空分布图像和轴向时空分布图像，将F帧惯性空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，得到惯性空化的横向时空分布图像和轴向时空分布图像。
上述步骤一，具体包括以下步骤：
1.1)采用聚焦超声辐照介质产生空化，采用线阵换能器被动接收聚焦超声辐照过程中的空化声散射信号，采用可编程全数字化超声成像设备的并行通道数据采集及存储模块对该线阵换能器接收的空化声散射信号进行采集；
1.2)利用巴特沃斯带通滤波器从线阵换能器第i(i＝1,2,...,N)个阵元接收的空化声散射信号中分别提取谐波、次谐波和超谐波成分，将这三种谐波成分相加得到稳态空化信号；从第i个阵元接收的空化声散射信号中减去稳态空化信号，然后利用巴特沃斯带阻滤波器滤除基波，得到惯性空化信号；
1.3)重复步骤1.2)，直至从线阵换能器的N个阵元接收的空化声散射信号中提取得到对应的稳态和/或惯性空化信号；
1.4)重复步骤1.1)～1.3)，直至采集到F帧空化声散射信号，并从每一帧空化声散射信号中提取得到稳态空化信号和/或惯性空化信号。
上述步骤二，具体包括以下步骤：
2.1)对线阵换能器第i(i＝1,2,...,N)个阵元接收的空化声散射信号在经过滤波之后进行延时和补偿，延时和补偿后的输出信号(即延时补偿信号)为：
其中，di(x,z)为成像位置(x,z)到第i个阵元(xi,0)的距离，η[di(x,z)]为补偿超声波球面传播衰减的接收阵列空间灵敏度补偿系数；pi(t)为第i个阵元接收的空化声散射信号在经过滤波之后得到的稳态或惯性空化信号；c为介质中声传播速度；
2.2)对步骤2.1)所得延时补偿信号进行Hilbert变换，得到第i个阵元的解析信号，计算第i个阵元解析信号的瞬时相位，根据瞬时相位标准差计算相位相干系数：
其中，γ为调整相位相干系数加权影响的控制参数，σ[Φ(x,z,t)]为信号瞬时相位标准差，Φ(x,z,t)为N个阵元解析信号的瞬时相位形成的矩阵，σ0为均匀分布[-π,π]的标准差；
2.3)利用步骤2.2)所得的相位相干系数对N个阵元对应延时补偿信号的加和信号进行加权，得波束合成输出qPCF(x,z,t)：
2.4)根据步骤2.3)所得波束合成输出，计算高分辨相干系数：
2.5)根据步骤2.4)所得的高分辨相干系数对N个阵元对应延时补偿信号的加和信号进行加权，得到波束合成输出qHRCF(x,z,t)：
2.6)在空化声散射信号采集时间区间[t0,t0+Δt]内对步骤2.5)所得的波束合成输出qHRCF(x,z,t)的平方进行积分，得到成像区域内每个成像位置处的空化能量I(x,z)：
其中，t0为空化声散射信号采集的起始时刻，Δt为空化声散射信号采集的时间长度；
2.7)对步骤一所得的与每一帧空化声散射信号对应的稳态空化信号和/或惯性空化信号均按照步骤2.1)～2.6)进行处理，得到F帧稳态空化图像和/或F帧惯性空化图像。
上述步骤三，具体包括以下步骤：
3.1)测量用于辐照介质的聚焦超声换能器的声场分布；
3.2)根据所述聚焦超声换能器的声场分布，计算得到所述聚焦超声换能器的焦域尺寸；分别寻找稳态空化图像时间序列和/或惯性空化图像时间序列中第k帧(k＝1,2,...,F)空化图像的空化能量最大值所在的成像位置，对于该成像位置位于焦域之外的对应帧空化图像，将该帧空化图像所有像素点赋值为0，从而完成对稳态空化图像时间序列和/或惯性空化图像时间序列的筛选；
3.3)经过步骤3.2)后，对稳态空化图像时间序列和/或惯性空化图像时间序列按照以下公式去除附加的空化能量：
其中，k＝2,3,...,F，Ik为经过步骤3.2)筛选后的稳态空化图像时间序列或惯性空化图像时间序列中的第k帧空化图像，为去除附加空化能量后的第k帧空化图像，为附加的空化能量：
其中，Ik-1为经过步骤3.2)筛选后的稳态空化图像时间序列或惯性空化图像时间序列中的第k-1帧空化图像，P为空化附加权重系数；
至此，得到去除了附加空化能量的F帧纯净稳态空化图像和/或F帧纯净惯性空化图像。
上述步骤四，具体包括以下步骤：
4.1)进行r次重复实验(同一参数条件下)，对每次重复实验获得的F帧空化声散射信号经空化信号提取后均按照步骤二和步骤三进行处理，得到每次重复实验对应的F帧纯净稳态空化图像和/或F帧纯净惯性空化图像；
4.2)将r次重复实验所得的第k帧纯净稳态空化图像都转换为l行×1列的列向量，l＝m×n，m和n分别为该图像的行数和列数，第k帧纯净稳态空化图像对应的r个列向量构成矩阵X，对矩阵X做标准化变换后构建协方差矩阵R，对R做特征值分解：
R＝UVUT
其中，V为特征值矩阵，该特征值矩阵的对角元素分别为λ1,λ2,...,λr，U＝[u1,u2,...,ur]为特征向量矩阵；
4.3)从步骤4.2)所得特征向量矩阵U中提取出前M个特征向量u1,u2,...,uM，并计算各个主成分分量：
其中i＝1,2,...,M，M为主成分数量；
4.4)计算步骤4.3)所得各主成分分量的方差，并利用该方差对主成分分量进行加权，得到加权后的主成分分量：
其中，为主成分分量中的第j个元素；
4.5)将步骤4.4)所得的加权后的主成分分量转换成m行×n列，得到第k帧稳态空化特征图像；
4.6)对F帧纯净稳态空化图像分别按照步骤4.2)～4.5)进行处理后，得到F帧稳态空化特征图像；
对r次重复实验所得纯净惯性空化图像，采用与步骤4.2)～4.6)相同的流程进行处理，得到F帧惯性空化特征图像。
上述步骤五，具体包括以下步骤：
5.1)对步骤四所得的F帧稳态空化特征图像A1,A2,...,AF，分别提取轴向最大能量分布曲线和横向最大能量分布曲线，根据曲线的半高宽，确定对应的轴向坐标z01和z02以及横向坐标x01和x02，根据轴向坐标和横向坐标对应计算横向平均能量分布和轴向平均能量分布k＝1,2,...,F；
5.2)对步骤5.1)所得的每一帧稳态空化特征图像的横向平均能量分布进行组合，得到稳态空化的横向时空分布图像；对步骤5.1)所得的每一帧稳态空化特征图像的轴向平均能量分布进行组合，得到稳态空化的轴向时空分布图像；
对步骤四所得的F帧惯性空化特征图像，采用与步骤5.1)～5.2)相同的流程进行处理，得到惯性空化的横向时空分布图像和轴向时空分布图像。
一种聚焦超声空化的实时高分辨时空分布成像系统，该成像系统包括空化发生装置以及空化信号检测装置，所述空化发生装置包括聚焦超声换能器、功率放大器以及控制聚焦超声换能器与功率放大器的时序同步的任意波形发生器，空化信号检测装置包括可编程全数字化超声成像设备，可编程全数字化超声成像设备包括用于被动接收在聚焦超声换能器辐照介质过程中产生的空化声散射信号的线阵换能器、用于对线阵换能器相应阵元接收的空化声散射信号进行采集和存储的模块(并行通道数据采集及存储模块)，以及空化实时高分辨时空分布成像模块；所述空化实时高分辨时空分布成像模块包括稳态和惯性空化信号提取子模块以及高分辨被动声学成像子模块；所述稳态和惯性空化信号提取子模块用于执行以上步骤一，所述高分辨被动声学成像子模块用于执行以上步骤二。
所述空化实时高分辨时空分布成像模块还包括稳态和惯性空化图像时间序列筛选及附加空化能量去除子模块、基于主成分分析的特征图像提取子模块，以及特征图像的各向平均能量分布组合子模块；所述稳态和惯性空化图像时间序列筛选及附加空化能量去除子模块用于执行以上步骤三，所述基于主成分分析的特征图像提取子模块用于执行以上步骤四，所述特征图像的各向平均能量分布组合子模块用于执行以上步骤五。
所述任意波形发生器一方面发出信号给功率放大器，经放大后输入给聚焦超声换能器，另一方面发出脉冲信号给可编程全数字化超声成像设备。聚焦超声换能器持续辐照介质以产生空化，由线阵换能器被动接收空化声散射信号，并通过可编程全数字化超声成像设备的并行通道数据采集及存储模块进行采集和存储，在辐照和采集过程进行F次后，对空化声散射信号采用空化实时高分辨时空分布成像模块进行处理。
所述成像系统还包括用于调整介质与聚焦超声辐照焦域的相对位置的三维移动装置。
本发明的有益效果体现在：
本发明针对传统发射/接收模式的主动超声成像无法对聚焦超声作用过程中空化活动的时空分布进行监控成像的缺陷，利用线阵换能器被动接收聚焦超声空化声散射信号并提取稳态、惯性空化信号；基于高分辨相干系数对现有被动声学成像算法进行修正(采用高分辨相干系数对线阵换能器阵元接收信号延时补偿后的加和信号进行加权)，从而得到高分辨的稳态、惯性空化图像，有效地提高了空化成像的空间分辨率，有助于获得空化的高分辨时空分布成像，本发明可应用于聚焦超声治疗中空化瞬态物理过程的监控和治疗效果评估。
进一步的，本发明通过主成分分析得到稳态、惯性空化特征图像序列，针对稳态空化、惯性空化分别通过提取多帧对应空化特征图像的横向和轴向平均能量分布，有效利用了每帧稳态、惯性空化特征图像中空化的空间分布信息，得到了稳态空化、惯性空化各自的横向时空分布和轴向时空分布。并且本发明通过筛选稳态、惯性空化时间序列图像，以及去除图像中的附加空化能量，可以得到更准确的空化时空分布。
附图说明
图1是本发明中聚焦超声空化发生及空化信号检测的系统框图；
图2是本发明中聚焦超声空化发生及空化信号检测的时序控制图；
图3是本发明中稳态和惯性空化信号的提取流程图；
图4是本发明中高分辨被动声学成像的流程图；
图5是本发明中根据高分辨被动声学成像和现有被动声学成像方法所得的稳态空化图像和惯性空化图像；其中：(a)和(b)分别为经过现有被动声学成像方法和本发明中高分辨被动声学成像方法所得的稳态空化图像，(c)和(d)分别为经过现有被动声学成像方法和本发明中高分辨被动声学成像方法所得的惯性空化图像；
图6是本发明中纯净的稳态和惯性空化图像的提取流程图；
图7是本发明中空化特征图像的提取流程图；
图8是本发明中空化横向时空分布和轴向时空分布的提取流程图。
具体实施方式
下面结合附图和实施例对本发明作进一步的详细说明。
本发明提供了一种基于被动声学成像的聚焦超声空化的实时高分辨时空分布成像方法，通过同步触发聚焦超声换能器和可编程全数字化超声成像设备来采集聚焦超声作用过程中的空化声散射信号，并设计滤波器提取出稳态和惯性空化信号；基于高分辨相干系数对现有被动声学成像方法进行改进，有效地提高成像空间分辨率，得到高分辨的稳态和惯性空化图像；接着通过测量聚焦超声换能器的焦域尺寸对稳态和惯性空化时间序列图像进行筛选并去除附加空化能量，得到纯净稳态和惯性空化图像；然后基于主成分分析得到稳态和惯性空化特征图像，最后提取多帧空化特征图像的横向和轴向平均能量分布得到稳态空化和惯性空化的横向时空分布和轴向时空分布，以下为具体步骤。
步骤一、搭建实验平台，聚焦超声辐照介质产生空化，利用线阵换能器被动接收聚焦超声辐照介质F次过程中对应产生的F帧空化声散射信号，对每一帧空化声散射信号通过滤波提取得到稳态空化信号和惯性空化信号；
参见图1，上述实验平台包括聚焦超声作用过程中空化声散射信号的高帧率被动采集装置，该装置包括：任意波形发生器1、功率放大器2、聚焦超声换能器3(例如，发射中心频率为1.2MHz)、可编程全数字化超声成像设备4。其中任意波形发生器1可实现任意幅度、频率、相位和波形的编写；功率放大器2为稳定可调节功率的信号放大设备，可将外部设备的输入信号进行不同程度的放大；聚焦超声换能器3为单阵元凹面球壳型换能器，可在介质6(例如，仿组织体模)内实现超声能量的有效聚集。可编程全数字化超声成像设备4采用线阵换能器5(例如，接收带宽为5～14MHz)接收空化声散射信号，在可编程全数字化超声成像设备4中可以通过程序编写来实现多种超声成像模式的切换，可编程全数字化超声成像设备4的并行通道数据采集及存储模块可实现对线阵换能器5接收到的原始超声数据的高帧率(例如，帧率≥5000Hz)采集和存储；线阵换能器5为诊断超声换能器，并与可编程全数字化超声成像设备4的接口相连接，在实验中通过编程将线阵换能器5所有阵元的发射变迹置为0，接收变迹置为1，以实现空化声散射信号的被动接收。介质6被固定在三维移动装置7，三维移动装置7可实现X、Y和Z三个方向的高精度移动，用以保证聚焦超声换能器1的焦域落在介质6中(焦域尺寸远小于介质6的尺寸)。实验操作在水槽8中进行，并在水槽8的侧壁和底部放置吸声材料9，以减少多次反射对实验的影响。
参见图2，被动接收空化声散射信号的时序控制由任意波形发生器1来实现，任意波形发生器的一个通道发出正弦信号为功率放大器2提供信号源，功率放大器2将正弦信号放大之后驱动聚焦超声换能器3工作，在聚焦超声能量达到空化阈值时辐照介质6并产生空化；任意波形发生器1的另一个通道发出脉冲触发信号用于触发可编程全数字化超声成像设备4，从而使得线阵换能器5被动接收空化声散射信号。需要注意的是，由于聚焦超声换能器3辐照介质6产生空化到空化声散射信号传播至线阵换能器5需要一定的时间，因此对于任意波形发生器1，上述脉冲触发信号应相对于上述正弦信号有一定延时(延时例如为80～100μs)。聚焦超声辐照介质6，使介质6产生空化，采用可编程全数字化超声成像设备的线阵换能器5被动接收聚焦超声作用过程中的空化声散射信号。
所述步骤一的具体流程如下：
(1.1)搭建好实验平台后，将介质6固定在三维移动装置7上，通过三维移动装置7对介质6位置的调整，使得聚焦超声换能器3的焦域位于介质6中；
(1.2)可编程全数字化超声成像设备4初始化后，聚焦超声换能器3辐照介质6产生空化，可编程全数字化超声成像设备4的并行通道数据采集及存储模块开始采集并存储线阵换能器阵元接收的空化声散射信号。
(1.3)对接收信号中的稳态和惯性空化信号进行提取
参见图3，从步骤(1.2)所得的空化声散射信号中提取第1个阵元接收的信号，设计多个阶数为4、带宽为300kHz的巴特沃斯带通滤波器，从第1个阵元接收的信号中提取谐波成分(中心频率为2f0,3f0,...,10f0,f0为聚焦超声换能器的发射中心频率)，并设计多个阶数为4、带宽为150kHz的巴特沃斯带通滤波器，从第1个阵元接收的信号中提取次谐波(中心频率为f0/2,f0/3,f0/4)和超谐波成分(中心频率为1.5f0,2.5f0,...,9.5f0)，将谐波、次谐波和超谐波三者相加即为稳态空化信号。从第1个阵元接收的信号中减去稳态空化信号，对剩余的信号利用阶数为4、带宽为300kHz的巴特沃斯带阻滤波器滤除掉基波(中心频率为f0)，滤波后的信号即为第1个阵元得到的惯性空化信号。
(1.4)重复(1.3)的过程，直到从线阵换能器5所有阵元(假设为N个，N例如为128)接收的信号中提取得到对应的稳态和惯性空化信号。
(1.5)重复步骤(1.2)～(1.4)，直至采集到F帧空化声散射信号，并从每一帧空化声散射信号中提取得到稳态空化信号和惯性空化信号。
步骤二、对由每一帧空化声散射信号中提取得到的稳态空化信号和惯性空化信号进行延时及补偿，得到稳态空化延时补偿信号和惯性空化延时补偿信号，通过Hilbert变换分别计算稳态空化延时补偿信号和惯性空化延时补偿信号的瞬时相位并根据对应瞬时相位的标准差计算相位相干系数；根据相位相干系数加权后所得的波束合成输出计算得到高分辨相干系数，利用高分辨相干系数对N个阵元对应的延时补偿信号的加和信号进行加权，对利用高分辨相干系数加权所得信号的平方在空化声散射信号采集时间区间(在上述脉冲触发信号触发后开始采集，例如，采样率为40～50MHz，每个通道的采样点数为4000～5000)内进行积分，得到成像区域内各成像位置处的空化能量，从而得到稳态空化和惯性空化的高分辨被动声学成像结果，由F帧空化声散射信号对应的稳态空化和惯性空化的高分辨被动声学成像结果构成稳态空化图像时间序列和惯性空化图像时间序列(图4)。
所述步骤二的具体流程如下：
(2.1)首先建立成像区域，一般来说，成像区域选择为聚焦超声焦点所在的区域。假设线阵换能器5有N个阵元并且聚焦在x-z平面上，则成像位置(x,z)到第i个阵元(xi,0)的距离为：
(2.2)对第i个阵元经过滤波之后得到的信号进行延时和补偿，延时和补偿后的输出信号为：
其中，η[di(x,z)]为超声波球面传播衰减的接收阵列空间灵敏度补偿系数，一般选择为pi(t)为第i个阵元接收的空化声散射信号经过滤波之后得到的信号，即步骤一提取到的稳态空化信号或惯性空化信号；c为介质中声传播速度；t表示时间；
(2.3)对步骤(2.2)所得信号进行Hilbert变换，得到第i个阵元的解析信号：
其中，j为虚数单位，*代表两个信号之间的卷积；
(2.4)根据步骤(2.3)所得解析信号的实部和虚部计算第i个阵元解析信号的瞬时相位：
其中imag[·]表示取虚部，real[·]表示取实部；
(2.5)根据步骤(2.4)得到N个阵元的解析信号的瞬时相位，由瞬时相位的标准差来计算相位相干系数：
其中，max{·}是为了保证相位相干系数的变化范围在0～1之间，γ为调整相位相干系数加权影响的控制参数(一般取为0.5～1)，σ[Φ(x,z,t)]为瞬时相位标准差，Φ(x,z,t)＝[φ1(x,z,t),φ2(x,z,t),...,φN(x,z,t)]为N个阵元解析信号的瞬时相位形成的矩阵，为均匀分布[-π,π]的标准差；
(2.6)根据步骤(2.5)所得的相位相干系数来计算相位相干系数加权后所得的波束合成输出：
(2.7)根据步骤(2.6)所得波束合成输出计算高分辨相干系数：
(2.8)根据步骤(2.7)所得的高分辨相干系数来计算最终的波束合成输出：
(2.9)在空化声散射信号采集时间区间[t0,t0+Δt]内对步骤(2.8)所得的qHRCF(x,z,t)的平方进行积分，则可得到成像区域内每个成像位置处的空化能量I(x,z)：
其中，t0为空化声散射信号采集的起始时刻(触发时刻)，Δt为空化声散射信号采集的时间长度；
(2.10)对步骤一所得的与每一帧空化声散射信号对应的稳态空化信号和惯性空化信号均按照步骤(2.1)～(2.9)进行处理，即进行稳态和惯性空化的高分辨被动声学成像，得到F帧稳态空化图像和F帧惯性空化图像；从图5可以看出，本发明中的高分辨被动声学成像方法可使得稳态空化图像和惯性空化图像的旁瓣水平得到有效降低、成像伪影得到有效抑制，从而使得图像的横向分辨率和轴向分辨率得到有效提高。
步骤三、搭建聚焦超声声场测量系统，得到聚焦超声换能器3的焦域尺寸。针对F帧稳态空化图像时间序列和F帧惯性空化图像时间序列，分别根据空化能量最大值所对应成像位置是否在焦域内，对空化图像时间序列进行筛选，然后去除由前一次辐照对后一次辐照的影响而形成的附加空化能量，得到F帧纯净稳态空化图像和F帧纯净惯性空化图像(图6)。
所述步骤三的具体流程如下：
(3.1)经过步骤二可以得到在特定声学参数下，随时间变化的稳态和惯性空化图像时间序列，但是由于空化的随机性、声学参数的改变及其他实验因素的影响，导致某些时刻在聚焦超声辐照下稳态或惯性空化并未发生，因此需要对不同时刻的空化成像结果进行甄别筛选。在筛选前首先需要测量聚焦超声换能器3的声场分布。
(3.2)搭建由聚焦超声换能器3、功率放大器2、针式水听器10、三维扫描装置11和超声信号接收器12构成的聚焦超声声场测量系统。由于功率较高有损坏水听器的风险，因此扫描过程中将功率放大器的输出功率设置为2W。首先将针式水听器10正对聚焦超声换能器3并固定于三维扫描装置11上，初步扫描声场得到声场分布参考图，将该参考图中声压最强点定为扫描平面中心点；然后设置扫描参数，开始进行声场扫描，由超声信号接收器12来接收焦域的超声信号。声场扫描在水槽8中进行，并在水槽8侧壁和底部放置吸声材料9。扫描过程中注意避免出现强反射物，并及时清除聚焦超声换能器3表面的气泡。
(3.3)根据步骤(3.2)得到聚焦超声换能器3的声场分布，计算出焦域的尺寸大小。寻找第1帧空化图像中能量最大值所在的坐标，并判断该坐标是否位于焦域，若不在，则将该帧空化图像所有像素点赋值为0(为了不影响时间上的连续性，该帧仍然参与图像时间序列的下一步处理)。根据这一判断方式则可完成对空化图像时间序列(步骤二得到的F帧稳态空化图像和F帧惯性空化图像)的筛选。
(3.4)由于聚焦超声辐照过程得到的每一帧空化图像既包括该次辐照产生的空化能量，也包括由于前一次辐照对后一次辐照的影响而造成的附加的空化能量。因此，为得到更准确的空化时空分布，需要将附加的空化能量去除：
其中，k＝2,3,...,F，F为帧数，Ik为经过步骤(3.3)筛选后的第k帧空化图像，为去除附加空化能量后的第k帧空化图像，为附加的空化能量；
(3.5)步骤(3.4)所述的附加空化能量为经过步骤(3.3)筛选后的第k-1帧空化图像Ik-1与空化附加权重系数P的乘积：
(3.6)步骤(3.5)所述的空化附加权重系数P为：
其中max(Ik)和max(Ik-1)分别为经过步骤(3.3)筛选后的第k帧和第k-1帧空化图像中的能量最大值；
(3.7)对经步骤(3.3)筛选所得的稳态和惯性空化图像均按照步骤(3.4)～(3.6)进行处理，则可得到去除了附加空化能量的F帧纯净稳态空化图像和F帧纯净惯性空化图像。
步骤四、在同一条件下按照步骤一至步骤三进行多次重复实验，得到多次重复实验下的纯净稳态空化图像和纯净惯性空化图像；对多次重复实验下的同一帧纯净稳态空化图像和纯净惯性空化图像分别进行主成分分析，并根据主成分分量的方差对主成分分量进行加权，得到F帧稳态空化特征图像和F帧惯性空化特征图像(图7)。
所述步骤四具体流程如下：
(4.1)同一参数下进行r次(例如，10～20次)重复实验，每次采集数据F帧，并对每帧数据经过稳态和惯性空化信号提取后均按照步骤二和步骤三处理，则每次重复实验可得F帧纯净稳态空化图像和F帧纯净惯性空化图像(图像均为m行×n列)，对纯净稳态空化图像和纯净惯性空化图像分别按照以下步骤处理。
(4.2)将r次重复实验所得的第1帧纯净空化图像都转换为列向量(l行×1列，l＝m×n)，则可得矩阵X：
X＝[E1,E2,...,Er]T (13)
其中，Ei为第i次重复实验所得纯净空化图像经过转换后得到的列向量，i＝1,2,...,r，[·]T代表矩阵的转置；
(4.3)对步骤(4.2)所得矩阵X做如下标准化变换，得到矩阵Z，矩阵Z的元素为：
其中，i＝1,2,...,r,j＝1,2,...,l，r和l分别为矩阵Z的行数和列数， xij为矩阵X中第i行第j列的元素；
(4.4)根据步骤(4.3)所得矩阵Z构建协方差矩阵：
(4.5)对步骤(4.4)所得协方差矩阵做特征值分解：
R＝UVUT (16)
其中，V为特征值矩阵，对角元素(特征值)分别为λ1,λ2,...,λr，且λ1≥λ2≥...≥λr，U＝[u1,u2,...,ur]为特征向量矩阵；
(4.6)根据来确定主成分数量M，从步骤(4.5)所得特征向量矩阵U中提取出前M个特征向量u1,u2,...,uM，并计算各个主成分分量：
其中，i＝1,2,...,M；
(4.7)计算步骤(4.6)所得各主成分分量的方差，并利用方差对主成分分量进行加权，得到加权后的主成分分量：
其中，为主成分分量中的第j个元素；
(4.8)将步骤(4.7)所得的加权后的主成分分量重新转换成m行×n列，则可得到r次重复实验所得第1帧空化特征图像；
(4.9)跳转至下一帧纯净空化图像，并重复步骤(4.2)～(4.8)，直到F帧纯净空化图像均被处理成空化特征图像。
(4.10)对步骤(4.1)所述的r次重复实验所得的纯净稳态空化图像和纯净惯性空化图像(每次重复实验均为F帧图像)分别按照步骤(4.2)～(4.9)进行处理后，得到F帧稳态空化特征图像和F帧惯性空化特征图像。
步骤五、针对步骤四所得的F帧稳态空化特征图像和F帧惯性空化特征图像，根据能量最大值所在坐标分别得到轴向和横向最大能量分布曲线，并根据其半高宽确定的两轴向坐标和横向坐标来计算横向和轴向平均能量分布，将F帧稳态空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，得到稳态空化的横向时空分布图像和轴向时空分布图像，将F帧惯性空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，得到惯性空化的横向时空分布图像和轴向时空分布图像(图8)。
(5.1)对步骤四所得的F帧空化特征图像记为A1,A2,...,AF，提取其中每一帧图像能量最大值所在坐标(x0,z0)，从而得到轴向最大能量分布曲线，根据该曲线的半高宽(曲线最大值下降到峰值一半时对应的宽度)，确定两轴向坐标z01和z02(z01<z02)，然后计算横向平均能量分布：
其中，k＝1,2,...,F，Ak,j为第k帧空化特征图像中轴向坐标为j时的横向能量分布；
(5.2)对上述F帧空化特征图像记为A1,A2,...,AF，提取其中每一帧图像能量最大值所在坐标(x0,z0)，从而得到横向最大能量分布曲线，根据该曲线的半高宽(曲线最大值下降到峰值一半时对应的宽度)，确定两横向坐标x01和x02(x01<x02)，然后计算轴向平均能量分布：
其中，k＝1,2,...,F，Ak,j为第k帧空化特征图像中横向坐标为j时的轴向能量分布；
(5.3)对步骤(5.1)和(5.2)所得的F帧空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，则可得到横向时空分布图像HT和轴向时空分布图像ZT：
(5.4)分别对步骤四得到的F帧稳态空化特征图像和F帧惯性空化特征图像按步骤(5.1)～(5.3)进行处理，则可分别得到稳态空化和惯性空化的横向时空分布和轴向时空分布。
本发明具有以下优点：
(1)本发明中线阵换能器只被动接收聚焦超声作用过程中产生的空化声散射信号，线阵换能器本身不发射信号，因此可以实现聚焦超声作用过程中空化活动的实时监控；与传统主动超声成像得到聚焦超声停止或作用间歇产生的空化微泡分布不同，本发明采用被动声学成像得到的是聚焦超声作用过程中空化活动的时空分布而非空化微泡分布。
(2)本发明中使用的被动声学成像算法在现有被动声学成像算法的基础上进一步改进，通过相位相干系数得到的高分辨相干系数，能有效地抑制旁瓣信号，从而大幅度地提高空化图像的空间分辨率。
(3)本发明中通过测量聚焦超声换能器的声场分布，将未产生空化的图像有效过滤，通过去除附加空化能量，避免了前一次聚焦超声辐照对后一次辐照的影响，得到纯净的空化图像，从而获得更准确的空化时空分布。本发明通过对多次重复实验下所得的同一帧空化图像进行主成分分析，可有效提取出稳态和惯性空化的特征图像，在此基础上得到稳态和惯性空化两种空化活动随时间和空间变化的规律。
(4)本发明一方面可以用于研究聚焦超声作用过程中空化活动的瞬态物理过程，另一方面也为组织热消融、组织毁损、超声溶栓、药物释放和血脑屏障开放等多种聚焦超声治疗中治疗区域变化和治疗效果评估提供了有效的手段。
总之，本发明提出了一种聚焦超声作用过程中空化活动的实时高分辨时空分布成像技术，可以克服主动超声成像只能得到空化微泡分布而不能得到空化活动时空分布，以及被动声学成像分辨率差且缺乏有效的空化时空分布计算方法的问题，为空化瞬态物理过程研究和聚焦超声治疗评价提供了一种有效手段，更为推进超声引导及监控的聚焦超声治疗系统在临床中的应用奠定了基础。
Claims (10)
1.一种聚焦超声空化的实时高分辨时空分布成像方法，该成像方法包括以下步骤：
步骤一：利用线阵换能器被动接收聚焦超声辐照介质F次过程中对应产生的F帧空化声散射信号，对其中每一帧空化声散射信号通过滤波提取得到稳态空化信号和/或惯性空化信号；其特征在于：
步骤二：对由每一帧空化声散射信号中提取得到的稳态空化信号和/或惯性空化信号进行延时及补偿，得到稳态空化延时补偿信号和/或惯性空化延时补偿信号，通过Hilbert变换分别计算稳态空化延时补偿信号和/或惯性空化延时补偿信号的瞬时相位并根据对应瞬时相位的标准差计算相位相干系数，根据利用相位相干系数加权后所得的波束合成输出计算得到高分辨相干系数，对利用高分辨相干系数加权后所得的波束合成输出的平方进行积分，得到稳态空化和/或惯性空化的高分辨被动声学成像结果，由F帧空化声散射信号对应的稳态空化和/或惯性空化的高分辨被动声学成像结果构成稳态空化图像时间序列和/或惯性空化图像时间序列。
2.根据权利要求1所述一种聚焦超声空化的实时高分辨时空分布成像方法，其特征在于：所述成像方法还包括以下步骤：
步骤三：针对所述稳态空化图像时间序列和/或惯性空化图像时间序列，分别根据其中每帧空化图像的空化能量最大值所对应成像位置是否在焦域内，对对应空化图像时间序列进行筛选，然后去除由前一次辐照对后一次辐照的影响而形成的附加空化能量，得到F帧纯净稳态空化图像和/或F帧纯净惯性空化图像；
步骤四：按照所述步骤一至步骤三进行多次重复实验，得到多次重复实验下的纯净稳态空化图像和/或纯净惯性空化图像；对多次重复实验下的同一帧纯净稳态空化图像和/或纯净惯性空化图像分别进行主成分分析，并根据主成分分量的方差对主成分分量进行加权，得到F帧稳态空化特征图像和/或F帧惯性空化特征图像；
步骤五：对F帧稳态空化特征图像和/或F帧惯性空化特征图像，分别根据对应F帧空化特征图像中各帧图像的空化能量最大值所对应成像位置得到轴向和横向最大能量分布曲线，并根据能量分布曲线半高宽确定的两轴向坐标和横向坐标计算得到横向和轴向平均能量分布；将F帧稳态空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，得到稳态空化的横向时空分布图像和轴向时空分布图像，将F帧惯性空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，得到惯性空化的横向时空分布图像和轴向时空分布图像。
3.根据权利要求1或2所述一种聚焦超声空化的实时高分辨时空分布成像方法，其特征在于：所述步骤一具体包括以下步骤：
1.1)采用聚焦超声辐照介质产生空化，采用线阵换能器被动接收聚焦超声辐照过程中的空化声散射信号，采用可编程全数字化超声成像设备的并行通道数据采集及存储模块对该线阵换能器接收的空化声散射信号进行采集；
1.2)利用巴特沃斯带通滤波器从线阵换能器第i(i＝1,2,...,N)个阵元接收的空化声散射信号中分别提取谐波、次谐波和超谐波成分，将这三种谐波成分相加得到稳态空化信号；从第i个阵元接收的空化声散射信号中减去稳态空化信号，然后利用巴特沃斯带阻滤波器滤除基波，得到惯性空化信号；
1.3)重复步骤1.2)，直至从线阵换能器的N个阵元接收的空化声散射信号中提取得到对应的稳态和/或惯性空化信号；
1.4)重复步骤1.1)～1.3)，直至采集到F帧空化声散射信号，并从每一帧空化声散射信号中提取得到稳态空化信号和/或惯性空化信号。
4.根据权利要求1或2所述一种聚焦超声空化的实时高分辨时空分布成像方法，其特征在于：所述步骤二具体包括以下步骤：
2.1)对线阵换能器第i(i＝1,2,...,N)个阵元接收的空化声散射信号在经过滤波之后进行延时和补偿，得延时补偿信号：
其中，di(x,z)为成像位置(x,z)到第i个阵元(xi,0)的距离，η[di(x,z)]为补偿超声波球面传播衰减的接收阵列空间灵敏度补偿系数；pi(t)为第i个阵元接收的空化声散射信号在经过滤波之后得到的稳态或惯性空化信号；c为介质中声传播速度；
2.2)对步骤2.1)所得延时补偿信号进行Hilbert变换，得到第i个阵元的解析信号，计算第i个阵元解析信号的瞬时相位，根据瞬时相位标准差计算相位相干系数：
其中，γ为调整相位相干系数加权影响的控制参数，σ[Φ(x,z,t)]为信号瞬时相位标准差，Φ(x,z,t)为N个阵元解析信号的瞬时相位形成的矩阵，σ0为均匀分布[-π,π]的标准差；
2.3)利用步骤2.2)所得的相位相干系数对N个阵元对应延时补偿信号的加和信号进行加权，得波束合成输出qPCF(x,z,t)：
2.4)根据步骤2.3)所得波束合成输出，计算高分辨相干系数：
2.5)根据步骤2.4)所得的高分辨相干系数对N个阵元对应延时补偿信号的加和信号进行加权，得到波束合成输出qHRCF(x,z,t)：
2.6)在空化声散射信号采集时间区间[t0,t0+Δt]内对步骤2.5)所得的波束合成输出qHRCF(x,z,t)的平方进行积分，得到成像区域内每个成像位置处的空化能量I(x,z)：
其中，t0为空化声散射信号采集的起始时刻，Δt为空化声散射信号采集的时间长度；
2.7)对步骤一所得的与每一帧空化声散射信号对应的稳态空化信号和/或惯性空化信号均按照步骤2.1)～2.6)进行处理，得到F帧稳态空化图像和/或F帧惯性空化图像。
5.根据权利要求2所述一种聚焦超声空化的实时高分辨时空分布成像方法，其特征在于：所述步骤三具体包括以下步骤：
3.1)测量用于辐照介质的聚焦超声换能器的声场分布；
3.2)根据所述聚焦超声换能器的声场分布，计算得到所述聚焦超声换能器的焦域尺寸；分别寻找稳态空化图像时间序列和/或惯性空化图像时间序列中第k帧(k＝1,2,...,F)空化图像的空化能量最大值所在的成像位置，对于该成像位置位于焦域之外的对应帧空化图像，将该帧空化图像所有像素点赋值为0，从而完成对稳态空化图像时间序列和/或惯性空化图像时间序列的筛选；
3.3)经过步骤3.2)后，对稳态空化图像时间序列和/或惯性空化图像时间序列按照以下公式去除附加的空化能量：
其中，k＝2,3,...,F，Ik为经过步骤3.2)筛选后的稳态空化图像时间序列或惯性空化图像时间序列中的第k帧空化图像，为去除附加空化能量后的第k帧空化图像，为附加的空化能量：
其中，Ik-1为经过步骤3.2)筛选后的稳态空化图像时间序列或惯性空化图像时间序列中的第k-1帧空化图像，P为空化附加权重系数。
6.根据权利要求2所述一种聚焦超声空化的实时高分辨时空分布成像方法，其特征在于：所述步骤四具体包括以下步骤：
4.1)进行r次重复实验，对每次重复实验获得的F帧空化声散射信号经空化信号提取后均按照步骤二和步骤三进行处理，得到每次重复实验对应的F帧纯净稳态空化图像和/或F帧纯净惯性空化图像；
4.2)将r次重复实验所得的第k帧纯净稳态空化图像都转换为l行×1列的列向量，l＝m×n，m和n分别为该图像的行数和列数，第k帧纯净稳态空化图像对应的r个列向量构成矩阵X，对矩阵X做标准化变换后构建协方差矩阵R，对R做特征值分解：
R＝UVUT
其中，V为特征值矩阵，该特征值矩阵的对角元素分别为λ1,λ2,...,λr，U＝[u1,u2,...,ur]为特征向量矩阵；
4.3)从步骤4.2)所得特征向量矩阵U中提取出前M个特征向量u1,u2,...,uM，并计算各个主成分分量：
其中i＝1,2,...,M，M为主成分数量；
4.4)计算步骤4.3)所得各主成分分量的方差，并利用该方差对主成分分量进行加权，得到加权后的主成分分量：
其中，为主成分分量中的第j个元素；
4.5)将步骤4.4)所得的加权后的主成分分量转换成m行×n列，得到第k帧稳态空化特征图像；
4.6)对F帧纯净稳态空化图像分别按照步骤4.2)～4.5)进行处理后，得到F帧稳态空化特征图像；
对r次重复实验所得纯净惯性空化图像，采用与步骤4.2)～4.6)相同的流程进行处理，得到F帧惯性空化特征图像。
7.根据权利要求2所述一种聚焦超声空化的实时高分辨时空分布成像方法，其特征在于：所述步骤五具体包括以下步骤：
5.1)对步骤四所得的F帧稳态空化特征图像A1,A2,...,AF，分别提取轴向最大能量分布曲线和横向最大能量分布曲线，根据曲线的半高宽，确定对应的轴向坐标z01和z02以及横向坐标x01和x02，根据轴向坐标和横向坐标对应计算横向平均能量分布和轴向平均能量分布k＝1,2,...,F；
5.2)对步骤5.1)所得的每一帧稳态空化特征图像的横向平均能量分布进行组合，得到稳态空化的横向时空分布图像；对步骤5.1)所得的每一帧稳态空化特征图像的轴向平均能量分布进行组合，得到稳态空化的轴向时空分布图像；
对步骤四所得的F帧惯性空化特征图像，采用与步骤5.1)～5.2)相同的流程进行处理，得到惯性空化的横向时空分布图像和轴向时空分布图像。
8.一种聚焦超声空化的实时高分辨时空分布成像系统，该成像系统包括空化发生装置以及空化信号检测装置，所述空化发生装置包括聚焦超声换能器(3)、与聚焦超声换能器(3)相连的功率放大器(2)以及控制聚焦超声换能器(3)与功率放大器(2)的时序同步的任意波形发生器(1)，空化信号检测装置包括可编程全数字化超声成像设备(4)，可编程全数字化超声成像设备(4)包括线阵换能器(5)以及空化实时高分辨时空分布成像模块；所述空化实时高分辨时空分布成像模块包括稳态和惯性空化信号提取子模块以及高分辨被动声学成像子模块；所述稳态和惯性空化信号提取子模块用于在线阵换能器(5)被动接收到聚焦超声辐照介质(6)过程中对应产生的F帧空化声散射信号后，对其中每一帧空化声散射信号通过滤波提取得到稳态空化信号和/或惯性空化信号；其特征在于：所述高分辨被动声学成像子模块用于对由每一帧空化声散射信号中提取得到的稳态空化信号和/或惯性空化信号进行延时及补偿、通过Hilbert变换分别计算稳态空化延时补偿信号和/或惯性空化延时补偿信号的瞬时相位、根据对应瞬时相位的标准差计算相位相干系数、根据利用相位相干系数加权后所得的波束合成输出计算高分辨相干系数，以及对利用高分辨相干系数加权后所得的波束合成输出的平方进行积分，从而得到由F帧空化声散射信号对应的稳态空化和/或惯性空化的高分辨被动声学成像结果构成的稳态空化图像时间序列和/或惯性空化图像时间序列。
9.根据权利要求8所述一种聚焦超声空化的实时高分辨时空分布成像系统，其特征在于：所述空化实时高分辨时空分布成像模块还包括稳态和惯性空化图像时间序列筛选及附加空化能量去除子模块、基于主成分分析的特征图像提取子模块，以及特征图像的各向平均能量分布组合子模块；所述稳态和惯性空化图像时间序列筛选及附加空化能量去除子模块用于针对所述稳态空化图像时间序列和/或惯性空化图像时间序列，分别根据其中每帧空化图像的空化能量最大值所对应成像位置是否在焦域内，对对应空化图像时间序列进行筛选，并在筛选后去除由前一次辐照对后一次辐照的影响而形成的附加空化能量，从而得到F帧纯净稳态空化图像和/或F帧纯净惯性空化图像；所述基于主成分分析的特征图像提取子模块用于对多次重复实验下的同一帧纯净稳态空化图像和/或纯净惯性空化图像分别进行主成分分析，并根据主成分分量的方差对主成分分量进行加权，从而得到F帧稳态空化特征图像和/或F帧惯性空化特征图像；所述特征图像的各向平均能量分布组合子模块用于将F帧稳态空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，从而得到稳态空化的横向时空分布图像和轴向时空分布图像，和/或将F帧惯性空化特征图像的横向平均能量分布和轴向平均能量分布分别进行组合，从而得到惯性空化的横向时空分布图像和轴向时空分布图像。
10.根据权利要求8所述一种聚焦超声空化的实时高分辨时空分布成像系统，其特征在于：所述任意波形发生器(1)一方面发出信号给功率放大器(2)，经放大后输入给聚焦超声换能器(3)，另一方面发出脉冲信号给可编程全数字化超声成像设备(4)，聚焦超声换能器(3)持续辐照介质(6)以产生空化，由线阵换能器(5)被动接收空化声散射信号，并通过可编程全数字化超声成像设备(4)的并行通道数据采集及存储模块进行采集及存储。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811083655.4A CN109431536B (zh) | 2018-09-17 | 2018-09-17 | 一种聚焦超声空化的实时高分辨时空分布成像方法与系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811083655.4A CN109431536B (zh) | 2018-09-17 | 2018-09-17 | 一种聚焦超声空化的实时高分辨时空分布成像方法与系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109431536A CN109431536A (zh) | 2019-03-08 |
| CN109431536B true CN109431536B (zh) | 2019-08-23 |
Family
ID=65532834
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811083655.4A Active CN109431536B (zh) | 2018-09-17 | 2018-09-17 | 一种聚焦超声空化的实时高分辨时空分布成像方法与系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109431536B (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110522992B (zh) * | 2019-07-22 | 2021-02-19 | 西安交通大学 | 基于空间非均匀聚焦涡旋声场的相变纳米液滴调控方法 |
| CN111220702B (zh) * | 2019-10-28 | 2023-01-13 | 大唐水电科学技术研究院有限公司 | 一种水轮机空蚀监测及评价方法 |
| CN111134719B (zh) * | 2019-12-19 | 2021-01-19 | 西安交通大学 | 一种聚焦超声辐照相变纳米液滴的主被动超声复合成像方法及系统 |
| CN113893468B (zh) * | 2020-06-22 | 2023-03-21 | 飞依诺科技(苏州)有限公司 | 治疗用超声波的调整方法、装置、计算机设备和存储介质 |
| CN112184879A (zh) * | 2020-08-25 | 2021-01-05 | 西安交通大学 | 一种反映空化泡尺寸时空分布的成像方法和系统 |
| CN112933435B (zh) * | 2021-02-05 | 2022-08-05 | 西安交通大学 | 低强度脉冲超声空化增效药物释放时的空化分类多参量分析与定征方法 |
| JP7137682B1 (ja) | 2021-11-29 | 2022-09-14 | ソニア・セラピューティクス株式会社 | 超音波治療装置 |
| CN114285494A (zh) * | 2022-01-06 | 2022-04-05 | 安徽省东超科技有限公司 | 多通道相控阵超声波发射系统 |
| CN116115260B (zh) * | 2023-01-03 | 2024-05-24 | 西安交通大学 | 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统 |
| CN118566852A (zh) * | 2024-06-06 | 2024-08-30 | 南京航空航天大学 | 一种抗间歇采样转发干扰的成像方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101642607A (zh) * | 2009-09-01 | 2010-02-10 | 西安交通大学 | 基于阵列换能器的低强度聚焦超声药物控释与监控装置 |
| CN102281918A (zh) * | 2008-11-05 | 2011-12-14 | 艾西斯创新有限公司 | 空化活动的绘制和表征 |
| CN103235041A (zh) * | 2013-04-26 | 2013-08-07 | 西安交通大学 | 基于超声主动空化成像的空化起始阈值分布重建方法 |
| CN104535645A (zh) * | 2014-12-27 | 2015-04-22 | 西安交通大学 | 微秒分辨空化时空分布的三维空化定量成像方法 |
| CN104887266A (zh) * | 2015-05-29 | 2015-09-09 | 西安交通大学 | 基于面阵的小区域三维被动空化成像及三维复合成像方法 |
| CN107260217A (zh) * | 2017-07-17 | 2017-10-20 | 西安交通大学 | 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统 |
| CN108392751A (zh) * | 2018-02-08 | 2018-08-14 | 浙江大学 | 一种实时监测高强聚焦超声治疗声空化的方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9669203B2 (en) * | 2011-03-01 | 2017-06-06 | University Of Cincinnati | Methods of enhancing delivery of drugs using ultrasonic waves and systems for performing the same |
| US9788811B2 (en) * | 2014-06-12 | 2017-10-17 | National Tsing Hua University | Imaging system of microbubble therapy and image evaluation method using the same |
-
2018
- 2018-09-17 CN CN201811083655.4A patent/CN109431536B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102281918A (zh) * | 2008-11-05 | 2011-12-14 | 艾西斯创新有限公司 | 空化活动的绘制和表征 |
| CN101642607A (zh) * | 2009-09-01 | 2010-02-10 | 西安交通大学 | 基于阵列换能器的低强度聚焦超声药物控释与监控装置 |
| CN103235041A (zh) * | 2013-04-26 | 2013-08-07 | 西安交通大学 | 基于超声主动空化成像的空化起始阈值分布重建方法 |
| CN104535645A (zh) * | 2014-12-27 | 2015-04-22 | 西安交通大学 | 微秒分辨空化时空分布的三维空化定量成像方法 |
| CN104887266A (zh) * | 2015-05-29 | 2015-09-09 | 西安交通大学 | 基于面阵的小区域三维被动空化成像及三维复合成像方法 |
| CN107260217A (zh) * | 2017-07-17 | 2017-10-20 | 西安交通大学 | 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统 |
| CN108392751A (zh) * | 2018-02-08 | 2018-08-14 | 浙江大学 | 一种实时监测高强聚焦超声治疗声空化的方法 |
Non-Patent Citations (1)
| Title |
|---|
| Spatial-temporal three-dimensional ultrasound plane-by-plane active cavitation mapping for high-intensity focused ultrasound in free field and pulsatile flow;Ting Ding et al，;《Ultrasonics》;20160412;第166-181页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109431536A (zh) | 2019-03-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109431536B (zh) | 一种聚焦超声空化的实时高分辨时空分布成像方法与系统 | |
| KR101651830B1 (ko) | 고강도 집속된 초음파를 위한 의료용 초음파 영상화에서의 피드백 | |
| JP6932192B2 (ja) | 超音波画像クラッタをフィルタリングする方法及びシステム | |
| Tong et al. | Multi-transmit beam forming for fast cardiac imaging-a simulation study | |
| US7101337B2 (en) | Method and non-invasive device for focusing acoustic waves | |
| CN111134719B (zh) | 一种聚焦超声辐照相变纳米液滴的主被动超声复合成像方法及系统 | |
| CN104535645B (zh) | 微秒分辨空化时空分布的三维空化定量成像方法 | |
| US11260247B2 (en) | Passive ultrasound imaging with sparse transducer arrays | |
| CN107260217B (zh) | 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统 | |
| Simpson et al. | Techniques for perfusion imaging with microbubble contrast agents | |
| CN104777484A (zh) | 压缩自适应波束合成的平面波超声成像和微泡成像的方法与系统 | |
| KR101610874B1 (ko) | 공간 일관성 기초 초음파 신호 처리 모듈 및 그에 의한 초음파 신호 처리 방법 | |
| CN109513123B (zh) | 一种基于半球阵的高分辨三维被动空化成像方法 | |
| CN111265245B (zh) | 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统 | |
| CN109864707B (zh) | 一种在有限视角情况下提高光声断层成像分辨率的方法 | |
| EP3665475B1 (en) | Shear wave elastography with ultrasound probe oscillation | |
| Prieur et al. | Observation of a cavitation cloud in tissue using correlation between ultrafast ultrasound images | |
| Gyöngy | Passive cavitation mapping for monitoring ultrasound therapy | |
| Park et al. | Photoacoustic imaging using array transducer | |
| Boulos et al. | Passive cavitation imaging using different advanced beamforming methods | |
| Wang et al. | High frame rate adaptive imaging using coherence factor weighting and the MVDR method | |
| Cheng et al. | Fourier based imaging method with steered plane waves and limited-diffraction array beams | |
| Li et al. | Evaluation of advanced Passive Acoustic Mapping (PAM) beamformers for high-duty-cycle HIFU ablated in Ex Vivo Tissue | |
| Tsysar et al. | Using a multi-rod waveguide system to create an ultrasound endoscope for imaging in aggressive liquids | |
| Guidi et al. | Experimental and simulation study of harmonic components generated by plane and focused waves |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |