CN116115260A - 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统 - Google Patents

高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统 Download PDF

Info

Publication number
CN116115260A
CN116115260A CN202310037068.6A CN202310037068A CN116115260A CN 116115260 A CN116115260 A CN 116115260A CN 202310037068 A CN202310037068 A CN 202310037068A CN 116115260 A CN116115260 A CN 116115260A
Authority
CN
China
Prior art keywords
transducer
signals
ultrasonic
signal
imaging
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
CN202310037068.6A
Other languages
English (en)
Other versions
CN116115260B (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
Publication of CN116115260A publication Critical patent/CN116115260A/zh
Application granted granted Critical
Publication of CN116115260B publication Critical patent/CN116115260B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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/44—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device
    • 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
    • A61B8/4488—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device characterised by features of the ultrasound transducer the transducer being a phased array
    • 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
    • A61B8/4494—Constructional features of the ultrasonic, sonic or infrasonic diagnostic device characterised by features of the ultrasound transducer characterised by the arrangement of the transducer elements
    • 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
    • A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • 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
    • A61B8/5215—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N7/00—Ultrasound therapy
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04—Analysing solids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N7/00—Ultrasound therapy
    • A61N2007/0039—Ultrasound therapy using microbubbles

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Animal Behavior & Ethology (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Public Health (AREA)
  • Biophysics (AREA)
  • Surgery (AREA)
  • Molecular Biology (AREA)
  • Medical Informatics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Gynecology & Obstetrics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Software Systems (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Acoustics & Sound (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明公开了一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统：同步聚焦超声换能器发射和超声成像换能器接收的时序并采集被动空化射频信号，根据延时叠加波束合成图获得聚焦超声换能器顶点定位坐标以准确地计算绝对渡越时间，根据实际中的参数和换能器空间设置构造被动阵列仿真信号，通过对根据被动阵列仿真信号所得的自适应权重进行主成分分析来构造维度转换矩阵，依据维度转换矩阵对被动空化射频信号进行快速特征空间自适应波束合成，得到发收时序同步超声被动空化成像结果。本发明通过对空化图像的高分辨率、高对比度和快速的重建，用于多种超声空化应用中空化时空动力学演变过程的监控。

Description

高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统
技术领域
本发明属于超声检测与超声成像技术领域，涉及一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统。
背景技术
超声波作用于介质产生的空化效应是组织消融、组织毁损、血栓溶解、血脑屏障开放等多种治疗超声应用的物理基础，空化的检测对于澄清空化物理机制和调控优化超声发射参数有重要意义。超声成像具有成像帧率高、可进行深部成像、成本低等优势，在空化检测方面发挥着重要作用，其可分为基于发射接收模式的主动成像和基于不发射只接收模式的被动成像(即超声被动空化成像)。超声被动空化成像避免了传统主动成像存在的聚焦超声换能器发射的信号干扰成像和超声成像换能器发射的信号干扰空化活动的问题，可在聚焦超声作用的过程中对空化活动进行实时监控且不会对空化活动造成干扰，在多种治疗超声应用中具有广泛的应用前景。
超声被动空化成像最初采用基于相对渡越时间信息的延时叠加积分方法来重建空化图像，然而在使用线阵换能器或相控阵换能器等小孔径超声换能器进行成像时，该方法会产生过大的轴向波束并伴随着严重的干扰伪迹，从而导致难以在轴向区分多个空化源。为解决此问题，研究人员先后发展了基于鲁棒波束合成、相位相干系数加权、频率耦合波束合成、延时相乘叠加波束合成等的超声被动空化成像方法，研究结果表明这些方法抑制干扰伪迹并改善成像分辨率；然而，由于小孔径的超声成像换能器的衍射模式的限制，成像结果的轴向分辨率仍然受到极大的限制。
发收时序同步超声被动空化成像技术是一种将聚焦超声换能器的发射时序与超声成像换能器的接收时序严格同步的成像技术，能够如超声主动成像一样通过绝对渡越时间信息对接收信号进行处理，从而使得成像结果的轴向分辨率不再受小孔径超声换能器衍射模式的限制，而是取决于聚焦超声换能器发射脉冲的长度，在短脉冲发射情境下获得高轴向分辨率。发收时序同步超声被动空化成像技术已用于检测生物组织中的空化成核、监控血脑屏障开放中的空化活动以及定征超声换能器的声场分布等方面，该技术也与脉冲逆转技术相结合来增强血脑屏障开放过程中弱空化信号的检测灵敏度。波束合成处理是实施发收时序同步超声被动空化成像技术的重要环节，目前主要采用的是传统的延时叠加波束合成方法；然而，该方法是非自适应的，其在波束合成过程中使用与接收信号无关的预设固定权重，会产生大的主瓣尺寸和高的旁瓣水平，导致成像结果的分辨率和对比度较差，从而不利于对空化活动的监控。
在发收时序同步超声被动空化成像中需要解决的难点问题是同时提高成像结果的分辨率和对比度以及以更快的计算速度达成此目的，这也是多种超声空化应用的重要需求。鉴于此，目前亟待提出一种高分辨高对比快速计算的发收时序同步超声被动空化成像技术。
发明内容
本发明的目的在于提供一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统。
为了实现上述目的，本发明采用了以下技术方案：
一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，包括以下步骤：
1)利用多通道波形发生器同步地触发脉冲发射/接收器和超声成像数据采集系统以对聚焦超声换能器发射和超声成像换能器接收的时序进行同步，利用超声成像数据采集系统采集由超声成像换能器被动接收的空化声辐射信号(即在对聚焦超声换能器发射和超声成像换能器接收的时序进行同步的情况下，对由超声成像换能器被动接收的空化声辐射信号进行采集)，得到被动空化射频信号并存储；
2)在成像坐标系下设置像素网格，计算聚焦超声换能器顶点坐标的遍历集合中每一顶点坐标对应的像素网格中各像素点处的绝对渡越时间，根据该绝对渡越时间对从步骤1采集得到的被动空化射频信号中选取的多帧被动空化射频信号分别进行延时叠加波束合成处理，得到所述遍历集合中每一顶点坐标对应的各像素点处的延时叠加波束合成信号，对所述遍历集合中每一顶点坐标对应的各像素点处的延时叠加波束合成信号进行希尔伯特变换，得到所述遍历集合中各顶点坐标对应的延时叠加波束合成图，根据所述遍历集合中各顶点坐标对应的延时叠加波束合成图寻找所述遍历集合中的最优顶点坐标，根据在选取的不同帧被动空化射频信号下分别寻找到的最优顶点坐标计算聚焦超声换能器顶点定位坐标；
3)设置声源辐射仿真信号和被动阵列仿真信号规模并在成像坐标系下设置仿真声源群中各声源的坐标，根据聚焦超声换能器顶点定位坐标计算仿真声源群中每一声源处的绝对渡越时间，根据该绝对渡越时间和声源辐射仿真信号构造超声成像换能器各阵元的仿真信号，根据超声成像换能器各阵元的仿真信号构造对应声源的阵列仿真信号，对各声源的阵列仿真信号进行叠加得到被动阵列仿真信号；
4)根据聚焦超声换能器顶点定位坐标计算所述像素网格中各像素点处的绝对渡越时间，根据该绝对渡越时间从被动阵列仿真信号中分别提取采样点，得到对应像素点处的全孔径采样信号，将对应像素点处的全孔径采样信号分割成多个重叠的子孔径采样信号并利用子孔径采样信号构造子孔径采样组合信号，计算子孔径采样组合信号的协方差矩阵并对该协方差矩阵进行对角加载，根据对角加载后的协方差矩阵计算对应像素点处的自适应权重，对各像素点处的自适应权重进行主成分分析，得到主成分向量，根据主成分向量构造维度转换矩阵；
5)根据步骤4中计算的绝对渡越时间从步骤1采集得到的任意一帧被动空化射频信号中提取并分割得到对应像素点处的多个重叠的子孔径采样信号，利用由该子孔径采样信号构造的子孔径采样组合信号计算相应的协方差矩阵并对该协方差矩阵进行对角加载，利用所述维度转换矩阵对对角加载后的协方差矩阵进行降维处理，得到对应像素点处的低维域协方差矩阵，根据低维域协方差矩阵计算对应像素点处的低维域自适应权重并将该低维域自适应权重投影到所述低维域协方差矩阵经特征分解所得的信号子空间上，得到对应像素点处的低维域特征空间自适应权重，根据所述维度转换矩阵和子孔径采样信号计算对应像素点处的低维域子孔径采样平均信号，根据所述低维域特征空间自适应权重和低维域子孔径采样平均信号计算对应像素点处的快速特征空间自适应波束合成信号，对各像素点处的快速特征空间自适应波束合成信号进行希尔伯特变换，得到快速特征空间自适应波束合成图，对快速特征空间自适应波束合成图进行归一化及对数化处理，得到对应帧被动空化射频信号的发收时序同步超声被动空化成像结果。
优选的，所述步骤1中，超声成像换能器为线阵换能器等诊断用超声换能器中的一种或多种；被动空化射频信号是使超声成像换能器工作在不发射只接收模式(即被动接收模式)、通过多通道波形发生器同步聚焦超声换能器发射和超声成像换能器接收的时序，并由超声成像数据采集系统采集而得到的。
优选的，所述步骤2中，所述遍历集合中聚焦超声换能器顶点坐标的x轴坐标、y轴坐标和z轴坐标的遍历区间的长度≥10mm；x轴坐标和y轴坐标的遍历步长为像素网格中x轴方向的像素间隔，z轴坐标的遍历步长为像素网格中z轴方向的像素间隔。
优选的，所述步骤2中，寻找所述遍历集合中的最优顶点坐标具体包括以下步骤：
2.1)统计所述遍历集合中每一顶点坐标对应的延时叠加波束合成图中像素值大于像素值阈值的像素点的数目，其中像素值阈值为该延时叠加波束合成图中最大像素值的0.1～0.5倍；
2.2)重复步骤2.1，得到所述遍历集合中各顶点坐标对应的像素值大于像素值阈值的像素点数目，寻找这些像素点数目中的最小像素点数目，该最小像素点数目对应的顶点坐标为最优顶点坐标。
优选的，所述步骤2中，聚焦超声换能器顶点定位坐标为随机选取的多帧(≥10)被动空化射频信号对应的最优顶点坐标的平均值。
优选的，所述步骤3中，声源辐射仿真信号为加Hanning窗的正弦信号，声源辐射仿真信号的时间长度与聚焦超声换能器发射脉冲信号的脉冲长度一致，声源辐射仿真信号的频率由被动空化射频信号的频谱分布的质心频率确定，声源辐射仿真信号的采样频率与被动空化射频信号的采样频率一致；被动阵列仿真信号的规模与被动空化射频信号的规模一致。
优选的，所述步骤3中，计算仿真声源群中每一声源处的绝对渡越时间所用的换能器空间设置与步骤1中聚焦超声换能器和超声成像换能器的空间设置相同。
优选的，所述步骤3中，构造的任意一声源的阵列仿真信号中超声成像换能器任意一阵元的仿真信号表示为：
Figure BDA0004049205960000041
其中，
Figure BDA0004049205960000044
为仿真声源群中第m个声源的阵列仿真信号中第i个阵元的仿真信号，i＝1,2,...,NE，NE为超声成像换能器的阵元数目，k＝1,2,...,NS，NS为被动空化射频信号的采样点数，m＝1,2,...,NA，NA为仿真声源群中声源的数目，round{·}表示四舍五入取整，
Figure BDA0004049205960000042
为根据聚焦超声换能器顶点定位坐标计算的第m个声源处的绝对渡越时间，τTrig为超声成像数据采集系统采集空化声辐射信号的触发延时，fs为被动空化射频信号的采样频率，sae为声源辐射仿真信号，
Figure BDA0004049205960000043
为声源辐射仿真信号的采样点数。
优选的，所述步骤4中，计算所述像素网格中各像素点处的绝对渡越时间所用的换能器空间设置与步骤1中聚焦超声换能器和超声成像换能器的空间设置相同。
优选的，所述步骤4中，对角加载后的协方差矩阵表示为：
Figure BDA0004049205960000051
Figure BDA0004049205960000052
Figure BDA0004049205960000053
Figure BDA0004049205960000054
其中，
Figure BDA00040492059600000510
和
Figure BDA0004049205960000055
分别为从被动阵列仿真信号中得到的第j个像素点处的全孔径采样信号中第l个、第l+1个和第l+L-1个阵元的信号，
Figure BDA0004049205960000056
和
Figure BDA0004049205960000057
分别为被动阵列仿真信号对应的第j个像素点处的子孔径采样信号、子孔径采样组合信号、子孔径采样组合信号的协方差矩阵和对角加载后的协方差矩阵，l＝1,2,...,NE-L+1，L为子孔径的长度，j＝1,2,...,NP，NP为像素点的数目，
Figure BDA0004049205960000058
为第j个像素点的坐标，[·]T表示转置，[·]H表示共轭转置，δ为对角加载因子，trace{·}表示求迹运算，I为单位矩阵。
优选的，所述步骤4中，子孔径的长度L为NE/2，对角加载因子δ为1/L。
优选的，所述步骤4中，对各像素点处的自适应权重进行主成分分析具体包括以下步骤：
将根据被动阵列仿真信号计算的所述像素网格中各像素点处的自适应权重作为样本，计算自适应权重的协方差矩阵，然后对该协方差矩阵进行特征分解，从而得到对应于L个降序排列的特征值的L个特征向量，即得到L个主成分向量。
优选的，所述步骤4中，维度转换矩阵的构造具体包括以下步骤：
从L个主成分向量中选择最大特征值对应的前Q个主成分向量v1,v2,...,vQ-1,vQ，并将主成分向量vQ替换为元素均为
Figure BDA0004049205960000059
的直流向量dc，所得维度转换矩阵表示为：
C＝[v1,v2,...,vQ-1,dc]。
优选的，所述步骤2、步骤3和步骤4中，绝对渡越时间的计算具体包括以下步骤：
S1以超声成像换能器的中心为原点、以超声成像换能器的横向和轴向分别为x轴和z轴，并以与超声成像换能器的成像平面垂直的方向为y轴，建立三维位置坐标系；
S2从所述三维位置坐标系中删除y轴，得到成像坐标系；
S3建立聚焦超声换能器顶点处切线或切面的方程：
对于换能器空间设置I，聚焦超声换能器顶点处切线的方程为：
sinα(x-xV)+cosα(z-zV)＝0
对于换能器空间设置II，聚焦超声换能器顶点处切面的方程为：
sinα(y-yV)+cosα(z-zV)＝0
其中，xV和zV分别为聚焦超声换能器顶点在所述三维位置坐标系中的x轴坐标和z轴坐标，yV为聚焦超声换能器顶点在所述三维位置坐标系中的y轴坐标；在换能器空间设置I中，聚焦超声换能器的中轴线在超声成像换能器的成像平面内，0≤α≤π/2；在换能器空间设置II中，聚焦超声换能器的中轴线与超声成像换能器的成像平面相交且与超声成像换能器的横向方向垂直，0＜α≤π/2；α为聚焦超声换能器的中轴线与超声成像换能器的中轴线之间的夹角；
S4建立成像坐标系中的点到聚焦超声换能器顶点处切线或切面的距离的计算公式：
对于换能器空间设置I，坐标为(x,z)的点到聚焦超声换能器顶点处切线的距离的计算公式为：
d1(x,z)＝|sinα(x-xV)+cosα(z-zV)|
其中，|·|表示取绝对值；
对于换能器空间设置II，坐标为(x,z)的点到聚焦超声换能器顶点处切面的距离的计算公式为：
d1(x,z)＝|-sinα(yV)+cosα(z-zV)|
S5建立成像坐标系中坐标为(x,z)的点到超声成像换能器各阵元的距离的计算公式：
Figure BDA0004049205960000061
其中，i＝1,2,...,NE，
Figure BDA0004049205960000062
为超声成像换能器第i个阵元的坐标；
S6根据步骤S4和S5中的距离的计算公式建立绝对渡越时间的计算公式：
τi(x,z)＝[d1(x,z)+di 2(x,z)]/c
其中，c为声传播速度。
对于步骤S4～S6中的计算公式：
令(x,z)为所述像素网格中任意一像素点的坐标
Figure BDA0004049205960000078
并令(xV,yV,zV)为聚焦超声换能器顶点坐标的遍历集合中任意一顶点坐标，则得到步骤2中任意一顶点坐标对应的任意一像素点处的绝对渡越时间；
令(x,z)为所述仿真声源群中任意一声源的坐标
Figure BDA0004049205960000079
并令(xV,yV,zV)为聚焦超声换能器顶点定位坐标(xVL,yVL,zVL)，则得到步骤3中任意一声源处的绝对渡越时间；
令(x,z)为所述像素网格中任意一像素点的坐标
Figure BDA00040492059600000710
并令(xV,yV,zV)为聚焦超声换能器顶点定位坐标(xVL,yVL,zVL)，则得到步骤4中任意一像素点处的绝对渡越时间。
优选的，所述步骤5中，低维域协方差矩阵的计算公式为：
Figure BDA0004049205960000071
其中，j＝1,2,...,NP，
Figure BDA0004049205960000072
Figure BDA0004049205960000073
和
Figure BDA0004049205960000074
分别为第p帧被动空化射频信号对应的第j个像素点处的低维域子孔径采样组合信号和子孔径采样组合信号，ξ为
Figure BDA0004049205960000075
中各元素的平方之和与子孔径采样信号的数目的比值。
优选的，所述步骤5中，低维域协方差矩阵经特征分解所得的信号子空间的维度根据互谱度量法确定，其中，互谱度量系数为1％～5％。
优选的，所述步骤5中，低维域子孔径采样平均信号的计算公式为：
Figure BDA0004049205960000076
其中，j＝1,2,...,NP，
Figure BDA0004049205960000077
为第p帧被动空化射频信号对应的第j个像素点处的第l个子孔径采样信号。
一种高分辨高对比快速计算的发收时序同步超声被动空化成像系统，该系统包括主控计算机、聚焦超声换能器、超声成像换能器和用于对聚焦超声换能器发射和超声成像换能器接收的时序进行同步的多通道波形发生器，所述主控计算机包括用于对通过采集超声成像换能器被动接收的空化声辐射信号得到的被动空化射频信号进行处理并生成发收时序同步超声被动空化成像结果的软件(所述主控计算机还用于存储所述被动空化射频信号)，所述软件包括聚焦超声换能器顶点定位坐标计算模块、被动阵列仿真信号构造模块、维度转换矩阵构造模块、快速特征空间自适应波束合成模块以及图像显示模块；
所述聚焦超声换能器顶点定位坐标计算模块用于执行上述步骤2)，主要是用于在成像坐标系下设置像素网格、计算聚焦超声换能器顶点坐标的遍历集合中每一顶点坐标对应的像素网格中各像素点处的绝对渡越时间、根据该绝对渡越时间对从采集得到的被动空化射频信号中选取的多帧被动空化射频信号分别进行延时叠加波束合成处理、对处理所得所述遍历集合中每一顶点坐标对应的各像素点处的延时叠加波束合成信号进行希尔伯特变换、根据变换所得所述遍历集合中各顶点坐标对应的延时叠加波束合成图寻找所述遍历集合中的最优顶点坐标，以及根据在选取的不同帧被动空化射频信号下分别寻找到的最优顶点坐标计算聚焦超声换能器顶点定位坐标；
所述被动阵列仿真信号构造模块用于执行上述步骤3)，主要是用于设置声源辐射仿真信号和被动阵列仿真信号规模并在成像坐标系下设置仿真声源群中各声源的坐标、根据所述聚焦超声换能器顶点定位坐标计算模块得到的聚焦超声换能器顶点定位坐标计算仿真声源群中每一声源处的绝对渡越时间、根据该绝对渡越时间和声源辐射仿真信号构造超声成像换能器各阵元的仿真信号、根据超声成像换能器各阵元的仿真信号构造对应声源的阵列仿真信号，以及通过对各声源的阵列仿真信号进行叠加得到被动阵列仿真信号；
所述维度转换矩阵构造模块用于执行上述步骤4)，主要是用于根据所述聚焦超声换能器顶点定位坐标计算模块得到的聚焦超声换能器顶点定位坐标计算所述像素网格中各像素点处的绝对渡越时间、根据该绝对渡越时间从所述被动阵列仿真信号构造模块得到的被动阵列仿真信号中分别提取采样点、将提取所得对应像素点处的全孔径采样信号分割成多个重叠的子孔径采样信号并利用子孔径采样信号构造子孔径采样组合信号、计算子孔径采样组合信号的协方差矩阵并对该协方差矩阵进行对角加载、根据对角加载后的协方差矩阵计算对应像素点处的自适应权重、对各像素点处的自适应权重进行主成分分析，以及根据分析所得主成分向量构造维度转换矩阵；
所述快速特征空间自适应波束合成模块用于执行上述步骤5)，主要是用于根据所述维度转换矩阵构造模块中计算的绝对渡越时间从采集得到的任意一帧被动空化射频信号中提取并分割得到对应像素点处的多个重叠的子孔径采样信号、利用由该子孔径采样信号构造的子孔径采样组合信号计算相应的协方差矩阵并对该协方差矩阵进行对角加载、利用所述维度转换矩阵构造模块得到的维度转换矩阵对对角加载后的协方差矩阵进行降维处理、根据降维处理所得对应像素点处的低维域协方差矩阵计算对应像素点处的低维域自适应权重并将该低维域自适应权重投影到所述低维域协方差矩阵经特征分解所得的信号子空间上、根据所述维度转换矩阵和子孔径采样信号计算对应像素点处的低维域子孔径采样平均信号、根据投影所得对应像素点处的低维域特征空间自适应权重和所述低维域子孔径采样平均信号计算对应像素点处的快速特征空间自适应波束合成信号，以及通过对各像素点处的快速特征空间自适应波束合成信号进行希尔伯特变换得到快速特征空间自适应波束合成图；
所述图像显示模块用于对所述快速特征空间自适应波束合成模块得到的快速特征空间自适应波束合成图进行归一化及对数化处理，对处理所得任意一帧被动空化射频信号对应的发收时序同步超声被动空化成像结果进行显示、或对多帧被动空化射频信号对应的发收时序同步超声被动空化成像结果进行动态的显示。
本发明的有益效果体现在：
本发明提出的高分辨高对比快速计算的发收时序同步超声被动空化成像方法，在同步发射接收时序并采集得到被动空化射频信号后，通过计算聚焦超声换能器顶点定位坐标、构造被动阵列仿真信号并据此构造维度转换矩阵以及进行快速特征空间自适应波束合成，有效提高发收时序同步超声被动空化成像的分辨率和对比度、同时提高成像的计算速度。本发明提出的成像方法对高强度组织毁损、低强度血脑屏障开放等多种超声空化应用中空化活动的时空动力学演变过程进行准确的实时动态监控，并为澄清这些应用背后的空化物理机制和相关的生物学机制以及调控和优化超声发射参数奠定坚实基础。
其中，本发明通过多通道波形发生器发射触发信号来同步聚焦超声换能器发射和超声成像换能器接收的时序，从而能够利用绝对渡越时间对被动空化射频信号进行处理，使得超声被动空化成像的轴向分辨率不再受限于换能器衍射模式而取决于发射脉冲的时间长度，可在短脉冲发射情境下获得高的轴向分辨率。本发明根据聚焦超声换能器顶点坐标的遍历集合中任意一顶点坐标对应的延时叠加波束合成图，从遍历集合中寻找最优顶点坐标并计算聚焦超声换能器顶点定位坐标，从而能够在三维位置坐标系中对聚焦超声换能器进行准确的定位，有助于提高绝对渡越时间的计算准确度。本发明通过数值仿真的方法构造被动阵列仿真信号并利用被动阵列仿真信号构造维度转换矩阵，然后利用维度转换矩阵对任意一帧被动空化射频信号进行处理，从而无需在处理每帧被动空化射频信号时每次构造相应的维度转换矩阵，也无需通过进行预实验的方法来构造维度转换矩阵；通过数值仿真方法得到的维度转换矩阵一经构造即可重复应用于同一情境下所产生的多帧被动空化射频信号，也可应用于不同的情境(如不同的脉冲重复频率等)，从而为数据处理提供了极大的便利性。本发明中对全孔径采样信号进行分割、利用子孔径采样信号构造子孔径采样组合信号、然后计算子孔径采样组合信号的协方差矩阵，从而对协方差矩阵进行良好的估计；同时对协方差矩阵进行对角加载，从而增强自适应波束合成的鲁棒性。本发明中采用的特征空间自适应波束合成集成了基于最小方差无失真响应的自适应波束合成方法和将自适应权重投影到特征分解所得信号子空间的方法，能够有效地减小主瓣尺寸并降低旁瓣水平，从而提高发收时序同步超声被动空化成像的分辨率和对比度。本发明通过对被动阵列仿真信号对应的自适应权重进行主成分分析来构造维度转换矩阵，并利用其对被动空化射频信号对应的协方差矩阵进行降维处理，将协方差矩阵的维度从子孔径的长度降低至一个更小的维度，使得协方差矩阵求逆和特征分解所需的巨大计算量得以减小，从而能够以更快的计算速度进行特征空间自适应波束合成，可有效地提高发收时序同步超声被动空化成像的计算速度。
进一步地，本发明中使超声成像换能器工作在不发射只接收模式来进行超声被动空化成像，避免了利用超声主动成像进行空化监控时聚焦超声换能器发射的信号会对成像造成干扰和超声成像换能器发射的信号会对空化活动造成干扰的问题，从而在聚焦超声作用过程中对空化活动进行实时监控且不会干扰空化活动。
进一步地，本发明根据延时叠加波束合成图中像素值大于像素值阈值的像素点的数目寻找聚焦超声换能器顶点坐标的遍历集合中的最优顶点坐标，该最优顶点坐标对应的最小像素点数目反映了最佳的成像分辨率，从而获得准确的聚焦超声换能器的位置信息，并提高发收时序同步超声被动空化成像的质量。
进一步地，本发明对随机选取的多帧被动空化射频信号所对应的最优顶点坐标进行平均，有助于降低误差，从而提高聚焦超声换能器顶点定位坐标的精度。
进一步地，本发明根据聚焦超声换能器发射脉冲信号的脉冲长度、被动空化射频信号的频谱分布的质心频率、信号采样频率、被动空化射频信号的规模这些实际参数以及与实际相同的聚焦超声换能器和超声成像换能器的空间设置来构造被动阵列仿真信号，该被动阵列仿真信号与实际所得被动空化射频信号有较高的契合度，根据该被动阵列仿真信号所得的自适应权重样本与实际中的自适应权重具有相似的统计特性，即样本中的主成分可以代表实际中的自适应权重，从而使得构造的维度转换矩阵更为可靠。
进一步地，本发明将根据被动阵列仿真信号计算的像素网格中各像素点处的自适应权重作为样本对自适应权重进行主成分分析，样本数(即像素点数目)远大于变量数(即子孔径的长度)，从而可满足统计学要求。
进一步地，本发明在构造维度转换矩阵时将最后一个主成分向量替换为直流向量，可避免由于主成分分析的均值中心化预处理而导致的低维域自适应波束合成的约束条件无法满足的问题。
进一步地，本发明提供了两种聚焦超声换能器和超声成像换能器的空间设置下绝对渡越时间的计算方法，可对目前超声空化研究中绝大多数常见的换能器空间设置下的绝对渡越时间进行计算，从而使得本发明提出的发收时序同步超声被动空化成像方法具有普遍适用性。
进一步地，本发明利用维度转换矩阵对被动空化射频信号对应的子孔径采样组合信号进行降维并通过低维域子孔径采样组合信号获得低维域协方差矩阵，该方法无需计算原始的高维域协方差矩阵就实现协方差矩阵的降维，在子孔径的长度较大时有效减小计算低维域协方差矩阵的计算量，从而提高发收时序同步超声被动空化成像的计算速度。
进一步地，本发明采用互谱度量法确定低维域协方差矩阵的信号子空间的维度，相比直接对特征值进行阈值化，更有益于成像质量的提高。
附图说明
图1为本发明实施例中实验系统的示意图；其中：1-透明水箱、2-吸声材料、3-聚焦超声换能器、4-脉冲发射/接收器、5-超声成像换能器、6-夹持装置、7-高精度三维位移台、8-超声成像数据采集系统、9-主控计算机、10-多通道波形发生器、11-多通道数字示波器。
图2为本发明实施例中得到的3帧发收时序同步超声被动空化成像结果，其中：(a)基于延时叠加波束合成的方法；(b)基于快速特征空间自适应波束合成的方法。
具体实施方式
下面结合附图和实施例对本发明作进一步的详细说明。所述实施例仅用于解释本发明，而非对本发明保护范围的限制。
参见图1，实验系统包括透明水箱1、吸声材料2、聚焦超声换能器3、脉冲发射/接收器4、超声成像换能器5、夹持装置6、高精度三维位移台7、超声成像数据采集系统8、主控计算机9、多通道波形发生器10和多通道数字示波器11。所述透明水箱1(例如由聚碳酸酯材料制成，长度、宽度和高度分别为50cm、28cm和28cm)的底部和侧壁放置吸声材料2，透明水箱1的上表面开放且一侧侧壁设有一可容纳聚焦超声换能器3的圆形孔洞，透明水箱1中充满除气水；聚焦超声换能器3(例如频率为1.6MHz、孔径为145mm)通过圆形孔洞固定于透明水箱1的侧壁，聚焦超声换能器3的接口端与脉冲发射/接收器4的输出端相连接；超声成像换能器5(例如阵元数目NE为128、阵元间距为0.3mm的线阵换能器)的探头端固定在夹持装置6(例如由聚碳酸酯材料制成)上，夹持装置6固定在高精度三维位移台7上，超声成像换能器5的接口端与超声成像数据采集系统8的换能器接口相连接；超声成像数据采集系统8的数据接口与主控计算机9相连接；多通道波形发生器10的两个通道在脉冲发射/接收器4和超声成像数据采集系统8工作在外部触发模式时分别与脉冲发射/接收器4的外部触发端口和超声成像数据采集系统8的外部触发端口相连接。
所述脉冲发射/接收器4根据设置的脉冲信号的频率、脉冲长度、脉冲重复频率和强度，在内部触发模式下激励聚焦超声换能器3向除气水中发射脉冲信号并在除气水中产生空化活动，或在外部触发模式下激励聚焦超声换能器3向介质(例如除气水)中发射脉冲信号并在介质中产生空化活动。
所述超声成像数据采集系统8在超声成像换能器5处于发射接收模式时采集B模式回波信号，得到B模式回波射频信号；超声成像数据采集系统8在超声成像换能器5处于不发射只接收模式(即被动接收模式)且聚焦超声换能器3的发射时序和超声成像换能器5的接收时序同步时采集空化声辐射信号，得到被动空化射频信号，此时超声成像数据采集系统8工作在外部触发模式下；其中超声成像换能器5的工作模式、超声成像数据采集系统8的触发模式、以及信号采样频率和信号采样点数在主控计算机9中设置，超声成像数据采集系统8采集得到的信号在主控计算机9中存储。
所述主控计算机9对超声成像数据采集系统8采集得到的B模式回波射频信号进行处理并生成B模式超声图像，主控计算机9对超声成像数据采集系统8采集得到的被动空化射频信号进行处理并生成高分辨率和高对比度的发收时序同步超声被动空化成像结果。
所述多通道波形发生器10利用其两个通道输出的触发信号同步地触发脉冲发射/接收器4和超声成像数据采集系统8，从而对聚焦超声换能器3的发射时序和超声成像换能器5的接收时序进行同步；在触发前将多通道波形发生器10的两个通道分别与多通道数字示波器11的两个通道相连接，从而利用多通道数字示波器11检查多通道波形发生器10的两个通道输出的触发信号。
本发明利用上述实验系统开展发收时序同步超声被动空化成像的流程如下：
(1)将脉冲发射/接收器4的触发模式设置为内部触发模式，将超声成像换能器5的工作模式设置为发射接收模式，根据B模式超声图像调整超声成像换能器5的位置，具体流程见下述步骤(1.1)～(1.4)。
(1.1)将脉冲发射/接收器4的触发模式设置为内部触发模式，并在脉冲发射/接收器4中设置激励聚焦超声换能器3的脉冲信号的频率(例如1.6MHz)、脉冲长度(例如5～10个周期)、脉冲重复频率(例如100～1000Hz)和强度(例如50～90)；
(1.2)在主控计算机9中将超声成像换能器5的工作模式设置为发射接收模式(例如聚焦扫描B模式)；
(1.3)通过脉冲发射/接收器4激励聚焦超声换能器3在除气水中产生空化活动，并通过B模式超声图像观察在除气水中产生的空化活动；
(1.4)根据B模式超声图像观察到的空化活动，通过高精度三维位移台7移动夹持装置6，将超声成像换能器5调整至一合适位置(例如在该位置处B模式超声图像中空化活动的信号最强)。
(2)在多通道波形发生器10的两个通道上分别设置触发信号并检查，分别在脉冲发射/接收器4和主控计算机9中设置参数，通过两个通道输出的触发信号同步地触发脉冲发射/接收器4和超声成像数据采集系统8，采集被动空化射频信号并存储，具体流程见下述步骤(2.1)～(2.6)。
(2.1)在多通道波形发生器10的两个通道上分别设置触发信号，其中第一个通道的触发信号(例如触发频率为1000Hz、触发次数为100，且无延时)用来触发脉冲发射/接收器4，第二个通道的触发信号(例如触发频率为1000Hz、触发次数为100，且触发延时τTrig为0～50μs)用来触发超声成像数据采集系统8；
(2.2)将多通道波形发生器10的两个通道分别与多通道数字示波器11的两个通道相连接，通过多通道数字示波器11检查多通道波形发生器10的两个通道输出的触发信号是否正确，若不正确，则对多通道波形发生器10中的参数设置进行检查、修改，直至触发信号正确无误，移除多通道数字示波器11；
(2.3)将多通道波形发生器10的两个通道分别与脉冲发射/接收器4和超声成像数据采集系统8的外部触发端口相连接，并将脉冲发射/接收器4和超声成像数据采集系统8的触发模式均设置为外部触发模式；
(2.4)在脉冲发射/接收器4中设置激励聚焦超声换能器3的脉冲信号的频率(例如1.6MHz)、脉冲长度(例如1个周期)、脉冲重复频率(例如1000Hz)和强度(例如90)；在主控计算机9中设置超声成像换能器5的工作模式为被动接收模式，设置信号采样频率fs(例如fs＝25MHz)和信号采样点数NS(例如NS＝5000)；
(2.5)利用多通道波形发生器10的两个通道输出的触发信号同步地触发脉冲发射/接收器4和超声成像数据采集系统8，以对聚焦超声换能器3发射和超声成像换能器5接收的时序进行同步；在第一个通道输出的触发信号的触发下，聚焦超声换能器3向介质(例如除气水)中发射脉冲信号并在该介质中产生空化活动，在第二个通道输出的触发信号的触发下，超声成像换能器5同步地对空化声辐射信号进行被动接收；
(2.6)利用超声成像数据采集系统8采集由超声成像换能器5被动接收的空化声辐射信号，得到NF(例如NF＝100)帧被动空化射频信号RFExpe,p(p＝1,2,...,NF)，其中第p帧被动空化射频信号RFExpe,p中包含NE个阵元的信号且信号采样点数为NS；将所得NF帧被动空化射频信号存储在主控计算机9中。
(3)建立聚焦超声换能器3顶点处切线或切面的方程，建立成像坐标系中的点到聚焦超声换能器3顶点处切线或切面的距离的计算公式，建立成像坐标系中的点到超声成像换能器5各阵元的距离的计算公式，建立绝对渡越时间的计算公式，具体流程见下述步骤(3.1)～(3.6)。
(3.1)以超声成像换能器5的中心为原点、以超声成像换能器5的横向和轴向分别为x轴和z轴，并以与超声成像换能器5的成像平面垂直的方向为y轴，建立三维位置坐标系；
(3.2)从步骤(3.1)所得三维位置坐标系中删除y轴，得到成像坐标系；
(3.3)建立聚焦超声换能器3顶点处切线或切面的方程；
对于换能器空间设置I，聚焦超声换能器3顶点处切线的方程为：
sinα(x-xV)+cosα(z-zV)＝0
对于换能器空间设置II，聚焦超声换能器3顶点处切面的方程为：
sinα(y-yV)+cosα(z-zV)＝0
其中，xV和zV分别为聚焦超声换能器3顶点在步骤(3.1)中所述三维位置坐标系或步骤(3.2)中所述成像坐标系中的x轴坐标和z轴坐标，yV为聚焦超声换能器3顶点在步骤(3.1)中所述三维位置坐标系中的y轴坐标；在换能器空间设置I中，聚焦超声换能器3的中轴线在超声成像换能器5的成像平面内，0≤α≤π/2；在换能器空间设置II中，聚焦超声换能器3的中轴线与超声成像换能器5的成像平面相交且与超声成像换能器5的横向方向垂直，0＜α≤π/2；α为聚焦超声换能器3的中轴线与超声成像换能器5的中轴线之间的夹角；
(3.4)根据步骤(3.3)所得方程建立步骤(3.2)中所述成像坐标系中的点到聚焦超声换能器3顶点处切线或切面的距离的计算公式：
对于换能器空间设置I，坐标为(x,z)的点到聚焦超声换能器3顶点处切线的距离的计算公式为：
Figure BDA0004049205960000151
由于
Figure BDA0004049205960000152
上式可化简为d1(x,z)＝|sinα(x-xV)+cosα(z-zV)|，其中|·|表示取绝对值；
对于换能器空间设置II，坐标为(x,z)的点到聚焦超声换能器3顶点处切面的距离的计算公式为：
Figure BDA0004049205960000153
由于
Figure BDA0004049205960000154
且成像坐标系中的点的y轴坐标为0，上式可化简为d1(x,z)＝|-sinα(yV)+cosα(z-zV)|；
(3.5)建立步骤(3.2)中所述成像坐标系中的点到超声成像换能器5各阵元的距离的计算公式：
Figure BDA0004049205960000155
其中，i＝1,2,...,NE，NE为超声成像换能器5的阵元数目，
Figure BDA0004049205960000169
为超声成像换能器5第i个阵元的坐标；
(3.6)根据步骤(3.4)和(3.5)所得计算公式建立绝对渡越时间的计算公式：
Figure BDA0004049205960000161
其中，c为声传播速度。
(4)根据聚焦超声换能器3顶点坐标的遍历集合中的顶点坐标对应的绝对渡越时间，对随机选取的被动空化射频信号进行延时叠加波束合成处理，根据延时叠加波束合成图从遍历集合中寻找最优顶点坐标，计算聚焦超声换能器3顶点定位坐标，具体流程见下述步骤(4.1)～(4.10)。
(4.1)在步骤(3.2)中所述成像坐标系下设置像素网格，该像素网格中像素点的坐标记为
Figure BDA0004049205960000162
其中，j＝1,2,...,NP，NP为像素点的数目；
(4.2)设置聚焦超声换能器3顶点坐标的遍历集合，其中，x轴坐标、y轴坐标和z轴坐标的遍历区间的长度≥10mm；x轴坐标和y轴坐标的遍历步长为步骤(4.1)中所述像素网格中x轴方向的像素间隔，z轴坐标的遍历步长为步骤(4.1)中所述像素网格中z轴方向的像素间隔；
(4.3)从步骤(2.6)所得NF帧被动空化射频信号中随机选取NCF(例如NCF＝10)帧被动空化射频信号，记为RFExpe,c(c＝1,2,...,NCF)；
(4.4)根据步骤(4.2)所得遍历集合中任意一顶点坐标计算所述像素网格中第j个像素点处的绝对渡越时间
Figure BDA0004049205960000163
其中，i＝1,2,...,NE，换能器空间设置如图1所示；即在步骤(3.4)和(3.5)所得计算公式中令(x,z)为
Figure BDA0004049205960000164
令
Figure BDA0004049205960000165
为所述遍历集合中任意一顶点坐标，并利用步骤(3.6)所得计算公式计算所述绝对渡越时间
Figure BDA0004049205960000166
(4.5)根据步骤(4.4)所得绝对渡越时间
Figure BDA0004049205960000167
对随机选取的第c帧被动空化射频信号RFExpe ,c进行延时叠加波束合成处理，得到第j个像素点处的延时叠加波束合成信号
Figure BDA0004049205960000168
Figure BDA0004049205960000171
其中，round{·}表示四舍五入取整，RFi Expe,c为RFExpe,c中第i个阵元的信号(i＝1,2,...,NE)，τTrig为步骤(2.1)中所述超声成像数据采集系统的触发延时；
(4.6)重复步骤(4.5)，计算得到步骤(4.1)中所述像素网格中各像素点处的延时叠加波束合成信号；
(4.7)对步骤(4.6)所得各像素点处的延时叠加波束合成信号进行希尔伯特变换，得到延时叠加波束合成图；
(4.8)设置像素值阈值并统计步骤(4.7)所得延时叠加波束合成图中像素值大于像素值阈值的像素点的数目，其中像素值阈值为该延时叠加波束合成图中最大像素值的0.1～0.5倍；
(4.9)重复步骤(4.4)～(4.8)，得到步骤(4.2)中所述遍历集合中所有顶点坐标各自对应的延时叠加波束合成图中像素值大于像素值阈值的像素点数目，从中寻找最小像素点数目并记录其对应的顶点坐标，得到随机选取的第c帧被动空化射频信号对应的最优顶点坐标；
(4.10)重复步骤(4.4)～(4.9)，得到随机选取的NCF帧被动空化射频信号对应的最优顶点坐标，对所得的NCF个最优顶点坐标分别在每个方向(x轴、y轴、z轴方向)上进行平均，得到聚焦超声换能器3顶点定位坐标(xVL,yVL,zVL)。
(5)根据聚焦超声换能器3顶点定位坐标计算仿真声源群中各声源处的绝对渡越时间，根据该绝对渡越时间和声源辐射仿真信号构造各声源的阵列仿真信号，叠加各声源的阵列仿真信号得到被动阵列仿真信号，具体流程见下述步骤(5.1)～(5.7)。
(5.1)在步骤(3.2)中所述成像坐标系下设置仿真声源群中各声源的坐标，其中第m个声源的坐标表示为
Figure BDA0004049205960000172
m＝1,2,...,NA，NA为声源的数目(例如10～100)，各声源坐标的分布为正态分布(例如x轴坐标和z轴坐标的中心分别为0mm和40mm，x轴坐标和z轴坐标的标准差分别为1mm和1mm)；
(5.2)设置声源辐射仿真信号为加Hanning窗的正弦信号，记为sae，该信号的采样点数为
Figure BDA0004049205960000173
该信号的时间长度与步骤(2.5)中聚焦超声换能器3发射脉冲信号的脉冲长度一致(例如1个周期)，该信号的频率由步骤(2.6)中所得被动空化射频信号的频谱分布的质心频率确定(例如为5MHz等所述质心频率附近的整数MHz级频率)，该信号的信号采样频率与步骤(2.4)中设置的信号采样频率fs一致(例如25MHz)；
(5.3)根据步骤(2.6)中所得被动空化射频信号的规模设置被动阵列仿真信号的规模，即被动阵列仿真信号中包含NE个阵元的仿真信号且各阵元的仿真信号的采样点数为NS；
(5.4)根据步骤(4.10)所得聚焦超声换能器3顶点定位坐标(xVL,yVL,zVL)计算步骤(5.1)中所述仿真声源群中第m个声源处的绝对渡越时间
Figure BDA0004049205960000181
其中，i＝1,2,...,NE，换能器空间设置参照图1；即在与实际相同的换能器空间设置下，在步骤(3.4)和(3.5)所得计算公式中令(x,z)为
Figure BDA0004049205960000182
令(xV,yV,zV)为(xVL,yVL,zVL)，并利用步骤(3.6)所得计算公式计算所述绝对渡越时间
Figure BDA0004049205960000183
(5.5)根据步骤(5.2)中所述声源辐射仿真信号sae和步骤(5.4)所得绝对渡越时间
Figure BDA0004049205960000184
构造超声成像换能器5第i个阵元的仿真信号
Figure BDA0004049205960000185
Figure BDA0004049205960000186
其中，round{·}表示四舍五入取整，τTrig为步骤(2.1)中所述超声成像数据采集系统的触发延时，k＝1,2,...,NS；
(5.6)重复步骤(5.5)，得到超声成像换能器5各阵元的仿真信号，根据超声成像换能器5各阵元的仿真信号构造第m个声源的阵列仿真信号RFSimu,m；
(5.7)重复步骤(5.4)～(5.6)，得到各声源的阵列仿真信号，将各声源的阵列仿真信号叠加得到被动阵列仿真信号RFSimu：
Figure BDA0004049205960000187
(6)根据聚焦超声换能器3顶点定位坐标计算像素网格中各像素点处的绝对渡越时间，从被动阵列仿真信号中提取采样点得到全孔径采样信号，通过分割全孔径采样信号和对角加载计算协方差矩阵并计算自适应权重，对各像素点处的自适应权重进行主成分分析，根据主成分向量构造维度转换矩阵，具体流程见下述步骤(6.1)～(6.10)。
(6.1)根据步骤(4.10)所得聚焦超声换能器3顶点定位坐标(xVL,yVL,zVL)计算步骤(4.1)中所述像素网格中第j个像素点处的绝对渡越时间
Figure BDA0004049205960000191
其中，i＝1,2,...,NE，j＝1,2,...,NP，换能器空间设置参照图1；即在与实际相同的换能器空间设置下，在步骤(3.4)和(3.5)所得计算公式中令(x,z)为
Figure BDA0004049205960000192
令(xV,yV,zV)为(xVL,yVL,zVL)，并利用步骤(3.6)所得计算公式计算所述绝对渡越时间
Figure BDA0004049205960000193
(6.2)根据步骤(6.1)所得绝对渡越时间
Figure BDA0004049205960000194
从步骤(5.7)所得被动阵列仿真信号RFSimu中提取采样点，得到全孔径采样信号
Figure BDA0004049205960000195
其中第i个阵元的信号为
Figure BDA0004049205960000196
Figure BDA0004049205960000197
其中，
Figure BDA0004049205960000198
为RFSimu中第i个阵元的信号(i＝1,2,...,NE)，round{·}表示四舍五入取整，τTrig为步骤(2.1)中所述超声成像数据采集系统的触发延时，fs为信号采样频率；
(6.3)将步骤(6.2)所得全孔径采样信号
Figure BDA0004049205960000199
分割成NE-L+1个重叠的子孔径采样信号
Figure BDA00040492059600001910
Figure BDA00040492059600001911
其中，l＝1,2,...,NE-L+1，L为子孔径的长度(例如NE/2)，[·]T表示转置；
(6.4)利用步骤(6.3)所得NE-L+1个子孔径采样信号
Figure BDA00040492059600001912
构造子孔径采样组合信号
Figure BDA00040492059600001913
Figure BDA00040492059600001914
(6.5)计算步骤(6.4)所得子孔径采样组合信号
Figure BDA00040492059600001915
的协方差矩阵
Figure BDA0004049205960000201
Figure BDA0004049205960000202
其中，[·]H表示共轭转置；
(6.6)对步骤(6.5)所得协方差矩阵
Figure BDA0004049205960000203
进行对角加载，得到对角加载后的协方差矩阵
Figure BDA0004049205960000204
Figure BDA0004049205960000205
其中，
Figure BDA0004049205960000206
为一L行L列的矩阵，δ为对角加载因子(例如1/L)，trace{·}表示求迹运算，I为单位矩阵；
(6.7)根据步骤(6.6)所得协方差矩阵
Figure BDA0004049205960000207
计算第j个像素点处的自适应权重
Figure BDA0004049205960000208
Figure BDA0004049205960000209
其中，inv[·]表示求逆运算，a为方向向量，该向量中的元素均为1；
(6.8)重复步骤(6.1)～(6.7)，得到像素网格中各像素点处的自适应权重；
(6.9)将步骤(6.8)所得各像素点处的自适应权重作为样本，对自适应权重进行主成分分析，得到L个主成分向量：
(6.9.1)计算自适应权重的协方差矩阵
Figure BDA00040492059600002010
Figure BDA00040492059600002011
其中，μSimu为各像素点处的自适应权重的平均，NP为像素点的数目；
(6.9.2)对步骤(6.9.1)所得协方差矩阵
Figure BDA00040492059600002012
进行特征分解：
Figure BDA00040492059600002013
其中，Λ的对角线元素为降序排列的特征值(L个)，V＝[v1,v2,...,vL]，V中包含了对应于L个降序排列的特征值的L个特征向量，即L个主成分向量；
(6.10)从步骤(6.9)所得L个主成分向量中选择最大特征值对应的前Q个主成分向量v1,v2,...,vQ-1,vQ(Q≤L，例如Q＝8)并将最后一个主成分向量vQ替换为元素均为
Figure BDA00040492059600002112
的直流向量dc，然后构造维度转换矩阵C：
C＝[v1,v2,...,vQ-1,dc]
其中，维度转换矩阵C为一L行Q列的矩阵。
(7)利用维度转换矩阵对任意一帧被动空化射频信号对应的各像素点处的对角加载后的协方差矩阵进行降维处理，计算低维域自适应权重并将其投影到低维域协方差矩阵经特征分解所得的信号子空间上，根据低维域特征空间自适应权重和低维域子孔径采样平均信号计算各像素点处的快速特征空间自适应波束合成信号，经希尔伯特变换、归一化及对数化处理得到发收时序同步超声被动空化成像结果，具体流程见下述步骤(7.1)～(7.9)。
(7.1)将步骤(6.2)中的被动阵列仿真信号RFSimu替换为步骤(2.6)所得的任意一帧被动空化射频信号，然后按步骤(6.1)～(6.6)依次计算步骤(4.1)中所述像素网格中第j个像素点处的绝对渡越时间
Figure BDA0004049205960000211
以及该帧被动空化射频信号对应的步骤(4.1)中所述像素网格中第j个像素点处的全孔径采样信号
Figure BDA0004049205960000212
NE-L+1个子孔径采样信号
Figure BDA0004049205960000213
子孔径采样组合信号
Figure BDA0004049205960000214
子孔径采样组合信号的协方差矩阵
Figure BDA0004049205960000215
对角加载后的协方差矩阵
Figure BDA0004049205960000216
其中，p＝1,2,...,NF，j＝1,2,...,NP：
(7.1.1)根据步骤(4.10)所得聚焦超声换能器3顶点定位坐标(xVL,yVL,zVL)计算步骤(4.1)中所述像素网格中第j个像素点处的绝对渡越时间
Figure BDA0004049205960000217
其中，i＝1,2,...,NE，j＝1,2,...,NP，换能器空间设置如图1所示；即在步骤(3.4)和(3.5)所得计算公式中令
Figure BDA0004049205960000218
令(xV,yV,zV)为(xVL,yVL,zVL)，并利用步骤(3.6)所得计算公式计算所述绝对渡越时间
Figure BDA0004049205960000219
(7.1.2)根据步骤(7.1.1)所得绝对渡越时间
Figure BDA00040492059600002110
从步骤(2.6)所得的第p帧被动空化射频信号RFExpe,p中提取采样点，得到全孔径采样信号
Figure BDA00040492059600002111
其中第i个阵元的信号为
Figure BDA0004049205960000221
Figure BDA0004049205960000222
其中，
Figure BDA0004049205960000223
为RFExpe,p中第i个阵元的信号(i＝1,2,...,NE)，round{·}表示四舍五入取整，τTrig为步骤(2.1)中所述超声成像数据采集系统的触发延时，fs为信号采样频率；
(7.1.3)将步骤(7.1.2)所得全孔径采样信号
Figure BDA0004049205960000224
分割成NE-L+1个重叠的子孔径采样信号
Figure BDA0004049205960000225
Figure BDA0004049205960000226
其中，l＝1,2,...,NE-L+1，L为子孔径的长度(例如NE/2)，[·]T表示转置；
(7.1.4)利用步骤(7.1.3)所得NE-L+1个子孔径采样信号
Figure BDA0004049205960000227
构造子孔径采样组合信号
Figure BDA0004049205960000228
Figure BDA0004049205960000229
(7.1.5)计算步骤(7.1.4)所得子孔径采样组合信号
Figure BDA00040492059600002210
的协方差矩阵
Figure BDA00040492059600002211
Figure BDA00040492059600002212
其中，[·]H表示共轭转置；
(7.1.6)对步骤(7.1.5)所得协方差矩阵
Figure BDA00040492059600002213
进行对角加载，得到对角加载后的协方差矩阵
Figure BDA00040492059600002214
Figure BDA00040492059600002215
其中，
Figure BDA00040492059600002216
为一L行L列的矩阵，δ为对角加载因子(例如1/L)，trace{·}表示求迹运算，I为单位矩阵；
(7.2)利用步骤(6.10)所得维度转换矩阵C对步骤(7.1)所得对角加载后的协方差矩阵
Figure BDA0004049205960000231
进行降维处理，得到低维域协方差矩阵
Figure BDA0004049205960000232
Figure BDA0004049205960000233
其中，[·]H表示共轭转置，
Figure BDA0004049205960000234
为一L行L列的矩阵，
Figure BDA0004049205960000235
为一Q行Q列的矩阵；通过上式协方差矩阵的维度从L行L列降低至Q行Q列；
上式可表示为：
Figure BDA0004049205960000236
其中，
Figure BDA0004049205960000237
Figure BDA0004049205960000238
为低维域子孔径采样组合信号，δ为步骤(7.1)中所述对角加载因子，ξ为步骤(7.1)所得子孔径采样组合信号
Figure BDA0004049205960000239
中各元素的平方之和与子孔径采样信号的数目(NE-L+1)的比值；
(7.3)根据步骤(7.2)所得低维域协方差矩阵
Figure BDA00040492059600002310
计算低维域自适应权重
Figure BDA00040492059600002311
Figure BDA00040492059600002312
其中，inv[·]表示求逆运算，b＝CHa，C为步骤(6.10)所得维度转换矩阵，a为步骤(6.7)中所述方向向量；
(7.4)对步骤(7.2)所得低维域协方差矩阵
Figure BDA00040492059600002313
进行特征分解，得到信号子空间
Figure BDA00040492059600002314
将步骤(7.3)所得低维域自适应权重
Figure BDA00040492059600002315
投影到该信号子空间上，得到低维域特征空间自适应权重
Figure BDA00040492059600002316
Figure BDA00040492059600002317
其中，信号子空间
Figure BDA00040492059600002318
的维度根据互谱度量法(例如互谱度量系数为1％～5％)确定；
(7.5)根据步骤(6.10)所得维度转换矩阵C和步骤(7.1)所得的NE-L+1个子孔径采样信号
Figure BDA0004049205960000243
计算低维域子孔径采样平均信号
Figure BDA0004049205960000244
Figure BDA0004049205960000241
(7.6)根据步骤(7.4)所得低维域特征空间自适应权重
Figure BDA0004049205960000247
和步骤(7.5)所得低维域子孔径采样平均信号
Figure BDA0004049205960000245
计算第j个像素点处的快速特征空间自适应波束合成信号
Figure BDA0004049205960000246
Figure BDA0004049205960000242
(7.7)重复步骤(7.1)～(7.6)，计算得到步骤(4.1)中所述像素网格中各像素点处的快速特征空间自适应波束合成信号；
(7.8)对步骤(7.7)所得各像素点处的快速特征空间自适应波束合成信号进行希尔伯特变换，得到快速特征空间自适应波束合成图；
(7.9)对步骤(7.8)所得快速特征空间自适应波束合成图进行归一化及对数化处理，得到任意一帧(即第p帧)被动空化射频信号对应的发收时序同步超声被动空化成像结果。
在图1所示实验系统中，聚焦超声换能器3与超声成像换能器5的空间设置属于上述换能器空间设置I，其中α＝π/2。对该实验系统所得的第30帧、第60帧和第90帧被动空化射频信号进行处理，得到如图2所示的第30帧、第60帧和第90帧发收时序同步超声被动空化成像结果，其中图2(a)和图2(b)分别为根据传统的基于延时叠加波束合成的方法和本发明所用的基于快速特征空间自适应波束合成的方法得到的成像结果，成像结果以40dB的动态范围进行显示。
可以看出，图2(a)所示的任意一帧成像结果中主瓣的尺寸较大且出现了明显的旁瓣干扰，表明成像分辨率和成像对比度较低(由于延时叠加波束合成的非自适应性，传统基于延时叠加波束合成的发收时序同步超声被动空化成像方法存在分辨率差且对比度低的缺陷)；而图2(b)所示的对应帧成像结果中的主瓣尺寸和旁瓣干扰水平明显降低，表明成像分辨率和成像对比度得到提升，同时也表明通过数值仿真方法构造的维度转换矩阵对于不同帧的被动空化射频信号具有普遍适用性。
采用成像结果中像素值大于-6dB的像素点的面积Area-6dB定量评价成像分辨率，Area-6dB越小，成像分辨率越高；采用CR(CR＝20lg(APVinside/APVoutside))定量评价成像对比度，CR越大，成像对比度越高，其中，APVinside和APVoutside分别为成像结果中像素值大于-6dB的像素点和像素值小于-6dB的像素点在未进行对数化处理的成像结果中对应的像素值的平均值。经过计算，图2(a)中3帧成像结果对应的Area-6dB分别为0.75mm2、1.51mm2和1.24mm2，对应的CR分别为39.81dB、38.11dB和41.48dB；图2(b)中3帧成像结果对应的Area-6dB分别为0.41mm2、0.80mm2和0.59mm2，对应的CR分别为47.69dB、44.37dB和49.20dB。上述分析表明，本发明能够有效地同时提高发收时序同步超声被动空化成像的分辨率和对比度。
特征空间自适应波束合成的计算负担主要由协方差矩阵的求逆和特征分解引起，对协方差矩阵进行求逆和特征分解所需要的浮点运算数为
Figure BDA0004049205960000251
其中，
Figure BDA0004049205960000252
为协方差矩阵的行数或列数。本发明利用通过数值仿真方法构造的维度转换矩阵将协方差矩阵的维度从L行L列(例如L＝64)降低至Q行Q列(例如Q＝8)，使得浮点运算数从2L3/3+21L3＝5679787降低至2Q3/3+21Q3＝11094，从而能够以更快的计算速度进行特征空间自适应波束合成，并因此提高发收时序同步超声被动空化成像的计算速度。
本发明具有以下优点：
(1)本发明通过对聚焦超声换能器发射和超声成像换能器接收的时序进行同步，允许利用绝对渡越时间信息对被动空化射频信号进行处理，从而可在短脉冲发射情境下使得超声被动空化成像的轴向分辨率得以提高。同时，本发明根据延时叠加波束合成图寻找最优顶点坐标并计算聚焦超声换能器顶点定位坐标，可为计算绝对渡越时间提供准确的聚焦超声换能器的位置信息，从而能够获得更好的成像质量。
(2)本发明通过数值仿真的方法构造维度转换矩阵，可避免在处理每帧被动空化射频信号时每次构造相应的维度转换矩阵的问题，同时无需通过预实验来构造维度转换矩阵；本发明构造的维度转换矩阵只需构造一次即可重复利用，从而为数据处理提供了极大的便利性。且本发明根据实际中的参数(聚焦超声换能器发射参数、被动空化射频信号采集参数、被动空化射频信号的频谱分布的质心频率)和与实际相同的换能器空间设置构造被动阵列仿真信号，根据该信号所得的自适应权重样本与实际中的自适应权重有相似的统计特性，从而增强了构造的维度转换矩阵的可靠性。且本发明在构造维度转换矩阵的过程中以像素网格中各像素点处的自适应权重作为主成分分析的样本并将最后一个主成分向量替换为直流向量，从而满足主成分分析的统计学要求并避免了低维域自适应波束合成的约束条件无法满足的问题。
(3)本发明采用的特征空间自适应波束合成，通过基于最小方差无失真响应的自适应波束合成和将自适应权重投影到特征分解所得信号子空间，可有效地减小主瓣尺寸并降低旁瓣水平，从而使得发收时序同步超声被动空化成像的分辨率和对比度得到同步的提升。且本发明通过子孔径平均和对角加载来计算协方差矩阵，可对协方差矩阵进行良好的估计并增强自适应波束合成的鲁棒性。且本发明采用互谱度量法来确定信号子空间的维度，也有益于成像质量的提高。
(4)本发明利用通过对被动阵列仿真信号对应的自适应权重进行主成分分析构造的维度转换矩阵，降低了协方差矩阵的维度，使得对协方差矩阵进行求逆和特征分解所需的巨大计算量得以降低，从而提高特征空间自适应波束合成的计算速度，并因此提高发收时序同步超声被动空化成像的计算速度。且本发明采用的对被动空化射频信号对应的子孔径采样组合信号进行降维、然后通过降维所得信号获得低维域协方差矩阵的方法，在子孔径的长度较大时可有效减小计算低维域协方差矩阵的计算量，也有助于发收时序同步超声被动空化成像计算速度的提高。
(5)本发明提供了两种聚焦超声换能器和超声成像换能器的空间设置下绝对渡越时间的计算方法，从而使得本发明提出的高分辨高对比快速计算的发收时序同步超声被动空化成像方法适用于目前超声空化研究中绝大多数常见的换能器空间设置。本发明可对超声场中空化活动的分布进行高分辨率、高对比度和快速的重建，为高强度组织毁损、低强度血脑屏障开放等多种超声空化应用中空化时空动力学演变过程的监控提供了有力的影像手段，并为相关生物物理机制的澄清以及超声发射参数的调控优化奠定了坚实基础。

Claims (10)

1.一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：包括以下步骤：
1)通过同步触发对聚焦超声换能器(3)发射和超声成像换能器(5)接收的时序进行同步，对由超声成像换能器被动接收的空化声辐射信号进行采集，得到被动空化射频信号；
2)在成像坐标系下设置像素网格，计算聚焦超声换能器(3)顶点坐标的遍历集合中每一顶点坐标对应的像素网格中各像素点处的绝对渡越时间，根据该绝对渡越时间对从步骤1采集得到的被动空化射频信号中选取的多帧被动空化射频信号分别进行延时叠加波束合成处理，得到所述遍历集合中每一顶点坐标对应的各像素点处的延时叠加波束合成信号，对所述遍历集合中每一顶点坐标对应的各像素点处的延时叠加波束合成信号进行希尔伯特变换，得到所述遍历集合中各顶点坐标对应的延时叠加波束合成图，根据所述遍历集合中各顶点坐标对应的延时叠加波束合成图寻找所述遍历集合中的最优顶点坐标，根据在选取的不同帧被动空化射频信号下分别寻找到的最优顶点坐标计算聚焦超声换能器(3)顶点定位坐标；
3)设置声源辐射仿真信号和被动阵列仿真信号规模并在成像坐标系下设置仿真声源群中各声源的坐标，根据聚焦超声换能器(3)顶点定位坐标计算仿真声源群中每一声源处的绝对渡越时间，根据该绝对渡越时间和声源辐射仿真信号构造超声成像换能器(5)各阵元的仿真信号，根据超声成像换能器(5)各阵元的仿真信号构造对应声源的阵列仿真信号，对各声源的阵列仿真信号进行叠加得到被动阵列仿真信号；
4)根据聚焦超声换能器(3)顶点定位坐标计算所述像素网格中各像素点处的绝对渡越时间，根据该绝对渡越时间从被动阵列仿真信号中分别提取采样点，得到对应像素点处的全孔径采样信号，将对应像素点处的全孔径采样信号分割成多个重叠的子孔径采样信号并利用子孔径采样信号构造子孔径采样组合信号，计算子孔径采样组合信号的协方差矩阵并对该协方差矩阵进行对角加载，根据对角加载后的协方差矩阵计算对应像素点处的自适应权重，对各像素点处的自适应权重进行主成分分析，得到主成分向量，根据主成分向量构造维度转换矩阵；
5)根据步骤4中计算的绝对渡越时间从步骤1采集得到的任意一帧被动空化射频信号中提取并分割得到对应像素点处的多个重叠的子孔径采样信号，利用由该子孔径采样信号构造的子孔径采样组合信号计算相应的协方差矩阵并对该协方差矩阵进行对角加载，利用所述维度转换矩阵对对角加载后的协方差矩阵进行降维处理，得到对应像素点处的低维域协方差矩阵，根据低维域协方差矩阵计算对应像素点处的低维域自适应权重并将该低维域自适应权重投影到所述低维域协方差矩阵经特征分解所得的信号子空间上，得到对应像素点处的低维域特征空间自适应权重，根据所述维度转换矩阵和子孔径采样信号计算对应像素点处的低维域子孔径采样平均信号，根据所述低维域特征空间自适应权重和低维域子孔径采样平均信号计算对应像素点处的快速特征空间自适应波束合成信号，对各像素点处的快速特征空间自适应波束合成信号进行希尔伯特变换，得到快速特征空间自适应波束合成图，对快速特征空间自适应波束合成图进行归一化及对数化处理，得到对应帧被动空化射频信号的发收时序同步超声被动空化成像结果。
2.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤2中，寻找所述遍历集合中的最优顶点坐标具体包括以下步骤：
2.1)统计所述遍历集合中每一顶点坐标对应的延时叠加波束合成图中像素值大于像素值阈值的像素点的数目，其中像素值阈值为该延时叠加波束合成图中最大像素值的0.1～0.5倍；
2.2)重复步骤2.1，得到所述遍历集合中各顶点坐标对应的像素值大于像素值阈值的像素点数目，寻找这些像素点数目中的最小像素点数目，该最小像素点数目对应的顶点坐标为最优顶点坐标。
3.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤3中，声源辐射仿真信号为加Hanning窗的正弦信号，声源辐射仿真信号的时间长度与聚焦超声换能器(3)发射脉冲信号的脉冲长度一致，声源辐射仿真信号的频率由被动空化射频信号的频谱分布的质心频率确定，声源辐射仿真信号的采样频率与被动空化射频信号的采样频率一致，被动阵列仿真信号的规模与被动空化射频信号的规模一致；所述步骤3中，计算仿真声源群中每一声源处的绝对渡越时间所用的换能器空间设置与步骤1中聚焦超声换能器(3)和超声成像换能器(5)的空间设置相同；所述步骤4中，计算所述像素网格中各像素点处的绝对渡越时间所用的换能器空间设置与步骤1中聚焦超声换能器(3)和超声成像换能器(5)的空间设置相同。
4.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤3中，构造的任意一声源的阵列仿真信号中超声成像换能器任意一阵元的仿真信号表示为：
Figure FDA0004049205950000031
其中，RFi Simu,m(k)为仿真声源群中第m个声源的阵列仿真信号中第i个阵元的仿真信号，i＝1,2,...,NE，NE为超声成像换能器(5)的阵元数目，k＝1,2,...,NS，NS为被动空化射频信号的采样点数，m＝1,2,...,NA，NA为仿真声源群中声源的数目，round{·}表示四舍五入取整，
Figure FDA0004049205950000032
为根据聚焦超声换能器(3)顶点定位坐标计算的第m个声源处的绝对渡越时间，τTrig为采集空化声辐射信号的触发延时，fs为被动空化射频信号的采样频率，sae为声源辐射仿真信号，
Figure FDA0004049205950000033
为声源辐射仿真信号的采样点数。
5.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤4中，对各像素点处的自适应权重进行主成分分析具体包括以下步骤：
将根据被动阵列仿真信号计算的所述像素网格中各像素点处的自适应权重作为样本，计算自适应权重的协方差矩阵，然后对该协方差矩阵进行特征分解，从而得到对应于L个降序排列的特征值的L个特征向量，即得到L个主成分向量；
所述步骤4中，维度转换矩阵的构造具体包括以下步骤：
从L个主成分向量中选择最大特征值对应的前Q个主成分向量v1,v2,...,vQ-1,vQ，并将主成分向量vQ替换为元素均为
Figure FDA0004049205950000034
的直流向量dc，所得维度转换矩阵表示为：
C＝[v1,v2,...,vQ-1,dc]。
6.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤2、步骤3和步骤4中，绝对渡越时间的计算具体包括以下步骤：
S1以超声成像换能器(5)的中心为原点、以超声成像换能器(5)的横向和轴向分别为x轴和z轴，并以与超声成像换能器(5)的成像平面垂直的方向为y轴，建立三维位置坐标系；
S2从所述三维位置坐标系中删除y轴，得到成像坐标系；
S3建立聚焦超声换能器(3)顶点处切线或切面的方程：
对于换能器空间设置I，聚焦超声换能器(3)顶点处切线的方程为：
sinα(x-xV)+cosα(z-zV)＝0
对于换能器空间设置II，聚焦超声换能器(3)顶点处切面的方程为：
sinα(y-yV)+cosα(z-zV)＝0
其中，xV和zV分别为聚焦超声换能器(3)顶点在所述三维位置坐标系中的x轴坐标和z轴坐标，yV为聚焦超声换能器(3)顶点在所述三维位置坐标系中的y轴坐标；在换能器空间设置I中，聚焦超声换能器(3)的中轴线在超声成像换能器(5)的成像平面内，0≤α≤π/2；在换能器空间设置II中，聚焦超声换能器(3)的中轴线与超声成像换能器(5)的成像平面相交且与超声成像换能器(5)的横向方向垂直，0＜α≤π/2；α为聚焦超声换能器(3)的中轴线与超声成像换能器(5)的中轴线之间的夹角；
S4建立成像坐标系中的点到聚焦超声换能器顶点处切线或切面的距离的计算公式：
对于换能器空间设置I，坐标为(x,z)的点到聚焦超声换能器(3)顶点处切线的距离的计算公式为：
d1(x,z)＝|sinα(x-xV)+cosα(z-zV)|
其中，|·|表示取绝对值；
对于换能器空间设置II，坐标为(x,z)的点到聚焦超声换能器(3)顶点处切面的距离的计算公式为：
d1(x,z)＝|-sinα(yV)+cosα(z-zV)|
S5建立成像坐标系中坐标为(x,z)的点到超声成像换能器(5)各阵元的距离的计算公式：
Figure FDA0004049205950000041
其中，i＝1,2,...,NE，
Figure FDA0004049205950000042
为超声成像换能器(5)第i个阵元的坐标；
S6根据步骤S4和S5中的距离的计算公式建立绝对渡越时间的计算公式：
Figure FDA0004049205950000051
其中，c为声传播速度。
7.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤5中，低维域协方差矩阵的计算公式为：
Figure FDA0004049205950000052
其中，j＝1,2,...,NP，NP为像素点的数目，NE为超声成像换能器(5)的阵元数目，L为子孔径的长度，
Figure FDA0004049205950000053
和
Figure FDA0004049205950000054
分别为第p帧被动空化射频信号对应的第j个像素点处的低维域子孔径采样组合信号和子孔径采样组合信号，δ为对角加载因子，ξ为
Figure FDA0004049205950000055
中各元素的平方之和与子孔径采样信号的数目的比值，C为维度转换矩阵。
8.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤5中，低维域协方差矩阵经特征分解所得的信号子空间的维度根据互谱度量法确定，其中，互谱度量系数为1％～5％。
9.根据权利要求1所述一种高分辨高对比快速计算的发收时序同步超声被动空化成像方法，其特征在于：所述步骤5中，低维域子孔径采样平均信号的计算公式为：
Figure FDA0004049205950000056
其中，j＝1,2,...,NP，NP为像素点的数目，NE为超声成像换能器(5)的阵元数目，L为子孔径的长度，C为维度转换矩阵，
Figure FDA0004049205950000057
为第p帧被动空化射频信号对应的第j个像素点处的第l个子孔径采样信号。
10.一种高分辨高对比快速计算的发收时序同步超声被动空化成像系统，该系统包括聚焦超声换能器顶点定位坐标计算模块、被动阵列仿真信号构造模块、维度转换矩阵构造模块、快速特征空间自适应波束合成模块以及图像显示模块；
所述聚焦超声换能器顶点定位坐标计算模块用于在成像坐标系下设置像素网格、计算聚焦超声换能器(3)顶点坐标的遍历集合中每一顶点坐标对应的像素网格中各像素点处的绝对渡越时间、根据该绝对渡越时间对从采集得到的被动空化射频信号中选取的多帧被动空化射频信号分别进行延时叠加波束合成处理、对处理所得所述遍历集合中每一顶点坐标对应的各像素点处的延时叠加波束合成信号进行希尔伯特变换、根据变换所得所述遍历集合中各顶点坐标对应的延时叠加波束合成图寻找所述遍历集合中的最优顶点坐标，以及根据在选取的不同帧被动空化射频信号下分别寻找到的最优顶点坐标计算聚焦超声换能器(3)顶点定位坐标；
所述被动阵列仿真信号构造模块用于设置声源辐射仿真信号和被动阵列仿真信号规模并在成像坐标系下设置仿真声源群中各声源的坐标、根据所述聚焦超声换能器顶点定位坐标计算模块得到的聚焦超声换能器(3)顶点定位坐标计算仿真声源群中每一声源处的绝对渡越时间、根据该绝对渡越时间和声源辐射仿真信号构造超声成像换能器(5)各阵元的仿真信号、根据超声成像换能器(5)各阵元的仿真信号构造对应声源的阵列仿真信号，以及通过对各声源的阵列仿真信号进行叠加得到被动阵列仿真信号；
所述维度转换矩阵构造模块用于根据所述聚焦超声换能器顶点定位坐标计算模块得到的聚焦超声换能器(3)顶点定位坐标计算所述像素网格中各像素点处的绝对渡越时间、根据该绝对渡越时间从所述被动阵列仿真信号构造模块得到的被动阵列仿真信号中分别提取采样点、将提取所得对应像素点处的全孔径采样信号分割成多个重叠的子孔径采样信号并利用子孔径采样信号构造子孔径采样组合信号、计算子孔径采样组合信号的协方差矩阵并对该协方差矩阵进行对角加载、根据对角加载后的协方差矩阵计算对应像素点处的自适应权重、对各像素点处的自适应权重进行主成分分析，以及根据分析所得主成分向量构造维度转换矩阵；
所述快速特征空间自适应波束合成模块用于根据所述维度转换矩阵构造模块中计算的绝对渡越时间从采集得到的任意一帧被动空化射频信号中提取并分割得到对应像素点处的多个重叠的子孔径采样信号、利用由该子孔径采样信号构造的子孔径采样组合信号计算相应的协方差矩阵并对该协方差矩阵进行对角加载、利用所述维度转换矩阵构造模块得到的维度转换矩阵对对角加载后的协方差矩阵进行降维处理、根据降维处理所得对应像素点处的低维域协方差矩阵计算对应像素点处的低维域自适应权重并将该低维域自适应权重投影到所述低维域协方差矩阵经特征分解所得的信号子空间上、根据所述维度转换矩阵和子孔径采样信号计算对应像素点处的低维域子孔径采样平均信号、根据投影所得对应像素点处的低维域特征空间自适应权重和所述低维域子孔径采样平均信号计算对应像素点处的快速特征空间自适应波束合成信号，以及通过对各像素点处的快速特征空间自适应波束合成信号进行希尔伯特变换得到快速特征空间自适应波束合成图；
所述图像显示模块用于对所述快速特征空间自适应波束合成模块得到的快速特征空间自适应波束合成图进行归一化及对数化处理，对处理所得任意一帧被动空化射频信号对应的发收时序同步超声被动空化成像结果进行显示、或对多帧被动空化射频信号对应的发收时序同步超声被动空化成像结果进行动态的显示。
CN202310037068.6A 2023-01-03 2023-01-10 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统 Active CN116115260B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202310004626 2023-01-03
CN2023100046269 2023-01-03

Publications (2)

Publication Number Publication Date
CN116115260A true CN116115260A (zh) 2023-05-16
CN116115260B CN116115260B (zh) 2024-05-24

Family

ID=86304250

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310037068.6A Active CN116115260B (zh) 2023-01-03 2023-01-10 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统

Country Status (1)

Country Link
CN (1) CN116115260B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104535645A (zh) * 2014-12-27 2015-04-22 西安交通大学 微秒分辨空化时空分布的三维空化定量成像方法
CN109431536A (zh) * 2018-09-17 2019-03-08 西安交通大学 一种聚焦超声空化的实时高分辨时空分布成像方法与系统
CN114098799A (zh) * 2021-10-27 2022-03-01 西安交通大学 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104535645A (zh) * 2014-12-27 2015-04-22 西安交通大学 微秒分辨空化时空分布的三维空化定量成像方法
WO2016101382A1 (zh) * 2014-12-27 2016-06-30 西安交通大学 微秒分辨空化时空分布的三维空化定量成像方法
CN109431536A (zh) * 2018-09-17 2019-03-08 西安交通大学 一种聚焦超声空化的实时高分辨时空分布成像方法与系统
CN114098799A (zh) * 2021-10-27 2022-03-01 西安交通大学 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
刘红飞;杜宏伟;: "一种改进的最小方差自适应波束形成超声成像方法", 北京生物医学工程, no. 02, 15 April 2014 (2014-04-15) *

Also Published As

Publication number Publication date
CN116115260B (zh) 2024-05-24

Similar Documents

Publication Publication Date Title
CN103536316B (zh) 一种空时平滑相干因子类自适应超声成像方法
EP0473959B1 (en) Method of transforming a multi-beam sonar image
JP6408297B2 (ja) ビームフォーミング方法、計測イメージング装置、及び、通信装置
CN104272134B (zh) 超声成像系统中的杂波抑制
US20150359512A1 (en) Synthetic aperture ultrasound system
CN109431536B (zh) 一种聚焦超声空化的实时高分辨时空分布成像方法与系统
Nair et al. Robust short-lag spatial coherence imaging
US11737733B2 (en) Method of, and apparatus for, determination of position in ultrasound imaging
WO2019145785A1 (en) Methods and instrumentation for estimation of wave propagation and scattering parmeters
JP2019535448A (ja) 超音波画像クラッタをフィルタリングする方法及びシステム
CN108836389B (zh) 平面波相关点相干自适应波束合成成像方法
KR20140012043A (ko) 영상 획득속도 최적화를 구비한 영상장치
US11857367B2 (en) Ultrasound method and apparatus
CN109513123B (zh) 一种基于半球阵的高分辨三维被动空化成像方法
CN112754527B (zh) 一种用于低频超声胸腔成像的数据处理方法
JP2014512923A (ja) オーバーラップする送信ビームにおける適格と評価された領域を使用する増強された超音波画像形成
CN111265245B (zh) 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
US20050154306A1 (en) Dort process-based method and system for adaptive beamforming in estimating the aberration in a medium
US11885917B2 (en) Methods and instrumentation for estimation of wave propagation and scattering parameters
Ramkumar et al. Compressed sensing approach for reducing the number of receive elements in synthetic transmit aperture imaging
CN116115260B (zh) 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统
CN117147694A (zh) 基于逆问题的超声全聚焦成像稀疏正则化重构方法及设备
EP4194893A1 (en) Method and system for processing beamformed data
CN114098799B (zh) 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统
EP4154031A1 (en) Methods and instrumentation for estimation of wave propagation and scattering parameters

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