CN112023283B - 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统 - Google Patents

基于高阶孔径自相关的环形多阵列超声被动成像方法及系统 Download PDF

Info

Publication number
CN112023283B
CN112023283B CN202010768521.7A CN202010768521A CN112023283B CN 112023283 B CN112023283 B CN 112023283B CN 202010768521 A CN202010768521 A CN 202010768521A CN 112023283 B CN112023283 B CN 112023283B
Authority
CN
China
Prior art keywords
array
ultrasonic
signal
imaging
autocorrelation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010768521.7A
Other languages
English (en)
Other versions
CN112023283A (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 CN202010768521.7A priority Critical patent/CN112023283B/zh
Publication of CN112023283A publication Critical patent/CN112023283A/zh
Application granted granted Critical
Publication of CN112023283B publication Critical patent/CN112023283B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N7/00—Ultrasound therapy
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61M—DEVICES FOR INTRODUCING MEDIA INTO, OR ONTO, THE BODY; DEVICES FOR TRANSDUCING BODY MEDIA OR FOR TAKING MEDIA FROM THE BODY; DEVICES FOR PRODUCING OR ENDING SLEEP OR STUPOR
    • A61M37/00—Other apparatus for introducing media into the body; Percutany, i.e. introducing medicines into the body by diffusion through the skin
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61M—DEVICES FOR INTRODUCING MEDIA INTO, OR ONTO, THE BODY; DEVICES FOR TRANSDUCING BODY MEDIA OR FOR TAKING MEDIA FROM THE BODY; DEVICES FOR PRODUCING OR ENDING SLEEP OR STUPOR
    • A61M37/00—Other apparatus for introducing media into the body; Percutany, i.e. introducing medicines into the body by diffusion through the skin
    • A61M37/0092—Other apparatus for introducing media into the body; Percutany, i.e. introducing medicines into the body by diffusion through the skin using ultrasonic, sonic or infrasonic vibrations, e.g. phonophoresis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H17/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the preceding groups
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/66—Analysis of geometric attributes of image moments or centre of gravity
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61M—DEVICES FOR INTRODUCING MEDIA INTO, OR ONTO, THE BODY; DEVICES FOR TRANSDUCING BODY MEDIA OR FOR TAKING MEDIA FROM THE BODY; DEVICES FOR PRODUCING OR ENDING SLEEP OR STUPOR
    • A61M37/00—Other apparatus for introducing media into the body; Percutany, i.e. introducing medicines into the body by diffusion through the skin
    • A61M2037/0007—Other apparatus for introducing media into the body; Percutany, i.e. introducing medicines into the body by diffusion through the skin having means for enhancing the permeation of substances through the epidermis, e.g. using suction or depression, electric or magnetic fields, sound waves or chemical agents
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10132—Ultrasound image

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Biomedical Technology (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Dermatology (AREA)
  • Medical Informatics (AREA)
  • Anesthesiology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Hematology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Geometry (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明公开了一种基于高阶孔径自相关的环形多阵列超声被动成像方法及系统：构建由K个超声阵列换能器组成的环形多阵列换能器组，建立成像坐标系并计算阵元坐标，然后划分成像网格；对其中任意一阵列换能器被动接收得到的超声被动射频信号进行时间移位，将被动时间移位信号转换为标准化信号后进行高阶孔径自相关处理，将阵元合成信号与高阶孔径自相关加权因子相乘，并将所得信号的模平方沿着信号采样点方向叠加，得到单阵列超声被动成像结果；对所有阵列换能器的单阵列超声被动成像结果进行累乘后计算K次算术方根，得到环形多阵列超声被动成像结果。本发明能够有效提高超声被动成像的分辨率和对比度，可用于超声空化效应的准确实时监控。

Description

基于高阶孔径自相关的环形多阵列超声被动成像方法及系统
技术领域
本发明属于超声检测与超声成像技术领域，涉及一种基于高阶孔径自相关的环形多阵列超声被动成像方法及系统，可用于超声空化效应的准确实时监控。
背景技术
超声空化效应的有效利用和合理控制依赖于影像监控技术的发展。相比磁共振成像技术，超声成像技术可实现动态实时显像且对空化声回波信号更加敏感，更适用于监控动态变化的空化效应。根据超声阵列换能器的工作模式，可将超声成像分为超声主动成像和超声被动成像，相比前者，超声被动成像具有更好的空化检测灵敏度，且能在超声辐照过程中对空化效应进行实时监控。
一般而言，超声被动成像技术使用单个一维超声阵列换能器(例如，线阵或相控阵超声换能器)来进行成像。影响超声被动成像质量的主要因素是阵列衍射模式的限制，即成像质量受到超声阵列换能器的长度、声源发射频率大小或换能器接收带宽大小以及声源与超声阵列换能器之间的距离这三个关键参数的制约。在该限制下单个声源的超声被动成像结果往往呈现出沿超声阵列换能器轴向延伸的“X”型分布；该分布中心的亮斑尺寸较大，造成了低的成像分辨率(由于亮斑的轴向长度相比横向长度更大，且对于上述三个关键参数更为敏感，因此阵列衍射模式造成的低分辨率主要体现在轴向分辨率上)；亮斑四周则为不存在真实声源的干扰伪影区域，造成了低的成像对比度。而在实际应用中，这三个参数很难进行大范围的调节，从而导致超声被动成像的性能难以得到大幅度的提升。
此外，在医学超声成像中，成像质量往往会受到组织非均一性、超声阵列换能器缺陷以及多个声散射源的相互作用这三大因素的影响。一般而言，这些因素对于超声主动成像所造成的影响并不是特别突出，尚在可接受范围内。然而，超声被动成像独有的阵列衍射模式限制会加剧这些因素对成像质量的影响，使得干扰伪影更严重也更容易产生，进一步降低了超声被动成像的分辨率和对比度，导致无法对声源位置、声源数目以及声源能量做出准确的判断，不利于对超声空化效应的监控。
有研究通过采用鲁棒自适应波束合成等技术来优化超声被动成像算法，使得成像质量得到一定程度的改善；然而传统基于单个超声阵列换能器的超声被动成像的轴向分辨率远低于横向分辨率这一根本问题并未得到有效解决。中国专利201810573085.0公开了一种基于改进DMAS算法的超声平面波成像方法，该专利中利用了延时乘累加算法(即二阶孔径自相关处理)来提高超声平面波成像质量。一方面，该专利所针对的超声平面波成像属于基于脉冲发射-回波接收的超声主动成像技术，而非超声被动成像技术；另一方面，该专利所采用的二阶孔径自相关处理在应用于超声被动成像时，成像分辨率和对比度仍有待改善。
超声被动成像技术的性能直接影响对超声空化效应监控的准确性，要得到空化源在空间中位置和分布的准确信息，就必须提高超声被动成像的分辨率和对比度；这是超声被动成像领域一直以来研究的热点问题和待攻克的难点问题。
发明内容
本发明的目的在于提供一种基于高阶孔径自相关的环形多阵列超声被动成像方法及系统。
为了实现上述目的，本发明采用了以下技术方案：
一种基于高阶孔径自相关的环形多阵列超声被动成像方法，包括以下步骤：
1)按照环状图形将K个超声阵列换能器排布成一环形多阵列换能器组，并建立环形多阵列成像坐标系，然后计算环形多阵列换能器组中各超声阵列换能器的各阵元坐标，并对环形多阵列成像的成像网格进行划分；
2)所述K个超声阵列换能器分别被动接收来自所述环状图形内部的声源的辐射信号，得到各自对应的超声被动射频信号；
3)计算任意一像素点到环形多阵列换能器组的第k个超声阵列换能器各阵元的超声传播时间，并将该超声传播时间转换为传播时间采样点数目，利用传播时间采样点数目对所述第k个超声阵列换能器被动接收得到的超声被动射频信号进行时间移位处理，得到被动时间移位信号，将被动时间移位信号以及该信号的模平方分别沿着第k个超声阵列换能器的阵元方向叠加，得到第k个超声阵列换能器在所述像素点处的阵元合成信号以及阵元信号总能量，其中，k＝1,2,...,K；
4)对所述被动时间移位信号进行标准化处理，得到标准化信号，对标准化信号进行高阶孔径自相关处理，得到第k个超声阵列换能器在所述像素点处的高阶孔径自相关合成信号，计算该高阶孔径自相关合成信号的模平方与所述阵元信号总能量的比值，得到第k个超声阵列换能器在所述像素点处的高阶孔径自相关加权因子；
5)将第k个超声阵列换能器在所述像素点处的阵元合成信号与相应的高阶孔径自相关加权因子相乘，得到高阶孔径自相关加权合成信号，对高阶孔径自相关加权合成信号的模平方沿着信号采样点方向进行叠加，得到第k个超声阵列换能器在所述像素点处的声能量值；
6)重复步骤3)～5)，直至得到第k个超声阵列换能器在所有像素点处的声能量值，即得到第k个超声阵列换能器的单阵列超声被动成像结果；
7)对环形多阵列换能器组中所有K个超声阵列换能器各自的单阵列超声被动成像结果进行累乘，并计算累乘结果的K次算术方根，得到环形多阵列超声被动成像结果，对环形多阵列超声被动成像结果进行归一化和对数化处理后进行显示。
优选的，所述K个超声阵列换能器选自超声线阵换能器或超声相控阵换能器等一维超声阵列换能器中的一种或多种。
优选的，所述K的取值为偶数，K最小为4；K的取值不宜过大(过大会增加环形多阵列换能器组的设计成本以及超声被动成像的计算量)；一般而言，K可选择为4～16之间的偶数。
优选的，所述步骤1)具体包括以下步骤：
1.1)将K个相同的超声阵列换能器沿一正K边形的各边布置，得到环形多阵列换能器组，将其中一超声阵列换能器设定为基准阵列换能器，并将其余超声阵列换能器设定为非基准阵列换能器，根据基准阵列换能器的位置建立环形多阵列成像坐标系；
1.2)在所述环形多阵列成像坐标系下，计算所述基准阵列换能器的各阵元坐标，根据基准阵列换能器的各阵元坐标以及各非基准阵列换能器的中心坐标，计算各非基准阵列换能器对应的各阵元坐标；
1.3)在所述环形多阵列成像坐标系下，确定环形多阵列成像的成像边界并设置环形多阵列成像的成像网格中的像素间距，然后根据所述成像边界以及成像网格中的像素间距计算成像网格中的像素数目。
优选的，所述步骤1.1)中，环形多阵列成像坐标系分别以基准阵列换能器的阵元方向、与该基准阵列换能器的阵元方向垂直的方向以及该基准阵列换能器的中心作为该坐标系的x轴、z轴和原点O。
优选的，所述步骤1.2)中，各非基准阵列换能器的中心坐标为正K边形对应边(除与基准阵列换能器对应的边之外)的中点的坐标，其中x轴坐标为与该中点相邻的正K边形两个顶点的x轴坐标之和的一半，z轴坐标为与该中点相邻的正K边形两个顶点的z轴坐标之和的一半，正K边形的各顶点坐标的计算公式表示为：
XVk＝RCcosαk+XC
ZVk＝-RCsinαk+ZC
其中，k＝1,2,...,K，(XVk,ZVk)为正K边形的第k个顶点的坐标，
Figure BDA0002615596410000041
为正K边形的外接圆的半径，LB为正K边形的边长，αk为正K边形的外接圆的圆心OC和点V0所形成的线段OCV0绕OC逆时针旋转到OC和正K边形的第k个顶点Vk所形成的线段OCVk的位置时转过的角度，点V0为经过OC且平行于x轴的直线与所述外接圆的两个交点中x轴坐标为正的交点，(XC,ZC)为正K边形的外接圆的圆心OC的坐标。
优选的，所述步骤1.2)中，K-1个非基准阵列换能器的各阵元坐标的计算公式表示为：
Figure BDA0002615596410000042
Figure BDA0002615596410000043
其中，
Figure BDA0002615596410000044
为第k个非基准阵列换能器的各阵元坐标，
Figure BDA0002615596410000045
为基准阵列换能器的各阵元坐标，i＝1,2,...,NE，NE为超声阵列换能器的阵元数目，βk为第k个非基准阵列换能器与环形多阵列成像坐标系的x轴方向之间的夹角，(XACk,ZACk)为第k个非基准阵列换能器的中心坐标。
优选的，所述步骤1.3)中，环形多阵列成像的成像边界视声源的位置而定，成像边界所确定的成像区域可覆盖正K边形内某一个或某几个感兴趣的声源位置，也可覆盖正K边形内所有的声源位置。最大成像区域对应的成像边界根据正K边形的所有顶点和所有边中点的x轴坐标的最小值和最大值，以及正K边形的所有顶点和所有边中点的z轴坐标的最小值和最大值来划定。
优选的，所述步骤1.3)中，成像网格沿x轴的像素数目为成像边界所确定的成像区域沿x轴的长度与像素沿x轴的间距的比值的向上取整结果，成像网格沿z轴的像素数目为成像边界所确定的成像区域沿z轴的长度与像素沿z轴的间距的比值的向上取整结果。
优选的，成像网格中像素沿x轴的间距和沿z轴的间距一般等于或小于超声阵列换能器的阵元间距。
优选的，所述步骤3)中，传播时间采样点数目的计算公式表示为：
Figure BDA0002615596410000051
其中，i＝1,2,...,NE，k＝1,2,...,K，round[·]表示四舍五入取整数，
Figure BDA0002615596410000052
为任意一像素点(坐标为(XP,ZP))到环形多阵列换能器组中第k个超声阵列换能器的第i个阵元的超声传播时间，
Figure BDA0002615596410000053
为所述像素点与环形多阵列换能器组中第k个超声阵列换能器的第i个阵元之间的距离，c为超声传播速度，Fs为超声被动射频信号采样率。
优选的，所述超声被动射频信号采样率一般设置为声源发射频率的10倍以上。
优选的，所述步骤3)中，阵元合成信号和阵元信号总能量的计算公式分别表示为：
Figure BDA0002615596410000054
Figure BDA0002615596410000055
其中，k＝1,2,...,K，EBk(XP,ZP,n)为阵元合成信号，ETEk(XP,ZP,n)为阵元信号总能量，
Figure BDA0002615596410000056
为被动时间移位信号：
Figure BDA0002615596410000057
其中，n＝1,2,...,NS，NS为超声被动射频信号采样点数目，
Figure BDA0002615596410000058
为环形多阵列换能器组中第k个超声阵列换能器各阵元(i＝1,2,...,NE)被动接收的超声射频信号(即上述第k个超声阵列换能器被动接收得到的超声被动射频信号)。
优选的，所述步骤4)中，标准化信号的计算公式表示为：
Figure BDA0002615596410000059
其中，
Figure BDA00026155964100000510
为对被动时间移位信号进行符号化处理后所得的符号信号，
Figure BDA00026155964100000511
为对被动时间移位信号进行绝对值化处理后所得的绝对值信号，P为高阶孔径自相关的阶数(P取大于等于2的正整数)。
优选的，所述步骤4)中，根据多项式定理对高阶孔径自相关合成信号的计算公式进行简化，简化的计算公式表示为：
Figure BDA0002615596410000061
其中，k＝1,2,...,K，ItmNum为将P拆分成若干个不超过P的正整数之和的拆分方案的数目，Coefb由多项式定理中累乘项的系数确定，P为高阶孔径自相关的阶数，
Figure BDA0002615596410000062
为将标准化信号
Figure BDA0002615596410000063
的p次方沿着阵元方向叠加所得的p次阵元标准化合成信号，
Figure BDA0002615596410000064
的计算公式表示为：
Figure BDA0002615596410000065
其中，p＝1,2,...,P，k＝1,2,...,K。
优选的，所述步骤4)中，高阶孔径自相关加权因子的计算公式表示为：
Figure BDA0002615596410000066
其中，k＝1,2,...,K，HEBk(XP,ZP,n)为高阶孔径自相关合成信号，ETEk(XP,ZP,n)为阵元信号总能量。
优选的，所述步骤5)中，高阶孔径自相关加权合成信号的计算公式表示为：
HCWEBk(XP,ZP,n)＝EBk(XP,ZP,n)HCWFk(XP,ZP,n)
其中，k＝1,2,...,K，HCWEBk(XP,ZP,n)为高阶孔径自相关加权合成信号，EBk(XP,ZP,n)为阵元合成信号，HCWFk(XP,ZP,n)为高阶孔径自相关加权因子。
一种基于高阶孔径自相关的环形多阵列超声被动成像系统，该系统包括阵元合成信号与阵元信号总能量计算模块、高阶孔径自相关加权因子计算模块、单阵列超声被动成像模块以及环形多阵列超声被动成像模块；
所述阵元合成信号与阵元信号总能量计算模块用于执行上述步骤3)，主要是利用成像区域中任意一像素点与第k个超声阵列换能器各阵元之间的超声传播时间对该超声阵列换能器被动接收得到的超声被动射频信号进行时间移位处理，并将处理所得被动时间移位信号以及该信号的模平方分别沿着第k个超声阵列换能器的阵元方向进行叠加；
所述高阶孔径自相关加权因子计算模块用于执行上述步骤4)，主要是对所述阵元合成信号与阵元信号总能量计算模块中的被动时间移位信号进行标准化处理以及对该标准化处理所得信号进行高阶孔径自相关处理，并计算高阶孔径自相关处理所得信号的模平方与所述被动时间移位信号的模平方沿阵元方向的叠加结果(即阵元信号总能量)的比值；
所述单阵列超声被动成像模块用于执行上述步骤5)，主要是针对成像区域各像素点，对所述阵元合成信号与阵元信号总能量计算模块中的被动时间移位信号沿阵元方向的叠加结果(即阵元合成信号)与所述高阶孔径自相关加权因子计算模块计算的比值(即高阶孔径自相关加权因子)进行相乘，并将相乘得到的高阶孔径自相关加权合成信号的模平方沿着信号采样点方向进行叠加；
所述环形多阵列超声被动成像模块用于执行上述步骤7)，主要是对所述单阵列超声被动成像模块在成像区域所有像素点处分别通过对高阶孔径自相关加权合成信号的模平方沿着信号采样点方向进行叠加而得到的针对每个超声阵列换能器的单阵列超声被动成像结果，进行累乘并计算累乘结果的K次算术方根。
优选的，所述高阶孔径自相关加权因子计算模块中，将标准化处理所得信号的p(p＝1,2,...,P，P为高阶孔径自相关的阶数)次方沿着阵元方向叠加，得到p次阵元标准化合成信号，并根据多项式定理利用p次标准化合成信号对高阶孔径自相关处理的计算过程进行简化。
优选的，所述系统还包括环形多阵列成像坐标系建立模块(执行步骤1.1))、环形多阵列阵元坐标计算模块(执行步骤1.2))以及环形多阵列成像网格划分模块(执行步骤1.3))。
优选的，所述系统具体包括上述环形多阵列换能器组以及与该环形多阵列换能器组相连的开放式数字化多通道超声成像平台，开放式数字化多通道超声成像平台用于设置环形多阵列换能器组中K个超声阵列换能器的发收参数，使得各超声阵列换能器工作在不发射只接收的被动模式，以及利用基于高阶孔径自相关的环形多阵列超声被动成像方法(上述步骤1)至步骤7)，并由对应模块执行)对各超声阵列换能器被动接收得到的超声被动射频信号进行处理，并对获得的环形多阵列超声被动成像结果进行显示。
优选的，所述系统还包括与所述环形多阵列换能器组相连的三维移动装置，该三维移动装置可带动环形多阵列换能器组沿着与环形多阵列成像平面垂直的方向移动，从而可实现三维全方位的超声被动成像。
本发明的有益效果体现在：
本发明利用多个超声阵列换能器构成的环形多阵列换能器组对声源进行成像，并对多个超声阵列换能器所得的单阵列超声被动成像结果进行复合(累乘及计算K次算术方根)，使得声源成像结果在多个不同的方向均具有较高的分辨率，同时使得单阵列超声被动成像中出现的严重干扰伪影得到抑制，从而有效地提高了超声被动成像的分辨率和对比度。
本发明中使用的高阶孔径自相关处理，使得接收孔径得到扩展且高频成分得到提高，能够大幅度地减小主瓣尺寸，且处理过程有效提取了接收信号之间的空间相干性，能够大幅度降低旁瓣干扰的影响；将所得的高阶孔径自相关加权因子与阵元合成信号相乘，能够有效增强对声源成像的聚焦能力(即减小主瓣尺寸并降低旁瓣水平)，从而提高了超声被动成像的分辨率和对比度。
本发明提出的超声被动成像方法可用于对超声辐照生物组织所产生的超声空化效应的位点进行高分辨的定位，并对多个超声空化位点在空间中的分布情况进行高分辨的成像，从而为超声空化诱导的高强度组织机械毁损和低强度血脑屏障开放等应用中空化效应的时空变化过程提供了有力的监控工具。
进一步地，本发明中，高阶孔径自相关处理的输入信号是对被动时间移位信号进行符号化和绝对值化处理后计算得到的标准化信号，而非被动时间移位信号，这使得高阶孔径自相关处理的输出信号(即高阶孔径自相关合成信号)与被动时间移位信号的物理量纲相一致，从而保证计算得到的高阶孔径自相关加权因子是无量纲的(一般而言，加权因子均应是无量纲数)。
进一步地，本发明针对高阶孔径自相关处理过程中大量的乘法和加法运算导致的计算量大的问题，首先计算p次阵元标准化合成信号，然后根据多项式定理对高阶孔径自相关合成信号的计算公式进行简化，显著降低了高阶孔径自相关处理中乘法和加法的运算次数，从而有效降低了计算量。
附图说明
图1为本发明实施例中环形多阵列换能器组(a)和环形多阵列成像坐标系(b)的示意图，其中：#1～#8为超声阵列换能器的编号。
图2为本发明实施例中环形多阵列阵元坐标的计算流程图。
图3为本发明实施例中正K边形外接圆(a)、正K边形顶点(b)以及正K边形各边中点(c)的示意图。
图4为本发明实施例中环形多阵列成像网格划分的示意图。
图5为本发明实施例中阵元合成信号与阵元信号总能量的计算流程图(a)、高阶孔径自相关加权因子的计算流程图(b)以及单阵列超声被动成像的流程图(c)。
图6为本发明实施例中环形多阵列超声被动成像的流程图。
图7为本发明实施例中三种情况下所得的单个声源的超声被动成像结果，其中：(a)使用单个阵列换能器且不使用高阶孔径自相关加权；(b)使用环形多阵列换能器组且不使用高阶孔径自相关加权；(c)使用环形多阵列换能器组且使用高阶(阶数等于4)孔径自相关加权。
图8为本发明实施例中三种情况下所得的多个声源的超声被动成像结果，其中：(a)使用单个阵列换能器且不使用高阶孔径自相关加权；(b)使用环形多阵列换能器组且不使用高阶孔径自相关加权；(c)使用环形多阵列换能器组且使用高阶(阶数等于4)孔径自相关加权。
图9为本发明实施例中三种情况下使用环形多阵列换能器组所得的多个声源的超声被动成像结果，其中：(a)使用二阶孔径自相关而不使用二阶孔径自相关加权；(b)使用二阶孔径自相关加权；(c)使用高阶(阶数等于4)孔径自相关加权。
具体实施方式
下面结合附图和实施例对本发明作进一步的详细说明。所述实施例仅用于解释本发明，而非对本发明保护范围的限制。
本发明提出了一种基于高阶孔径自相关的环形多阵列超声被动成像方法，能够有效提高超声被动成像的分辨率并提高成像对比度，从而为超声诊疗等应用中声源的精准定位提供了一种行之有效的解决手段。
参见图1，首先构建由多个超声阵列换能器(例如，一维线阵换能器)组成的环形多阵列换能器组，并建立环形多阵列成像坐标系，具体流程如下述步骤(1.1)～(1.3)所示。
(1.1)参见图1(a)，将K(例如，K＝8)个完全相同的超声阵列换能器(以下将超声阵列换能器简称为阵列换能器)沿着一正K边形的各边长放置，构成环形多阵列换能器组；K个超声阵列换能器分别被动接收来自正K边形内部的声源的辐射信号(声源由位于正K边形所在平面外的超声激励换能器发射穿过正K边形内部的波束，从而在正K边形内部产生)，得到各自对应的超声被动射频信号；
环形多阵列换能器组中各个阵列换能器的中心与正K边形的各边中点重合，并且各阵列换能器的阵元朝向正K边形内部；由于阵列换能器除了阵元主体外还包括外壳体，因此正K边形的边长LB应大于阵列换能器的长度(例如，LB＝55mm)；
(1.2)将环形多阵列换能器组中的任意一阵列换能器设定为基准阵列换能器(例如，选择与水平方向平行的阵列换能器的其中一个)，将该基准阵列换能器标记为1号；将除基准阵列换能器以外的其余K-1个阵列换能器设定为非基准阵列换能器，并按顺时针顺序分别标记为2号至K号；
(1.3)参见图1(b)，根据步骤(1.2)设定的1号基准阵列换能器的位置建立环形多阵列成像坐标系，其中1号基准阵列换能器的阵元方向为坐标系的x轴，与1号基准阵列换能器的阵元方向垂直的方向为坐标系的z轴，1号基准阵列换能器的中心为坐标系的原点O，即1号基准阵列换能器的中心坐标(XAC1,ZAC1)为(0,0)。
参见图2，在环形多阵列成像坐标系下，计算环形多阵列换能器组中各阵列换能器的各阵元坐标，即计算环形多阵列阵元坐标，具体流程如下述步骤(2.1)～(2.5)所示；
(2.1)在步骤(1.3)建立的环形多阵列成像坐标系下，计算步骤(1.2)设定的1号基准阵列换能器的各阵元坐标
Figure BDA0002615596410000101
Figure BDA0002615596410000102
Figure BDA0002615596410000103
其中，i＝1,2,...,NE，NE为阵列换能器的阵元数目(例如，128)，EP为阵列换能器的阵元间距(例如，0.3mm)；
(2.2)计算正K边形的外接圆的半径RC和外接圆的圆心OC的坐标(XC,ZC)(参见图3(a))：
Figure BDA0002615596410000104
XC＝0
Figure BDA0002615596410000105
(2.3)根据步骤(2.2)所得外接圆半径RC和外接圆的圆心坐标(XC,ZC)计算正K边形的各顶点坐标(XVk,ZVk)(参见图3(b)中V1至V8)：
XVk＝RCcosαk+XC
ZVk＝-RCsinαk+ZC
其中，k＝1,2,...,K，αk为线段OCV0绕OC逆时针旋转到线段OCVk的位置时转过的角度(
Figure BDA0002615596410000111
)，V0为经过OC且平行于x轴的直线与外接圆的两个交点中x轴坐标为正的交点(参见图3(a))，Vk为正K边形的第k个顶点；
(2.4)根据步骤(2.3)所得正K边形的各顶点坐标(XVk,ZVk)计算正K边形各边(除与1号基准阵列换能器对应的边之外)的中点的坐标(XACk,ZACk)(k＝2,3,...,K，参见图3(c)中M2至M8)，即得到K-1个非基准阵列换能器的中心坐标：
Figure BDA0002615596410000112
Figure BDA0002615596410000113
(2.5)根据步骤(2.1)所得1号基准阵列换能器的各阵元坐标
Figure BDA0002615596410000114
以及步骤(2.4)所得各非基准阵列换能器(2号至8号非基准阵列换能器)的中心坐标(XACk,ZACk)，计算环形多阵列成像坐标系下各非基准阵列换能器的各阵元坐标
Figure BDA0002615596410000115
Figure BDA0002615596410000116
Figure BDA0002615596410000117
其中，i＝1,2,...,NE，k＝2,3,...,K，βk为k号非基准阵列换能器与环形多阵列成像坐标系中x轴正方向之间的夹角(
Figure BDA0002615596410000118
Figure BDA0002615596410000119
)。
参见图4，在环形多阵列成像坐标系下，通过确定环形多阵列成像边界和设置成像网格中的像素间距，对环形多阵列成像网格进行划分，具体流程如下述步骤(3.1)～(3.4)所示；
(3.1)根据环形多阵列换能器组中阵元方向与z轴平行的两个非基准阵列换能器(7号和3号非基准阵列换能器)的中心坐标确定环形多阵列成像的x轴边界；其中，x轴边界的左边界XAC7和右边界XAC3分别为步骤(2.4)所得的7号非基准阵列换能器和3号非基准阵列换能器的中心坐标的x轴坐标；
(3.2)根据环形多阵列换能器组中的1号基准阵列换能器以及阵元方向与x轴平行的另一个非基准阵列换能器(5号非基准阵列换能器)的中心坐标确定环形多阵列成像的z轴边界；其中，z轴边界的上边界ZAC1和下边界ZAC5分别为步骤(1.3)和(2.4)所得的1号基准阵列换能器和5号非基准阵列换能器的中心坐标的z轴坐标；
(3.3)分别设置环形多阵列成像网格中的像素沿x轴的间距SPx(例如，0.2mm)和沿z轴的间距SPz(例如，0.2mm)；
(3.4)根据步骤(3.1)和(3.2)所得环形多阵列成像的x轴边界和z轴边界以及步骤(3.3)设置的像素间距，计算环形多阵列成像网格中的像素数目：
Figure BDA0002615596410000121
Figure BDA0002615596410000122
其中，ceil(·)表示向上取整数，NPx和NPz分别为沿x轴和沿z轴的像素数目，即环形多阵列成像网格中包含了NPx×NPz个像素。
参见图5(a)，对环形多阵列换能器组中任意一阵列换能器被动接收得到的超声被动射频信号进行时间移位处理，在此基础上即可计算得到对应单个阵列换能器的阵元合成信号和阵元信号总能量，具体流程如下述步骤(4.1)～(4.5)所示；
(4.1)针对所述环形多阵列成像网格中任意一像素点坐标(XP,ZP)，分别计算该像素点坐标(XP,ZP)与环形多阵列换能器组中k号阵列换能器(此处及以下的流程描述中不再区分基准阵列换能器和非基准阵列换能器)的各阵元坐标之间的距离
Figure BDA0002615596410000123
Figure BDA0002615596410000124
其中，
Figure BDA0002615596410000131
为步骤(2.1)和(2.5)所得的环形多阵列换能器组中k号阵列换能器的第i个阵元的坐标；
(4.2)根据步骤(4.1)所得距离
Figure BDA0002615596410000132
计算超声传播时间，并根据超声被动射频信号采样率将该超声传播时间转换为传播时间采样点数目
Figure BDA0002615596410000133
Figure BDA0002615596410000134
其中，i＝1,2,...,NE，k＝1,2,...,K，round[·]表示四舍五入取整数，c为超声传播速度(例如，1480m/s)，Fs为超声被动射频信号采样率(例如，20MHz)；
(4.3)根据步骤(4.2)所得传播时间采样点数目
Figure BDA0002615596410000135
对环形多阵列换能器组中k号阵列换能器被动接收得到的超声被动射频信号
Figure BDA0002615596410000136
进行时间移位处理，得到k号阵列换能器的被动时间移位信号
Figure BDA0002615596410000137
Figure BDA0002615596410000138
其中，i＝1,2,...,NE，k＝1,2,...,K，n＝1,2,...,NS，NS为超声被动射频信号采样点数目(例如，2000)；
(4.4)将步骤(4.3)所得被动时间移位信号
Figure BDA0002615596410000139
沿着阵元方向叠加得到k号阵列换能器在像素点坐标(XP,ZP)处的阵元合成信号EBk(XP,ZP,n)：
Figure BDA00026155964100001310
其中，k＝1,2,...,K；
(4.5)将步骤(4.3)所得被动时间移位信号
Figure BDA00026155964100001311
的模平方沿着阵元方向叠加得到k号阵列换能器在像素点坐标(XP,ZP)处的阵元信号总能量ETEk(XP,ZP,n)：
Figure BDA00026155964100001312
其中，k＝1,2,...,K。
参见图5(b)，将被动时间移位信号转换为标准化信号后进行高阶孔径自相关处理，并根据所得高阶孔径自相关合成信号计算高阶孔径自相关加权因子，具体流程如下述步骤(5.1)～(5.4)所示；
(5.1)对步骤(4.3)所得k号阵列换能器的被动时间移位信号
Figure BDA0002615596410000141
分别进行符号化处理和绝对值化处理，得到符号信号
Figure BDA0002615596410000142
和绝对值信号
Figure BDA0002615596410000143
Figure BDA0002615596410000144
Figure BDA0002615596410000145
其中，i＝1,2,...,NE，k＝1,2,...,K，sign[·]表示对信号做符号化处理，即将信号中大于0的值置为1，小于0的值置为-1，等于0的值保持不变，|·|表示对信号做绝对值化处理；
(5.2)将步骤(5.1)所得的符号信号
Figure BDA0002615596410000146
与绝对值信号
Figure BDA0002615596410000147
的
Figure BDA0002615596410000148
次方相乘，得到k号阵列换能器的标准化信号
Figure BDA0002615596410000149
Figure BDA00026155964100001410
其中，i＝1,2,...,NE，k＝1,2,...,K，P为高阶孔径自相关的阶数(例如，P为大于2的正整数)；
(5.3)对步骤(5.2)所得标准化信号
Figure BDA00026155964100001411
进行阶数为P的高阶孔径自相关处理，得到k号阵列换能器在像素点坐标(XP,ZP)处的高阶孔径自相关合成信号HEBk(XP,ZP,n)：
Figure BDA00026155964100001412
其中，k＝1,2,...,K；
上述高阶孔径自相关合成信号HEBk(XP,ZP,n)的计算公式需要大量的乘法和加法运算，计算量过大；为此，根据下述步骤(5.3.1)和(5.3.2)对其进行简化以减小计算量：
(5.3.1)将步骤(5.2)所得标准化信号
Figure BDA0002615596410000151
的p次方沿着阵元方向叠加，得到p次阵元标准化合成信号
Figure BDA0002615596410000152
Figure BDA0002615596410000153
其中，p＝1,2,...,P，k＝1,2,...,K；
(5.3.2)根据多项式定理
Figure BDA0002615596410000154
利用步骤(5.3.1)所得p次阵元标准化合成信号
Figure BDA0002615596410000155
对步骤(5.3)所述的高阶孔径自相关合成信号HEBk(XP,ZP,n)的计算公式进行简化，得到简化公式：
Figure BDA0002615596410000156
其中，k＝1,2,...,K，上式中共有ItmNum项，ItmNum为将P拆分成若干个不超过P的正整数之和的拆分方案的数目，每一项的系数Coefb由多项式定理中累乘项的系数来确定；
以P＝4为例，简化后的高阶孔径自相关合成信号HEBk(XP,ZP,n)的计算公式表示为：
Figure BDA0002615596410000157
在步骤(5.3)中，计算高阶孔径自相关合成信号HEBk(XP,ZP,n)所需的乘法和加法运算的总次数为
Figure BDA0002615596410000158
(
Figure BDA0002615596410000159
表示从NE个不同阵元中取出P个的所有组合的个数)；步骤(5.3.2)所得简化公式的计算量主要是由p次阵元标准化合成信号
Figure BDA00026155964100001510
的计算贡献，P个p次阵元标准化合成信号
Figure BDA0002615596410000161
所需的乘法和加法运算的总次数为
Figure BDA0002615596410000162
这远远小于简化前的计算量。
(5.4)计算步骤(5.3)或步骤(5.3.2)所得高阶孔径自相关合成信号HEBk(XP,ZP,n)的模平方与步骤(4.5)所得阵元信号总能量ETEk(XP,ZP,n)的比值，得到k号阵列换能器在像素点坐标(XP,ZP)处的高阶孔径自相关加权因子HCWFk(XP,ZP,n)：
Figure BDA0002615596410000163
其中，k＝1,2,...,K。
参见图5(c)，通过将环形多阵列成像网格中各像素点坐标处的阵元合成信号与高阶孔径自相关加权因子相乘，并将所得信号的模平方沿着信号采样点方向进行叠加，得到环形多阵列成像网格中所有像素点坐标处的声能量值，即得到上述单个阵列换能器对应的单阵列超声被动成像结果，具体流程如下述步骤(6.1)～(6.3)所示；
(6.1)将步骤(4.4)所得阵元合成信号EBk(XP,ZP,n)与步骤(5.4)所得高阶孔径自相关加权因子HCWFk(XP,ZP,n)相乘，得到k号阵列换能器在像素点坐标(XP,ZP)处的高阶孔径自相关加权合成信号HCWEBk(XP,ZP,n)：
HCWEBk(XP,ZP,n)＝EBk(XP,ZP,n)HCWFk(XP,ZP,n)
其中，k＝1,2,...,K；
(6.2)对步骤(6.1)所得高阶孔径自相关加权合成信号HCWEBk(XP,ZP,n)的模平方沿着信号采样点方向进行叠加，得到k号阵列换能器在像素点坐标(XP,ZP)处的声能量值Energyk(XP,ZP)：
Figure BDA0002615596410000164
其中，k＝1,2,...,K；
(6.3)重复步骤(4.1)～(4.5)、步骤(5.1)～(5.4)以及步骤(6.1)～(6.2)，直至得到环形多阵列成像网格中所有NPx×NPz个像素点对应的声能量值，从而得到k号阵列换能器的单阵列超声被动成像结果，记为UPIk，该结果为一个NPz行NPx列的二维矩阵。
参见图6，对所有阵列换能器的单阵列超声被动成像结果进行累乘并计算累乘结果的K次算术方根，然后进行显示，即完成环形多阵列超声被动成像，具体流程如下述步骤(7.1)～(7.2)所示；
(7.1)根据步骤(6.3)可分别得到1号至K号阵列换能器的单阵列超声被动成像结果UPIk(k＝1,2,...,K)，对该结果进行累乘，然后计算累乘结果的K次算术方根，得到环形多阵列超声被动成像结果MUPI(该结果为一个NPz行NPx列的二维矩阵)：
Figure BDA0002615596410000171
(7.2)超声成像一般以对数形式进行显示，因此以步骤(7.1)所得环形多阵列超声被动成像结果MUPI中的最大值为基准对其进行归一化处理，然后进行对数化处理，得到对数化的环形多阵列超声被动成像结果MUPILog：
Figure BDA0002615596410000172
其中，max(·)表示取最大值。
图7和图8所示的单个声源和多个声源的超声被动成像结果(图像显示的动态范围设置为-80～0dB)中，单个声源的坐标为(0mm,66.4mm)；多个声源中包括9个声源，坐标(从左至右、从上至下)依次为(-10mm,56.4mm)、(0mm,56.4mm)、(10mm,56.4mm)、(-10mm,66.4mm)、(0mm,66.4mm)、(10mm,66.4mm)、(-10mm,76.4mm)、(0mm,76.4mm)和(10mm,76.4mm)。图7和图8中x轴对应于步骤(1.2)所述1号基准阵列换能器的横向(即阵元方向)，z轴对应于1号基准阵列换能器的轴向(即与阵元方向垂直的方向)。
图7和图8中(a)～(c)分别为以下三种情况下所得的超声被动成像结果：
(a)使用单个阵列换能器且不使用高阶孔径自相关加权(即目前最常用的超声被动成像方法)：对应结果为步骤(6.1)中的高阶孔径自相关加权因子置为1时，步骤(6.3)所得的1号阵列换能器的单阵列超声被动成像结果UPI1经过归一化和对数化处理后所得的结果；
(b)使用环形多阵列换能器组且不使用高阶孔径自相关加权(即在目前最常用的超声被动成像方法的基础上使用环形多阵列)：对应结果为步骤(6.1)中的高阶孔径自相关加权因子置为1时，步骤(7.2)所得的对数化的环形多阵列超声被动成像结果MUPILog；
(c)使用环形多阵列换能器组且使用高阶孔径自相关加权(即本发明提出的超声被动成像方法)：对应结果为步骤(5.2)中的高阶孔径自相关的阶数置为4时，步骤(7.2)所得的对数化的环形多阵列超声被动成像结果MUPILog。
从单个声源的超声被动成像结果可以看出，图7(a)中像素呈现沿z轴方向延伸的“X”型分布，中心亮斑尺寸较大且干扰伪影很强，表明该图像具有低的分辨率和对比度。相比图7(a)，图7(b)的中心亮斑明显减小且干扰伪影得到明显抑制，表明成像分辨率和对比度得到提升。相比图7(b)，图7(c)的中心亮斑尺寸进一步减小且旁瓣水平明显降低，表明成像分辨率和对比度得到进一步提升。
从多个声源的超声被动成像结果可以看出，图8(a)中无法区分出9个声源的空间位置，沿z轴分布的每列的三个声源粘连在一起，这主要是由于单阵列超声被动成像极低的轴向分辨率引起的；且图8(a)的干扰伪影相比图7(a)更为严重，这主要是由于多个声源之间的相互干扰造成的，表明该图像具有极低的分辨率和对比度。相比图8(a)，图8(b)中各声源被清晰地区分开来且干扰伪影得到大幅度降低，表明成像分辨率和对比度得到大幅度提升。相比图8(b)，图8(c)中各声源对应的亮斑尺寸进一步减小且旁瓣水平得到大幅度降低，表明成像分辨率和对比度得到进一步提升。
为了定量地评价本发明提出的基于高阶孔径自相关的环形多阵列超声被动成像方法的性能，采用A-6dB和ICR两个指标分别评价成像分辨率和对比度；其中，A-6dB定义为对数化的成像结果中小于-6dB的像素所对应的面积，该值越小，成像分辨率越高；ICR定义为20lg(MPV1/MPV2)，该值越大，成像对比度越高，MPV1和MPV2分别为未经对数化处理的成像结果中大于最大像素值一半的像素值的平均值和小于最大像素值一半的像素值的平均值。
单个声源的定量评价结果表明，图7(a)、图7(b)和图7(c)所对应的A-6dB分别为17.9mm2、2.4mm2和0.8mm2，对应的ICR分别为36.0dB、49.7dB和93.4dB。相比图7(a)，图7(b)的A-6dB减小了15.5mm2且ICR提高了13.7dB；相比图7(b)，图7(c)的A-6dB减小了1.6mm2且ICR提高了43.7dB。多个声源的定量评价结果表明，图8(a)、图8(b)和图8(c)所对应的A-6dB分别为128.7mm2、24.2mm2和4.4mm2，对应的ICR分别为21.5dB、27.9dB和72.0dB。相比图8(a)，图8(b)的A-6dB减小了104.5mm2且ICR提高了6.4dB；相比图8(b)，图8(c)的A-6dB减小了19.8mm2且ICR提高了44.1dB。以上结果定量地证明了相比传统的单阵列超声被动成像方法，环形多阵列超声被动成像方法能够有效地提高成像分辨率和对比度；在环形多阵列超声被动成像的实施过程中引入高阶孔径自相关加权后，能够进一步提高成像分辨率和对比度。
图9中(a)～(c)分别为以下三种情况下使用环形多阵列换能器组所得的多个声源(9个声源，各个声源的坐标与图8一致)的超声被动成像结果(图像显示的动态范围设置为-80～0dB，图中x轴对应于步骤(1.2)所述1号基准阵列换能器的横向(即阵元方向)，z轴对应于1号基准阵列换能器的轴向(即与阵元方向垂直的方向))：
(a)使用二阶孔径自相关而不使用二阶孔径自相关加权：对应结果为将步骤(5.2)中的高阶孔径自相关的阶数置为2，并利用步骤(5.3)或步骤(5.3.2)所得二阶孔径自相关合成信号替换掉步骤(6.2)中的高阶孔径自相关加权合成信号时，步骤(7.2)所得的对数化的环形多阵列超声被动成像结果MUPILog；
(b)使用二阶孔径自相关加权：对应结果为将步骤(5.2)中的高阶孔径自相关的阶数置为2时，步骤(7.2)所得的对数化的环形多阵列超声被动成像结果MUPILog；
(c)使用高阶孔径自相关加权：对应结果为将步骤(5.2)中的高阶孔径自相关的阶数置为4时，步骤(7.2)所得的对数化的环形多阵列超声被动成像结果MUPILog。
从超声被动成像结果可以看出，相比图9(a)，图9(b)中各声源对应的亮斑尺寸减小且旁瓣水平降低，表明成像分辨率和对比度得到提升。相比图9(b)，图9(c)中各声源对应的亮斑尺寸减小且旁瓣水平降低，表明成像分辨率和对比度得到提升。根据前述A-6dB和ICR两个指标分别对成像分辨率和对比度进行定量评价，结果表明图9(a)、图9(b)和图9(c)所对应的A-6dB分别为10.7mm2、7.2mm2和4.4mm2，对应的ICR分别为46.5dB、66.9dB和72.0dB。相比图9(a)，图9(b)的A-6dB减小了3.5mm2且ICR提高了20.4dB；相比图9(b)，图9(c)的A-6dB减小了2.8mm2且ICR提高了5.1dB。以上结果定量地证明，相比使用二阶孔径自相关而不使用二阶孔径自相关加权的情况，使用二阶孔径自相关加权能够有效地提高成像分辨率和对比度；相比使用二阶孔径自相关加权的情况，使用高阶孔径自相关加权能够有效地提高成像分辨率和对比度。
本发明具有以下优点：
(1)传统基于单个超声阵列换能器的超声被动成像方法由于受到阵列衍射模式的限制导致轴向分辨率远低于横向分辨率且容易出现严重的干扰伪影，而本发明则采用由多个超声阵列换能器组成的环形多阵列换能器组来进行超声被动成像。通过对在不同角度放置的多个超声阵列换能器对应的单阵列超声被动成像结果进行复合(累乘及计算K次算术方根)，充分利用了各个角度的单阵列超声被动成像的良好的横向分辨率，有效增强了对声源的分辨能力以及对干扰伪影的抑制程度，从而提高了超声被动成像的分辨率和对比度。
(2)本发明在传统超声被动成像算法的基础上引入了高阶孔径自相关处理过程，该处理过程一方面可提高接收孔径和高频成分，另一方面使得接收信号之间的空间相干性得到非常充分的利用。该处理过程获得的高阶孔径自相关加权因子能够更好地评价成像聚焦性能，从而更好地对主瓣信号和旁瓣信号进行区分。将该加权因子与阵元合成信号相乘后，一方面可以有效克服阵列衍射模式对成像分辨率的限制，另一方面可以有效地降低旁瓣对成像结果的影响，从而进一步提高超声被动成像的分辨率和对比度。
(3)本发明首先对被动时间移位信号进行符号化和绝对值化处理，然后对标准化信号进行高阶孔径自相关处理，使得高阶孔径自相关合成信号具有与被动时间移位信号相一致的物理量纲，从而保证高阶孔径自相关加权因子是无量纲的。在本发明涉及的高阶孔径自相关处理过程中，需要进行大量的乘法和加法运算，从而导致了极大的计算量。本发明首先根据标准化信号计算p次阵元标准化合成信号，然后根据多项式定理对高阶孔径自相关处理过程进行简化，使得计算所需的乘法运算和加法运算的次数大大降低，有效降低了高阶孔径自相关处理过程的计算复杂度，从而减小了环形多阵列超声被动成像的计算时间。
(4)本发明可对不同形式的超声场(例如，聚焦波声场、非聚焦声波场和驻波声场等)作用下、多种超声场联合作用下以及超声场与其他物理场联合作用下多个空化源以及单个空化源内多个空化微泡的空间分布进行准确监控，进一步可用于研究超声诱导空化背后的复杂物理机制，从而为组织消融、组织机械毁损、血管栓块溶解和血脑屏障开放及药物控释等应用的超声激励参数的控制和优化奠定基础。
(5)本发明提出的环形多阵列超声被动成像的主要目的是为了克服传统单阵列超声被动成像在分辨率和对比度上的限制，但在此基础上通过一外部的三维移动装置带动环形多阵列换能器组沿着与x-z成像平面垂直的y轴方向移动，即可实现高性能的三维超声被动成像，从而为超声空化提供三维全方位的监控手段。

Claims (8)

1.一种基于高阶孔径自相关的环形多阵列超声被动成像方法，其特征在于：包括以下步骤：
1)按照环状图形将K个超声阵列换能器排布成环形多阵列换能器组，并建立环形多阵列成像坐标系，然后计算环形多阵列换能器组中各超声阵列换能器的各阵元坐标，并对环形多阵列成像的成像网格进行划分，所述成像网格中任意一像素点坐标为(XP,ZP)；
2)所述K个超声阵列换能器分别被动接收来自所述环状图形内部的声源的辐射信号，得到各自对应的超声被动射频信号；
3)计算任意一像素点到环形多阵列换能器组的第k个超声阵列换能器各阵元的超声传播时间，并将该超声传播时间转换为传播时间采样点数目，利用传播时间采样点数目对所述第k个超声阵列换能器被动接收得到的超声被动射频信号进行时间移位处理，得到被动时间移位信号，将被动时间移位信号以及该信号的模平方分别沿着第k个超声阵列换能器的阵元方向叠加，得到第k个超声阵列换能器在所述像素点处的阵元合成信号以及阵元信号总能量，其中，k＝1,2,...,K；
4)对所述被动时间移位信号进行标准化处理，得到标准化信号，对标准化信号进行高阶孔径自相关处理，得到第k个超声阵列换能器在所述像素点处的高阶孔径自相关合成信号，计算该高阶孔径自相关合成信号的模平方与所述阵元信号总能量的比值，得到第k个超声阵列换能器在所述像素点处的高阶孔径自相关加权因子；
所述步骤4)中，标准化信号的计算公式表示为：
Figure FDA0003166941730000011
其中，i＝1,2,...,NE，NE为超声阵列换能器的阵元数目，k＝1,2,...,K，
Figure FDA0003166941730000012
为对被动时间移位信号进行符号化处理后所得的符号信号，
Figure FDA0003166941730000013
为对被动时间移位信号进行绝对值化处理后所得的绝对值信号，n＝1,2,...,NS，NS为超声被动射频信号采样点数目，P为高阶孔径自相关的阶数，P取大于等于2的正整数；
所述步骤4)中，高阶孔径自相关合成信号的计算公式表示为：
Figure FDA0003166941730000014
上式中共有ItmNum项，ItmNum为将P拆分成若干个不超过P的正整数之和的拆分方案的数目，每一项的系数Coefb由多项式定理中累乘项的系数确定，P为高阶孔径自相关的阶数，P取大于等于2的正整数，k＝1,2,...,K，
Figure FDA0003166941730000021
为将标准化信号
Figure FDA0003166941730000022
的p次方沿着阵元方向叠加所得的p次阵元标准化合成信号，
Figure FDA0003166941730000023
的计算公式表示为：
Figure FDA0003166941730000024
其中，p＝1,2,...,P，k＝1,2,...,K；
5)将第k个超声阵列换能器在所述像素点处的阵元合成信号与相应的高阶孔径自相关加权因子相乘，得到高阶孔径自相关加权合成信号，对高阶孔径自相关加权合成信号的模平方沿着信号采样点方向进行叠加，得到第k个超声阵列换能器在所述像素点处的声能量值；
6)重复步骤3)～5)，直至得到第k个超声阵列换能器在所有像素点处的声能量值，即得到第k个超声阵列换能器的单阵列超声被动成像结果；
7)对环形多阵列换能器组中所有K个超声阵列换能器各自的单阵列超声被动成像结果进行累乘，并计算累乘结果的K次算术方根，得到环形多阵列超声被动成像结果，对环形多阵列超声被动成像结果进行归一化和对数化处理后进行显示。
2.根据权利要求1所述一种基于高阶孔径自相关的环形多阵列超声被动成像方法，其特征在于：所述K个超声阵列换能器选自一维超声阵列换能器中的一种或多种；所述K为4～16之间的偶数。
3.根据权利要求1所述一种基于高阶孔径自相关的环形多阵列超声被动成像方法，其特征在于：所述步骤1)具体包括以下步骤：
1.1)将K个相同的超声阵列换能器沿正K边形的各边布置，得到环形多阵列换能器组，将其中一超声阵列换能器设定为基准阵列换能器，并将其余超声阵列换能器设定为非基准阵列换能器，根据基准阵列换能器的位置建立环形多阵列成像坐标系；
1.2)在所述环形多阵列成像坐标系下，计算所述基准阵列换能器的各阵元坐标，根据基准阵列换能器的各阵元坐标以及各非基准阵列换能器的中心坐标，计算各非基准阵列换能器对应的各阵元坐标；
1.3)在所述环形多阵列成像坐标系下，确定环形多阵列成像的成像边界并设置环形多阵列成像的成像网格中的像素间距，然后根据所述成像边界以及成像网格中的像素间距计算成像网格中的像素数目。
4.根据权利要求1所述一种基于高阶孔径自相关的环形多阵列超声被动成像方法，其特征在于：所述步骤3)中，传播时间采样点数目的计算公式表示为：
Figure FDA0003166941730000031
其中，i＝1,2,...,NE，NE为超声阵列换能器的阵元数目，k＝1,2,...,K，round[·]表示四舍五入取整数，
Figure FDA0003166941730000032
为任意一坐标为(XP,ZP)的像素点到环形多阵列换能器组中第k个超声阵列换能器的第i个阵元的超声传播时间，
Figure FDA0003166941730000033
为所述像素点与环形多阵列换能器组中第k个超声阵列换能器的第i个阵元之间的距离，c为超声传播速度，Fs为超声被动射频信号采样率。
5.根据权利要求1所述一种基于高阶孔径自相关的环形多阵列超声被动成像方法，其特征在于：所述步骤3)中，阵元合成信号和阵元信号总能量的计算公式分别表示为：
Figure FDA0003166941730000034
Figure FDA0003166941730000035
其中，i＝1,2,...,NE，NE为超声阵列换能器的阵元数目，k＝1,2,...,K，EBk(XP,ZP,n)为阵元合成信号，ETEk(XP,ZP,n)为阵元信号总能量，
Figure FDA0003166941730000036
为被动时间移位信号：
Figure FDA0003166941730000037
其中，n＝1,2,...,NS，NS为超声被动射频信号采样点数目，RFi k(n)为环形多阵列换能器组中第k个超声阵列换能器被动接收得到的超声被动射频信号，
Figure FDA0003166941730000038
为传播时间采样点数目。
6.根据权利要求1所述一种基于高阶孔径自相关的环形多阵列超声被动成像方法，其特征在于：所述步骤4)中，高阶孔径自相关加权因子的计算公式表示为：
Figure FDA0003166941730000041
其中，k＝1,2,...,K，HEBk(XP,ZP,n)为高阶孔径自相关合成信号，ETEk(XP,ZP,n)为阵元信号总能量。
7.根据权利要求1所述一种基于高阶孔径自相关的环形多阵列超声被动成像方法，其特征在于：所述步骤5)中，高阶孔径自相关加权合成信号的计算公式表示为：
HCWEBk(XP,ZP,n)＝EBk(XP,ZP,n)HCWFk(XP,ZP,n)
其中，k＝1,2,...,K，HCWEBk(XP,ZP,n)为高阶孔径自相关加权合成信号，EBk(XP,ZP,n)为阵元合成信号，HCWFk(XP,ZP,n)为高阶孔径自相关加权因子。
8.一种基于高阶孔径自相关的环形多阵列超声被动成像系统，其特征在于：该系统包括阵元合成信号与阵元信号总能量计算模块、高阶孔径自相关加权因子计算模块、单阵列超声被动成像模块以及环形多阵列超声被动成像模块；
所述阵元合成信号与阵元信号总能量计算模块用于利用成像区域中坐标为(XP,ZP)的任意一像素点与第k个超声阵列换能器各阵元之间的超声传播时间对该超声阵列换能器被动接收得到的超声被动射频信号进行时间移位处理，并将处理所得被动时间移位信号以及该信号的模平方分别沿着第k个超声阵列换能器的阵元方向进行叠加；
所述高阶孔径自相关加权因子计算模块用于对所述阵元合成信号与阵元信号总能量计算模块中的被动时间移位信号进行标准化处理以及对该标准化处理所得信号进行高阶孔径自相关处理，并计算高阶孔径自相关处理所得信号的模平方与所述被动时间移位信号的模平方沿阵元方向的叠加结果的比值；
标准化信号的计算公式表示为：
Figure FDA0003166941730000042
其中，i＝1,2,...,NE，NE为超声阵列换能器的阵元数目，k＝1,2,...,K，
Figure FDA0003166941730000043
为对被动时间移位信号进行符号化处理后所得的符号信号，
Figure FDA0003166941730000044
为对被动时间移位信号进行绝对值化处理后所得的绝对值信号，n＝1,2,...,NS，NS为超声被动射频信号采样点数目，P为高阶孔径自相关的阶数，P取大于等于2的正整数；
高阶孔径自相关合成信号的计算公式表示为：
Figure FDA0003166941730000051
上式中共有ItmNum项，ItmNum为将P拆分成若干个不超过P的正整数之和的拆分方案的数目，每一项的系数Coefb由多项式定理中累乘项的系数确定，P为高阶孔径自相关的阶数，P取大于等于2的正整数，k＝1,2,...,K，
Figure FDA0003166941730000052
为将标准化信号
Figure FDA0003166941730000053
的p次方沿着阵元方向叠加所得的p次阵元标准化合成信号，
Figure FDA0003166941730000054
的计算公式表示为：
Figure FDA0003166941730000055
其中，p＝1,2,...,P，k＝1,2,...,K；
所述单阵列超声被动成像模块用于对所述阵元合成信号与阵元信号总能量计算模块中的被动时间移位信号沿阵元方向的叠加结果与所述高阶孔径自相关加权因子计算模块计算的比值进行相乘，并将相乘得到的高阶孔径自相关加权合成信号的模平方沿着信号采样点方向进行叠加；
所述环形多阵列超声被动成像模块用于对所述单阵列超声被动成像模块在成像区域所有像素点处分别通过对高阶孔径自相关加权合成信号的模平方沿着信号采样点方向进行叠加而得到的针对每个超声阵列换能器的单阵列超声被动成像结果，进行累乘并计算累乘结果的K次算术方根。
CN202010768521.7A 2020-08-03 2020-08-03 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统 Active CN112023283B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010768521.7A CN112023283B (zh) 2020-08-03 2020-08-03 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010768521.7A CN112023283B (zh) 2020-08-03 2020-08-03 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统

Publications (2)

Publication Number Publication Date
CN112023283A CN112023283A (zh) 2020-12-04
CN112023283B true CN112023283B (zh) 2021-09-07

Family

ID=73583793

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010768521.7A Active CN112023283B (zh) 2020-08-03 2020-08-03 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统

Country Status (1)

Country Link
CN (1) CN112023283B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112754527B (zh) * 2020-12-28 2023-10-20 沈阳工业大学 一种用于低频超声胸腔成像的数据处理方法
CN114098799B (zh) * 2021-10-27 2023-06-27 西安交通大学 一种单脉冲内超声空化的快速低伪影实时动态成像方法与系统
CN113768541B (zh) * 2021-10-27 2024-02-13 之江实验室 一种复杂曲面超声阵列换能器阵元位置误差校正方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1279844A (zh) * 1997-10-06 2001-01-10 趣点公司 带有△∑反馈控制的成束超声成像器
CN108670304A (zh) * 2018-06-06 2018-10-19 东北大学 一种基于改进dmas算法的超声平面波成像方法
CN109513123A (zh) * 2018-12-28 2019-03-26 西安交通大学 一种基于半球阵的高分辨三维被动空化成像方法
CN111050655A (zh) * 2017-05-28 2020-04-21 利兰斯坦福初级大学董事会 通过非线性定位进行的超声成像
CN111265245A (zh) * 2020-01-20 2020-06-12 西安交通大学 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102279354B (zh) * 2011-07-01 2013-07-31 西安交通大学 用于变压器局部放电定位的十字形超声阵列传感器及方法
US9844359B2 (en) * 2013-09-13 2017-12-19 Decision Sciences Medical Company, LLC Coherent spread-spectrum coded waveforms in synthetic aperture image formation
US10786227B2 (en) * 2014-12-01 2020-09-29 National Institute Of Advanced Industrial Science And Technology System and method for ultrasound examination
CN107049252A (zh) * 2017-03-29 2017-08-18 华北电力大学（保定） 一种生物磁光声联合内窥成像方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1279844A (zh) * 1997-10-06 2001-01-10 趣点公司 带有△∑反馈控制的成束超声成像器
CN111050655A (zh) * 2017-05-28 2020-04-21 利兰斯坦福初级大学董事会 通过非线性定位进行的超声成像
CN108670304A (zh) * 2018-06-06 2018-10-19 东北大学 一种基于改进dmas算法的超声平面波成像方法
CN109513123A (zh) * 2018-12-28 2019-03-26 西安交通大学 一种基于半球阵的高分辨三维被动空化成像方法
CN111265245A (zh) * 2020-01-20 2020-06-12 西安交通大学 基于双约束鲁棒Capon波束合成和多重变迹互相关的被动空化成像方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
《基于线调频频散变换的多模式超声兰姆波时延估计》;行鸿彦 等;《声学学报》;20020831;全文 *
《眼科高频环阵列超声波数字化成像技术的研究》;王立伟 等;《北京生物医学工程》;20111030;全文 *

Also Published As

Publication number Publication date
CN112023283A (zh) 2020-12-04

Similar Documents

Publication Publication Date Title
CN112023283B (zh) 基于高阶孔径自相关的环形多阵列超声被动成像方法及系统
Cohen et al. Sparse convolutional beamforming for ultrasound imaging
JP6030207B2 (ja) 超音波合成イメージングの装置と方法
EP3043715B1 (en) Ultrasound imaging using apparent point-source transmit transducer
Wang et al. MVDR-based coherence weighting for high-frame-rate adaptive imaging
JP2013503681A (ja) 対側アレイベースの経頭蓋超音波収差補正
CN103536316A (zh) 一种空时平滑相干因子类自适应超声成像方法
US11650300B2 (en) Ultrasound system and method for suppressing noise using per-channel weighting
CN109513123B (zh) 一种基于半球阵的高分辨三维被动空化成像方法
US20180092627A1 (en) Ultrasound signal processing device, ultrasound signal processing method, and ultrasound diagnostic device
US11408987B2 (en) Ultrasonic imaging with multi-scale processing for grating lobe suppression
WO2019057592A1 (en) METHODS AND SYSTEMS FOR ULTRASONIC CONTRAST STRENGTHENING
JP2022023982A (ja) 超音波画像を処理するための方法及びシステム
US8394027B2 (en) Multi-plane/multi-slice processing for 2-D flow imaging in medical diagnostic ultrasound
Cohen et al. Sparse convolutional beamforming for 3-D ultrafast ultrasound imaging
Jørgensen et al. Row–column beamformer for fast volumetric imaging
Hazard et al. Effects of motion on a synthetic aperture beamformer for real-time 3D ultrasound
Ponnle et al. Suppression of grating lobe artifacts in ultrasound images formed from diverging transmitting beams by modulation of receiving beams
US11883240B2 (en) Sparse convolutional beamforming for ultrasound imaging
Xu et al. Adaptive minimum variance beamforming combined with phase coherence imaging for ultrasound imaging
CN110869799A (zh) 用于处理超声图像的方法和系统
Guo et al. Fourier-Domain Beamforming and Sub-Nyquist Sampling for Coherent Pixel-Based Ultrasound Imaging
Wang et al. A mixed transmitting–receiving beamformer with a robust generalized coherence factor: Enhanced resolution and contrast
Wang et al. Multiline acquisition beamforming for ultrasound computed tomography
Peralta et al. Optimization of transducer distribution and transmit sequence in coherent-multi transducer ultrasound (CoMTUS) imaging

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