CN114098799A - 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统 - Google Patents

一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统 Download PDF

Info

Publication number
CN114098799A
CN114098799A CN202111263927.0A CN202111263927A CN114098799A CN 114098799 A CN114098799 A CN 114098799A CN 202111263927 A CN202111263927 A CN 202111263927A CN 114098799 A CN114098799 A CN 114098799A
Authority
CN
China
Prior art keywords
matrix
cavitation
array element
array
elements
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
CN202111263927.0A
Other languages
English (en)
Other versions
CN114098799B (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 CN202111263927.0A priority Critical patent/CN114098799B/zh
Publication of CN114098799A publication Critical patent/CN114098799A/zh
Application granted granted Critical
Publication of CN114098799B publication Critical patent/CN114098799B/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
    • 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
    • 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)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • General Health & Medical Sciences (AREA)
  • Biomedical Technology (AREA)
  • Animal Behavior & Ethology (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • General Physics & Mathematics (AREA)
  • Gynecology & Obstetrics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明公开了一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统：通过间隔抽取阵元分割全阵元孔径并分别针对每个M阵元孔径构造稀疏延时信息矩阵和被动空化射频信号扩展矩阵，将两个矩阵相乘得到每个M阵元孔径的波束合成信号矩阵并对其进行标准化处理，利用所有M阵元孔径的标准化波束合成信号矩阵计算互相关系数向量并对其元素进行阈值化及归一化处理后对全阵元孔径波束合成信号矩阵进行按行加权，根据相关系数加权矩阵计算得到由空化起始和终止时刻确定的空化起止时间段内多个时间窗的被动空化成像结果。本发明可提高被动空化成像的计算速度、抑制被动空化成像的干扰伪影、并对单个脉冲内的超声空化进行实时动态成像。

Description

一种单脉冲内超声空化的快速低伪影实时动态成像方法与 系统
技术领域
本发明属于超声检测与超声成像技术领域，具体涉及一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统。
背景技术
超声空化是指在聚焦或非聚焦超声换能器(即空化源激励超声换能器)发射声波激励下介质中微泡成核、膨胀、收缩以及坍塌的一系列物理过程，在高强度组织机械毁损、低强度药物控释等多种超声治疗应用中发挥着极为重要的作用。由于超声空化本身是随机不可控的，因此需要发展有效的成像手段监控空化过程，以避免超声空化对组织产生的意外损伤。被动空化成像是一种通过多阵元超声换能器各个阵元同步被动接收空化源辐射信号并对其进行处理从而获得空化空间分布的成像技术，该成像技术可在空化源激励超声换能器作用的过程中对空化进行实时成像，具有广泛的临床应用前景。
在被动空化成像过程中，首先在某一像素坐标处对多阵元超声换能器接收的被动空化射频信号进行波束合成处理，然后经过时间积分得到该像素坐标处的能量值，最后遍历成像区域内所有像素坐标得到一幅图像；即该成像技术是基于逐像素波束合成而实现的。当成像区域内像素数目较少时，逐像素波束合成方法所需的计算时间尚可接受。而超声空化往往需要在一个像素间距较小的较大成像区域中进行观察，此时逐像素波束合成方法会耗费大量计算时间，使得计算速度受到极大限制。此外，当被动空化射频信号的采样点数较多或信号采样频率较高时，被动空化成像的计算速度也会进一步下降。
目前，被动空化成像技术的实现普遍使用的是一种基于延时叠加的波束合成方法，其性能主要取决于空化源位置、空化源频率以及多阵元超声换能器孔径。对于一些大孔径超声换能器，如用于脑部超声成像的半球阵超声换能器，该方法可获得较好的成像效果。然而，为便于临床应用，多阵元超声换能器往往采用线阵换能器等小孔径超声换能器，此时该方法会沿着换能器轴向产生高水平的“X”状干扰伪影。此外，多阵元超声换能器缺陷、生物组织异质性以及空化源相互作用等因素会进一步加重干扰伪影，从而导致被动空化成像技术监控空化过程的准确性和可靠性难以保证。
在空化源激励超声换能器作用的过程中，除了相邻激励脉冲间超声空化的相互影响外，还应关注单个激励脉冲内部超声空化在时间和空间上的动力学演变过程，这有助于研究多种超声空化介导的应用背后的更深层次的物理机制。在连续的多个脉冲作用下，通过对各个脉冲下多阵元超声换能器接收的被动空化射频信号进行处理，可获得各个脉冲下的空化图像，从而对多个脉冲作用下超声空化的变化过程进行观察。然而，每个脉冲下的空化图像只能反映该脉冲整个作用过程中超声空化的累计空间分布，而无法对该脉冲作用过程中超声空化的时空变化过程进行动态成像。
尽管目前已发展出一些改进的被动空化成像方法，但大多数方法采用的依然是逐像素计算模式，在抑制干扰伪影的同时不可避免地进一步降低了计算速度；少数方法能够提高计算速度，但对于抑制干扰伪影并无贡献；同时，目前并无有效的方法对单脉冲内超声空化进行准确的动态成像。鉴于此，亟待提出一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统。
发明内容
本发明的目的在于提供一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统。
为了实现上述目的，本发明采用了以下技术方案：
一种单脉冲内超声空化的快速低伪影实时动态成像方法，包括以下步骤：
1)根据多阵元超声换能器所在位置建立被动空化成像坐标系，在该成像坐标系下规划被动空化成像区域，根据由该成像区域内任意一像素到多阵元超声换能器任意一阵元的距离构成的像素阵元距离集合中的最大距离和最小距离计算临界采样点数；
2)通过间隔抽取阵元的方式对多阵元超声换能器的全阵元孔径进行分割，得到多个M阵元孔径，针对每个M阵元孔径初始化延时信息矩阵，将根据所述被动空化成像区域内任意一像素到M阵元孔径各阵元的距离计算的相对延时转化为延时采样点数并对延时采样点数进行基准化处理，得到基准化延时采样点数，分别计算在所述延时信息矩阵中像素坐标的行索引和基准化延时采样点数的列索引并依据该行、列索引构造延时信息矩阵，对延时信息矩阵进行稀疏化处理，得到每个M阵元孔径的稀疏延时信息矩阵；
3)根据所述临界采样点数对多阵元超声换能器被动接收所得的被动空化射频信号矩阵进行前后向补零处理，得到全阵元孔径补零信号矩阵，对全阵元孔径补零信号矩阵按照所述间隔抽取阵元的方式进行分割，得到每个M阵元孔径的补零信号矩阵，将该补零信号矩阵转化为列向量并由该列向量通过截取元素生成多个子列向量，利用生成的子列向量构造每个M阵元孔径的被动空化射频信号扩展矩阵；
4)将所述每个M阵元孔径的稀疏延时信息矩阵和被动空化射频信号扩展矩阵相乘得到波束合成信号矩阵并对该波束合成信号矩阵进行标准化处理，得到对应M阵元孔径的标准化波束合成信号矩阵，将所有M阵元孔径的标准化波束合成信号矩阵分为两组并分别进行叠加，得到两个叠加矩阵，对两个叠加矩阵进行点乘，得到相关性矩阵并计算互相关系数向量，对互相关系数向量中的元素进行阈值化及归一化处理并利用处理所得结果对由所有M阵元孔径的波束合成信号矩阵叠加所得的全阵元孔径波束合成信号矩阵进行按行加权，得到相关系数加权矩阵；
5)对所述相关系数加权矩阵的各元素进行平方，得到功率矩阵，对功率矩阵中的元素沿着行方向在由空化起始时刻和空化终止时刻确定的空化起止时间段经分段所得的时间窗上进行叠加，得到对应时间窗的能量向量，将该能量向量转化为能量矩阵并进行转置，得到对应时间窗的被动空化成像结果。
优选的，所述步骤1)中，被动空化成像坐标系的x轴为多阵元超声换能器的阵元方向，被动空化成像坐标系的z轴为与阵元方向垂直的方向；像素阵元距离集合中的最大距离为所述被动空化成像区域中x轴坐标最小且z轴坐标最大的像素到多阵元超声换能器末阵元的距离和所述被动空化成像区域中x轴坐标最大且z轴坐标最大的像素到多阵元超声换能器首阵元的距离中的最大值，像素阵元距离集合中的最小距离为所述被动空化成像区域中x轴坐标最小且z轴坐标最小的像素的z轴坐标。
优选的，所述步骤1)中，临界采样点数的计算公式表示为：
NCritSamp＝ceil(Dmax/c×Fs)-floor(Dmin/c×Fs)+1
其中，NCritSamp为临界采样点数，ceil(·)表示向上取整，floor(·)表示向下取整，Dmax和Dmin分别为像素阵元距离集合中的最大距离和最小距离，c为超声波在介质中的传播速度，Fs为被动空化射频信号采样频率。
优选的，所述步骤2)中，间隔抽取阵元的方式为：自多阵元超声换能器的第s个阵元开始，每隔NE/M个阵元抽取1个阵元，构成第s个M阵元孔径；其中s＝1,2,...,NE/M，NE为多阵元超声换能器的阵元数目，M为M阵元孔径的阵元数目。
优选的，所述步骤2)中，延时采样点数的计算公式表示为：
Figure BDA0003324784670000031
其中，
Figure BDA0003324784670000041
为根据被动空化成像区域内坐标为
Figure BDA0003324784670000042
的像素到第s个M阵元孔径中坐标为
Figure BDA0003324784670000043
的第is个阵元的距离计算的相对延时所对应的延时采样点数，round[·]表示四舍五入取整，is＝1,2,...,M，s＝1,2,...,NE/M，j＝1,2,...,NPx，k＝1,2,...,NPz，NPx和NPz分别为被动空化成像区域沿x轴和沿z轴的像素数目。
优选的，所述步骤2)中，基准化延时采样点数利用根据像素阵元距离集合中的最小距离计算的最小相对延时对应的最小延时采样点数来计算，基准化延时采样点数的计算公式表示为：
Figure BDA0003324784670000044
其中，is＝1,2,...,M，s＝1,2,...,NE/M，j＝1,2,...,NPx，k＝1,2,...,NPz。
优选的，所述步骤2)中，延时信息矩阵的构造包括以下步骤：当第s个M阵元孔径的延时信息矩阵As的第m行的列索引n为nD时，第m行第n列对应的元素值
Figure BDA0003324784670000045
为1，否则为0；其中，m为坐标为
Figure BDA0003324784670000046
的像素在延时信息矩阵As中的行索引，nD为第s个M阵元孔径的第is个阵元的基准化延时采样点数在延时信息矩阵As中的列索引：
m＝j+(k-1)NPx
Figure BDA0003324784670000047
其中，j＝1,2,...,NPx，k＝1,2,...,NPz，is＝1,2,...,M，s＝1,2,...,NE/M。
优选的，所述步骤2)中，对延时信息矩阵的稀疏化处理包括以下步骤：将延时信息矩阵中的所有零元素去除而仅存储所有的非零元素及其对应的索引。
优选的，所述步骤3)中，被动空化射频信号矩阵是通过使多阵元超声换能器各个阵元同步被动接收(即此时各个阵元不发射信号)超声激励所产生的空化源的辐射信号而得到的，多阵元超声换能器一般为线阵超声换能器等诊断用超声换能器中的一种或多种。
优选的，所述步骤3)中，按照下述公式对被动空化射频信号矩阵进行前后向补零处理：
Figure BDA0003324784670000051
其中，
Figure BDA0003324784670000052
为全阵元孔径补零信号矩阵
Figure BDA0003324784670000053
中的元素，
Figure BDA0003324784670000054
为NE行(NAcqSamp+2NCritSamp-2)列的矩阵，RFi,u为被动空化射频信号矩阵RF中的元素，RF为NE行NAcqSamp列的矩阵，i＝1,2,...,NE，v＝1,2,...,NAcqSamp+2NCritSamp-2，NAcqSamp为被动空化射频信号采样点数。
优选的，所述步骤3)中，截取元素包括以下步骤：对由第s个M阵元孔径的补零信号矩阵转化得到的M×(NAcqSamp+2NCritSamp-2)行1列的列向量bs，自列向量bs第1个元素开始截取M×NCritSamp个元素得到第1个子列向量
Figure BDA0003324784670000055
自第M+1个元素开始截取M×NCritSamp个元素得到第2个子列向量
Figure BDA0003324784670000056
直到自第M×(NAcqSamp+NCritSamp-2)+1个元素开始截取M×NCritSamp个元素得到第NAcqSamp+NCritSamp-1个子列向量
Figure BDA0003324784670000057
从而得到第s个M阵元孔径对应的NAcqSamp+NCritSamp-1个M×NCritSamp行1列的子列向量
Figure BDA0003324784670000058
该子列向量中的元素
Figure BDA0003324784670000059
为列向量bs中的元素
Figure BDA00033247846700000510
其中，s＝1,2,...,NE/M，ti＝1,2,...,NAcqSamp+NCritSamp-1，n＝1,2,...,M×NCritSamp。
优选的，所述步骤3)中，M阵元孔径的被动空化射频信号扩展矩阵为：
Figure BDA00033247846700000511
其中，
Figure BDA00033247846700000512
为第s个M阵元孔径的被动空化射频信号扩展矩阵，该矩阵的行数为M×NCritSamp，列数为NAcqSamp+NCritSamp-1；
Figure BDA00033247846700000513
分别为从第s个M阵元孔径的补零信号矩阵转化得到的列向量中通过截取元素得到的第1个、第2个、...、第NAcqSamp+NCritSamp-1个子列向量。
优选的，所述步骤4)中，对波束合成信号矩阵的标准化处理包括以下步骤：对于第s个M阵元孔径的波束合成信号矩阵Qs的每一行，首先将该行元素减去该行元素的均值得到该行的去均值化元素，然后将该行的去均值化元素除以该行的去均值化元素的平方和的平方根，得到第s个M阵元孔径的标准化波束合成信号矩阵
Figure BDA0003324784670000061
即
Figure BDA0003324784670000062
中的元素
Figure BDA0003324784670000063
的计算公式表示为：
Figure BDA0003324784670000064
其中，
Figure BDA0003324784670000065
为第s个M阵元孔径的波束合成信号矩阵Qs中的元素，s＝1,2,...,NE/M，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1。
优选的，所述步骤4)中，将所有M阵元孔径的标准化波束合成信号矩阵分为前半组和后半组两组，前半组由
Figure BDA0003324784670000066
组成，后半组由
Figure BDA0003324784670000067
组成。
优选的，所述步骤4)中，互相关系数向量中的元素的计算公式表示为：
Figure BDA0003324784670000068
其中，
Figure BDA0003324784670000069
为相关性矩阵R中的元素，m＝1,2,..,NPx×NPz。
优选的，所述步骤4)中，对互相关系数向量中的元素的阈值化处理包括以下步骤：当互相关系数向量中的某一元素小于阈值时，将该元素置为阈值，否则保持不变；阈值一般小于等于0.001。
优选的，所述步骤4)中，相关系数加权矩阵中的元素的计算公式表示为：
Figure BDA00033247846700000610
Figure BDA00033247846700000611
其中，
Figure BDA00033247846700000612
为全阵元孔径波束合成信号矩阵Q中的元素，
Figure BDA00033247846700000613
为归一化互相关系数向量
Figure BDA00033247846700000614
中的元素，
Figure BDA00033247846700000615
为对互相关系数向量中的元素进行阈值化处理所得阈值化互相关系数向量
Figure BDA00033247846700000616
中的元素，max{·}表示取最大值，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1。
优选的，所述步骤5)中，空化起始时刻和空化终止时刻的计算包括以下步骤：
5.1)对所述相关系数加权矩阵进行绝对值化处理，将处理所得结果转化为列向量；
5.2)将步骤5.1)所得列向量中的元素降序排列，然后将所得各元素的索引按照相关系数加权矩阵的行数和列数转化为矩阵行索引和矩阵列索引；
5.3)对由步骤5.2)中所得矩阵行索引构成的矩阵行索引向量中的元素进行去重复处理，即对于重复出现的某一元素只保留第一次出现的该元素，从去重复处理后得到的矩阵行索引向量中根据多空化源像素数目提取多空化源行索引；
5.4)对每个多空化源行索引在所述相关系数加权矩阵中对应的行向量进行希尔伯特变换及绝对值化处理，得到对应的包络行向量；
5.5)对所得包络行向量中的元素进行阈值化处理，得到阈值化包络行向量，寻找阈值化包络行向量的首非零元素索引(即阈值化包络行向量中第一个非零元素的索引)和末非零元素索引(即阈值化包络行向量中最后一个非零元素的索引)；
5.6)根据所有多空化源行索引对应的首非零元素索引的最小值确定空化起始时刻，根据所有多空化源行索引对应的末非零元素索引的最大值确定空化终止时刻。
优选的，所述步骤5.5)中，对包络行向量中的元素的阈值化处理包括以下步骤：当包络行向量中的某一元素小于阈值时，将该元素置为零，否则保持不变；阈值一般小于包络行向量中所有元素的最大值的0.5倍。
一种单脉冲内超声空化的快速低伪影实时动态成像系统，该系统包括临界采样点数计算模块、M阵元孔径稀疏延时信息矩阵构造模块、M阵元孔径被动空化射频信号扩展矩阵构造模块、M阵元孔径波束合成及互相关加权模块和单脉冲内超声空化实时动态成像模块；
所述临界采样点数计算模块用于执行上述步骤1)，主要是用于根据多阵元超声换能器所在位置建立被动空化成像坐标系并在该成像坐标系下规划被动空化成像区域，以及根据由该成像区域内任意一像素到多阵元超声换能器任意一阵元的距离构成的像素阵元距离集合中的最大距离和最小距离计算临界采样点数；
所述M阵元孔径稀疏延时信息矩阵构造模块用于执行上述步骤2)，主要是用于通过间隔抽取阵元的方式对多阵元超声换能器的全阵元孔径进行分割、针对分割所得的每个M阵元孔径初始化延时信息矩阵、将根据所述临界采样点数计算模块规划的被动空化成像区域内任意一像素到M阵元孔径各阵元的距离计算的相对延时转化为延时采样点数并对延时采样点数进行基准化处理、分别计算在所述延时信息矩阵中像素坐标的行索引和基准化延时采样点数的列索引，以及依据该行、列索引构造延时信息矩阵并对每个M阵元孔径的延时信息矩阵进行稀疏化处理；
所述M阵元孔径被动空化射频信号扩展矩阵构造模块用于执行上述步骤3)，主要是用于根据所述临界采样点数计算模块得到的临界采样点数对多阵元超声换能器被动接收所得的被动空化射频信号矩阵进行前后向补零处理、对处理所得全阵元孔径补零信号矩阵按照所述M阵元孔径稀疏延时信息矩阵构造模块中的间隔抽取阵元的方式进行分割、将分割所得每个M阵元孔径的补零信号矩阵转化为列向量并由该列向量通过截取元素生成多个子列向量，以及利用生成的子列向量构造每个M阵元孔径的被动空化射频信号扩展矩阵；
所述M阵元孔径波束合成及互相关加权模块用于执行上述步骤4)，主要是用于将所述M阵元孔径稀疏延时信息矩阵构造模块经稀疏化处理得到的每个M阵元孔径的稀疏延时信息矩阵和所述M阵元孔径被动空化射频信号扩展矩阵构造模块得到的每个M阵元孔径的被动空化射频信号扩展矩阵对应相乘、对相乘所得每个M阵元孔径的波束合成信号矩阵进行标准化处理、将处理所得所有M阵元孔径的标准化波束合成信号矩阵分为两组并分别进行叠加、对叠加所得两个叠加矩阵进行点乘并计算互相关系数向量，以及对互相关系数向量中的元素进行阈值化及归一化处理并利用处理所得结果对由所有M阵元孔径的波束合成信号矩阵叠加所得的全阵元孔径波束合成信号矩阵进行按行加权；
所述单脉冲内超声空化实时动态成像模块用于执行上述步骤5)，主要是用于对所述M阵元孔径波束合成及互相关加权模块经对全阵元孔径波束合成信号矩阵进行按行加权得到的相关系数加权矩阵的各元素进行平方、对平方后所得功率矩阵中的元素沿着行方向在由空化起始时刻和空化终止时刻确定的空化起止时间段经分段所得的时间窗上进行叠加，以及将叠加所得能量向量转化为能量矩阵并进行转置，从而得到对应时间窗的被动空化成像结果。
优选的，所述系统还包括空化起始和终止时刻自适应计算模块，用于执行上述步骤5.1)～5.6)。
本发明的有益效果体现在：
本发明提出的单脉冲内超声空化的快速低伪影实时动态成像方法，通过间隔抽取阵元分割全阵元孔径并分别将针对每个M阵元孔径构造的稀疏延时信息矩阵和被动空化射频信号扩展矩阵相乘得到对应的波束合成信号矩阵，利用所有M阵元孔径的标准化波束合成信号矩阵计算互相关系数向量并经阈值化及归一化处理后对全阵元孔径波束合成信号矩阵进行按行加权，根据相关系数加权矩阵计算得到空化起止时间段内多个时间窗的被动空化成像结果，从而可提高被动空化成像的计算速度、抑制被动空化成像的干扰伪影、并对单个脉冲内的超声空化进行实时动态成像。本发明可对单个超声激励脉冲内空化的时空动力学演变过程进行快速、低伪影、实时动态的成像，为监控高强度组织机械毁损、低强度药物控释等多种超声空化介导的应用、以及在单脉冲内研究超声空化瞬态物理机制和控制超声空化瞬态物理过程提供了高性能的影像手段。
本发明根据像素阵元距离集合中的最大距离和最小距离计算临界采样点数并据此分别确定稀疏延时信息矩阵和被动空化射频信号扩展矩阵的规模，并对延时采样点数进行基准化处理，可构造出最小的稀疏延时信息矩阵和被动空化射频信号扩展矩阵，将矩阵相乘的计算复杂度优化至最小，从而可在不影响被动空化成像结果的前提下有效地降低计算量，使得被动空化成像的计算速度得到提高。
本发明在构造被动空化射频信号扩展矩阵的过程中根据临界采样点数对被动空化射频信号矩阵进行前后向补零处理，一方面避免了不进行补零或补零不足时不能充分利用被动空化射频信号以至于被动空化成像结果错误的问题，另一方面避免了补零过多时被动空化射频信号扩展矩阵过大以至于计算量增加的问题。
本发明通过构造稀疏延时信息矩阵和被动空化射频信号扩展矩阵并将两个矩阵相乘的方法对每个M阵元孔径进行波束合成处理，可直接一次性得到包含了所有像素坐标的波束合成信号的波束合成信号矩阵，使得波束合成处理所需的时间相比传统的逐像素波束合成方法大幅度减小，从而能够以更快的计算速度进行被动空化成像。
本发明通过对每个M阵元孔径的波束合成信号矩阵进行标准化处理和对分组所得两个叠加矩阵进行点乘的方法得到相关性矩阵并进一步得到互相关系数向量，其中基于矩阵的运算过程可一次性得到所有像素坐标的互相关系数、分组叠加则避免了对M阵元孔径的波束合成信号彼此之间进行互相关处理的遍历循环操作，这使得计算量有效降低，从而提高了被动空化成像的计算速度。
本发明利用M阵元孔径的波束合成信号之间的相关性在空化源处高而在非空化源处低的特点，通过间隔抽取阵元分割全阵元孔径并利用互相关系数向量对全阵元孔径波束合成信号矩阵加权，有效地突显出空化源与非空化源间的能量差异，从而可抑制被动空化成像中的干扰伪影；且在计算互相关系数向量的过程中利用了所有M阵元孔径的波束合成信号矩阵，充分利用了这一相关性特点，使得干扰伪影的抑制更为鲁棒。
进一步地，本发明对延时信息矩阵进行稀疏化处理，一方面仅存储延时信息矩阵中非零元素及其索引，从而节省了数据存储空间；另一方面避免了与延时信息矩阵中零值元素对应的运算过程，使得计算量有效降低，从而提高了被动空化成像的计算速度。
进一步地，本发明中稀疏延时信息矩阵的构造取决于空化成像区域、多阵元超声换能器参数和信号采样频率，而这些参数在一个特定的应用情境中是固定的，因此在特定应用情境中该矩阵只需构造一次即可重复利用，为不同超声空化源激励条件下(如不同的超声激励声压、超声激励脉冲长度及重复频率等)的被动空化成像提供了极大的便利。
进一步地，本发明中多阵元超声换能器的各个阵元工作在同步被动接收模式，避免了空化源激励超声换能器发射信号的干扰，从而可在空化发生的过程中对空化进行实时成像。
进一步地，本发明在构造被动空化射频信号扩展矩阵的过程中将M阵元孔径的补零信号矩阵转化为列向量并从中截取M阵元孔径的阵元数与临界采样点数的乘积个元素，有效地减少了构造被动空化射频信号扩展矩阵所需的时间，从而提高了被动空化成像的计算速度。
进一步地，本发明对互相关系数向量中的元素进行阈值化处理，消除了互相关系数向量中的负值元素向被动空化成像结果中引入的干扰伪影，从而可使得被动空化成像中的干扰伪影得到进一步的抑制。
进一步地，本发明在对互相关系数向量中的元素进行阈值化处理后进行归一化处理，可有效避免因换能器缺陷、组织异质性和空化源相互作用等复杂因素造成空化源坐标处波束合成信号相关性降低、进而导致空化能量估计错误的问题。
进一步地，本发明通过相关系数加权矩阵的绝对值化处理及列向量转化、列向量中元素的降序排列、矩阵行索引向量中元素的去重复处理、相关系数加权矩阵中多空化源行索引对应行向量的希尔伯特变换、包络行向量中元素的阈值化处理等步骤，自适应地计算出空化起始和终止时刻并据此确定空化起止时间段，通过对该时间段进行分段得到多个时间窗，从而可对单脉冲内的超声空化进行动态成像。
附图说明
图1为本发明实施例中被动空化成像坐标系和被动空化成像区域的示意图。
图2为本发明实施例中构造M阵元孔径的稀疏延时信息矩阵的流程图。
图3为本发明实施例中构造M阵元孔径的被动空化射频信号扩展矩阵的流程图。
图4为本发明实施例中计算相关系数加权矩阵的流程图。
图5为本发明实施例中计算空化起始和终止时刻的流程图。
图6为本发明实施例中计算空化起止时间段内时间窗的被动空化成像结果的流程图。
图7为本发明实施例中计算相关系数加权矩阵所需的计算时间与传统方法所需的计算时间的比较，其中：传统方法是指传统的逐像素波束合成方法，本发明所用方法是指本发明中所用的基于稀疏延时信息矩阵和被动空化射频信号扩展矩阵相乘的波束合成方法。
图8为本发明实施例中分别使用两种方法得到的8个时间窗内的被动空化成像结果，其中：(a)传统的延时叠加波束合成方法；(b)本发明中所用的基于M阵元孔径波束合成及互相关加权的方法；8个时间窗分别为0～5μs、5～10μs、10～15μs、15～20μs、20～25μs、25～30μs、30～35μs和35～40μs。
具体实施方式
下面结合附图和实施例对本发明作进一步的详细说明。
(1.1)参见图1，根据多阵元超声换能器(例如，线阵换能器)所在位置建立被动空化成像坐标系(以下简称为成像坐标系)，多阵元超声换能器的阵元方向为成像坐标系的x轴，与阵元方向垂直的方向为成像坐标系的z轴，多阵元超声换能器的中心为成像坐标系的原点；
(1.2)在步骤(1.1)所述成像坐标系下多阵元超声换能器的各阵元坐标为
Figure BDA0003324784670000111
其中，i＝1,2,...,NE，NE为多阵元超声换能器的阵元数目(例如，128)；将i＝1和i＝NE时对应的阵元分别定义为首阵元和末阵元，其坐标分别为
Figure BDA0003324784670000112
和
Figure BDA0003324784670000113
(1.3)在步骤(1.1)所述成像坐标系下规划被动空化成像区域(以下简称为成像区域)，该区域中共包含NPx×NPz个像素，像素坐标表示为
Figure BDA0003324784670000114
其中j＝1,2,...,NPx，k＝1,2,...,NPz，NPx和NPz分别为沿x轴和沿z轴的像素数目；
将所述成像区域中x轴坐标最小且z轴坐标最小的像素定义为左上边界像素，其坐标为
Figure BDA0003324784670000115
将所述成像区域中x轴坐标最小且z轴坐标最大的像素定义为左下边界像素，其坐标为
Figure BDA0003324784670000116
将所述成像区域中x轴坐标最大且z轴坐标最大的像素定义为右下边界像素，其坐标为
Figure BDA0003324784670000117
(1.4)步骤(1.3)所述成像区域内任意一像素到多阵元超声换能器任意一阵元的距离构成了像素阵元距离集合，该集合中共有NPx×NPz×NE个元素；
该集合中的最小距离Dmin为步骤(1.3)所述左上边界像素的z轴坐标
Figure BDA0003324784670000121
该集合中的最大距离Dmax为步骤(1.3)所述左下边界像素到步骤(1.2)所述末阵元的距离DLD-L和步骤(1.3)所述右下边界像素到步骤(1.2)所述首阵元的距离DRD-F的最大值：
Dmax＝max{DLD-L,DRD-F}
Figure BDA0003324784670000122
Figure BDA0003324784670000123
其中，max{·}表示取最大值；
(1.5)根据步骤(1.4)所述像素阵元距离集合中的最大距离Dmax与最小距离Dmin之差计算临界采样点数NCritSamp：
NCritSamp＝ceil(Dmax/c×Fs)-floor(Dmin/c×Fs)+1
其中，ceil(·)表示向上取整，floor(·)表示向下取整，c为超声波在介质中的传播速度(例如，1480m/s)，Fs为被动空化射频信号采样频率(例如，50MHz)。
参见图2，通过间隔抽取阵元分割全阵元孔径并针对所得的每个M阵元孔径初始化延时信息矩阵，计算延时采样点数并做基准化处理，根据基准化延时采样点数构造延时信息矩阵并进行稀疏化处理，得到每个M阵元孔径的稀疏延时信息矩阵，具体流程见下述步骤(2.1)～(2.9)；
(2.1)通过间隔抽取阵元的方式对全阵元孔径(含NE个阵元)进行分割，得到NE/M个具有M个阵元(例如，M＝8)的子孔径(以下称为M阵元孔径)；具体的间隔抽取阵元的方式为：自第s个阵元开始，每隔NE/M个阵元抽取1个阵元，构成第s个M阵元孔径；将第s个M阵元孔径对应的阵元索引记为is＝1,2,...,M，其中s＝1,2,...,NE/M；
(2.2)对于步骤(2.1)所述第s个M阵元孔径，初始化延时信息矩阵As；该矩阵的行数为步骤(1.1)所述x轴像素数目和z轴像素数目的乘积NPx×NPz，该矩阵的列数为步骤(2.1)所述M阵元孔径的阵元数目与步骤(1.5)所述临界采样点数的乘积M×NCritSamp；该矩阵中的元素为
Figure BDA0003324784670000131
其中s＝1,2,...,NE/M，m＝1,2,...,NPx×NPz，n＝1,2,...,M×NCritSamp，元素值均为零；
(2.3)根据步骤(1.3)所述成像区域内坐标为
Figure BDA0003324784670000132
的像素到步骤(2.1)所述第s个M阵元孔径中坐标为
Figure BDA0003324784670000133
的第is个阵元的距离计算对应的相对延时，并利用步骤(1.5)所述信号采样频率将其转化为延时采样点数
Figure BDA0003324784670000134
Figure BDA0003324784670000135
其中，round[·]表示四舍五入取整，is＝1,2,...,M，s＝1,2,...,NE/M，j＝1,2,...,NPx，k＝1,2,...,NPz；
(2.4)根据步骤(1.4)所得像素阵元距离集合中的最小距离Dmin计算最小相对延时，并根据步骤(1.5)所述信号采样频率Fs将其转化为最小延时采样点数；然后根据所得最小延时采样点数对步骤(2.3)所得延时采样点数
Figure BDA0003324784670000136
进行基准化处理，得到步骤(2.1)所述第s个M阵元孔径的第is个阵元的基准化延时采样点数
Figure BDA0003324784670000137
Figure BDA0003324784670000138
其中，is＝1,2,...,M，s＝1,2,...,NE/M，j＝1,2,...,NPx，k＝1,2,...,NPz；
(2.5)分别计算在步骤(2.2)所述延时信息矩阵As中步骤(1.3)所述坐标为
Figure BDA0003324784670000139
的像素的行索引m，和步骤(2.4)所得第is个阵元的基准化延时采样点数
Figure BDA00033247846700001310
的列索引nD：
m＝j+(k-1)NPx
Figure BDA00033247846700001311
其中，j＝1,2,...,NPx，k＝1,2,...,NPz，is＝1,2,...,M，s＝1,2,...,NE/M；
(2.6)对于步骤(2.2)所述延时信息矩阵As，当第m行的列索引n为nD时，第m行第n列对应的元素值
Figure BDA0003324784670000141
为1，否则为0，即：
Figure BDA0003324784670000142
(2.7)重复步骤(2.3)～(2.6)，完成对延时信息矩阵As的构造，该矩阵为非零即一的矩阵；
(2.8)对步骤(2.7)所得延时信息矩阵As进行稀疏化处理，即将As中的所有零元素去除而仅存储所有的非零元素及其对应的索引，得到第s个M阵元孔径的稀疏延时信息矩阵并记为
Figure BDA0003324784670000143
(2.9)重复步骤(2.2)～(2.8)，得到NE/M个M阵元孔径对应的稀疏延时信息矩阵。
参见图3，根据临界采样点数对被动空化射频信号矩阵进行前后向补零处理，对全阵元孔径补零信号矩阵进行分割并将每个M阵元孔径的补零信号矩阵转化为列向量，通过截取元素构造每个M阵元孔径的被动空化射频信号扩展矩阵，具体流程见下述步骤(3.1)～(3.6)；
(3.1)使多阵元超声换能器各个阵元同步被动接收(即各个阵元不发射信号)超声激励所产生的空化源的辐射信号，得到NE行NAcqSamp列的被动空化射频信号矩阵RF，该矩阵中的元素为RFi,u，其中，i＝1,2,...,NE，u＝1,2,...,NAcqSamp，NAcqSamp为被动空化射频信号采样点数；基于多阵元超声换能器各个阵元同步被动接收的被动空化成像可在空化发生过程中对空化进行实时成像；
根据步骤(1.5)所得临界采样点数NCritSamp对RF进行前后向补零处理，得到NE行(NAcqSamp+2NCritSamp-2)列的全阵元孔径补零信号矩阵
Figure BDA0003324784670000144
该矩阵中的元素为：
Figure BDA0003324784670000145
其中，i＝1,2,...,NE，v＝1,2,...,NAcqSamp+2NCritSamp-2；
(3.2)按照步骤(2.1)所述间隔抽取阵元的方式将步骤(3.1)所得全阵元孔径补零信号矩阵
Figure BDA0003324784670000151
分割为NE/M个M阵元孔径的补零信号矩阵
Figure BDA0003324784670000152
该矩阵的行数为M，列数为NAcqSamp+2NCritSamp-2，该矩阵中的元素为
Figure BDA0003324784670000153
其中s＝1,2,...,NE/M，is＝1,2,...,M，v＝1,2,...,NAcqSamp+2NCritSamp-2；
(3.3)将步骤(3.2)所得第s个M阵元孔径的补零信号矩阵
Figure BDA0003324784670000154
转化为M×(NAcqSamp+2NCritSamp-2)行1列的列向量bs，该列向量中的元素为：
Figure BDA0003324784670000155
其中，p＝is+(v-1)M，p＝1,2,...,M×(NAcqSamp+2NCritSamp-2)，is＝1,2,...,M，s＝1,2,...,NE/M，v＝1,2,...,NAcqSamp+2NCritSamp-2；
(3.4)对于步骤(3.3)所得的列向量bs，通过截取元素生成多个子列向量：自第1个元素开始截取M×NCritSamp个元素得到第1个子列向量
Figure BDA0003324784670000156
自第M+1个元素开始截取M×NCritSamp个元素得到第2个子列向量
Figure BDA0003324784670000157
直到自第M×(NAcqSamp+NCritSamp-2)+1个元素开始截取M×NCritSamp个元素得到第NAcqSamp+NCritSamp-1个子列向量
Figure BDA0003324784670000158
从而得到NAcqSamp+NCritSamp-1个M×NCritSamp行1列的子列向量
Figure BDA0003324784670000159
该子列向量中的元素
Figure BDA00033247846700001510
为列向量bs中的元素
Figure BDA00033247846700001511
Figure BDA00033247846700001512
其中，s＝1,2,...,NE/M，ti＝1,2,...,NAcqSamp+NCritSamp-1，n＝1,2,...,M×NCritSamp；
(3.5)根据步骤(3.4)所得的NAcqSamp+NCritSamp-1个子列向量
Figure BDA00033247846700001513
来构造第s个M阵元孔径的被动空化射频信号扩展矩阵
Figure BDA00033247846700001514
Figure BDA00033247846700001515
所得矩阵的行数为M×NCritSamp，列数为NAcqSamp+NCritSamp-1；
(3.6)重复步骤(3.3)～(3.5)，得到NE/M个M阵元孔径对应的被动空化射频信号扩展矩阵。
参见图4，将每个M阵元孔径的稀疏延时信息矩阵与被动空化射频信号扩展矩阵相乘得到波束合成信号矩阵，对其进行标准化处理，计算互相关系数向量，对该向量中元素进行阈值化及归一化处理后对全阵元孔径波束合成信号矩阵进行按行加权，得到相关系数加权矩阵，具体流程见下述步骤(4.1)～(4.10)；
(4.1)将步骤(2.8)所得第s个M阵元孔径的稀疏延时信息矩阵
Figure BDA0003324784670000161
与步骤(3.6)所得第s个M阵元孔径的被动空化射频信号扩展矩阵
Figure BDA0003324784670000162
相乘(即波束合成处理)，得到第s个M阵元孔径的波束合成信号矩阵Qs：
Figure BDA0003324784670000163
该矩阵中的元素为：
Figure BDA0003324784670000164
其中，s＝1,2,...,NE/M，m＝1,2,...,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1；
(4.2)将步骤(4.1)所得第s个M阵元孔径的波束合成信号矩阵Qs进行标准化处理(即对于矩阵Qs的每一行，首先将该行元素减去该行元素的均值得到该行的去均值化元素，然后将该行的去均值化元素除以该行的去均值化元素的平方和的平方根)，从而得到第s个M阵元孔径的标准化波束合成信号矩阵
Figure BDA0003324784670000165
该矩阵中的元素为：
Figure BDA0003324784670000166
其中，s＝1,2,...,NE/M，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1；
(4.3)重复步骤(4.1)～(4.2)，得到NE/M个M阵元孔径的波束合成信号矩阵Qs以及NE/M个M阵元孔径的标准化波束合成信号矩阵
Figure BDA0003324784670000167
(4.4)将步骤(4.3)所得NE/M个M阵元孔径的标准化波束合成信号矩阵
Figure BDA0003324784670000168
分为前半组和后半组两组，组内矩阵数目均为NE/2M，前半组由
Figure BDA0003324784670000169
组成，后半组由
Figure BDA0003324784670000171
组成，分别对前半组和后半组的标准化波束合成信号矩阵进行叠加，得到前半组叠加矩阵
Figure BDA0003324784670000172
和后半组叠加矩阵
Figure BDA0003324784670000173
Figure BDA0003324784670000174
矩阵
Figure BDA0003324784670000175
和
Figure BDA0003324784670000176
中的元素分别为：
Figure BDA0003324784670000177
其中，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1；
(4.5)对步骤(4.4)所得前半组叠加矩阵
Figure BDA0003324784670000178
和后半组叠加矩阵
Figure BDA0003324784670000179
进行点乘(即各元素对应相乘)，得到相关性矩阵R，该矩阵中的元素为：
Figure BDA00033247846700001710
其中，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1；
(4.6)将步骤(4.5)所得相关性矩阵R沿着行方向进行叠加得到相关性向量，该向量为NPx×NPz行1列的列向量；将所得相关性向量除以步骤(4.4)所述组内矩阵数目NE/2M的平方，从而得到互相关系数向量r，该向量中的元素为：
Figure BDA00033247846700001711
其中，m＝1,2,..,NPx×NPz；
(4.7)对步骤(4.6)所得互相关系数向量r中的元素进行阈值化处理，得到阈值化互相关系数向量
Figure BDA00033247846700001712
该向量中的元素为：
Figure BDA00033247846700001713
其中，m＝1,2,..,NPx×NPz，ε为阈值(例如，0.001)；
(4.8)对步骤(4.7)所得阈值化互相关系数向量
Figure BDA00033247846700001714
中的元素以
Figure BDA00033247846700001715
中所有元素的最大值为基准进行归一化处理，得到归一化互相关系数向量
Figure BDA00033247846700001716
该向量中的元素为：
Figure BDA0003324784670000181
其中，m＝1,2,..,NPx×NPz；
(4.9)将步骤(4.3)所得NE/M个M阵元孔径的波束合成信号矩阵Qs叠加得到全阵元孔径波束合成信号矩阵Q：
Figure BDA0003324784670000182
该矩阵中的元素为：
Figure BDA0003324784670000183
其中，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1；
(4.10)利用步骤(4.8)所得归一化互相关系数向量
Figure BDA0003324784670000184
对步骤(4.9)所得全阵元孔径波束合成信号矩阵Q进行按行加权(即将归一化互相关系数向量
Figure BDA0003324784670000185
每一行的元素与全阵元孔径波束合成信号矩阵Q对应行的元素相乘)，得到相关系数加权矩阵
Figure BDA0003324784670000186
该矩阵中的元素为：
Figure BDA0003324784670000187
其中，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1。
参见图5，对相关系数加权矩阵进行绝对值化处理后转化为列向量并对其元素进行降序排列，对矩阵行索引向量中的元素进行去重复处理并根据计算的多空化源像素数目得到多空化源行索引，基于希尔伯特变换计算其对应的包络行向量并对该向量中的元素进行阈值化处理，通过寻找首非零元素索引和末非零元素索引得到空化起始和终止时刻，具体流程见下述步骤(5.1)～(5.9)；
(5.1)对步骤(4.10)所得相关系数加权矩阵
Figure BDA0003324784670000188
进行绝对值化处理，并将所得结果转化为(NPx×NPz)×(NAcqSamp+NCritSamp-1)行1列的列向量w，该列向量中的元素为：
Figure BDA0003324784670000189
其中，|·|表示取绝对值，l＝m+(ti-1)×(NPx×NPz)，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1，l＝1,2,...,(NPx×NPz)×(NAcqSamp+NCritSamp-1)；
(5.2)对步骤(5.1)所得列向量w中的元素进行降序排列得到列向量
Figure BDA0003324784670000191
然后将所得列向量
Figure BDA0003324784670000192
中各元素的索引按照相关系数加权矩阵
Figure BDA0003324784670000193
的行数NPx×NPz和列数NAcqSamp+NCritSamp-1转化为矩阵行索引和矩阵列索引，
Figure BDA0003324784670000194
中所有元素对应的矩阵行索引构成了(NPx×NPz)×(NAcqSamp+NCritSamp-1)行1列的矩阵行索引向量IRow；
(5.3)对步骤(5.2)所得的矩阵行索引向量IRow中的元素进行去重复处理，即若某一元素重复出现，则只保留第一次出现的该元素而去除后续重复出现的该元素，从而得到(NPx×NPz)行1列的矩阵行索引向量
Figure BDA0003324784670000195
(5.4)计算多空化源像素数目NPMCAV，从步骤(5.3)得到的(NPx×NPz)行1列的矩阵行索引向量
Figure BDA0003324784670000196
中提取前NPMCAV个元素，得到NPMCAV个多空化源行索引，即
Figure BDA0003324784670000197
步骤(5.4)所述多空化源像素数目NPMCAV根据下述步骤(5.4.1)～(5.4.4)计算：
(5.4.1)计算被动空化成像的点扩散函数在给定z轴坐标处的x轴半高宽HWx和z轴半高宽HWz：
Figure BDA0003324784670000198
Figure BDA0003324784670000199
其中，zG为给定z轴坐标(例如，为z轴坐标最小值与最大值之和的一半，即
Figure BDA00033247846700001910
f0为空化源频率(例如，5MHz)，pitch为步骤(1.1)所述多阵元超声换能器的阵元间距(例如，0.3mm)；
(5.4.2)根据步骤(5.4.1)所得x轴半高宽HWx和z轴半高宽HWz计算点扩散面积AreaPS：
AreaPS＝π×HWx×HWz/4
(5.4.3)根据步骤(5.4.2)所得点扩散面积AreaPS计算点扩散面积对应的像素数目NPPS：
NPPS＝ceil[AreaPS/(Spcx×Spcz)]
其中，ceil[·]表示向上取整，Spcx为x轴像素间距(例如，0.15mm)，Spcz为z轴像素间距(例如，0.5mm)；
(5.4.4)将步骤(5.4.3)所得点扩散面积对应的像素数目NPPS与空化源数目NCAV(例如，2)相乘，得到多空化源像素数目NPMCAV：
NPMCAV＝NPPS×NCAV
(5.5)对步骤(5.4)所得的第mC个多空化源行索引在步骤(4.10)所得相关系数加权矩阵
Figure BDA0003324784670000201
中对应的1行(NAcqSamp+NCritSamp-1)列的行向量
Figure BDA0003324784670000202
进行希尔伯特变换及绝对值化处理，得到第mC个多空化源行索引的包络行向量
Figure BDA0003324784670000203
Figure BDA0003324784670000204
其中，IU为虚数单位，hilbert(·)表示希尔伯特变换，mC＝1,2,...,NPMCAV；
(5.6)对步骤(5.5)所得包络行向量
Figure BDA0003324784670000205
中的元素
Figure BDA0003324784670000206
进行阈值化处理，得到阈值化包络行向量
Figure BDA0003324784670000207
该行向量中的元素为：
Figure BDA0003324784670000208
其中，δ为阈值(例如，包络行向量
Figure BDA0003324784670000209
中所有元素的最大值的0.3倍)，ti＝1,2,...,NAcqSamp+NCritSamp-1；
(5.7)寻找步骤(5.6)所得阈值化包络行向量
Figure BDA00033247846700002010
中第一个非零元素的索引
Figure BDA00033247846700002011
(即首非零元素索引)和最后一个非零元素的索引
Figure BDA00033247846700002012
(即末非零元素索引)；
(5.8)重复步骤(5.5)～(5.7)，得到NPMCAV个多空化源行索引对应的首非零元素索引和末非零元素索引；
(5.9)由步骤(5.8)所得的NPMCAV个多空化源行索引对应的首非零元素索引的最小值确定空化起始时刻TStart，由步骤(5.8)所得的NPMCAV个多空化源行索引对应的末非零元素索引的最大值确定空化终止时刻Tend。
参见图6，对空化起止时间段进行分段，将相关系数加权矩阵各元素平方所得功率矩阵中的元素沿着行方向在各个时间窗上进行叠加得到能量向量，将该向量转化为矩阵并进行转置得到对应时间窗的被动空化成像结果，具体流程见下述步骤(6.1)～(6.6)；
(6.1)按照一定的时间窗长Δt对由步骤(5.9)所得的空化起始时刻TStart和空化终止时刻Tend确定的空化起止时间段ti＝TStart,TStart+1,...,TEnd进行分段，得到
Figure BDA0003324784670000211
个时间窗，其中
Figure BDA0003324784670000212
floor[·]表示向下取整；第nSub个时间窗的时刻为ti＝TStart+(nSub-1)Δt,TStart+(nSub-1)Δt+1,...,TStart+nSubΔt-1，
Figure BDA0003324784670000213
(6.2)对步骤(4.10)所得相关系数加权矩阵
Figure BDA0003324784670000214
的各元素进行平方得到功率矩阵P，该矩阵中的元素为：
Figure BDA0003324784670000215
其中，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1；
(6.3)对步骤(6.2)所得功率矩阵P中的元素沿着行方向在步骤(6.1)所得的第nSub个时间窗上进行叠加，得到第nSub个时间窗的能量向量
Figure BDA0003324784670000216
该向量中的元素为：
Figure BDA0003324784670000217
其中，
Figure BDA0003324784670000218
m＝1,2,..,NPx×NPz；
(6.4)将步骤(6.3)所得第nSub个时间窗的能量向量
Figure BDA0003324784670000219
转化为能量矩阵
Figure BDA00033247846700002110
该矩阵中的元素为：
Figure BDA00033247846700002111
其中，
Figure BDA00033247846700002112
j＝1,2,...,NPx，k＝1,2,...,NPz，m＝j+(k-1)NPx；
(6.5)对步骤(6.4)所得能量矩阵
Figure BDA00033247846700002113
进行转置，得到第nSub个时间窗的被动空化成像结果；
(6.6)重复步骤(6.3)～(6.5)，得到
Figure BDA0003324784670000221
个时间窗的被动空化成像结果，从而实现单脉冲内超声空化的实时动态成像。
采用数值仿真方法验证本发明提出的单脉冲内超声空化的快速低伪影实时动态成像方法的性能，仿真参数如下所述：空化源频率为5MHz，脉冲长度为30μs；空化源数目为2，空化源的坐标分别为(-5,35mm)和(5,45mm)，位于(5,45mm)处的空化源相对位于(-5,35mm)处的空化源有10μs的延时；线阵换能器的阵元数目为128，阵元间距为0.3mm；超声波传播速度为1480m/s；信号采样频率为50MHz。根据上述参数可分别得到线阵换能器各个阵元对应的信号，其中有用信号包含在中间的信号段中(采样点数目为2700)。按照10dB的信噪比分别向各个阵元对应的信号段中加入高斯白噪声，则可得到前述被动空化射频信号矩阵。计算6次重复仿真所得的下述计算时间和伪影抑制率的均值和标准差，并表示为mean±std的形式，其中mean代表均值，std代表标准差。
下面分别对本发明提出的单脉冲内超声空化的快速低伪影实时动态成像方法的快速计算性能、伪影抑制性能以及单脉冲内实时动态成像功能进行分析。
(1)在计算相关系数加权矩阵的过程中，传统的逐像素波束合成方法和本发明中所用的基于稀疏延时信息矩阵和被动空化射频信号扩展矩阵相乘的波束合成方法所需的计算时间如图7所示。可以看出，传统方法所需的计算时间为148.9±0.7s；而本发明所用方法所需的计算时间仅为12.7±0.1s，计算速度相比传统方法提高10倍以上。上述分析表明本发明中所用的基于稀疏延时信息矩阵和被动空化射频信号扩展矩阵相乘的波束合成方法能够以更快的计算速度进行被动空化成像。
(2)传统的延时叠加波束合成方法和本发明中所用的基于M阵元孔径波束合成及互相关加权的方法得到的8个时间窗内的被动空化成像结果如图8(a)和图8(b)分别所示，8个时间窗分别为0～5μs，5～10μs，10～15μs，15～20μs，20～25μs，25～30μs，30～35μs和35～40μs，被动空化成像结果经过归一化和对数化处理后进行显示，显示范围为-30～0dB。可以看出，图8(a)中两个空化源周围均存在沿z轴延伸的高水平“X”状干扰伪影，严重影响图像质量；相比图8(a)，图8(b)中的“X”状干扰伪影得到明显抑制，图像质量明显提升。
利用伪影抑制率ASR定量评价对干扰伪影的抑制能力，ASR＝10lg(EC/EA)，其中EC和EA分别为被动空化成像结果中大于最大像素值一半的像素值的平均值和小于最大像素值一半的像素值的平均值。计算结果表明，图8(a)中8个时间窗对应的ASR分别为17.06±0.07dB，15.67±0.05dB，14.54±0.06dB，14.55±0.01dB，14.54±0.01dB，14.73±0.01dB，15.58±0.04dB和17.48±0.11dB，图8(b)中8个时间窗对应的ASR分别为29.42±0.04dB，28.09±0.18dB，25.26±0.03dB，24.99±0.01dB，24.96±0.01dB，25.12±0.02dB，26.44±0.05dB和27.72±0.05dB。相比图8(a)，图8(b)中8个时间窗对应的ASR分别提高了12.35±0.04dB，12.42±0.15dB，10.72±0.05dB，10.44±0.01dB，10.42±0.01dB，10.40±0.01dB，10.86±0.01dB和10.24±0.06dB。上述分析表明本发明中所用的基于M阵元孔径波束合成及互相关加权的方法能够对被动空化成像中的干扰伪影进行有效抑制。
(3)一方面，由于在被动空化成像过程中多阵元超声换能器的各个阵元同步被动接收空化源的辐射信号，因此在空化发生的过程中就可对空化进行实时成像。另一方面，从图8(a)和图8(b)可以看出，位于(-5,35mm)处的空化源在0～5μs的时间窗内开始出现，而位于(5,45mm)处的空化源在10～15μs的时间窗内开始出现；两个空化源出现的时间差为10μs，这与前述的仿真参数设置一致。上述分析表明本发明中被动空化成像具有实时性，且空化起始和终止时刻自适应计算方法能够准确地确定空化起止时间段，因此能够对单脉冲内的超声空化进行准确的实时动态成像。
本发明具有以下优点：
(1)在大多数被动空化成像方法的实施过程中采用的是一种逐像素波束合成的处理方法，该方法在实际应用中通常会耗费大量的计算时间，从而使得被动空化成像的计算速度受到极大限制。而本发明在波束合成处理中则采用了构造稀疏延时信息矩阵和被动空化射频信号扩展矩阵并将两个矩阵相乘的方法，该方法可一次性得到所有像素坐标的波束合成结果，相比逐像素波束合成方法计算时间呈几何数量级的减少，从而大幅度提高了被动空化成像的计算速度。
(2)本发明在基于矩阵相乘的波束合成处理中，通过计算临界采样点数将稀疏延时信息矩阵和被动空化射频信号扩展矩阵相乘的计算复杂度优化至最小；稀疏延时信息矩阵构造一次后即可多次重复利用的特点为被动空化成像提供了极大便利；对延时信息矩阵的稀疏化处理可节省数据存储空间并降低计算量；基于临界采样点数的前后向补零可避免不补零或补零不足导致成像结果错误以及补零过多导致计算量增加的问题；通过将补零信号矩阵转化为列向量及截取元素可减少构造被动空化射频信号扩展矩阵所需的时间；上述方法均有助于进一步提高被动空化成像的计算速度。
(3)传统被动空化成像常采用的是延时叠加波束合成方法，该方法难以有效对抗换能器缺陷、组织异质性和空化源相互作用等复杂因素造成的干扰伪影。而本发明采用间隔抽取阵元的方式将全阵元孔径分割为多个M阵元孔径，然后计算互相关系数向量并利用其对全阵元孔径波束合成信号矩阵进行加权，充分利用了M阵元孔径的波束合成信号在空化源处具有高相关性而在非空化源处具有低相关性的特点，从而可鲁棒地抑制被动空化成像中的干扰伪影。
(4)本发明通过标准化处理M阵元孔径的波束合成信号矩阵及点乘分组所得两个叠加矩阵的方法来计算互相关系数向量，可一次性得到所有像素坐标处的互相关系数并避免了M阵元孔径间互相关处理所需的大量遍历循环操作，从而提高了被动空化成像的计算速度。本发明对互相关系数向量中的元素进行了阈值化及归一化处理，其中，阈值化处理可进一步抑制互相关系数向量中的负值元素向被动空化成像结果中引入的干扰伪影，归一化处理可有效避免空化源坐标处波束合成信号相关性降低导致的空化能量估计错误问题。
(5)本发明中在对空化成像时多阵元超声换能器的各个阵元工作在同步被动接收模式，从而可在空化发生的过程中对空化进行实时成像。本发明中采用的基于相关系数加权矩阵绝对值化处理及列向量转化、列向量元素降序排列、矩阵行索引向量元素去重复处理、多空化源像素数目计算、相关系数加权矩阵行向量希尔伯特变换、包络行向量元素阈值化处理等的空化起始和终止时刻自适应计算方法，能够准确地确定出空化起止时间段，在此基础上对其进行分段得到多个时间窗的被动空化成像结果，从而可实现单脉冲内超声空化的实时动态成像。
(6)本发明提出的单脉冲内超声空化的快速低伪影实时动态成像方法兼顾了超声空化成像在快速计算性能、伪影抑制性能以及单脉冲内实时动态成像功能三方面的需求，其性能和功能不受限于超声空化源激励形式(如不同的超声激励换能器、不同的超声激励脉冲参数等)，可用于观察多种超声空化介导的应用中单脉冲内空化的时空瞬态动力学过程，并进一步研究和分析单脉冲内超声空化诱发的组织生物学效应。
(7)本发明提出的方法具有较高的扩展性，首先，除本发明实施例中所述的线阵换能器外，本发明提出的方法也适用于相控阵换能器等其他一维阵列超声换能器；其次，在本发明提出的方法的基础上，通过将一维阵列超声换能器与一外部的三维扫描装置相连接，可实现基于逐面扫描的空化三维成像；再次，本发明提出的方法经过拓展后也可适用于面阵换能器等二维阵列超声换能器，从而可直接对空化进行三维成像。

Claims (10)

1.一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：包括以下步骤：
1)根据多阵元超声换能器所在位置建立被动空化成像坐标系，在该成像坐标系下规划被动空化成像区域，根据由该成像区域内任意一像素到多阵元超声换能器任意一阵元的距离构成的像素阵元距离集合中的最大距离和最小距离计算临界采样点数；
2)通过间隔抽取阵元的方式对多阵元超声换能器的全阵元孔径进行分割，得到多个M阵元孔径，针对每个M阵元孔径初始化延时信息矩阵，将根据所述被动空化成像区域内任意一像素到M阵元孔径各阵元的距离计算的相对延时转化为延时采样点数并对延时采样点数进行基准化处理，得到基准化延时采样点数，分别计算在所述延时信息矩阵中像素坐标的行索引和基准化延时采样点数的列索引并依据该行、列索引构造延时信息矩阵，对延时信息矩阵进行稀疏化处理，得到每个M阵元孔径的稀疏延时信息矩阵；
3)根据所述临界采样点数对多阵元超声换能器被动接收所得的被动空化射频信号矩阵进行前后向补零处理，得到全阵元孔径补零信号矩阵，对全阵元孔径补零信号矩阵按照所述间隔抽取阵元的方式进行分割，得到每个M阵元孔径的补零信号矩阵，将该补零信号矩阵转化为列向量并由该列向量通过截取元素生成多个子列向量，利用生成的子列向量构造每个M阵元孔径的被动空化射频信号扩展矩阵；
4)将所述每个M阵元孔径的稀疏延时信息矩阵和被动空化射频信号扩展矩阵相乘得到波束合成信号矩阵并对该波束合成信号矩阵进行标准化处理，得到对应M阵元孔径的标准化波束合成信号矩阵，将所有M阵元孔径的标准化波束合成信号矩阵分为两组并分别进行叠加，得到两个叠加矩阵，对两个叠加矩阵进行点乘，得到相关性矩阵并计算互相关系数向量，对互相关系数向量中的元素进行阈值化及归一化处理并利用处理所得结果对由所有M阵元孔径的波束合成信号矩阵叠加所得的全阵元孔径波束合成信号矩阵进行按行加权，得到相关系数加权矩阵；
5)对所述相关系数加权矩阵的各元素进行平方，得到功率矩阵，对功率矩阵中的元素沿着行方向在由空化起始时刻和空化终止时刻确定的空化起止时间段经分段所得的时间窗上进行叠加，得到对应时间窗的能量向量，将该能量向量转化为能量矩阵并进行转置，得到对应时间窗的被动空化成像结果。
2.根据权利要求1所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤1)中，临界采样点数的计算公式表示为：
NCritSamp＝ceil(Dmax/c×Fs)-floor(Dmin/c×Fs)+1
其中，NCritSamp为临界采样点数，ceil(·)表示向上取整，floor(·)表示向下取整，Dmax和Dmin分别为像素阵元距离集合中的最大距离和最小距离，c为超声波在介质中的传播速度，Fs为被动空化射频信号采样频率。
3.根据权利要求1所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤2)中，间隔抽取阵元的方式为：自多阵元超声换能器的第s个阵元开始，每隔NE/M个阵元抽取1个阵元，构成第s个M阵元孔径；其中s＝1,2,...,NE/M，NE为多阵元超声换能器的阵元数目，M为M阵元孔径的阵元数目。
4.根据权利要求1所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤2)中，延时采样点数的计算公式表示为：
Figure FDA0003324784660000021
其中，
Figure FDA0003324784660000022
为根据被动空化成像区域内坐标为
Figure FDA0003324784660000023
的像素到第s个M阵元孔径中坐标为
Figure FDA0003324784660000024
的第is个阵元的距离计算的相对延时所对应的延时采样点数，round[·]表示四舍五入取整，is＝1,2,...,M，s＝1,2,...,NE/M，j＝1,2,...,NPx，k＝1,2,...,NPz，NPx和NPz分别为被动空化成像区域沿x轴和沿z轴的像素数目，c为超声波在介质中的传播速度，Fs为被动空化射频信号采样频率；
基准化延时采样点数利用根据像素阵元距离集合中的最小距离计算的最小相对延时对应的最小延时采样点数来计算，基准化延时采样点数的计算公式表示为：
Figure FDA0003324784660000025
其中，is＝1,2,...,M，s＝1,2,...,NE/M，j＝1,2,...,NPx，k＝1,2,...,NPz，floor(·)表示向下取整；
延时信息矩阵的构造包括以下步骤：当第s个M阵元孔径的延时信息矩阵As的第m行的列索引n为nD时，第m行第n列对应的元素值
Figure FDA0003324784660000026
为1，否则为0；其中，m为坐标为
Figure FDA0003324784660000027
的像素在延时信息矩阵As中的行索引，nD为第s个M阵元孔径的第is个阵元的基准化延时采样点数在延时信息矩阵As中的列索引：
m＝j+(k-1)NPx
Figure FDA0003324784660000031
其中，j＝1,2,...,NPx，k＝1,2,...,NPz，is＝1,2,...,M，s＝1,2,...,NE/M；
对延时信息矩阵的稀疏化处理包括以下步骤：将延时信息矩阵中的所有零元素去除而仅存储所有的非零元素及其对应的索引。
5.根据权利要求1所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤3)中，被动空化射频信号矩阵是通过使多阵元超声换能器各个阵元同步被动接收超声激励所产生的空化源的辐射信号而得到的；按照下述公式对被动空化射频信号矩阵进行前后向补零处理：
Figure FDA0003324784660000032
其中，
Figure FDA0003324784660000033
为全阵元孔径补零信号矩阵
Figure FDA0003324784660000034
中的元素，
Figure FDA0003324784660000035
为NE行(NAcqSamp+2NCritSamp-2)列的矩阵，RFi,u为被动空化射频信号矩阵RF中的元素，RF为NE行NAcqSamp列的矩阵，i＝1,2,...,NE，v＝1,2,...,NAcqSamp+2NCritSamp-2，NAcqSamp为被动空化射频信号采样点数，NCritSamp为临界采样点数。
6.根据权利要求1所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤3)中，截取元素包括以下步骤：对由第s个M阵元孔径的补零信号矩阵转化得到的M×(NAcqSamp+2NCritSamp-2)行1列的列向量bs，自列向量bs第1个元素开始截取M×NCritSamp个元素得到第1个子列向量
Figure FDA0003324784660000036
自第M+1个元素开始截取M×NCritSamp个元素得到第2个子列向量
Figure FDA0003324784660000037
......；直到自第M×(NAcqSamp+NCritSamp-2)+1个元素开始截取M×NCritSamp个元素得到第NAcqSamp+NCritSamp-1个子列向量
Figure FDA0003324784660000038
从而得到第s个M阵元孔径对应的NAcqSamp+NCritSamp-1个M×NCritSamp行1列的子列向量
Figure FDA0003324784660000041
该子列向量中的元素
Figure FDA0003324784660000042
为列向量bs中的元素
Figure FDA0003324784660000043
其中，s＝1,2,...,NE/M，ti＝1,2,...,NAcqSamp+NCritSamp-1，n＝1,2,...,M×NCritSamp，NAcqSamp为被动空化射频信号采样点数，NCritSamp为临界采样点数。
7.根据权利要求1所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤4)中，对波束合成信号矩阵的标准化处理包括以下步骤：对于第s个M阵元孔径的波束合成信号矩阵Qs的每一行，首先将该行元素减去该行元素的均值得到该行的去均值化元素，然后将该行的去均值化元素除以该行的去均值化元素的平方和的平方根，得到第s个M阵元孔径的标准化波束合成信号矩阵
Figure FDA0003324784660000044
即
Figure FDA0003324784660000045
中的元素
Figure FDA0003324784660000046
的计算公式表示为：
Figure FDA0003324784660000047
其中，
Figure FDA0003324784660000048
为第s个M阵元孔径的波束合成信号矩阵Qs中的元素，s＝1,2,...,NE/M，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1，，NAcqSamp为被动空化射频信号采样点数，NCritSamp为临界采样点数；
互相关系数向量中的元素的计算公式表示为：
Figure FDA0003324784660000049
其中，
Figure FDA00033247846600000410
为相关性矩阵R中的元素，m＝1,2,..,NPx×NPz；
对互相关系数向量中的元素的阈值化处理包括以下步骤：当互相关系数向量中的某一元素小于阈值时，将该元素置为阈值，否则保持不变；阈值一般小于等于0.001；
相关系数加权矩阵中的元素的计算公式表示为：
Figure FDA00033247846600000411
Figure FDA00033247846600000412
其中，
Figure FDA00033247846600000413
为全阵元孔径波束合成信号矩阵Q中的元素，
Figure FDA00033247846600000414
为归一化互相关系数向量
Figure FDA0003324784660000051
中的元素，
Figure FDA0003324784660000052
为对互相关系数向量中的元素进行阈值化处理所得阈值化互相关系数向量
Figure FDA0003324784660000053
中的元素，max{·}表示取最大值，m＝1,2,..,NPx×NPz，ti＝1,2,...,NAcqSamp+NCritSamp-1，，NAcqSamp为被动空化射频信号采样点数，NCritSamp为临界采样点数。
8.根据权利要求1所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤5)中，空化起始时刻和空化终止时刻的计算包括以下步骤：
5.1)对所述相关系数加权矩阵进行绝对值化处理，将处理所得结果转化为列向量；
5.2)将步骤5.1)所得列向量中的元素降序排列，然后将所得各元素的索引按照相关系数加权矩阵的行数和列数转化为矩阵行索引和矩阵列索引；
5.3)对由步骤5.2)中所得矩阵行索引构成的矩阵行索引向量中的元素进行去重复处理，即对于重复出现的某一元素只保留第一次出现的该元素，从去重复处理后得到的矩阵行索引向量中根据多空化源像素数目提取多空化源行索引；
5.4)对每个多空化源行索引在所述相关系数加权矩阵中对应的行向量进行希尔伯特变换及绝对值化处理，得到对应的包络行向量；
5.5)对所得包络行向量中的元素进行阈值化处理，得到阈值化包络行向量，寻找阈值化包络行向量的首非零元素索引和末非零元素索引；
5.6)根据所有多空化源行索引对应的首非零元素索引的最小值确定空化起始时刻，根据所有多空化源行索引对应的末非零元素索引的最大值确定空化终止时刻。
9.根据权利要求8所述一种单脉冲内超声空化的快速低伪影实时动态成像方法，其特征在于：所述步骤5.5)中，对包络行向量中的元素的阈值化处理包括以下步骤：当包络行向量中的某一元素小于阈值时，将该元素置为零，否则保持不变；阈值一般小于包络行向量中所有元素的最大值的0.5倍。
10.一种单脉冲内超声空化的快速低伪影实时动态成像系统，其特征在于：该系统包括临界采样点数计算模块、M阵元孔径稀疏延时信息矩阵构造模块、M阵元孔径被动空化射频信号扩展矩阵构造模块、M阵元孔径波束合成及互相关加权模块和单脉冲内超声空化实时动态成像模块；
所述临界采样点数计算模块用于根据多阵元超声换能器所在位置建立被动空化成像坐标系并在该成像坐标系下规划被动空化成像区域，以及根据由该成像区域内任意一像素到多阵元超声换能器任意一阵元的距离构成的像素阵元距离集合中的最大距离和最小距离计算临界采样点数；
所述M阵元孔径稀疏延时信息矩阵构造模块用于通过间隔抽取阵元的方式对多阵元超声换能器的全阵元孔径进行分割、针对分割所得的每个M阵元孔径初始化延时信息矩阵、将根据所述临界采样点数计算模块规划的被动空化成像区域内任意一像素到M阵元孔径各阵元的距离计算的相对延时转化为延时采样点数并对延时采样点数进行基准化处理、分别计算在所述延时信息矩阵中像素坐标的行索引和基准化延时采样点数的列索引，以及依据该行、列索引构造延时信息矩阵并对每个M阵元孔径的延时信息矩阵进行稀疏化处理；
所述M阵元孔径被动空化射频信号扩展矩阵构造模块用于根据所述临界采样点数计算模块得到的临界采样点数对多阵元超声换能器被动接收所得的被动空化射频信号矩阵进行前后向补零处理、对处理所得全阵元孔径补零信号矩阵按照所述M阵元孔径稀疏延时信息矩阵构造模块中的间隔抽取阵元的方式进行分割、将分割所得每个M阵元孔径的补零信号矩阵转化为列向量并由该列向量通过截取元素生成多个子列向量，以及利用生成的子列向量构造每个M阵元孔径的被动空化射频信号扩展矩阵；
所述M阵元孔径波束合成及互相关加权模块用于将所述M阵元孔径稀疏延时信息矩阵构造模块经稀疏化处理得到的每个M阵元孔径的稀疏延时信息矩阵和所述M阵元孔径被动空化射频信号扩展矩阵构造模块得到的每个M阵元孔径的被动空化射频信号扩展矩阵对应相乘、对相乘所得每个M阵元孔径的波束合成信号矩阵进行标准化处理、将处理所得所有M阵元孔径的标准化波束合成信号矩阵分为两组并分别进行叠加、对叠加所得两个叠加矩阵进行点乘并计算互相关系数向量，以及对互相关系数向量中的元素进行阈值化及归一化处理并利用处理所得结果对由所有M阵元孔径的波束合成信号矩阵叠加所得的全阵元孔径波束合成信号矩阵进行按行加权；
所述单脉冲内超声空化实时动态成像模块用于对所述M阵元孔径波束合成及互相关加权模块经对全阵元孔径波束合成信号矩阵进行按行加权得到的相关系数加权矩阵的各元素进行平方、对平方后所得功率矩阵中的元素沿着行方向在由空化起始时刻和空化终止时刻确定的空化起止时间段经分段所得的时间窗上进行叠加，以及将叠加所得能量向量转化为能量矩阵并进行转置，从而得到对应时间窗的被动空化成像结果。
CN202111263927.0A 2021-10-27 2021-10-27 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统 Active CN114098799B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111263927.0A CN114098799B (zh) 2021-10-27 2021-10-27 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111263927.0A CN114098799B (zh) 2021-10-27 2021-10-27 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统

Publications (2)

Publication Number Publication Date
CN114098799A true CN114098799A (zh) 2022-03-01
CN114098799B CN114098799B (zh) 2023-06-27

Family

ID=80377545

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111263927.0A Active CN114098799B (zh) 2021-10-27 2021-10-27 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统

Country Status (1)

Country Link
CN (1) CN114098799B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116115260A (zh) * 2023-01-03 2023-05-16 西安交通大学 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111265245A (zh) * 2020-01-20 2020-06-12 西安交通大学 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
CN112023283A (zh) * 2020-08-03 2020-12-04 西安交通大学 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统
WO2021004076A1 (zh) * 2019-07-05 2021-01-14 山东大学 基于人工智能芯片的适形穿戴式生物信息监测设备及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2021004076A1 (zh) * 2019-07-05 2021-01-14 山东大学 基于人工智能芯片的适形穿戴式生物信息监测设备及系统
CN111265245A (zh) * 2020-01-20 2020-06-12 西安交通大学 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
CN112023283A (zh) * 2020-08-03 2020-12-04 西安交通大学 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
齐雁;谭冠政;范必双;: "基于FPGA的医学超声成像数字波束合成器设计", 计算机测量与控制, no. 04 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116115260A (zh) * 2023-01-03 2023-05-16 西安交通大学 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统
CN116115260B (zh) * 2023-01-03 2024-05-24 西安交通大学 高分辨高对比快速计算的发收时序同步超声被动空化成像方法及系统

Also Published As

Publication number Publication date
CN114098799B (zh) 2023-06-27

Similar Documents

Publication Publication Date Title
CN110023782B (zh) 用于对超声图像杂波滤波的方法和系统
CN111134719B (zh) 一种聚焦超声辐照相变纳米液滴的主被动超声复合成像方法及系统
US6695778B2 (en) Methods and systems for construction of ultrasound images
CN104688271B (zh) 合成聚焦超声成像方法和装置
US20210132223A1 (en) Method and Apparatus for Ultrasound Imaging with Improved Beamforming
CN104272134B (zh) 超声成像系统中的杂波抑制
US11737733B2 (en) Method of, and apparatus for, determination of position in ultrasound imaging
CN108836389B (zh) 平面波相关点相干自适应波束合成成像方法
US20120163691A1 (en) Time-domain estimator for image reconstruction
CN104739448A (zh) 一种超声成像方法及装置
Rasmussen et al. 3D ultrasound imaging performance of a row-column addressed 2D array transducer: a simulation study
Wei et al. High frame rate volumetric imaging of microbubbles using a sparse array and spatial coherence beamforming
CN111265245B (zh) 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统
CN114098799A (zh) 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统
CN106940883B (zh) 基于超声系统点扩散函数仿真和压缩感知的超声成像方法
Noda et al. Ultrasound imaging with a flexible probe based on element array geometry estimation using deep neural network
CN105266847A (zh) 基于压缩感知自适应波束合成的脉冲逆转谐波平面波快速造影成像方法
Varray et al. Nonlinear radio frequency image simulation for harmonic imaging: Creanuis
CN101961251B (zh) 一种医学超声诊断系统中实时计算变迹曲线的方法及装置
CN107204021B (zh) 基于高斯函数探头响应模型和压缩感知的超声成像方法
Luo et al. Fundamental performance assessment of 2-D myocardial elastography in a phased-array configuration
CN109766646B (zh) 一种基于稀疏通道回波数据重建的超声成像方法及装置
Wang et al. An easily-achieved time-domain beamformer for ultrafast ultrasound imaging based on compressive sensing
US20220155439A1 (en) Method and apparatus for adaptive beamforming
CN106997045B (zh) 基于超声系统点扩散函数测量和压缩感知的超声成像方法

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