一种高分辨快速计算的被动超声成像方法及系统

1.本发明属于超声检测与超声成像技术领域，具体涉及一种高分辨快速计算的被动超声成像方法及系统。


背景技术：

2.为保证超声治疗的精准性和安全性并促进其临床应用，需要发展可靠的医学影像监控技术。超声成像在超声治疗监控中具有广阔的应用前景，其可分为两类，第一类是换能器工作在发射接收模式的主动超声成像，第二类是换能器工作在不发射只接收模式的被动超声成像。主动超声成像由于其发射的探测脉冲与超声治疗信号互相干扰，因此只能和超声治疗过程交错穿插进行，不能实现真正的实时监控。而被动超声成像由于不发射探测脉冲，因此不存在信号干扰的问题，从而可在超声治疗的过程中进行实时监控；且为了临床应用推广，被动超声成像常使用线阵换能器等用于医学诊断的阵列超声换能器进行成像。被动超声成像在超声治疗中的应用主要体现在对空化的检测上，近年来被动超声成像在空化动力学研究以及空化介导的多种超声治疗应用的监控等方面发挥着越来越重要的作用。
3.然而，被动超声成像还存在两个需要解决的问题，即分辨性能差和计算速度慢。被动超声成像的分辨性能主要取决于阵列超声换能器的长度、接收信号的频率以及声源与换能器之间的距离，阵列超声换能器长度越大、接收信号频率越大或声源距离换能器越近，则成像分辨性能越好。而由于阵列超声换能器的长度会受到人体某些部位(例如心脏)的限制、接收信号频率会受到换能器带宽或组织高频衰减的限制、且声源深度在实际应用中可调整的空间并不是很大，因此成像分辨性能往往不能得到有效的改善；另外，组织异质性、换能器缺陷以及散射子之间的相互作用干扰也会造成被动超声成像分辨性能的下降。被动超声成像的计算速度取决于信号采样点数目、成像视野尺寸和像素尺寸。被动超声成像时首先需要对成像区域进行规划，然后通过延时叠加积分方法处理接收信号并计算成像区域内每个像素的像素值，最后由成像区域内的所有像素构成一幅图像，即通过逐个像素点计算来实现成像。在信号采样点数目给定的条件下，成像视野越大或者像素越小，计算速度越慢。而在实际应用中，往往需要选择较大的成像视野尺寸和较小的像素尺寸来进行成像，这就会增加计算的时间，从而不利于实时成像。
4.目前有研究人员提出了一些改进的被动超声成像算法，例如基于自适应波束合成、相干系数加权或孔径自相关的时域算法，以及基于频率叠加、自适应波束合成或角谱法的频域算法；然而这些算法要么只能改善图像质量而无法提高计算速度，要么只能提高计算速度而无法改善图像质量。如何同时提高被动超声成像的分辨性能和计算速度，一直是被动超声成像领域存在的一个技术难题，同时也是超声治疗监控的重要需求。鉴于此，亟待提出一种高分辨快速计算的被动超声成像方法及系统，从而为超声治疗过程中空化的准确实时监控提供强有力的工具，并为超声精准诊疗奠定坚实基础。


技术实现要素：

5.本发明的目的在于提供一种高分辨快速计算的被动超声成像方法及系统，可有效提高被动超声成像的分辨性能和计算速度。
6.为了实现上述目的，本发明采用了以下技术方案：
7.一种高分辨快速计算的被动超声成像方法，包括以下步骤：
8.1)利用阵列超声换能器依次被动地接收经连续的多个等间隔相位的超声脉冲信号作用于介质后产生的多个声辐射信号，对多个声辐射信号进行采样，得到多个单相位信号，将多个单相位信号叠加，得到相位叠加信号；
9.2)对相位叠加信号进行绝对值化处理，得到绝对值信号，对绝对值信号进行二值化处理，得到二值化信号，根据二值化信号分别建立以采样点索引为元素的阵元最小采样点集合和阵元最大采样点集合，对阵元最小采样点集合中元素的最小值进行下调，得到最小有效采样点的索引，对阵元最大采样点集合中元素的最大值进行上调，得到最大有效采样点的索引，根据最小有效采样点和最大有效采样点的索引在信号采样点方向上对相位叠加信号进行截断，得到全阵元有效信号；
10.3)分别对全阵元有效信号进行前向阵元截断处理和后向阵元截断处理，得到前向阵元截断有效信号和后向阵元截断有效信号；
11.4)根据扩展阵元数目和扩展信号采样点数目，分别沿阵元方向和信号采样点方向对全阵元有效信号、前向阵元截断有效信号和后向阵元截断有效信号进行补零，得到全阵元扩展信号、前向阵元截断扩展信号和后向阵元截断扩展信号，对全阵元扩展信号、前向阵元截断扩展信号和后向阵元截断扩展信号分别进行二维快速傅里叶变换，得到全阵元时空初始谱、前向阵元截断时空初始谱和后向阵元截断时空初始谱；
12.5)根据信号采样率和扩展信号采样点数目计算时间频率，根据阵列超声换能器的阵元间距和扩展阵元数目计算空间频率，根据时间频率和空间频率计算时空谱传递因子；
13.6)根据时空谱传递因子计算一个深度下的时空谱传递函数，将全阵元时空初始谱、前向阵元截断时空初始谱和后向阵元截断时空初始谱分别与计算得到的一个深度下的时空谱传递函数相乘，得到全阵元时空传递谱、前向阵元截断时空传递谱和后向阵元截断时空传递谱，沿阵元方向分别对全阵元时空传递谱、前向阵元截断时空传递谱和后向阵元截断时空传递谱进行快速傅里叶逆变换并将逆变换所得结果的模平方沿信号采样点方向叠加，得到全阵元能量、前向阵元截断能量和后向阵元截断能量，分别从全阵元能量、前向阵元截断能量和后向阵元截断能量中截取有效元素，得到全阵元有效能量、前向阵元截断有效能量和后向阵元截断有效能量，计算前向阵元截断有效能量和后向阵元截断有效能量的相关系数并对所得相关系数进行阈值化处理，得到阈值化相关系数，将全阵元有效能量与阈值化相关系数相乘，得到一个深度下的像素向量；
14.7)按照设置的深度下限、深度间隔和深度数目改变深度并重复步骤6)，直至得到所有深度下的像素向量，将所有深度下的像素向量构成像素矩阵，对像素矩阵进行归一化和对数化处理，得到被动超声成像结果。
15.优选的，所述步骤1)中，阵列超声换能器选自线阵换能器等诊断用阵列超声换能器中的一种或多种；被动接收是通过将阵列超声换能器的工作模式设置为不发射只接收来实现的。
16.优选的，所述步骤2)中，二值化处理所用的阈值为绝对值信号最大值的0.1～0.3倍。
17.优选的，所述步骤2)中，阵元最小采样点集合和阵元最大采样点集合的建立包括以下步骤：
18.2.1)分别初始化阵元最小采样点集合和阵元最大采样点集合为空集；
19.2.2)针对二值化信号中的一个二值化阵元信号进行如下处理：
20.若存在值等于1的采样点且其数量大于1，则将所有采样点的索引中的最小值存入阵元最小采样点集合minsampset，并将所有采样点的索引中的最大值存入阵元最大采样点集合maxsampset；
21.若存在值等于1的采样点且其数量等于1，则将该采样点的索引同时存入阵元最小采样点集合minsampset和阵元最大采样点集合maxsampset；
22.2.3)重复步骤2.2)，直至针对二值化信号中的所有二值化阵元信号的处理均完毕，从而得到阵元最小采样点集合minsampset和阵元最大采样点集合maxsampset。
23.优选的，所述步骤2)中，最小有效采样点和最大有效采样点的索引的计算公式表示为：
24.mineffecsamp＝floor(min{minsampset}
×
dar)
25.maxeffecsamp＝ceil(max{maxsampset}
×
uar)
26.其中，mineffecsamp和maxeffecsamp分别为最小有效采样点的索引和最大有效采样点的索引，floor(
·
)和ceil(
·
)分别表示向下取整和向上取整，min{
·
}和max{
·
}分别表示求最小值和求最大值，dar和uar分别为下调率和上调率，dar＝1-ar，uar＝1 ar，ar为0.1～0.2。
27.优选的，所述步骤3)中，前向阵元和后向阵元分别为阵列超声换能器的前一半阵元和后一半阵元。
28.优选的，所述步骤3)中，前向阵元截断有效信号和后向阵元截断有效信号表示为：
[0029][0030][0031]
其中，pmef(n,i)和pmeb(n,i)分别为前向阵元截断有效信号和后向阵元截断有效信号，pmea(n,i)为全阵元有效信号，n＝1,2,...,ns
ef
，i＝1,2,...,ne，ns
ef
为有效信号采样点数目，ne为阵列超声换能器的阵元数目，ne
cut
为阵列超声换能器的阵元数目的一半。
[0032]
优选的，所述步骤4)中，扩展阵元数目和扩展信号采样点数目的计算公式表示为：
[0033][0034][0035]
其中，ne
zc
和ns
zc
分别为扩展阵元数目和扩展信号采样点数目。
[0036]
优选的，所述步骤4)中，全阵元扩展信号、前向阵元截断扩展信号和后向阵元截断扩展信号表示为：
[0037][0038][0039][0040]
其中，pmza(n,i)、pmzf(n,i)和pmzb(n,i)分别为全阵元扩展信号、前向阵元截断扩展信号和后向阵元截断扩展信号，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
。
[0041]
优选的，所述步骤4)中，二维快速傅里叶变换是指沿阵元方向进行快速傅里叶变换后再沿信号采样点方向进行快速傅里叶变换。
[0042]
优选的，所述步骤5)中，时空谱传递因子的计算公式表示为：
[0043][0044]
其中，c为超声波传播速度，f
t
(n)和fs(i)分别为时间频率和空间频率，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
，r
samp
为信号采样率，pitch为阵列超声换能器的阵元间距。
[0045]
优选的，所述步骤6)中，时空谱传递函数的计算公式表示为：
[0046]
trsffunz(n,i)＝exp[j
·z·
trsffac(n,i)]
[0047]
其中，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
，j为虚数单位，z为深度，trsffac(n,i)为时空谱传递因子。
[0048]
优选的，所述步骤6)中，截取的有效元素为第1到第ne个元素。
[0049]
优选的，所述步骤6)中，前向阵元截断有效能量和后向阵元截断有效能量的相关系数的计算公式表示为：
[0050][0051]
其中，和分别为前向阵元截断有效能量和后向阵元截断有效能量。
[0052]
优选的，所述步骤6)中，相关系数的阈值化处理按照以下公式进行：
[0053][0054]
其中，为阈值化相关系数，ξ为阈值(为小于10-3
的正数)。
[0055]
一种高分辨快速计算的被动超声成像系统，该系统包括单相位信号叠加模块、有效采样点截断模块、前后向阵元截断模块、时空初始谱计算模块、时空谱传递因子计算模块、像素向量计算模块和像素矩阵计算及处理模块；
[0056]
所述单相位信号叠加模块用于执行上述步骤1)，主要是用于利用阵列超声换能器依次被动地接收经连续的多个等间隔相位的超声脉冲信号作用于介质后产生的多个声辐射信号，对多个声辐射信号进行采样，以及将采样所得多个单相位信号叠加；
[0057]
所述有效采样点截断模块用于执行上述步骤2)，主要是用于对单相位信号叠加模块所得相位叠加信号进行绝对值化处理，对所得绝对值信号进行二值化处理，根据所得二值化信号分别建立以采样点索引为元素的阵元最小采样点集合和阵元最大采样点集合，对阵元最小采样点集合中元素的最小值和阵元最大采样点集合中元素的最大值进行下调和上调，以及根据下调所得最小有效采样点的索引和上调所得最大有效采样点的索引在信号采样点方向上对相位叠加信号进行截断；
[0058]
所述前后向阵元截断模块用于执行上述步骤3)，主要是用于分别对有效采样点截断模块所得全阵元有效信号进行前向阵元截断处理和后向阵元截断处理；
[0059]
所述时空初始谱计算模块用于执行上述步骤4)，主要是用于根据扩展阵元数目和扩展信号采样点数目，分别沿阵元方向和信号采样点方向对有效采样点截断模块所得全阵元有效信号、前后向阵元截断模块所得前向阵元截断有效信号和后向阵元截断有效信号进行补零，以及对补零所得全阵元扩展信号、前向阵元截断扩展信号和后向阵元截断扩展信号进行二维快速傅里叶变换；
[0060]
所述时空谱传递因子计算模块用于执行上述步骤5)，主要是用于根据信号采样率和扩展信号采样点数目计算时间频率，根据阵列超声换能器的阵元间距和扩展阵元数目计算空间频率，以及根据时间频率和空间频率计算时空谱传递因子；
[0061]
所述像素向量计算模块用于执行上述步骤6)，主要是用于根据时空谱传递因子计算模块所得时空谱传递因子计算一个深度下的时空谱传递函数，将时空初始谱计算模块所得全阵元时空初始谱、前向阵元截断时空初始谱和后向阵元截断时空初始谱分别与计算得到的一个深度下的时空谱传递函数相乘，沿阵元方向分别对相乘所得全阵元时空传递谱、前向阵元截断时空传递谱和后向阵元截断时空传递谱进行快速傅里叶逆变换并将逆变换所得结果的模平方沿信号采样点方向叠加，分别从叠加所得全阵元能量、前向阵元截断能量和后向阵元截断能量中截取有效元素，计算通过从前向阵元截断能量和后向阵元截断能量中截取有效元素所得前向阵元截断有效能量和后向阵元截断有效能量的相关系数并对所得相关系数进行阈值化处理，以及将通过从全阵元能量中截取有效元素所得全阵元有效能量与阈值化处理所得阈值化相关系数相乘；
[0062]
所述像素矩阵计算及处理模块用于执行上述步骤7)，主要是用于按照设置的深度下限、深度间隔和深度数目改变像素向量计算模块中的深度，将像素向量计算模块所得所有深度下的像素向量构成像素矩阵，以及对像素矩阵进行归一化和对数化处理，从而得到
被动超声成像结果。
[0063]
本发明的有益效果体现在：
[0064]
本发明通过以下几个方面解决了现有被动超声成像方法存在分辨性能差且计算速度慢的问题：
[0065]
1)通过叠加k个单相位信号得到相位叠加信号；相比单相位信号，相位叠加信号中的k次谐波和k的整数次谐波得到增强同时其他次谐波得到抑制，即高频成分得到有效增强；而由于被动超声成像的分辨性能取决于信号频率的高低(频率越高，分辨性能越好)，因此单相位信号的叠加可有效提高被动超声成像的分辨性能。
[0066]
2)根据最小有效采样点和最大有效采样点的索引在信号采样点方向对相位叠加信号进行截断，使得无效采样点不参与后续的处理，从而可有效地减少计算时间，提高被动超声成像的计算速度；且最小有效采样点和最大有效采样点的索引是通过下调阵元最小采样点集合中元素的最小值和上调阵元最大采样点集合中元素的最大值分别得到的，可提高对相位叠加信号进行截断的准确性。
[0067]
3)对全阵元、前向阵元截断和后向阵元截断扩展信号分别进行的快速傅里叶变换，以及对全阵元、前向阵元截断和后向阵元截断时空传递谱分别进行的快速傅里叶逆变换，有助于提高被动超声成像的计算速度。
[0068]
4)通过将全阵元、前向阵元截断和后向阵元截断时空初始谱与时空谱传递函数相乘、然后沿阵元方向分别对各自的时空传递谱进行快速傅里叶逆变换并将各自所得结果的模平方沿信号采样点方向叠加，能够一次性获得任意一深度下所有阵元位置处对应的像素值，即采用的是逐行计算模式；相比传统被动超声成像算法中采用的逐点计算模式，逐行计算模式可提高被动超声成像的计算速度。
[0069]
5)利用前向阵元截断和后向阵元截断的有效能量在声源位置处相关性高而在非声源位置处相关性低的特点，通过计算前向阵元截断和后向阵元截断的有效能量的相关系数并进行阈值化处理，提供了一种沿着深度方向的自适应权重，将全阵元有效能量与阈值化相关系数相乘，可有效地抑制旁瓣伪影，从而可有效地提高被动超声成像的分辨性能。
[0070]
进一步地，本发明中在对全阵元有效信号进行前向阵元截断处理和后向阵元截断处理时采用的前向阵元和后向阵元分别为阵列超声换能器的前一半阵元和后一半阵元，可有效减小主瓣尺寸并滤除旁瓣伪影，从而可获得更好的成像分辨性能。
[0071]
进一步地，本发明中根据扩展阵元数目分别对全阵元、前向阵元截断和后向阵元截断有效信号沿阵元方向补零，有助于降低被动超声成像的背景水平，从而提高被动超声成像的分辨性能；本发明中根据扩展信号采样点数目分别对全阵元、前向阵元截断和后向阵元截断有效信号沿信号采样点方向补零至2的幂次，有助于提高快速傅里叶变换的运算性能，从而提高被动超声成像的计算速度。
[0072]
进一步地，本发明中对前向阵元截断和后向阵元截断的有效能量的相关系数进行阈值化处理时所采用的阈值为一较小的正数，一方面可避免对数化处理时产生无意义的像素值，另一方面也有助于更彻底地滤除旁瓣伪影。
附图说明
[0073]
图1为本发明实施例中生成相位叠加信号的示意图。
[0074]
图2为本发明实施例中对相位叠加信号进行截断的流程图。
[0075]
图3为本发明实施例中全阵元有效信号(a)、前向阵元截断有效信号(b)和后向阵元截断有效信号(c)的示意图。
[0076]
图4为本发明实施例中计算全阵元时空初始谱、前向阵元截断时空初始谱和后向阵元截断时空初始谱的流程图。
[0077]
图5为本发明实施例中计算时空谱传递因子的流程图。
[0078]
图6为本发明实施例中计算任意一深度下的像素向量的流程图。
[0079]
图7为本发明实施例中计算像素矩阵并对其进行处理的流程图。
[0080]
图8为三种情况下的被动超声成像结果：(a)单相位超声脉冲作用且不利用前后向阵元截断有效能量的阈值化相关系数；(b)多相位超声脉冲作用且不利用前后向阵元截断有效能量的阈值化相关系数；(c)多相位超声脉冲作用且利用前后向阵元截断有效能量的阈值化相关系数。
[0081]
图9为传统的被动超声成像方法(a)和本发明提出的被动超声成像方法(b)所得的成像结果。
具体实施方式
[0082]
下面结合附图和实施例对本发明作进一步的详细说明。
[0083]
本发明提出了一种高分辨快速计算的被动超声成像方法，并通过数值仿真验证其性能。
[0084]
(一)被动超声成像方法
[0085]
参见图1，对在多相位超声脉冲连续作用过程中阵列超声换能器被动接收的声辐射信号进行采样，以及对多个单相位信号进行叠加，具体流程见下述步骤(1.1)～(1.3)；
[0086]
(1.1)以下式所示的连续的k个等间隔相位的超声脉冲信号作用于介质(例如，水、生物组织、流动微泡或液滴)：
[0087][0088]
其中，k＝1,2,...,k，k为相位数目(例如，6)，p0为峰值声压(例如，1mpa)，f0为超声脉冲信号的频率(例如，1.2mhz)，t为时间(例如，作用时长为12个周期)，为初始相位(例如，0)；
[0089]
(1.2)由阵列超声换能器(例如，阵元数目ne＝128、阵元间距pitch＝0.3mm的线阵换能器)依次被动(即不向外部发射探测脉冲)地接收上述k个等间隔相位的超声脉冲信号作用于介质后产生的k个声辐射信号，按照一定的信号采样率r
samp
(例如，50mhz)和信号采样点数目ns(例如，4000)对这k个声辐射信号采样后得到k个单相位信号，记为psk(n,i)，其中k＝1,2,...,k，n＝1,2,...,ns，i＝1,2,...,ne；单相位信号psk(n,i)中包含着对应于超声脉冲信号频率f0的基波、二次谐波和高次谐波等；
[0090]
(1.3)将步骤(1.2)所得的与k个等间隔相位对应的k个单相位信号psk(n,i)叠加，得到相位叠加信号pm(n,i)：
[0091][0092]
其中，n＝1,2,...,ns，i＝1,2,...,ne。
[0093]
参见图2，对相位叠加信号进行绝对值化处理及二值化处理，建立阵元最小采样点集合和阵元最大采样点集合，根据下调率和上调率计算最小有效采样点和最大有效采样点，然后对相位叠加信号进行截断得到全阵元有效信号，具体流程见下述步骤(2.1)～(2.4)；
[0094]
(2.1)对步骤(1.3)所得相位叠加信号pm(n,i)进行绝对值化处理，得到绝对值信号然后对绝对值信号进行二值化处理，得到二值化信号
[0095][0096]
其中，n＝1,2,...,ns，i＝1,2,...,ne，thre为二值化阈值，thre设置为绝对值信号最大值的0.1～0.3倍；
[0097]
(2.2)建立阵元最小采样点集合和阵元最大采样点集合
[0098]
(2.2.1)分别初始化阵元最小采样点集合minsampset和阵元最大采样点集合maxsampset为空集；
[0099]
(2.2.2)针对步骤(2.1)所得二值化信号中的第i个二值化阵元信号：
[0100]
若存在值等于1的采样点且其数量大于1，则将所有这些值等于1的采样点的索引中的最小值存入阵元最小采样点集合minsampset，并将所有这些值等于1的采样点的索引中的最大值存入阵元最大采样点集合maxsampset；
[0101]
若存在值等于1的采样点且其数量等于1，则将该采样点的索引同时存入阵元最小采样点集合minsampset和阵元最大采样点集合maxsampset；
[0102]
若不存在值等于1的采样点，则不执行任何操作；
[0103]
(2.2.3)重复步骤(2.2.2)，直至将ne个二值化阵元信号均处理完毕，从而得到阵元最小采样点集合minsampset和阵元最大采样点集合maxsampset；
[0104]
(2.3)根据下调率对步骤(2.2)所得阵元最小采样点集合minsampset中元素的最小值进行下调，得到最小有效采样点的索引mineffecsamp；根据上调率对阵元最大采样点集合maxsampset中元素的最大值进行上调，得到最大有效采样点的索引maxeffecsamp：
[0105]
mineffecsamp＝floor(min{minsampset}
×
dar)
[0106]
maxeffecsamp＝ceil(max{maxsampset}
×
uar)
[0107]
其中，floor(
·
)和ceil(
·
)分别表示向下取整和向上取整，min{
·
}和max{
·
}分别表示求最小值和求最大值，dar和uar分别为下调率和上调率，dar＝1-ar，uar＝1 ar，ar设置为0.1～0.2；
[0108]
(2.4)根据步骤(2.3)所得最小有效采样点的索引mineffecsamp和最大有效采样点的索引maxeffecsamp在信号采样点方向上对步骤(1.3)所得相位叠加信号pm(n,i)进行截断，得到全阵元有效信号pmea(n,i)，其中，n＝1,2,...,ns
ef
，i＝1,2,...,ne，ns
ef
为有效信号采样点数目，ns
ef
＝maxeffecsamp-mineffecsamp 1。
[0109]
参见图3，根据从阵列超声换能器各阵元中选出的截断阵元所对应的索引，分别对全阵元有效信号进行前向阵元截断处理和后向阵元截断处理，得到前向阵元截断有效信号和后向阵元截断有效信号，具体流程见下述步骤(3.1)～(3.3)；
[0110]
(3.1)设置截断阵元对应的索引ne
cut
(例如，ne
cut
为阵元数目ne的一半)；
[0111]
(3.2)对步骤(2.4)所得全阵元有效信号pmea(n,i)进行前向阵元截断处理，即保持前ne
cut
个阵元对应的信号不变并将后ne-ne
cut
个阵元对应的信号置为零，得到前向阵元截断有效信号pmef(n,i)：
[0112][0113]
其中，n＝1,2,...,ns
ef
，i＝1,2,...,ne；
[0114]
(3.3)对步骤(2.4)所得全阵元有效信号pmea(n,i)进行后向阵元截断处理，即将前ne
cut
个阵元对应的信号置为零并保持后ne-ne
cut
个阵元对应的信号不变，得到后向阵元截断有效信号pmeb(n,i)：
[0115][0116]
其中，n＝1,2,...,ns
ef
，i＝1,2,...,ne。
[0117]
参见图4，分别沿阵元方向和信号采样点方向对全阵元有效信号、前向阵元截断有效信号和后向阵元截断有效信号进行补零，然后沿阵元方向和信号采样点方向进行快速傅里叶变换，得到全阵元时空初始谱、前向阵元截断时空初始谱和后向阵元截断时空初始谱，具体流程见下述步骤(4.1)～(4.5)；
[0118]
(4.1)设置扩展阵元数目ne
zc
：
[0119][0120]
其中，ne为阵列超声换能器的的阵元数目，ceil(
·
)表示向上取整；
[0121]
(4.2)设置扩展信号采样点数目ns
zc
：
[0122][0123]
其中，ns
ef
为步骤(2.4)所述的有效信号采样点数目；
[0124]
(4.3)根据步骤(4.1)和(4.2)设置的扩展阵元数目ne
zc
和扩展信号采样点数目ns
zc
，分别沿阵元方向和信号采样点方向对步骤(2.4)所得全阵元有效信号pmea(n,i)、步骤(3.2)所得前向阵元截断有效信号pmef(n,i)和步骤(3.3)所得后向阵元截断有效信号pmeb(n,i)进行补零，得到全阵元扩展信号pmza(n,i)、前向阵元截断扩展信号pmzf(n,i)和后向阵元截断扩展信号pmzb(n,i)：
[0125]
[0126][0127][0128]
其中，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
；
[0129]
(4.4)沿阵元方向对步骤(4.3)所得全阵元扩展信号pmza(n,i)、前向阵元截断扩展信号pmzf(n,i)和后向阵元截断扩展信号pmzb(n,i)进行快速傅里叶变换，分别得到全阵元空间谱sa(n,i)、前向阵元截断空间谱sf(n,i)和后向阵元截断空间谱sb(n,i)：
[0130]
sa(n,i)＝f
elem
[pmza(n,i)]
[0131]
sf(n,i)＝f
elem
[pmzf(n,i)]
[0132]
sb(n,i)＝f
elem
[pmzb(n,i)]
[0133]
其中，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
，f
elem
代表沿阵元方向的快速傅里叶变换；
[0134]
(4.5)沿信号采样点方向对步骤(4.4)所得全阵元空间谱sa(n,i)、前向阵元截断空间谱sf(n,i)和后向阵元截断空间谱sb(n,i)进行快速傅里叶变换，分别得到全阵元时空初始谱tsa0(n,i)、前向阵元截断时空初始谱tsf0(n,i)和后向阵元截断时空初始谱tsb0(n,i)：
[0135]
tsa0(n,i)＝f
samp
[sa(n,i)]
[0136]
tsf0(n,i)＝f
samp
[sf(n,i)]
[0137]
tsb0(n,i)＝f
samp
[sb(n,i)]
[0138]
其中，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
，f
samp
代表沿信号采样点方向的快速傅里叶变换。
[0139]
参见图5，计算时间频率和空间频率，然后计算时空谱传递因子，具体流程见下述步骤(5.1)～(5.3)；
[0140]
(5.1)根据步骤(1.2)所述信号采样率r
samp
和步骤(4.2)设置的扩展信号采样点数目ns
zc
计算时间频率f
t
(n)：
[0141][0142]
其中，n＝1,2,...,ns
zc
；
[0143]
(5.2)根据阵元间距pitch和步骤(4.1)设置的扩展阵元数目ne
zc
计算空间频率fs(i)：
[0144][0145]
其中，i＝1,2,...,ne
zc
；
[0146]
(5.3)根据步骤(5.1)所得时间频率f
t
(n)和(5.2)所得空间频率fs(i)计算时空谱
传递因子trsffac(n,i)：
[0147][0148]
其中，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
，c为超声波传播速度(例如，1480m/s)。
[0149]
参见图6，计算时空谱传递函数并将全阵元时空初始谱、前向阵元截断时空初始谱和后向阵元截断时空初始谱分别与时空谱传递函数相乘，沿阵元方向分别对全阵元时空传递谱、前向阵元截断时空传递谱和后向阵元截断时空传递谱进行快速傅里叶逆变换并将逆变换所得结果的模平方沿信号采样点方向叠加，分别从全阵元能量、前向阵元截断能量和后向阵元截断能量中截取有效元素，计算前向阵元截断有效能量和后向阵元截断有效能量的相关系数并对所得相关系数进行阈值化处理，将全阵元有效能量与阈值化相关系数相乘，得到像素向量，具体流程见下述步骤(6.1)～(6.7)；
[0150]
(6.1)根据步骤(5.3)所得时空谱传递因子trsffac(n,i)计算任意一深度z下的时空谱传递函数trsffunz(n,i)：
[0151]
trsffunz(n,i)＝exp[j
·z·
trsffac(n,i)]
[0152]
其中，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
，j为虚数单位；
[0153]
(6.2)将步骤(4.5)所得全阵元时空初始谱tsa0(n,i)、前向阵元截断时空初始谱tsf0(n,i)和后向阵元截断时空初始谱tsb0(n,i)分别与步骤(6.1)所得时空谱传递函数trsffunz(n,i)相乘，得到全阵元时空传递谱tsaz(n,i)、前向阵元截断时空传递谱tsfz(n,i)和后向阵元截断时空传递谱tsbz(n,i)：
[0154]
tsaz(n,i)＝tsa0(n,i)
·
trsffunz(n,i)
[0155]
tsfz(n,i)＝tsf0(n,i)
·
trsffunz(n,i)
[0156]
tsbz(n,i)＝tsb0(n,i)
·
trsffunz(n,i)
[0157]
其中，n＝1,2,...,ns
zc
，i＝1,2,...,ne
zc
；
[0158]
(6.3)沿阵元方向分别对步骤(6.2)所得全阵元时空传递谱tsaz(n,i)、前向阵元截断时空传递谱tsfz(n,i)和后向阵元截断时空传递谱tsbz(n,i)进行快速傅里叶逆变换，并将快速傅里叶逆变换所得结果的模平方沿信号采样点方向叠加，得到全阵元能量eaz(i)、前向阵元截断能量efz(i)和后向阵元截断能量ebz(i)：
[0159][0160][0161][0162]
其中，i＝1,2,...,ne
zc
，代表沿阵元方向的快速傅里叶逆变换；
[0163]
(6.4)分别从步骤(6.3)所得全阵元能量eaz(i)、前向阵元截断能量efz(i)和后向阵元截断能量ebz(i)中截取有效的第1～ne个元素，得到全阵元有效能量前向阵元截断有效能量和后向阵元截断有效能量其中，i＝1,2,...,ne；
[0164]
(6.5)计算步骤(6.4)所得前向阵元截断有效能量和后向阵元截断有效能量的相关系数γz：
[0165][0166]
(6.6)对步骤(6.5)所得相关系数γz进行阈值化处理，得到阈值化相关系数
[0167][0168]
其中，ξ为阈值(为一较小的正实数，例如，10-4
)；
[0169]
(6.7)将步骤(6.4)所得全阵元有效能量与步骤(6.6)所得阈值化相关系数相乘，得到任意一深度z下的像素向量pvz(i)：
[0170][0171]
其中，i＝1,2,...,ne。
[0172]
参见图7，设置被动超声成像的深度下限、深度间隔和深度数目，计算所有深度下的像素向量并构成像素矩阵，对像素矩阵进行归一化和对数化处理，得到被动超声成像结果，具体流程见下述步骤(7.1)～(7.3)；
[0173]
(7.1)设置被动超声成像的深度下限z
min
(例如，20mm)、深度间隔z
step
(例如，0.5mm)和深度数目nz(例如，121)；
[0174]
(7.2)从步骤(7.1)所述深度下限z
min
开始、以步骤(7.1)所述深度间隔z
step
为步长，重复步骤(6.1)～(6.7)，得到nz个深度下的像素向量；由nz个深度下的像素向量构成像素矩阵pm(z,i)：
[0175][0176]
其中，z＝z
min
,z
min
z
step
,...,z
min
(nz-1)z
step
，i＝1,2,...,ne；
[0177]
(7.3)对步骤(7.2)所得像素矩阵pm(z,i)进行归一化和对数化处理，得到标准像素矩阵npm(z,i)，即被动超声成像结果：
[0178]
npm(z,i)＝10lg{pm(z,i)/max[pm(z,i)]}
[0179]
其中，z＝z
min
,z
min
z
step
,...,z
min
(nz-1)z
step
，i＝1,2,...,ne，max[
·
]表示求最大值。
[0180]
(二)仿真及结果
[0181]
所采用的仿真模型为yang-church模型，声源所在位置为(-2.5,50mm)。为定量评价本发明提出的被动超声成像方法的分辨性能，采用主瓣面积和旁瓣水平两个指标，其中主瓣面积定义为步骤(7.3)所得被动超声成像结果中像素值大于-3db的所有像素的面积，旁瓣水平定义为步骤(7.3)所得被动超声成像结果中小于-3db的所有像素值之和；主瓣面积越小且旁瓣水平越低，说明成像分辨性能越好。为定量评价本发明提出的被动超声成像
方法的计算速度，分别记录实现成像方法的程序运行的开始时间和结束时间，程序运行时间越短，说明计算速度越快。
[0182]
参见图8，在三种情况下进行被动超声成像(动态范围设置为20db)：
[0183]
(a)单相位超声脉冲作用，此时，步骤(1.1)中的相位数目k为1，步骤(1.3)中的相位叠加信号pm(n,i)为单相位信号ps1(n,i)；且不利用前后向阵元截断有效能量的阈值化相关系数，此时，步骤(6.7)中的像素向量pvz(i)直接为全阵元有效能量
[0184]
(b)多相位超声脉冲作用，此时，步骤(1.1)中的相位数目k为6；且不利用前后向阵元截断有效能量的阈值化相关系数，此时，步骤(6.7)中的像素向量pvz(i)直接为全阵元有效能量
[0185]
(c)多相位超声脉冲作用，此时，步骤(1.1)中的相位数目k为6；且利用前后向阵元截断有效能量的阈值化相关系数，此时，步骤(6.7)中的像素向量pvz(i)为全阵元有效能量与阈值化相关系数的乘积。
[0186]
情况(a)的成像结果(图8a)中，中心亮斑较大且x型旁瓣伪影较大，主瓣面积和旁瓣水平分别为1.80mm2和-2.79
×
105db；相比情况(a)，情况(b)的成像结果(图8b)中，中心亮斑减小且x型旁瓣伪影减小，主瓣面积和旁瓣水平分别减小0.90mm2和0.30
×
105db；相比情况(b)，情况(c)的成像结果(图8c)中，中心亮斑减小且x型旁瓣伪影减小，主瓣面积和旁瓣水平分别减小0.30mm2和4.90
×
105db。以上不同情况的成像结果说明，在本发明提出的被动超声成像方法中，多相位超声脉冲作用和前后向阵元截断有效能量的阈值化相关系数有利于减小主瓣面积和旁瓣水平。
[0187]
参见图9，图9a为传统的被动超声成像方法(即使用延时叠加积分算法处理单相位超声脉冲作用下所得的单相位信号)所得的成像结果；图9b为本发明提出的被动超声成像方法所得的成像结果(与图8c相同)；动态范围设置为20db。结果表明，传统的被动超声成像方法所得成像结果(图9a)的中心亮斑较大且x型旁瓣伪影较大，而本发明提出的被动超声成像方法所得成像结果(图9b)的中心亮斑减小且x型旁瓣伪影减小(图9a和图9b的主瓣面积分别为2.25mm2和0.60mm2，旁瓣水平分别为-3.59
×
105db和-7.99
×
105db)；另外，传统的被动超声成像方法的运行时间为80.46秒，而本发明提出的被动超声成像方法的运行时间仅为2.37秒。以上结果说明，本发明提出的被动超声成像方法相比传统的被动超声成像方法具有更好的分辨性能同时具有更快的计算速度。
[0188]
(三)本发明提出的被动超声成像方法的优点：
[0189]
(1)传统被动超声成像中处理的是单相位超声脉冲作用下得到的单相位信号，而本发明则对多相位超声脉冲作用下的多个单相位信号进行叠加，这使得谐波成分得到有效增强，从而可增强高频成分；而由于被动超声成像的分辨性能依赖于信号频率的大小，因此通过叠加多个单相位信号并利用所得相位叠加信号进行成像可有效地提高被动超声成像的分辨性能。
[0190]
(2)本发明利用了前向阵元截断有效能量和后向阵元截断有效能量在声源位置处相关性高而在非声源位置处相关性低的特点，二者之间的相关系数经阈值化处理后提供了一种深度方向的自适应权重，将全阵元有效能量与阈值化相关系数相乘能够有效降低非声源位置处的像素值，从而可减小主瓣尺寸并抑制旁瓣伪影，因此进一步提高了被动超声成
像的分辨性能。
[0191]
(3)本发明仅对最小有效采样点和最大有效采样点之间的有效采样段进行处理，避免了有效采样段之外的无效采样点增加计算量的问题，从而可提高被动超声成像的计算速度；另外，本发明中对全阵元、前向阵元截断和后向阵元截断扩展信号和时空传递谱分别进行的快速傅里叶变换和快速傅里叶逆变换，有利于提高被动超声成像的计算速度。
[0192]
(4)本发明摒弃了传统被动超声成像方法所采用的逐点计算模式，而是采用了逐行计算模式，即通过将全阵元、前向阵元截断和后向阵元截断时空初始谱与时空谱传递函数相乘、然后沿阵元方向分别对各自的时空传递谱进行快速傅里叶逆变换并将各自所得结果的模平方沿信号采样点方向叠加，然后一次性获得任意一深度下所有阵元位置处对应的像素值，可大幅度减小计算量、提高计算效率，从而可提高被动超声成像的计算速度。
[0193]
(5)本发明可同时提高被动超声成像的分辨性能和计算速度，可对组织消融、组织毁损、血脑屏障开放等多种超声治疗应用中空化的空间信息进行高分辨且快速的重建，一方面为超声治疗过程中空化的准确监控提供了有力工具并为精准超声治疗奠定了基础，另一方面为推进被动超声成像在超声成像系统上的实时实施及其在超声治疗中的应用奠定了基础。此外，本发明提出的被动超声成像方法也为其他被动声学领域(例如，水下声源被动定位)的技术发展提供了参考和借鉴。