基于特征空间自适应波束合成的频域被动空化成像及频率复合成像方法与流程

本发明属于超声检测与超声成像技术领域，具体涉及一种基于特征空间自适应波束合成的频域被动空化成像及频率复合成像方法。

背景技术：

在生物医学超声领域，超声空化一直是研究的热点，在肿瘤热消融、体外碎石、超声溶栓以及药物释放等聚焦超声治疗中空化都起了至关重要的作用。在肿瘤热消融中，空化微泡会对聚焦超声波束形成遮挡从而影响聚焦性能，导致损伤区域不可预测，降低治疗的准确性；而另一方面空化可提高组织的声热转换效率，从而降低治疗能量并缩短治疗时间。在碎石中，空化微泡瞬间坍塌时释放的能量在局部区域内产生的高温、高压、冲击波和射流等物理现象加快了结石的破碎。在超声溶栓中，空化微泡振动以及破裂产生的高速射流以及冲击波等可有效地分解血栓，从而提高溶栓效率。然而由于空化的不可控性和随机性，可能会对正常组织形成不可逆转的损伤，达不到治疗目的。此外，研究表明，即使是在诊断超声声压下，空化微泡形成的剪切力和射流也可对周围细胞或组织产生生物效应，从而产生血管损伤和破裂的风险。因此有效地控制和利用空化对于聚焦超声治疗方案的调控与优化有着重要意义，而要实现这一点则必须要依赖于空化检测技术的发展。

近年来，空化检测技术从一维发展到了二维，利用临床上常用的线阵和相控阵等超声诊断换能器可实现对空化空间分布的定征。目前已发展了传统b模式超声成像、平面波超快速超声成像以及基于时序控制的超声逐次扫描成像等主动空化成像方法，通过调节发射模式同时改进成像算法性能可实现高对比度和高分辨率的空化成像。然而，由于在聚焦超声治疗过程中会对主动空化成像产生强烈的干扰，因此只能选择在聚焦超声完全停止作用之后或是在下一次聚焦超声脉冲还未到来前进行成像，而正由于这一点，主动空化成像只能检测聚焦超声不作用的时候产生的空化微泡，而不能检测聚焦超声作用过程中产生的空化行为，也就不能实现聚焦超声治疗的实时监控成像。为了弥补主动空化成像的不足，近年来有学者提出了一种被动空化成像方法，该方法由于不发射只被动接收空化信号，因此不受聚焦超声信号的干扰，从而可以实现聚焦超声治疗过程中空化的有效定征。被动空化成像已被用于监控组织热消融、组织毁损、药物释放以及血脑屏障开放等方面。目前最流行也是最常用的被动空化成像方法是在时域中进行，之后有学者提出了频域的被动空化成像。相比时域被动空化成像，频域被动空化成像可直接在频域选择成像对象是稳态空化或是惯性空化，只有被选择的频带信息做波束合成，因此相比时域方法可降低计算量；另外，频域被动空化成像直接通过傅里叶变换实现信号的时间延时，可以避免由于低采样率的离散时间采样造成的延时误差。然而传统的频域被动空化成像方法从本质上来讲是一种不依赖于数据本身特性的盲波束合成方法，会造成极为严重的成像伪影，成像质量不佳，不能提高空化空间定位的精准度，因此不利于聚焦超声治疗过程的实时监控。因此如何实现高分辨率的频域被动空化成像是本领域的一大挑战，也是亟需克服的一个难点。

技术实现要素：

本发明的目的在于提供一种基于特征空间自适应波束合成的频域被动空化成像及频率复合成像方法。

为了实现上述目的，本发明采用了以下技术方案：

一种基于特征空间自适应波束合成的频域被动空化成像方法，包括以下步骤：

步骤一：由超声成像换能器被动接收空化射频信号，对空化射频信号进行傅里叶变换得到频域空化信号，对频域空化信号作相位偏移处理后利用归一化余弦变迹函数进行变迹处理，得到超声成像换能器所有阵元的经过相位偏移和阵元变迹处理后的频域空化信号；

步骤二：根据步骤一所得的所有阵元的经过相位偏移和阵元变迹处理后的频域空化信号构造协方差矩阵，利用该协方差矩阵及导向矢量不确定集参数构建自适应波束合成的代价函数，通过对所述代价函数的约束条件引入拉格朗日乘数构造格朗日约束函数，对拉格朗日约束函数求导并令导函数为零，得到导向矢量表达式，将导向矢量表达式代入所述代价函数的约束条件中并利用牛顿迭代法求解得到拉格朗日乘数，根据求解得到的拉格朗日乘数计算导向矢量并将导向矢量作归一化处理，根据归一化处理后的导向矢量和所述协方差矩阵得到自适应波束合成的加权矢量；

步骤三：对所述协方差矩阵进行特征空间分解，得到信号子空间，将步骤二所得的自适应波束合成的加权矢量投影到信号子空间上，利用投影后的加权矢量对步骤一所得的所有阵元的经过相位偏移和阵元变迹处理后的频域空化信号进行加权求和，然后去除直流分量，得到特征空间自适应波束合成的输出信号；对特征空间自适应波束合成的输出信号进行平方并在整个频域上进行叠加，得到全频域的被动空化成像结果；

步骤四：根据空化伪影比选择最优的导向矢量不确定集参数，根据最优的导向矢量不确定集参数确定特征空间自适应波束合成的最优输出信号；对该最优输出信号分别在整个频域中的不同子频域上进行叠加，得到不同子频域的被动空化成像结果。

一种基于特征空间自适应波束合成的被动空化频率复合成像方法，包括以下步骤：

步骤一、步骤二、步骤三及步骤四与上述频域被动空化成像方法相同；

步骤五：对所述不同子频域的被动空化成像结果做标准化处理，然后将其中两个以上子频域的经过标准化处理的被动空化成像结果进行叠加，得到被动空化频率复合成像结果。

上述步骤一，具体包括以下步骤：

1.1)利用超声成像换能器被动接收时域空化信号，所述时域空化信号是指聚焦超声换能器焦域处散射的空化射频信号；

1.2)对超声成像换能器第i个阵元接收到的时域空化信号进行傅里叶变换，得到频域空化信号：

其中，i＝1,2,...,n，n为超声成像换能器阵元数，t1是采集时域空化信号的开始时刻，t2是采集时域空化信号的停止时刻，j是虚数单位，pi(t,x,z)是第i个阵元接收到的时域空化信号，f表示频率；

1.3)对步骤1.2)所得频域空化信号进行相位偏移处理：

其中，pofi(f,x,z)为第i个阵元的相位偏移；

1.4)对步骤1.3)所得的相位偏移处理后的频域空化信号进行阵元变迹处理：

其中，apodi(f,x,z)为第i个阵元的归一化余弦变迹函数；

1.5)重复步骤1.2)～1.4)，直到得到超声成像换能器的n个阵元的经过相位偏移和阵元变迹处理后的频域空化信号。

上述步骤二，具体包括以下步骤：

2.1)利用步骤一所得的所有阵元的经过相位偏移和阵元变迹处理后的频域空化信号构建协方差矩阵r：

其中，s(f,x,z)＝[s1(f,x,z)；s2(f,x,z)；...；sn(f,x,z)]为超声成像换能器的n个阵元的经过相位偏移和阵元变迹处理后的频域空化信号矩阵，表示对于任意的频率f；

2.2)利用步骤2.1)所得协方差矩阵构建自适应波束合成的代价函数：

其中，a为要求解的导向矢量，为假定的导向矢量，||·||代表欧几里得范数，为约束条件，ε为导向矢量不确定集参数；

2.3)经过步骤2.2)后，引入拉格朗日乘数并构造拉格朗日约束函数：

其中，λ为拉格朗日乘数；

2.4)将步骤2.3)所得的拉格朗日约束函数对导向矢量a求导并令求导所得导函数为零，然后利用矩阵求逆引理，得到导向矢量表达式：

其中，i为单位矩阵；

2.5)将步骤2.4)所得的导向矢量表达式代入步骤2.2)所得的代价函数的约束条件中，则约束条件改写为：

2.6)通过步骤2.1)所得的协方差矩阵r的特征空间分解对步骤2.5)所得的约束条件进行转化，然后采用牛顿迭代法来求解转化后的约束条件中的拉格朗日乘数λ，然后通过步骤2.4)所得的导向矢量表达式计算出导向矢量，并做归一化处理，得到归一化导向矢量

2.7)根据步骤2.6)所得的归一化导向矢量和步骤2.1)所得的协方差矩阵r计算自适应波束合成的加权矢量w：

上述步骤三，具体包括以下步骤：

3.1)对根据步骤一得到的所有阵元的经过相位偏移和阵元变迹处理后的频域空化信号所构建的协方差矩阵进行特征空间分解：

r＝uλu^h＝rs+ri

其中，u＝[u1,u2,...,un]为特征向量矩阵，λ＝diag[γ1,γ2,...,γn]为对角矩阵，该矩阵对角元素为特征值，且γ1≥γ2≥...≥γn，rs为信号子空间的协方差矩阵，ri为噪声子空间的协方差矩阵；

3.2)将自适应波束合成的加权矢量投影到信号子空间，得到投影后的加权矢量

其中，为归一化导向矢量，us＝[u1,u2,...,ul]为前l个较大特征值对应的特征向量构成的信号子空间，λs＝diag[γ1,γ2,...,γl]为信号子空间的特征值矩阵，diag[·]表示矢量的对角化；l为特征值中大于最大特征值δ倍的所有特征值的数目，δ为0.2～0.5；

3.3)利用超声成像换能器每个阵元的经过相位偏移和阵元变迹处理后的频域空化信号和步骤3.2)所得的投影后的加权矢量计算所述每个阵元的加权频域空化信号：

其中，为步骤3.2)所得的投影后的加权矢量的第i个元素，si(f,x,z)为超声成像换能器第i个阵元的经过相位偏移和阵元变迹处理后的频域空化信号，其中，i＝1,2,...,n，n为超声成像换能器阵元数；

3.4)对超声成像换能器n个阵元的加权频域空化信号进行叠加，得到信号q(f,x,z)：

3.5)从步骤3.4)所得信号q(f,x,z)中去除掉直流分量，则得到特征空间自适应波束合成的输出信号对进行平方并在整个频域上进行叠加，得到每个成像位置(x,z)的能量：

其中，表示对于任意的频率f。

上述步骤四，具体包括以下步骤：

4.1)构造基于平均能量的空化伪影比指标car：

car＝20log10(imeancav/imeanart)

其中，imeancav为空化区域或空化区域选定部分的平均能量，imeanart为伪影区域或伪影区域选定部分的平均能量，car的单位为db；

4.2)根据空化区域和伪影区域的空化伪影比或空化区域和伪影区域不同选定部分的空化伪影比的均值，评价不同ε下的全频域被动空化成像质量，ε为导向矢量不确定集参数，并选取当car或car均值最高时对应的ε作为最优的导向矢量不确定集参数，根据最优的导向矢量不确定集参数得到特征空间自适应波束合成的最优输出信号其中，ε的遍历范围为0.01～n，遍历步长为0.01～1；

4.3)选择不同的子频域，对步骤4.2)所得的特征空间自适应波束合成的最优输出信号在不同子频域分别进行叠加，则得到不同子频域的被动空化成像：

其中，fr表示给定的某一子频域，f1为该子频域的频率下限，f2为该子频域的频率上限。

上述步骤五，具体包括以下步骤：

5.1)对步骤四所得的不同子频域的被动空化成像结果进行标准化处理，然后进行叠加，得到被动空化频率复合成像结果：

其中，nfr为参与频率复合成像的子频域数目，为参与频率复合成像的经过标准化处理后的第i个子频域的被动空化成像结果。

一种基于特征空间自适应波束合成的频域被动空化成像系统，包括相位偏移和阵元变迹处理模块、自适应波束合成加权矢量计算模块、全频域的被动空化成像模块以及子频域的被动空化成像模块；所述相位偏移和阵元变迹处理模块用于执行上述步骤一；所述自适应波束合成加权矢量计算模块用于执行上述步骤二；所述全频域的被动空化成像模块用于执行上述步骤三；所述子频域的被动空化成像模块用于执行上述步骤四。

一种基于特征空间自适应波束合成的被动空化频率复合成像系统，包括上述相位偏移和阵元变迹处理模块、上述自适应波束合成加权矢量计算模块、上述全频域的被动空化成像模块、上述子频域的被动空化成像模块以及被动空化频率复合成像模块；所述被动空化频率复合成像模块用于执行上述步骤五。

本发明的有益效果体现在：

本发明针对传统频域被动空化成像方法在监控聚焦超声治疗方面存在的缺陷，通过对频域空化信号进行相位偏移和阵元变迹处理以及构造频域的协方差矩阵，构建频域自适应波束合成的代价函数，在此基础上求得自适应波束合成的加权矢量；然后将加权矢量投影到频域的协方差矩阵经特征空间分解所得的信号子空间上，并采用投影后的加权矢量进行波束合成，对成像伪影进行了有效的抑制，提高了不同子频域被动空化成像的空间分辨性能，对于聚焦超声治疗过程中空化的实时监控有着重要的意义。

本发明基于高分辨性能的子频域被动空化成像，通过对子频域被动空化成像结果进行标准化处理和叠加实现了被动空化频率复合成像，进一步提高了成像质量。

附图说明

图1是本发明实施例中时域空化信号转换到频域并对频域空化信号进行相位偏移和阵元变迹的流程图；

图2为某超声线阵换能器获取到的时域空化信号和经过相位偏移和阵元变迹处理后的频域空化信号；

图3为本发明实施例中自适应波束合成的加权矢量的计算流程图；

图4为自适应波束合成中协方差矩阵、拉格朗日乘数迭代值、导向矢量和加权矢量的结果；

图5为本发明实施例中基于特征空间自适应波束合成的全频域被动空化成像流程图；

图6为特征空间自适应波束合成中信号和噪声子空间的协方差矩阵、投影后的加权矢量和全频域被动空化成像的结果；

图7为本发明实施例中不同子频域的被动空化成像及被动空化频率复合成像的流程图；

图8为不同子频域的被动空化成像及被动空化频率复合成像结果。

具体实施方式

下面结合附图和实施例对本发明作进一步的详细说明。

本发明提供一种基于特征空间自适应波束合成的频域被动空化成像及频率复合成像方法，包括以下步骤：

步骤一：超声成像换能器被动接收时域空化信号，对时域空化信号进行傅里叶变换得到频域空化信号，并对频域空化信号作相位偏移处理，然后利用归一化余弦变迹函数得到经过相位偏移和阵元变迹处理后的频域空化信号。具体步骤如下(图1)：

(1.1)将超声成像设备设置为只被动接收信号的模式，通过任意波形发生器编写波形来同步触发聚焦超声换能器和超声成像设备，由超声成像设备的超声成像换能器被动接收聚焦超声换能器作用过程中焦域处散射的空化射频信号，将该空化射频信号称为时域空化信号，该时域空化信号由超声成像设备的并行通道数据采集及存储模块采集；

(1.2)假设超声成像换能器(例如，超声线阵换能器)有n个阵元，对第i个阵元采集到的时域空化信号进行傅里叶变换，得到频域空化信号pi(f,x,z)：

其中i＝1,2,...,n，t1是采集时域空化信号的开始时刻，t2是采集时域空化信号的停止时刻，j是虚数单位，pi(t,x,z)是第i个阵元采集到的时域空化信号，t表示时间；

(1.3)计算从成像位置(x,z)到第i个阵元(xi,0)的声传播时间，即第i个阵元的时间延时：

其中，c是超声波的声传播速度；

(1.4)根据步骤(1.3)所得的时间延时计算第i个阵元的相位偏移：

其中，τi为τi(x,z)的简称；

(1.5)根据步骤(1.2)所得频域空化信号和步骤(1.4)所得的相位偏移得到相位偏移处理后的频域空化信号

(1.6)根据步骤(1.2)所得频域空化信号计算归一化变迹函数：

其中，表示对于任意的频率f；

(1.7)根据成像位置(x,z)和第i个阵元的位置(xi,0)来计算余弦变迹函数：

(1.8)根据步骤(1.6)所得归一化变迹函数和步骤(1.7)所得余弦变迹函数计算归一化余弦变迹函数：

apodi(f,x,z)＝apod1i(f,x,z)·apod2i(f,x,z)

(1.9)利用步骤(1.8)所得归一化余弦变迹函数对步骤(1.5)所得的相位偏移处理后的频域空化信号进行阵元变迹处理：

(1.10)跳转至下一个阵元，重复步骤(1.2)～(1.9)，直到得到n个阵元的经过相位偏移和阵元变迹处理后的频域空化信号。

参见图2，其中(a)为超声成像换能器阵列示意图，阵元总数为128，孔径大小为38mm，(b)为第64个阵元采集到的时域空化信号，(c)为第64个阵元的经过相位偏移和阵元变迹处理后的频域空化信号。

步骤二：根据步骤一所得的所有阵元的经过相位偏移和阵元变迹处理后的频域空化信号来构造协方差矩阵，利用该协方差矩阵及导向矢量不确定集参数构建自适应波束合成的代价函数，对代价函数的约束条件引入拉格朗日乘数并获得导向矢量表达式，利用牛顿迭代法求得约束条件中的拉格朗日乘数，在此基础上计算导向矢量并对其作归一化处理，最终得到自适应波束合成的加权矢量。具体步骤如下(图3)：

(2.1)利用步骤一所得的n个阵元的经过相位偏移和阵元变迹处理后的频域空化信号来构建协方差矩阵：

其中，s(f,x,z)＝[s1(f,x,z)；s2(f,x,z)；...；sn(f,x,z)]为n个阵元的经过相位偏移和阵元变迹处理后的频域空化信号矩阵，[·]^h代表矩阵的共轭转置；

(2.2)利用步骤(2.1)所得协方差矩阵构建自适应波束合成的代价函数：

其中，a为要求解的导向矢量，为假定的导向矢量，||·||代表欧几里得范数，为约束条件，ε为描述导向矢量不确定集的参数，简称为不确定集参数，[·]-¹代表矩阵求逆；

(2.3)为了求解步骤(2.2)所述的代价函数，引入拉格朗日乘数并构造拉格朗日约束函数：

其中，λ为拉格朗日乘数；

(2.4)将步骤(2.3)所述的拉格朗日约束函数对导向矢量a求导并令求导所得的导函数为零，则得到导向矢量表达式：

其中，r为步骤(2.1)所得的协方差矩阵r(f,x,z)，以下简写为r，i为单位矩阵；

(2.5)根据矩阵求逆引理将步骤(2.4)所述的导向矢量表达式重写为：

(2.6)将步骤(2.5)所得的导向矢量表达式代入步骤(2.2)的约束条件中，则约束条件改写为：

(2.7)为了方便计算，对协方差矩阵r做特征空间分解r＝uλu^h，其中u为特征向量矩阵，λ为特征值矩阵，并令则步骤(2.6)所得的约束条件可转化为：

其中，bi为b的第i个元素，γi为特征值矩阵λ的第i个对角元素；

(2.8)采用牛顿迭代法来求解经步骤(2.7)转化后的约束条件中的拉格朗日乘数λ，将求解得到的λ代入步骤(2.5)所得的导向矢量表达式中，则可计算出来导向矢量，对计算出来的导向矢量做归一化处理：

步骤(2.8)所述的牛顿迭代法中，首先设置收敛条件和迭代初始值，当收敛条件满足时，迭代停止，得到拉格朗日乘数。牛顿迭代法的迭代关系式为：

其中k＝1,2,...，λk为拉格朗日乘数的第k次迭代值，λk+1为拉格朗日乘数的第k+1次迭代值，f(λk)为第k次迭代的函数，f'(λk)为f(λk)的导数；拉格朗日乘数迭代初始值λ1设置为其中γ1为步骤(2.7)所述特征值矩阵λ的第1个对角元素；

(2.9)根据步骤(2.8)所得的归一化导向矢量和步骤(2.1)所得的协方差矩阵计算自适应波束合成的加权矢量：

参见图4，其中(a)为根据128个阵元的经过相位偏移和阵元变迹处理后的频域空化信号得到的协方差矩阵；(b)拉格朗日乘数随着迭代次数增多的变化曲线，可看出第6次迭代时已开始收敛；(c)为归一化的导向矢量，(d)为自适应波束合成的加权矢量。

步骤三：对步骤二所得协方差矩阵进行特征空间分解，得到信号子空间，并将步骤二所得的自适应波束合成的加权矢量投影到信号子空间上，利用投影后的加权矢量对步骤一所得的所有阵元的经过相位偏移和阵元变迹处理后的频域空化信号进行加权求和，再去除掉直流分量，得到特征空间自适应波束合成的输出信号；对特征空间自适应波束合成的输出信号进行平方并在整个频域上进行叠加，得到全频域的被动空化成像结果。具体步骤如下(图5)：

(3.1)自适应波束合成算法的性能下降主要是由于噪声子空间的扰动引起的，其中信号子空间包含着有用的空化信息，而噪声子空间则包含影响成像质量的伪影信息。为了将伪影有效抑制，对步骤(2.1)所得的协方差矩阵r进行特征空间分解：

r＝uλu^h＝rs+ri

其中，u＝[u1,u2,...,un]为特征向量矩阵，λ＝diag[γ1,γ2,...,γn]为对角矩阵，对角元素为特征值，γ1≥γ2≥...≥γn，rs为信号子空间的协方差矩阵，ri为噪声子空间的协方差矩阵；

步骤(3.1)所得的信号子空间的协方差矩阵具体为：

其中，us＝[u1,u2,...,ul]为前l个较大特征值对应的特征向量构成的特征向量矩阵，λs＝diag[γ1,γ2,...,γl]为信号子空间的特征值矩阵，diag[·]表示矢量的对角化；

步骤(3.1)所得的噪声子空间的协方差矩阵具体为：

其中，ui＝[ul+1,ul+2,...,un]为余下的较小特征值对应的特征向量构成的特征向量矩阵，λi＝diag[γl+1,γl+2,...,γn]为噪声子空间的特征值矩阵；

参数l为特征值γ1,γ2,...,γn中大于最大特征值(γ1)δ倍的所有特征值的数目，δ一般选择为0.2～0.5；

(3.2)令则步骤(2.9)所得的自适应波束合成的加权矢量可表示为：

其中，为步骤(2.8)所得的归一化导向矢量；

由于理想情况下，信号子空间正交于噪声子空间，因此将自适应波束合成的加权矢量投影到信号子空间，依据步骤(3.5)可得投影后的加权矢量为：

(3.3)利用步骤一所得的每个阵元的经过相位偏移和阵元变迹处理后的频域空化信号和步骤(3.2)所得的投影后的加权矢量来计算每个阵元的加权频域空化信号：

其中，为步骤(3.2)所得投影后的加权矢量的第i个元素，si(f,x,z)为步骤一所得的第i个阵元的经过相位偏移和阵元变迹处理后的频域空化信号；

(3.4)对超声成像换能器n个阵元的加权频域空化信号进行叠加，则可得到：

(3.5)从步骤(3.4)所得信号q(f,x,z)中去除掉直流分量，则得到特征空间自适应波束合成的输出信号：

(3.6)对步骤(3.5)所得的特征空间自适应波束合成的输出信号进行平方并在整个频域上进行叠加，得到每个成像位置(x,z)的能量，即全频域的被动空化成像结果：

参见图6，其中(a)为信号子空间的协方差矩阵，(b)噪声子空间的协方差矩阵，(c)为自适应波束合成的加权矢量投影到信号子空间后得到的加权矢量，(d)为全频域被动空化成像的结果。

步骤四：构造空化伪影比指标，根据空化伪影比计算结果来选择最优的导向矢量不确定集参数，从而确定特征空间自适应波束合成的最优输出信号；在此基础上获得子频域的被动空化成像结果，对不同子频域的被动空化成像结果做标准化处理并进行叠加，得到被动空化频率复合成像结果。具体步骤如下(图7)：

(4.1)构造基于平均能量的空化伪影比指标：

car＝20log10(imeancav/imeanart)

其中imeancav为空化区域的平均能量，imeanart为伪影区域的平均能量，car的单位为db；

(4.2)设置不确定集参数ε的遍历范围0.01～n，遍历步长一般可设置为0.01～1，根据步骤一、二和三所述的方法，得到不同的不确定集参数ε下的全频域被动空化成像结果；

(4.3)依据步骤(4.1)所述的空化伪影比指标car，计算不同的不确定集参数ε下的全频域被动空化成像结果的car，以此来评价不同不确定集参数ε下的全频域被动空化成像质量，并选取当car最高时对应的ε作为最优的不确定集参数；

空化区域与伪影区域的边界根据全频域被动空化成像结果中能量最大值在横向和轴向下降到能量最大值一半时所对应的横向坐标和轴向坐标来确定，边界内为空化区域，边界外为伪影区域；

(4.4)在最优不确定集参数ε的基础上进行步骤一、步骤二和步骤三中步骤(3.1)～(3.5)所述的处理，得到特征空间自适应波束合成的最优输出信号

(4.5)在整个频域中选择不同的子频域，得到不同的频率成分，如次谐波(f0/2±f0/8)、谐波(nf0±f0/4,n＝2,3,...,10)、超谐波[(2n+1)f0/2±f0/8,n＝1,2,...,10]和宽带噪声[nf0/2+f0/8～(n+1)f0/2-f0/8,n＝1,2,...,20]等，其中f0为产生空化的聚集超声换能器的发射频率；

(4.6)对步骤(4.4)所述的特征空间自适应波束合成的最优输出信号在步骤(4.5)所述的不同子频域进行叠加，则可得到不同子频域的被动空化成像：

其中，f1为子频域的频率下限，f2为子频域的频率上限，fr表示一给定的子频域；

(4.7)对步骤(4.6)所得的不同子频域的被动空化成像结果进行如下标准化处理：

其中，表示对于任意的成像位置(x,z)，m为成像位置的总数目；

(4.8)对步骤(4.7)所得的标准化处理后的不同子频域的被动空化成像结果进行叠加，则可得到被动空化频率复合成像结果：

其中，nfr为参与频率复合成像的子频域数目；

参见图8，其中(a)、(b)、(c)、(d)为根据传统频域被动空化成像方法所得到的不同子频域(次谐波、谐波、超谐波、宽带噪声)的成像结果；(e)、(f)、(g)、(h)为根据本发明提出的基于特征空间自适应波束合成的频域被动空化成像方法所得到的不同子频域(次谐波、谐波、超谐波、宽带噪声)的成像结果；比较(a)～(d)和(e)～(h)可知，本发明提出的方法可有效抑制传统频域被动空化成像的伪影，因而提高了空化的空间定位精度。另外，图8(i)为被动空化频率复合成像结果，经过频率(上述次谐波、谐波、超谐波、宽带噪声四种频率成分)复合之后，可以看出，相对于各个子频域的成像结果，成像质量进一步提高。

对于图8(i)的结果，本发明提出以下分析：聚焦超声作用过程中，前一次空化散射的多个频率成分之间互相干扰，进而影响后续的空化行为，而传统频域被动空化成像只针对某一个子频域进行成像，无法实现空化活动在多个子频域的有效定征，同时这也是造成成像质量下降的一个原因。

本发明具有以下优点：

(1)本发明将采集到的空化散射信号在频域进行波束合成，可直接在频域选择频率成分，有效地提高了计算效率；波束合成过程中采用频域的相位偏移而非时域的时间延时，因此降低了由于实际超声成像设备采样率受限产生延时误差而造成的成像伪影；

(2)本发明提出的方法是在频域进行自适应波束合成，解决了传统频域被动空化成像方法不依赖于数据本身特性的问题，并将加权矢量在信号子空间上进行投影，在对旁瓣伪影进行大幅度抑制的同时保留了空化信息，从而实现了高分辨率的频域被动空化成像；

(3)本发明利用空化伪影比这一指标选择出了最优的不确定集参数，在此基础上得到了次谐波、谐波、超谐波和宽带噪声等不同频率成分对应的子频域的高分辨率被动空化成像，从而为聚焦超声作用过程中稳态空化和惯性空化的实时检测和定量分析奠定了基础；

(4)本发明将不同子频域的被动空化成像进行复合，得到了高分辨率的频率复合成像，成像质量进一步提高，且提出的频率复合成像方法可对不同子频域宽度和不同子频域数目下的多个频率成分进行综合定征；

(5)本发明提出的频域被动空化成像及频率复合成像方法可对不同聚焦超声换能器、不同超声成像换能器、不同聚焦超声发射参数和不同实验设置下产生的空化进行成像，同时适用于不同治疗部位和不同的聚焦超声治疗应用，此外也可用于超声空化瞬态物理过程的分析和研究。