CN107260217A - 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统 - Google Patents

用于脑部聚焦超声空化实时监控的三维无源成像方法及系统 Download PDF

Info

Publication number
CN107260217A
CN107260217A CN201710582424.7A CN201710582424A CN107260217A CN 107260217 A CN107260217 A CN 107260217A CN 201710582424 A CN201710582424 A CN 201710582424A CN 107260217 A CN107260217 A CN 107260217A
Authority
CN
China
Prior art keywords
mrow
msub
cavitation
ultrasonic
msup
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.)
Granted
Application number
CN201710582424.7A
Other languages
English (en)
Other versions
CN107260217B (zh
Inventor
万明习
路舒宽
余先波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201710582424.7A priority Critical patent/CN107260217B/zh
Publication of CN107260217A publication Critical patent/CN107260217A/zh
Application granted granted Critical
Publication of CN107260217B publication Critical patent/CN107260217B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/08—Detecting organic movements or changes, e.g. tumours, cysts, swellings
    • A61B8/0808—Detecting organic movements or changes, e.g. tumours, cysts, swellings for diagnosis of the brain
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/44—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device
    • A61B8/4483—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device characterised by features of the ultrasound transducer
    • 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/483—Diagnostic techniques involving the acquisition of a 3D volume of data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/58—Testing, adjusting or calibrating the diagnostic device
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Neurology (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Gynecology & Obstetrics (AREA)
  • Surgical Instruments (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明提供了一种用于脑部聚焦超声空化实时监控的三维无源成像方法及系统。在对颅骨遮挡导致的空化信号失真进行校准的基础上，对面阵两个方向的阵元做鲁棒的Capon波束合成以大幅度抑制来自其他方向的干扰和空化微泡间的相互作用；并基于空化信号相位间的差异引入相位相干系数对波束合成算法进行修正，从而在抑制成像伪影的同时提高成像分辨率；最后对三维空化体数据阈值化、平滑等处理。本发明可以克服传统磁共振和主动超声成像监控方法对空化的检测灵敏度不足、主动超声成像监控无法实现实时监控以及传统无源成像方法性能受限等问题，为临床提供了一种脑部聚焦超声治疗过程中空化效应的实时监控手段，使得脑部治疗的实时反馈与控制成为可能。

Description

用于脑部聚焦超声空化实时监控的三维无源成像方法及系统
技术领域
本发明属于超声检测与超声成像技术领域，特别涉及一种适用于脑部聚焦超声空化实时监控的三维无源成像方法及系统。
背景技术
空化效应是指介质中的空化核在一些外加物理场(例如超声场、激光场、微波场)的作用下产生的膨胀、收缩、塌陷的动力学过程。脑部聚焦超声治疗主要是利用包膜微泡在超声激励下的空化效应来达到治疗脑部疾病的目的。已有研究表明，聚焦超声联合微泡开放血脑屏障的主要机制是空化效应，在较低声压下微泡发生稳态空化，此时与大脑毛细血管尺寸相当的微泡在振动过程中与毛细血管相互作用，将构成血脑屏障的内皮细胞分开，从而有利于血脑屏障的开放。另外，空化效应已被证实为超声溶栓的主要机理，超声空化效应使栓塞部位的包膜微泡发生振动及瞬间破裂，血栓表面变软，通过损坏栓块中的纤维蛋白结构，有利于血栓的崩溃瓦解，增加溶栓药物与栓块的结合位点，从而加速溶栓。作为一种人工空化核，包膜微泡可以极大地降低空化阈值从而减少脑组织的额外损伤。然而，微泡在诊断超声声压下的惯性收缩力也可能会导致微泡坍塌，从而对脑部细胞结构产生破坏或者对血管壁形成不必要的损伤。由于空化效应瞬时、随机、不可控，因此对脑部聚焦超声治疗中的空化效应进行实时监控是非常必要的。
近年来，国内外研究表明，磁共振引导的经颅聚焦超声技术成为脑部超声治疗的新手段。磁共振监控的优势主要体现在较高的图像分辨率、对组织温度实时变化的敏感性以及与聚焦超声治疗过程互不干扰等方面，然而常规的磁共振成像主要关注组织结构的变化，难以实现对空化瞬态物理过程的实时检测与监控。传统以聚焦波为发射方式的B模式超声成像在脑部监控存在超声信号衰减、探测深度以及检测灵敏度等相互制约的矛盾。且为了保证B模式成像与聚焦超声信号互不影响，B模式成像只能在聚焦超声脉冲发射之后或者不连续的聚焦超声激励下提供反馈图像信息，无法对聚焦超声作用过程中的空化进行监控，因此B模式成像也无法实现对脑部超声治疗过程中空化的实时监控成像。以平面波为发射方式的超声成像极大地提高了成像帧频，但是也只能对聚焦超声治疗完成之后的空化进行快速的监控，依然无法与聚焦超声治疗信号做到同步。
因此，为了克服上述问题，近年来兴起了基于诊断超声换能器的无源成像技术。无源成像技术是使诊断超声换能器不发射脉冲信号只接收聚焦超声换能器焦点处的空化信号，因此可以对聚焦超声作用过程中的空化效应进行实时监控且相比主动超声成像技术具有更高的检测灵敏度。无源成像技术根据无源波束合成方法获得空化源的空间分布，以此来对聚焦超声治疗进行实时监控。聚焦超声治疗监控的准确性取决于无源波束合成方法的性能。当前，脑部聚焦超声空化的无源成像技术主要是根据CT数据获取成像区域有效声速，然后使用基于延时-叠加-积分的时间曝光声学(Time Exposure Acoustic,TEA)方法进行成像。由于受成像深度及换能器孔径的限制，此方法的横向/轴向分辨率较差，且微泡间的相互作用会形成严重的成像伪影。在非脑部聚焦超声空化的无源成像中，鲁棒的Capon波束合成(Robust Capon Beamformer,RCB)方法通过对每个成像位置的导向矢量重新估计并计算加权矢量，可以有效地消除微泡相互作用产生的伪影并同时提高横向/轴向分辨率，但成像分辨率依然较差且成像伪影依然不能被完全消除。且当前无源成像技术均基于一维线阵超声换能器，只能提供脑部治疗过程中空化效应的平面分布，不能得到空化的三维空间分布，无法满足脑部聚焦超声治疗实时监控的临床需求。
鉴于以上内容，很有必要提出一种针对脑部聚焦超声空化实时监控的三维无源成像方法。
发明内容
针对磁共振成像对脑部空化微泡信号的检测不灵敏、传统B模式和平面波主动成像的探测深度和检测灵敏度互相制约且无法监控聚焦超声作用过程中的空化效应、以及现有二维无源成像技术在分辨率低、伪影高和无法提供空化三维空间信息等方面的问题，本发明的目的在于提供一种用于脑部聚焦超声空化实时监控的三维无源成像方法及系统。
为实现上述目的，本发明采用了以下技术方案：
一种用于脑部聚焦超声空化实时监控的三维无源成像方法，包括以下步骤：
步骤一：计算超声面阵换能器每个阵元方向上超声波的有效声速；
步骤二：利用拉格朗日插值法对超声面阵换能器接收的聚焦超声换能器焦点区域的空化信号进行插值，然后根据步骤一所得有效声速对空化信号进行延时；
步骤三：构造第j列阵元经插值及延时后的空化信号的协方差矩阵，然后基于鲁棒的Capon波束合成方法计算最优加权系数；对第j列阵元经插值及延时后的空化信号进行二值化并计算相位相干系数；根据该最优加权系数以及相位相干系数，计算第j列阵元的空化源强度波束合成信号；
步骤四：根据步骤三计算的超声面阵换能器各列阵元的空化源强度波束合成信号构造协方差矩阵，然后基于鲁棒的Capon波束合成方法计算最优加权系数；对各列阵元的空化源强度波束合成信号进行二值化并计算相位相干系数；
步骤五：根据各列阵元的空化源强度波束合成信号以及步骤四所得最优加权系数和相位相干系数计算超声面阵换能器所有阵元的空化源强度波束合成信号和每个成像位置的空化源能量，得到空化源三维体数据，然后进行脑部聚焦超声空化的三维显示。
上述步骤一中，每个阵元方向上超声波的有效声速按照以下公式计算：
其中K为除颅骨以外不同介质的数目，Li,j为空化信号到达超声面阵换能器阵元(i,j)的距离，Di,j,k和ci,j,k分别为空化信号到达阵元(i,j)的声传播路径上除颅骨以外不同介质的厚度和声速，Di,j,s为空化信号到达阵元(i,j)的声传播路径上颅骨对应位置的厚度，ci,j,s为超声波在所述颅骨对应位置的声速；
超声波在所述颅骨对应位置的声速ci,j,s根据以下公式计算：
ci,j,s＝c1Φi,j,s+c2(1-Φi,j,s)
其中Φi,j,s为所述颅骨对应位置的孔隙度，c1为超声波在水中的传播速度，c2为超声波在颅骨中的最大传播速度。
上述步骤二，具体包括以下步骤：
(2.1)设置超声面阵换能器的阵元变迹，使换能器工作在不发射脉冲信号只接收信号的模式，同步触发聚焦超声换能器和超声面阵换能器，由超声面阵换能器采集聚焦超声换能器作用过程中的焦点处空化的三维射频信号pi,j，其中i＝1,2,...,M；j＝1,2,...,N；M和N分别为超声面阵换能器y和x方向的阵元数目；
(2.2)根据以下拉格朗日插值多项式对每个阵元采集到的射频信号进行插值：
其中，为插值后的空化信号，ki,j(tl,t)为拉格朗日基本多项式，pi,j(tl)为空化信号采样点，P为采样点个数；
(2.3)根据三维成像区域和成像精度，对步骤(2.2)中所得信号进行延时处理，得到延时后的信号si,j(x,y,z,t)：
其中di,j(x,y,z)为成像位置(x,y,z)到超声面阵换能器阵元位置的空间距离，ci,j为阵元方向上超声波的有效声速。
上述步骤三，具体包括以下步骤：
(3.1)在一定的时间区间[T0,T0+ΔT]内构造协方差矩阵：
其中sj(x,y,z,t)为第j列阵元经插值及延时后的空化信号矩阵，[·]T代表矩阵的转置；
(3.2)对步骤(3.1)所得协方差矩阵做特征值分解：
其中Vj为特征值矩阵，Uj中第i列为Vj中对角线上第i个元素对应的特征向量；
(3.3)根据步骤(3.2)中所得特征值矩阵Vj和特征向量矩阵Uj以及拉格朗日迭代法计算每个成像位置的导向矢量
其中为M行1列(M×1)预设导向矢量，λ为拉格朗日乘数，I为M行M列(M×M)单位矩阵；
(3.4)基于鲁棒的Capon波束合成方法并利用步骤(3.1)所得协方差矩阵Rj和步骤(3.3)所得导向矢量计算最优加权系数wj：
其中||·||代表欧几里得范数，[·]-1代表矩阵的逆；
(3.5)将第j列阵元经插值及延时后的空化信号si,j(x,y,z,t)中大于零的点置为1，小于等于零的点置为-1，得到二值化信号bi,j(x,y,z,t)，然后计算相位相干系数
其中p为可调参数；
(3.6)利用步骤(3.4)所得最优加权系数以及步骤(3.5)所得相位相干系数计算第j列阵元的空化源强度波束合成信号：
上述步骤四，具体包括以下步骤：
(4.1)计算超声面阵换能器各列阵元的空化源强度波束合成信号矩阵：
q(x,y,z,t)＝[q1(x,y,z,t)；q2(x,y,z,t)；...；qN(x,y,z,t)]
其中qj(x,y,z,t)为第j列阵元的空化源强度波束合成信号；
(4.2)在一定的时间区间[T0,T0+ΔT]内，根据步骤(4.1)所得空化源强度波束合成信号矩阵q(x,y,z,t)计算协方差矩阵：
其中[·]T代表矩阵的转置；
(4.3)对步骤(4.2)所得协方差矩阵做特征值分解：
R(x,y,z)＝UVUT
其中V为特征值矩阵，U中第j列为V中对角线上第j个元素对应的特征向量，j＝1,2,...,N；
(4.4)根据步骤(4.3)中所得特征值矩阵V和特征向量矩阵U以及拉格朗日迭代法计算每个成像位置的导向矢量
其中为N行1列(N×1)预设导向矢量，λ为拉格朗日乘数，I为N行N列(N×N)单位矩阵；
(4.5)基于鲁棒的Capon波束合成方法并利用步骤(4.2)所得协方差矩阵R和步骤(4.4)所得导向矢量计算最优加权系数w：
其中||·||代表欧几里得范数，[·]-1代表矩阵的逆；
(4.6)将每列阵元的空化源强度波束合成信号qj(x,y,z,t)中大于零的点置为1，小于等于零的点置为-1，得到二值化信号bj(x,y,z,t)，然后计算相位相干系数SCFp(x,y,z,t)：
其中p为可调参数。
上述步骤五中，所有阵元的空化源强度波束合成信号按照以下公式计算：
Q(x,y,z,t)＝wTq(x,y,z,t)SCFp(x,y,z,t)
在一定的时间区间[T0,T0+ΔT]内对Q(x,y,z,t)的平方进行积分，得到每个成像位置处的空化源能量I(x,y,z)：
上述步骤五中，三维显示具体步骤为：
(5.1)对所得的三维体数据进行阈值化处理，阈值一般取体数据最大值的0.1～0.2倍；
(5.2)选择滤波器对步骤(5.1)所得阈值化后的三维体数据进行平滑处理，然后提取出体素等值面；
(5.3)根据步骤(5.2)所得等值面构造三维曲面，然后设置颜色、光照、视角等进行三维显示。
一种用于脑部聚焦超声空化实时监控的三维无源成像系统，包括超声面阵换能器、声速计算模块、空化信号校准模块、自适应波束合成模块以及三维显示模块；
所述声速计算模块计算超声面阵换能器每个阵元方向上超声波的有效声速(参见上述步骤一)；
所述空化信号校准模块对超声面阵换能器接收的聚焦超声换能器焦点区域的空化信号进行插值以及延时处理(参见上述步骤二)；
所述自适应波束合成模块对超声面阵换能器两个方向的阵元做鲁棒的Capon波束合成，并基于空化信号相位间的差异引入相位相干系数对波束合成算法进行修正；
所述三维显示模块利用自适应波束合成模块得到的空化源三维体数据进行脑部聚焦超声空化的三维显示。
本发明的有益效果体现在：
本发明是一种针对脑部聚焦超声空化实时监控的三维无源成像技术，针对磁共振成像、传统B模式和平面波主动成像以及现有无源成像方法在监控脑部聚焦超声治疗方面存在的缺陷，通过同步触发聚焦超声治疗换能器和超声面阵换能器采集治疗过程中空化的三维射频信号并根据CT扫描数据对信号延时进行校准。基于鲁棒的Capon波束合成方法对面阵阵元的延时信号进行波束合成，对来自其他方向的信号进行了大幅度的抑制；根据通道信号相位间的差异引入了相位相干系数对波束合成算法进行改进，提高了成像分辨率，从而准确地反映了空化在三维空间中的分布，对于脑部聚焦超声治疗过程中空化效应的实时监控有着重要的意义。
附图说明
图1为脑部CT扫描数据及用于三维无源成像方法的颅骨内声速获取示意图；
图2为超声面阵换能器获取的空化信号的插值及延时流程图；
图3为某超声面阵换能器有无颅骨遮挡时空化信号的时频域对比图；
图4为超声面阵换能器第j列阵元的空化源强度波束合成的流程图；
图5为超声面阵换能器每列阵元波束合成中所得的协方差矩阵、特征向量矩阵、特征值、导向矢量、最优加权系数和相位相干系数结果；
图6为超声面阵换能器所有阵元的空化源强度波束合成的流程图；
图7为超声面阵换能器所有阵元波束合成中所得的协方差矩阵、特征向量矩阵、特征值、导向矢量、最优加权系数和相位相干系数结果；
图8为基于不同无源成像算法的三维空化成像x-z切面结果；
图9为三维无源成像显示流程图；
图10为不同条件下的三维无源成像显示结果图。
具体实施方式
下面结合附图和实施例对本发明做详细描述。
本发明涉及到的用于脑部聚焦超声空化实时监控的三维无源成像方法包括以下具体步骤：
步骤一：根据颅骨CT扫描数据和基于颅骨孔隙度的经验模型计算超声波在颅骨中的传播速度，计算超声面阵换能器每个阵元方向上超声波的有效声速。其具体步骤如下：
(1.1)获取颅骨CT扫描切片数据，图1中(a)为获取的颅骨CT扫描切片数据，(b)为利用(a)数据得到的颅骨三维重建结果。根据图1(a)所示的颅骨CT扫描切片数据计算空化信号到达阵元(i,j)的声传播路径上颅骨对应位置的孔隙度：
其中Hi,j为亨氏单位；
(1.2)在颅骨孔隙度经验模型的基础上计算超声波在颅骨对应位置的声速ci,j,s：
ci,j,s＝c1Φi,j,s+c2(1-Φi,j,s)
其中c1为超声波在水中的传播速度，c2为超声波在颅骨中的最大传播速度；
(1.3)图1中(c)为颅骨孔隙度计算结果，(d)和(e)分别为颅骨内部密度和超声波传播速度的分布结果，密度和声速分布图像基本一致，证实超声波在密度较大的地方传播较快；且两幅图像均显示颅骨内侧和外层密度和声速较大，而中间层减小；这是因为颅骨内外两侧为密质层而中间层为松质层，有静脉通过。(f)为(e)的局部放大结果，用于三维无源成像中有效声速的计算；成像区域选择颅骨较薄的颞窗部位，以尽可能减小空化信号的衰减；
(1.4)在实际实验过程中，通过测量空化信号到超声面阵换能器传播路径上不同介质(例如颅骨、液体、脑组织、脑组织仿体等)的厚度和声速，计算每个阵元方向上的有效声速
其中K为除颅骨以外不同介质的数目，Li,j为空化信号到达超声面阵换能器阵元(i,j)的距离，Di,j,k和ci,j,k分别为空化信号到达阵元(i,j)声传播路径上除颅骨以外不同介质的厚度和声速，Di,j,s为颅骨对应位置的厚度，ci,j,s为根据步骤(1.2)计算的声速。
步骤二：超声面阵换能器接收聚焦超声换能器焦点区域的空化信号，利用拉格朗日插值法对空化信号进行插值，并根据步骤一所得有效声速对空化信号进行延时，得到延时信号。其具体步骤如下(图2)：
(2.1)设置超声面阵换能器的阵元变迹，使换能器工作在不发射脉冲信号只接收信号的模式，通过波形发生器编写触发波形来同步触发聚焦超声换能器和面阵换能器，采集聚焦超声换能器作用过程中的焦点处空化的三维射频信号pi,j，其中i＝1,2,...,M；j＝1,2,...,N；M和N分别为面阵换能器y和x方向的阵元数目(即面阵中阵元的行数和列数)；
(2.2)由于实际实验中信号采样率的限制，可能会造成信号延时的精度不够，因此在信号延时之前对每个阵元采集到的空化信号进行插值；假设原信号采样点为pi,j(t0),pi,j(t1),...,pi,j(tP-1)，对P个采样点进行拉格朗日插值，计算当前列当前阵元的拉格朗日基本多项式：
(2.3)设置插值精度，一般在每两个采样点之间插入4～8个点，根据拉格朗日插值多项式计算原采样点之外的插值点的函数值，当前列当前阵元的拉格朗日插值多项式为：
其中ki,j(tl,t)为拉格朗日基本多项式，为插值后的空化信号；
(2.4)跳转至当前列下一个阵元重复步骤(2.2)和(2.3)，直到当前列所有阵元采集到的空化信号均插值完毕；
(2.5)跳转至下一列阵元并重复步骤(2.2)～(2.4)，直到面阵所有阵元采集到的空化信号均插值完毕；
(2.6)规划三维成像区域和成像精度，一般三维成像区域横向选择为面阵换能器的尺寸，轴向选择为成像区域与面阵之间距离的2～3倍；当需要观察局部放大的空化分布时，可适当地减小成像区域。成像精度一般为所选择的成像区域中有(100～200)×(100～200)×(100～200)个体素，成像精度根据需要可以适当增加或减小；
(2.7)计算每个成像位置(x,y,z)到面阵换能器当前列当前阵元位置(xi,yj,0)的空间距离：
(2.8)结合步骤(1.4)计算得到的有效声速ci,j对插值后的空化信号做延时处理，延时后的信号为：
(2.9)跳转至当前列下一个阵元并重复步骤(2.7)和(2.8)，直到当前列所有阵元采集到的信号均延时完毕；
(2.10)跳转至下一列阵元并重复步骤(2.7)～(2.9)，直到面阵所有阵元采集到的信号均延时完毕；
如图3所示，其中(a)为超声面阵换能器阵列示意图，由M×N个阵元组成，阵元均匀分布在32行32列栅格上，阵元总数为1024，孔径大小为38.4×38.4mm。(b)和(c)为无颅骨遮挡条件下经插值及延时的空化信号和相应的频谱分布；(d)和(e)为有颅骨遮挡条件下经插值及延时的空化信号和相应的频谱分布。结果显示有颅骨遮挡的空化信号强度明显低于无颅骨遮挡情况，说明颅骨对空化信号有一定的衰减；无颅骨遮挡下的空化信号频谱成分尤其高频成分丰富而有颅骨遮挡下空化信号中约7MHz以上的高频成分被抑制，说明颅骨对于空化信号相当于一个低通滤波器。
步骤三：构造每列阵元经插值及延时后的空化信号的协方差矩阵，然后基于鲁棒的Capon波束合成方法(Li J,Stoica P,Wang Z.On robust Capon beamforming anddiagonal loading[J].Signal Processing IEEE Transactions on,2003,51(7):1702-1715.)计算最优加权系数；对每列阵元经插值及延时后的空化信号进行二值化并计算相位相干系数；根据该最优加权系数以及相位相干系数，计算每列阵元的空化源强度波束合成信号。其具体步骤如下(图4)：
(3.1)提取第j列M个阵元经过拉格朗日插值及延时之后的空化信号；
(3.2)在一定的时间区间[T0,T0+ΔT]内构造协方差矩阵：
其中sj(x,y,z,t)＝[s1,j(x,y,z,t)；s2,j(x,y,z,t)；...；sM,j(x,y,z,t)]为第j列M个阵元经插值及延时后的空化信号矩阵，[·]T代表矩阵的转置；时间长度ΔT取决于聚焦超声每次照射的时间；
(3.3)对步骤(3.2)所得协方差矩阵做特征值分解：
其中Vj为特征值矩阵，Uj中第i列为Vj中对角线上第i个元素(特征值)对应的特征向量，i＝1,2,...,M；
(3.4)设置预设导向矢量和拉格朗日乘数λ，结合步骤(3.3)中所得特征值矩阵Vj和特征向量矩阵Uj计算每个成像位置的导向矢量
其中，I为单位矩阵(M×M)；
(3.5)基于鲁棒的Capon波束合成方法并利用步骤(3.2)所得协方差矩阵Rj和步骤(3.4)所得导向矢量计算最优加权系数wj：
其中||·||代表欧几里得范数，[·]-1代表矩阵的逆；
(3.6)将第j列当前阵元经插值及延时后的空化信号中大于零的点置为1，小于等于零的点置为-1，得到二值化信号：
(3.7)跳转至下一个阵元并重复步骤(3.6)，直到第j列所有阵元经插值及延时后的空化信号均已被二值化；
(3.8)根据步骤(3.7)所得第j列M个阵元的二值化信号计算相位相干系数
其中p为可调参数以改善图像质量，p的范围一般为0～2；
(3.9)利用步骤(3.5)所得最优加权系数以及步骤(3.8)所得相位相干系数计算第j列M个阵元的空化源强度波束合成信号：
如图5所示，(a)为每列阵元经插值及延时后的空化信号矩阵sj(x,y,z,t)所得到的协方差矩阵，(b)和(c)分别为(a)中协方差矩阵经过特征值分解得到的特征向量矩阵和特征值，(d)为根据(b)中特征向量矩阵和(c)中特征值计算得到的导向矢量，(e)为利用(a)中协方差矩阵和(d)中导向矢量得到的最优加权系数，(f)为利用步骤(3.8)得到的相位相干系数。
步骤四：根据步骤三计算的超声面阵换能器各列阵元的空化源强度波束合成信号构造协方差矩阵，然后基于鲁棒的Capon波束合成方法计算最优加权系数；对各列阵元的空化源强度波束合成信号进行二值化并计算相位相干系数。其具体步骤如下(图6)：
(4.1)重复步骤(3.1)～(3.9)，直到第1～N列阵元的空化源强度波束合成信号均被计算完毕；
(4.2)构造第1～N列阵元的空化源强度波束合成信号矩阵：
q(x,y,z,t)＝[q1(x,y,z,t)；q2(x,y,z,t)；...；qN(x,y,z,t)]
(4.3)根据步骤(4.2)所得第1～N列阵元的空化源强度波束合成信号矩阵q(x,y,z,t)计算协方差矩阵：
其中[·]T代表矩阵的转置；
(4.4)对步骤(4.3)所得协方差矩阵做特征值分解：
R(x,y,z)＝UVUT
其中V为特征值矩阵，U中第j列为V中对角线上第j个元素(特征值)对应的特征向量，j＝1,2,...,N；
(4.5)设置预设导向矢量和拉格朗日乘数λ，结合步骤(4.4)中所得特征值矩阵V和特征向量矩阵U计算每个成像位置的导向矢量
其中，I为单位矩阵(N×N)；
(4.6)基于鲁棒的Capon波束合成方法并利用步骤(4.3)所得协方差矩阵R和步骤(4.5)所得导向矢量计算最优加权系数w：
其中||·||代表欧几里得范数，[·]-1代表矩阵的逆；
(4.7)将第j列阵元的空化源强度波束合成信号中大于零的点置为1，小于等于零的点置为-1，得到二值化信号bj(x,y,z,t)：
(4.8)跳转至下一列阵元并重复步骤(4.7)，直到N列阵元的空化源强度波束合成信号均已被二值化；
(4.9)根据步骤(4.8)所得N个二值化信号计算相位相干系数SCFp(x,y,z,t)：
其中p为可调参数以改善图像质量，p的范围一般为0～2；
如图7所示，(a)为根据N列阵元的空化源强度波束合成信号矩阵q(x,y,z,t)所得到的协方差矩阵，(b)和(c)为(a)中协方差矩阵经过特征值分解得到的特征向量矩阵和特征值，(d)为根据(b)中特征向量矩阵和(c)中特征值计算得到的导向矢量，(e)为利用(a)中协方差矩阵和(d)中导向矢量得到的最优加权系数，(f)为利用步骤(4.9)得到的相位相干系数。
步骤五：根据步骤四所得每列阵元的空化源强度波束合成信号、最优加权系数和相位相干系数计算所有阵元的空化源强度波束合成信号和每个成像位置的空化源能量，得到空化源三维体数据，对三维体数据进行阈值化、平滑等处理，然后进行脑部聚焦超声空化的三维显示。其具体步骤如下：
(5.1)利用步骤(4.2)所得第1～N列阵元的空化源强度波束合成信号矩阵q(x,y,z,t)、步骤(4.6)所得最优加权系数以及步骤(4.9)所得相位相干系数计算所有阵元的空化源强度波束合成信号：
Q(x,y,z,t)＝wTq(x,y,z,t)SCFp(x,y,z,t)
(5.2)在一定的时间区间[T0,T0+ΔT]内对步骤(5.1)所得Q(x,y,z,t)的平方进行积分，得到每个成像位置处的空化源能量(空化源三维体数据)：
如图8所示，(a)～(c)分别为使用传统的TEA和RCB算法以及本发明所述方法得到的三维无源成像的x-z切面，结果显示传统TEA算法分辨率较差，成像伪影严重干扰了空化区域，RCB算法在一定程度上提高了分辨率并抑制了伪影，但在换能器远场仍存在较大的伪影，本发明则能在提高分辨率同时对远场的伪影进行大幅度抑制。(d)和(e)分别为x、z方向上的分辨率曲线(虚线为TEA，点虚线为RCB，实线为本发明方法)，(f)为伪影抑制程度的定量化(1～3分别代表TEA，RCB和本发明方法)，定量结果显示本发明所述方法具有良好的分辨性能，使得空化被局限在一个小区域内，从而实现空化在三维空间的精准定位。
如图9所示，三维无源成像的显示流程为：
(5.3)从步骤(5.2)获取空化源三维体数据之后，设置阈值为所有体数据最大值的0.1～0.2倍，将体数据值小于此阈值的置为零；
(5.4)选择滤波器为Savitzky-Golay滤波器，设置Savitzky-Golay滤波器的窗宽、多项式模型的阶数以及平滑点数对步骤(5.3)所得阈值化后的体数据进行平滑处理；
(5.5)计算步骤(5.4)所得平滑后的体数据的等值面的表面和顶点，提取出所有体素的等值面；
(5.6)根据步骤(5.5)所得等值面构造三维曲面，计算等值面的顶点法线，最后设置颜色、光照、视角等进行三维显示。
如图10所示，(a)～(d)分别为对不同空化微泡数目(20,50,100,200)的三维无源成像显示结果，结果说明本发明所述方法得到的三维无源成像具有足够高的精度。
本发明以相位相干系数修正的鲁棒的Capon波束合成算法为核心，利用超声面阵换能器接收脑部聚焦超声治疗过程中的空化信号，一方面可以弥补磁共振监控和传统主动超声成像方法的缺陷，另一方面可以解决现有无源成像方法在成像分辨率受限、成像伪影严重和三维空间信息丢失等方面的难题；通过改变聚焦超声治疗参数可以得到任意治疗模式下的空化图像，改变信号采集触发参数可以得到任意帧频的时间序列空化图像，二者结合起来可反应治疗程度和治疗效果，为临床提供了一种脑部聚焦超声治疗过程中空化效应的实时监控手段，使得脑部治疗的实时反馈与控制成为可能。
本发明具有以下优点：
(1)针对磁共振成像与传统主动超声成像无法与脑部聚焦超声治疗过程同步的问题，通过设置面阵阵元变迹函数使其工作在不发射只接收模式，从而实现了脑部聚焦超声治疗中空化效应的实时监控，为脑部聚焦超声治疗的临床应用提供了一种实时反馈的手段。
(2)利用颅骨CT扫描数据对面阵换能器采集到的空化信号进行校准，在此基础上对空化信号做沿着面阵两个方向的鲁棒的Capon波束合成，可以克服颅骨遮挡引起的信号失真问题，同时可以减弱来自于其他方向的干扰信号和空化微泡间的相互作用，使得三维无源成像的伪影被大幅度抑制。
(3)根据面阵阵元采集到的空化信号之间的相位差异引入了相位相干系数对波束合成算法进行修正，在抑制了成像伪影的同时可以提高成像分辨率，使得对空化的三维空间定位更加精准。
(4)高分辨及低伪影的三维无源成像方法使得对临床治疗的实时反馈更加直观，更有利于对脑部治疗过程的控制；提出的三维无源成像方法不仅可以用于脑部聚焦超声治疗的实时监控，也可用于其他部位(如肝脏、肾脏等)和其他治疗方式(如激光、微波等)的实时监控。

Claims (10)

1.一种用于脑部聚焦超声空化实时监控的三维无源成像方法，其特征在于：包括以下步骤：
步骤一：计算超声面阵换能器每个阵元方向上超声波的有效声速；
步骤二：利用拉格朗日插值法对超声面阵换能器接收的聚焦超声换能器焦点区域的空化信号进行插值，然后根据步骤一所得有效声速对空化信号进行延时；
步骤三：构造超声面阵换能器每列阵元经插值及延时后的空化信号的协方差矩阵，然后基于鲁棒的Capon波束合成方法计算最优加权系数；对所述每列阵元经插值及延时后的空化信号进行二值化并计算位相干系数；根据该最优加权系数以及相位相干系数，计算所述每列阵元的空化源强度波束合成信号；
步骤四：根据步骤三计算的超声面阵换能器各列阵元的空化源强度波束合成信号构造协方差矩阵，然后基于鲁棒的Capon波束合成方法计算最优加权系数；对所述各列阵元的空化源强度波束合成信号进行二值化并计算相位相干系数；
步骤五：根据所述各列阵元的空化源强度波束合成信号以及步骤四所得最优加权系数和相位相干系数计算超声面阵换能器所有阵元的空化源强度波束合成信号和每个成像位置的空化源能量，得到空化源三维体数据，然后进行脑部聚焦超声空化的三维显示。
2.根据权利要求1所述的方法，其特征在于：所述步骤一中，每个阵元方向上超声波的有效声速按照以下公式计算：
<mrow> <msub> <mi>c</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>K</mi> </munderover> <mfrac> <msub> <mi>D</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <msub> <mi>L</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mfrac> <msub> <mi>c</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mo>+</mo> <mfrac> <msub> <mi>D</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>s</mi> </mrow> </msub> <msub> <mi>L</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> </mfrac> <msub> <mi>c</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>s</mi> </mrow> </msub> </mrow>
其中K为除颅骨以外不同介质的数目，Li,j为空化信号到达超声面阵换能器阵元(i,j)的距离，Di,j,k和ci,j,k分别为空化信号到达阵元(i,j)的声传播路径上除颅骨以外不同介质的厚度和声速，Di,j,s为空化信号到达阵元(i,j)的声传播路径上颅骨对应位置的厚度，ci,j,s为空化信号到达阵元(i,j)的声传播路径上颅骨对应位置的声速。
3.根据权利要求2所述的方法，其特征在于：在颅骨孔隙度经验模型的基础上计算ci,j,s：
ci,j,s＝c1Φi,j,s+c2(1-Φi,j,s)
其中Φi,j,s为颅骨对应位置的孔隙度，c1为超声波在水中的传播速度，c2为超声波在颅骨中的最大传播速度。
4.根据权利要求1所述的方法，其特征在于：所述步骤二具体包括以下步骤：
(2.1)设置超声面阵换能器的阵元变迹，使超声面阵换能器工作在不发射脉冲信号只接收信号的模式，同步触发聚焦超声换能器和超声面阵换能器，由超声面阵换能器采集聚焦超声换能器作用过程中的焦点处空化的三维射频信号pi,j，其中i＝1,2,...,M；j＝1,2,...,N；M和N分别为超声面阵换能器y和x方向的阵元数目；
(2.2)根据以下拉格朗日插值多项式对每个阵元采集到的射频信号进行插值：
<mrow> <msubsup> <mi>p</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> <mi>L</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>P</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>p</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msub> <mi>k</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>l</mi> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow>
其中，为插值后的空化信号，ki,j(tl,t)为拉格朗日基本多项式，pi,j(tl)为空化信号采样点，P为采样点个数；
(2.3)根据三维成像区域和成像精度，对步骤(2.2)中所得信号进行延时处理，得到延时后的信号si,j(x,y,z,t)：
<mrow> <msub> <mi>s</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>p</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> <mi>L</mi> </msubsup> <mo>&amp;lsqb;</mo> <mi>t</mi> <mo>+</mo> <msub> <mi>d</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>c</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mo>&amp;rsqb;</mo> </mrow>
其中di,j(x,y,z)为成像位置(x,y,z)到超声面阵换能器阵元位置的空间距离，ci,j为阵元方向上超声波的有效声速。
5.根据权利要求1所述的方法，其特征在于：所述步骤三具体包括以下步骤：
(3.1)在一定的时间区间[T0,T0+ΔT]内构造协方差矩阵Rj：
<mrow> <msub> <mi>R</mi> <mi>j</mi> </msub> <mo>=</mo> <msubsup> <mo>&amp;Integral;</mo> <msub> <mi>T</mi> <mn>0</mn> </msub> <mrow> <msub> <mi>T</mi> <mn>0</mn> </msub> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>T</mi> </mrow> </msubsup> <msub> <mi>s</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msub> <mi>s</mi> <mi>j</mi> </msub> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>d</mi> <mi>t</mi> </mrow>
其中sj(x,y,z,t)为第j列M个阵元经插值及延时后的空化信号矩阵，[·]T代表矩阵的转置；
(3.2)对步骤(3.1)所得协方差矩阵做特征值分解：
<mrow> <msub> <mi>R</mi> <mi>j</mi> </msub> <mo>=</mo> <msub> <mi>U</mi> <mi>j</mi> </msub> <msub> <mi>V</mi> <mi>j</mi> </msub> <msubsup> <mi>U</mi> <mi>j</mi> <mi>T</mi> </msubsup> </mrow>
其中Vj为特征值矩阵，Uj中第i列为Vj中对角线上第i个元素对应的特征向量，i＝1,2,...,M；
(3.3)根据步骤(3.2)中所得特征值矩阵Vj和特征向量矩阵Uj以及拉格朗日迭代法计算每个成像位置的导向矢量
<mrow> <msub> <mover> <mi>a</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> <mo>=</mo> <mover> <mi>a</mi> <mo>&amp;OverBar;</mo> </mover> <mo>-</mo> <msub> <mi>U</mi> <mi>j</mi> </msub> <msup> <mrow> <mo>(</mo> <mi>I</mi> <mo>+</mo> <msub> <mi>&amp;lambda;V</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msubsup> <mi>U</mi> <mi>j</mi> <mi>T</mi> </msubsup> <mover> <mi>a</mi> <mo>&amp;OverBar;</mo> </mover> </mrow>
其中为M行1列预设导向矢量，λ为拉格朗日乘数，I为M行M列单位矩阵；
(3.4)基于鲁棒的Capon波束合成方法并利用步骤(3.1)所得协方差矩阵Rj和步骤(3.3)所得导向矢量计算最优加权系数wj：
<mrow> <msub> <mi>w</mi> <mi>j</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msub> <mover> <mi>a</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> <mo>|</mo> <mo>|</mo> </mrow> <mi>M</mi> </mfrac> <mfrac> <mrow> <msubsup> <mi>R</mi> <mi>j</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <msub> <mover> <mi>a</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> </mrow> <mrow> <msubsup> <mover> <mi>a</mi> <mo>^</mo> </mover> <mi>j</mi> <mi>T</mi> </msubsup> <msubsup> <mi>R</mi> <mi>j</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <msub> <mover> <mi>a</mi> <mo>^</mo> </mover> <mi>j</mi> </msub> </mrow> </mfrac> </mrow>
其中||·||代表欧几里得范数，[·]-1代表矩阵的逆；
(3.5)将第j列阵元经插值及延时后的空化信号si,j(x,y,z,t)中大于零的点置为1，小于等于零的点置为-1，得到二值化信号bi,j(x,y,z,t)，然后计算相位相干系数
<mrow> <msubsup> <mi>SCF</mi> <mi>j</mi> <mi>p</mi> </msubsup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mn>1</mn> <mo>-</mo> <msqrt> <mrow> <mn>1</mn> <mo>-</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <mfrac> <mn>1</mn> <mi>M</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>b</mi> <mrow> <mi>i</mi> <mo>,</mo> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <msup> <mo>|</mo> <mi>p</mi> </msup> </mrow>
其中p为可调参数；
(3.6)利用步骤(3.4)所得最优加权系数以及步骤(3.5)所得相位相干系数计算第j列阵元的空化源强度波束合成信号：
<mrow> <msub> <mi>q</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>w</mi> <mi>j</mi> <mi>T</mi> </msubsup> <msub> <mi>s</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>SCF</mi> <mi>j</mi> <mi>p</mi> </msubsup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>.</mo> </mrow>
6.根据权利要求1所述的方法，其特征在于：所述步骤四具体包括以下步骤：
(4.1)计算超声面阵换能器各列阵元的空化源强度波束合成信号矩阵：
q(x,y,z,t)＝[q1(x,y,z,t)；q2(x,y,z,t)；...；qN(x,y,z,t)]
其中qj(x,y,z,t)为第j列阵元的空化源强度波束合成信号，j＝1,2,...,N；
(4.2)在一定的时间区间[T0,T0+ΔT]内，根据步骤(4.1)所得空化源强度波束合成信号矩阵q(x,y,z,t)计算协方差矩阵R：
<mrow> <mi>R</mi> <mo>=</mo> <msubsup> <mo>&amp;Integral;</mo> <msub> <mi>T</mi> <mn>0</mn> </msub> <mrow> <msub> <mi>T</mi> <mn>0</mn> </msub> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>T</mi> </mrow> </msubsup> <mi>q</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>q</mi> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mi>d</mi> <mi>t</mi> </mrow>
其中[·]T代表矩阵的转置；
(4.3)对步骤(4.2)所得协方差矩阵做特征值分解：
R＝UVUT
其中V为特征值矩阵，U中第j列为V中对角线上第j个元素对应的特征向量，j＝1,2,...,N；
(4.4)根据步骤(4.3)中所得特征值矩阵V和特征向量矩阵U以及拉格朗日迭代法计算每个成像位置的导向矢量
<mrow> <mover> <mi>a</mi> <mo>^</mo> </mover> <mo>=</mo> <mover> <mi>a</mi> <mo>&amp;OverBar;</mo> </mover> <mo>-</mo> <mi>U</mi> <msup> <mrow> <mo>(</mo> <mi>I</mi> <mo>+</mo> <mi>&amp;lambda;</mi> <mi>V</mi> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>U</mi> <mi>T</mi> </msup> <mover> <mi>a</mi> <mo>&amp;OverBar;</mo> </mover> </mrow>
其中为N行1列预设导向矢量，λ为拉格朗日乘数，I为N行N列单位矩阵；
(4.5)基于鲁棒的Capon波束合成方法并利用步骤(4.2)所得协方差矩阵R和步骤(4.4)所得导向矢量计算最优加权系数w：
<mrow> <mi>w</mi> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <mover> <mi>a</mi> <mo>^</mo> </mover> <mo>|</mo> <mo>|</mo> </mrow> <mi>N</mi> </mfrac> <mfrac> <mrow> <msup> <mi>R</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mover> <mi>a</mi> <mo>^</mo> </mover> </mrow> <mrow> <msup> <mover> <mi>a</mi> <mo>^</mo> </mover> <mi>T</mi> </msup> <msup> <mi>R</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mover> <mi>a</mi> <mo>^</mo> </mover> </mrow> </mfrac> </mrow>
其中||·||代表欧几里得范数，[·]-1代表矩阵的逆；
(4.6)将每列阵元的空化源强度波束合成信号qj(x,y,z,t)中大于零的点置为1，小于等于零的点置为-1，得到二值化信号bj(x,y,z,t)，然后计算相位相干系数SCFp(x,y,z,t)：
<mrow> <msup> <mi>SCF</mi> <mi>p</mi> </msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mn>1</mn> <mo>-</mo> <msqrt> <mrow> <mn>1</mn> <mo>-</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mi>b</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <msup> <mo>|</mo> <mi>p</mi> </msup> </mrow>
其中p为可调参数。
7.根据权利要求1所述的方法，其特征在于：所述步骤五中，所有阵元的空化源强度波束合成信号Q(x,y,z,t)按照以下公式计算：
Q(x,y,z,t)＝wTq(x,y,z,t)SCFp(x,y,z,t)
其中，[·]T代表矩阵的转置；q(x,y,z,t)为超声面阵换能器各列阵元的空化源强度波束合成信号矩阵；w为步骤四中得到的最优加权系数；SCFp(x,y,z,t)为步骤四中得到的相位相干系数；
在一定的时间区间[T0,T0+ΔT]内对Q(x,y,z,t)的平方进行积分，得到每个成像位置处的空化源能量I(x,y,z)：
<mrow> <mi>I</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mi>&amp;Delta;</mi> <mi>T</mi> </mrow> </mfrac> <msubsup> <mo>&amp;Integral;</mo> <msub> <mi>T</mi> <mn>0</mn> </msub> <mrow> <msub> <mi>T</mi> <mn>0</mn> </msub> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>T</mi> </mrow> </msubsup> <mi>Q</mi> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>,</mo> <mi>z</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mi>d</mi> <mi>t</mi> <mo>.</mo> </mrow>
8.根据权利要求1所述的方法，其特征在于：所述步骤五中，三维显示的具体步骤为：
(5.1)对所得的三维体数据进行阈值化处理，阈值取体数据最大值的0.1～0.2倍；
(5.2)对步骤(5.1)所得阈值化后的三维体数据进行平滑处理，然后提取出体素等值面；
(5.3)根据步骤(5.2)所得等值面构造三维曲面，然后设置颜色、光照及视角后进行三维显示。
9.一种用于脑部聚焦超声空化实时监控的三维无源成像系统，其特征在于：包括声速计算模块、空化信号校准模块、自适应波束合成模块以及三维显示模块；
所述声速计算模块计算超声面阵换能器每个阵元方向上超声波的有效声速；
所述空化信号校准模块对超声面阵换能器接收的聚焦超声换能器焦点区域的空化信号进行插值以及延时处理；
所述自适应波束合成模块对超声面阵换能器两个方向的阵元做鲁棒的Capon波束合成，并基于空化信号相位间的差异引入相位相干系数对波束合成算法进行修正；
所述三维显示模块利用自适应波束合成模块得到的空化源三维体数据进行脑部聚焦超声空化的三维显示。
10.根据权利要求9所述的系统，其特征在于：所述自适应波束合成具体包括以下步骤：
构造超声面阵换能器每列阵元经插值及延时后的空化信号的协方差矩阵，然后基于鲁棒的Capon波束合成方法计算最优加权系数；对所述每列阵元经插值及延时后的空化信号进行二值化并计算相位相干系数；根据该最优加权系数以及相位相干系数，计算所述每列阵元的空化源强度波束合成信号；
根据超声面阵换能器各列阵元的空化源强度波束合成信号构造协方差矩阵，然后基于鲁棒的Capon波束合成方法计算最优加权系数；对各列阵元的空化源强度波束合成信号进行二值化并计算相位相干系数；然后计算所有阵元的空化源强度波束合成信号。
CN201710582424.7A 2017-07-17 2017-07-17 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统 Active CN107260217B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710582424.7A CN107260217B (zh) 2017-07-17 2017-07-17 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710582424.7A CN107260217B (zh) 2017-07-17 2017-07-17 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统

Publications (2)

Publication Number Publication Date
CN107260217A true CN107260217A (zh) 2017-10-20
CN107260217B CN107260217B (zh) 2018-07-17

Family

ID=60072897

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710582424.7A Active CN107260217B (zh) 2017-07-17 2017-07-17 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统

Country Status (1)

Country Link
CN (1) CN107260217B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108078590A (zh) * 2018-01-03 2018-05-29 声泰特（成都）科技有限公司 基于超声频谱多普勒的血流动力学可视化方法与系统
CN108175439A (zh) * 2017-12-08 2018-06-19 西安交通大学 低频经颅专用超声换能器及其应用
CN109171811A (zh) * 2018-09-25 2019-01-11 西安交通大学 基于特征空间自适应波束合成的频域被动空化成像及频率复合成像方法
CN109431536A (zh) * 2018-09-17 2019-03-08 西安交通大学 一种聚焦超声空化的实时高分辨时空分布成像方法与系统
CN109513123A (zh) * 2018-12-28 2019-03-26 西安交通大学 一种基于半球阵的高分辨三维被动空化成像方法
CN109938771A (zh) * 2019-03-30 2019-06-28 河南省省立医院有限公司 三维盆底超声图像处理控制系统
CN110507355A (zh) * 2019-09-20 2019-11-29 青岛海信医疗设备股份有限公司 一种超声成像系统、方法、设备及介质
CN111265245A (zh) * 2020-01-20 2020-06-12 西安交通大学 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
CN111466949A (zh) * 2020-04-13 2020-07-31 剑桥大学南京科技创新中心有限公司 一种mmse波束形成器、mmse波束形成方法、计算机可读存储介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1215616A (zh) * 1998-01-25 1999-05-05 重庆医科大学附属第二医院 高强度聚焦超声肿瘤扫描治疗系统
US20050283079A1 (en) * 2004-06-22 2005-12-22 Steen Erik N Method and apparatus for medical ultrasound navigation user interface
US20100317971A1 (en) * 2009-05-04 2010-12-16 Siemens Medical Solutions Usa, Inc. Feedback in medical ultrasound imaging for high intensity focused ultrasound
CN104535645A (zh) * 2014-12-27 2015-04-22 西安交通大学 微秒分辨空化时空分布的三维空化定量成像方法
CN104777485A (zh) * 2015-04-20 2015-07-15 西安交通大学 超声二维面阵的三维宽波束小区域快速空化成像方法
CN104887266A (zh) * 2015-05-29 2015-09-09 西安交通大学 基于面阵的小区域三维被动空化成像及三维复合成像方法
CN105188559A (zh) * 2013-02-28 2015-12-23 爱飞纽医疗机械贸易有限公司 检测空化的方法及超声医疗设备
CN103235041B (zh) * 2013-04-26 2016-03-02 西安交通大学 基于超声主动空化成像的空化起始阈值分布重建方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1215616A (zh) * 1998-01-25 1999-05-05 重庆医科大学附属第二医院 高强度聚焦超声肿瘤扫描治疗系统
US20050283079A1 (en) * 2004-06-22 2005-12-22 Steen Erik N Method and apparatus for medical ultrasound navigation user interface
US20100317971A1 (en) * 2009-05-04 2010-12-16 Siemens Medical Solutions Usa, Inc. Feedback in medical ultrasound imaging for high intensity focused ultrasound
CN105188559A (zh) * 2013-02-28 2015-12-23 爱飞纽医疗机械贸易有限公司 检测空化的方法及超声医疗设备
CN103235041B (zh) * 2013-04-26 2016-03-02 西安交通大学 基于超声主动空化成像的空化起始阈值分布重建方法
CN104535645A (zh) * 2014-12-27 2015-04-22 西安交通大学 微秒分辨空化时空分布的三维空化定量成像方法
CN104777485A (zh) * 2015-04-20 2015-07-15 西安交通大学 超声二维面阵的三维宽波束小区域快速空化成像方法
CN104887266A (zh) * 2015-05-29 2015-09-09 西安交通大学 基于面阵的小区域三维被动空化成像及三维复合成像方法

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108175439A (zh) * 2017-12-08 2018-06-19 西安交通大学 低频经颅专用超声换能器及其应用
CN108078590A (zh) * 2018-01-03 2018-05-29 声泰特（成都）科技有限公司 基于超声频谱多普勒的血流动力学可视化方法与系统
CN109431536A (zh) * 2018-09-17 2019-03-08 西安交通大学 一种聚焦超声空化的实时高分辨时空分布成像方法与系统
CN109431536B (zh) * 2018-09-17 2019-08-23 西安交通大学 一种聚焦超声空化的实时高分辨时空分布成像方法与系统
CN109171811A (zh) * 2018-09-25 2019-01-11 西安交通大学 基于特征空间自适应波束合成的频域被动空化成像及频率复合成像方法
CN109171811B (zh) * 2018-09-25 2019-08-23 西安交通大学 基于特征空间自适应波束合成的频域被动空化成像及频率复合成像方法
CN109513123A (zh) * 2018-12-28 2019-03-26 西安交通大学 一种基于半球阵的高分辨三维被动空化成像方法
CN109938771A (zh) * 2019-03-30 2019-06-28 河南省省立医院有限公司 三维盆底超声图像处理控制系统
CN110507355A (zh) * 2019-09-20 2019-11-29 青岛海信医疗设备股份有限公司 一种超声成像系统、方法、设备及介质
CN111265245A (zh) * 2020-01-20 2020-06-12 西安交通大学 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
CN111265245B (zh) * 2020-01-20 2021-02-09 西安交通大学 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
CN111466949A (zh) * 2020-04-13 2020-07-31 剑桥大学南京科技创新中心有限公司 一种mmse波束形成器、mmse波束形成方法、计算机可读存储介质

Also Published As

Publication number Publication date
CN107260217B (zh) 2018-07-17

Similar Documents

Publication Publication Date Title
CN107260217B (zh) 用于脑部聚焦超声空化实时监控的三维无源成像方法及系统
US9730676B2 (en) Ultrasound imaging system using beamforming techniques for phase coherence grating lobe suppression
JP6030207B2 (ja) 超音波合成イメージングの装置と方法
CN109431536B (zh) 一种聚焦超声空化的实时高分辨时空分布成像方法与系统
US8734352B2 (en) Spatially-fine shear wave dispersion ultrasound vibrometry sampling
EP2849651B1 (en) System and method for measurement of shear wave speed from multi-directional wave fields
CN104887266B (zh) 基于面阵的小区域三维被动空化成像及三维复合成像方法
CN111134719B (zh) 一种聚焦超声辐照相变纳米液滴的主被动超声复合成像方法及系统
US8818064B2 (en) Time-domain estimator for image reconstruction
CN110248606A (zh) 气穴定位
US9366753B2 (en) Systems and methods for ultrasound retrospective transmit focus beamforming
JP2014506523A (ja) 非合焦超音波による超音波振動法
CN103356189A (zh) 磁共振和超声参数化图像融合
CN108836389B (zh) 平面波相关点相干自适应波束合成成像方法
WO2016101382A1 (zh) 微秒分辨空化时空分布的三维空化定量成像方法
KR101610874B1 (ko) 공간 일관성 기초 초음파 신호 처리 모듈 및 그에 의한 초음파 신호 처리 방법
CN111265245B (zh) 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
CN109513123B (zh) 一种基于半球阵的高分辨三维被动空化成像方法
Li et al. Time-resolved passive cavitation mapping using the transient angular spectrum approach
Hazard et al. Effects of motion on a synthetic aperture beamformer for real-time 3D ultrasound
JP2003180688A (ja) 幅広ビーム映像化
US11432804B2 (en) Methods and systems for processing an unltrasound image
Nordenfur Comparison of pushing sequences for shear wave elastography
CN112764040B (zh) 一种基于射线理论相位修正的合成孔径波束形成方法
Afrakhteh et al. 2d/3d echocardiography frame rate enhancement by means of a novel spatio-temporal reconstruction technique

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