CN105266847B - 基于压缩感知自适应波束合成的脉冲逆转谐波平面波快速造影成像方法 - Google Patents

Abstract

本发明提供了一种基于压缩感知自适应波束合成的脉冲逆转谐波平面波快速造影成像方法，使得高分辨率、造影组织比的平面波造影图像实时成像成为可能。本方法仅对超声阵列探头阵元接收到的造影回波的2阶谐波分量，进行压缩随机采样处理，大幅度剔除冗余数据，减少波束合成运算量；进一步采用压缩感知理论对自适应波束合成算法性能优化，在保证图像分辨率不降低的前提下，实现快速波束合成运算；联合脉冲逆转谐波造影成像方式，换能器交替发射一对相位相反的平面波，高灵敏地检测微泡的非线性谐波分量、大幅度改善图像造影组织比的同时，使得快速合成一幅平面波造影图像只需2次平面波发射，避免了微泡多次辐射造成的破坏，准确地扑捉微泡的瞬态分布。

Description

基于压缩感知自适应波束合成的脉冲逆转谐波平面波快速造影成像方法

技术领域

本发明属于超声造影成像技术领域，具体涉及一种基于压缩感知自适应波束合成的脉冲逆转谐波(PIH-CSBF)平面波快速造影成像方法。

背景技术

超声造影成像是利用包膜微泡在声场激励下的强非线性响应，以增强超声对血流的检测能力，有效地从周围组织中鉴别微血管分布的成像方式，可得到高分辨性能的微血管造影图像。它在实时描绘造影灌注的动态分布过程，同时可从中提取与生理病理相关的血流动力学参量，使得心血管疾病与恶性肿瘤新生血管的早期诊断和血流灌注时相分析成为可能，并在临床广泛应用。传统的超声造影成像技术采取的聚焦波逐条扫描成像方式，以获得对比度、分辨率较佳的造影图像。但该扫描成像方式存在扫描线间延时导致的伪实时成像、成像帧率难以提高以及在焦点区域对造影微泡破坏率较高的问题。虽然，可以通过降低成像发射声压以减少造影微泡的破坏率、减少扫描线次数并降低成像深度以提高图像帧率；尽管如此，其改善依然是局限性的，未能从技术上突破，且以损失细节分辨率、信噪比为代价。

因此，为克服上述问题，近年兴起了以平面波为发射方式的造影成像技术。但因平面波发射声压较低缺乏焦点的问题，其成像造影组织比(Contrast-to-tissueratio,CTR)和分辨率受限，为此发展了多角度发射的相干复合平面波成像方法、脉冲逆转微泡小波平面波造影成像等方法，以提高造影图像的CTR和分辨率，并在血流仿体、兔肾脏造影团注成像中应用。遗憾的是，上述平面波造影成像方法均为离线处理方法，它们均采用延时叠加(DelayandSum,DAS)波束合成法。但DAS没有考虑通道数据的差异，加权系数固定，其主瓣宽度过宽，旁瓣高度过高，所得图像质量较差。为此就有上述脉冲逆转微泡小波变换与多角度发射相干复合平面造影成像技术的出现，以改善DAS所得造影图像的图像质量。近期有文献显示已有学者为了提高横向分辨率和对比度，通过自适应波束合成方法以改善DAS的技术缺陷。

尽管自适应波束合成相比DAS可得到高分辨率的造影图像，但其较复杂的算法引入更多的计算量。在实际应用中，为获得高CTR与分辨率造影图像而采用的多角度发射相干复合与脉冲逆转微泡小波变换造影成像方法联合自适应波束合成技术均极为复杂，在面对成像系统高采样率、平面波快速发射采集的极为庞大的数据量时，使得超快平面波造影成像难以实时实现。因此，如何在保证造影图像CTR与分辨率的前提下，实现平面波造影实时成像是本领域一大挑战。

发明内容

本发明的目的在于提供一种基于压缩感知自适应波束合成的脉冲逆转谐波平面波快速造影成像方法。

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

步骤一超声阵列换能器按脉冲逆转模式发射平面波，其交替发射相位分别为0°相与180°相，计算机采集超声阵列换能器阵元的0°相与180°相射频回波信号X₀(t)与X₁₈₀(t)；

步骤二将0°相与180°相射频回波信号加和运算所得射频回波信号X(t)，依据第i条待合成扫描线上的第d个目标点空间位置与动态孔径N，对N个有效阵元通道数据延时后得到加和延时时域信号X_id(t)；

步骤三对加和延时时域信号X_id(t)通过傅立叶变换转换为频域信号X_id(ω)，并选取造影回波的谐波分量内的频点X_id(ω_k)，进行按比率的压缩随机采样抽取，得到压缩采样后的造影谐波频域信号

步骤四通过构造所选频点的稳健自相关矩阵，计算该频点所对应的最小方差自适应波束合成的最优加权系数w_id(ω_k)，从而得到该频点的最优波束合成的频域输出

步骤五按照步骤四遍历所有抽取的造影谐波频点并由多点欠定系统聚焦求解(RM-FOCUSS)算法重建谐波分量的频域信息，然后重建第i条扫描线的频域信息y_i(ω)；

步骤六按扫描线波束合成次序，重复步骤二～步骤五得到所有扫描线的频域信息Y(ω)，并将加和射频数据的频域信息Y(ω)通过逆傅立叶变换变换到时域，得到Y(t)，对Y(t)进行包络检波、对数压缩、坐标变换等图像后处理得到高质量的平面波造影图像。

具体步骤如下：

一、本发明所采用的技术方案步骤一，具体步骤如下：

1)超声阵列换能器按脉冲逆转模式交替发射中心频率为f₀的一对相位相反的平面波，其相位分别为0°相与180°相，探头工作频率f₀范围：2.0MHz～5.0MHz；

2)计算机采集超声阵列换能器阵元的0°相与180°相射频回波信号X₀(t)与X₁₈₀(t)；

3)其成像特征为，合成一副平面波造影图像，仅需发射2次，即：一对相位相反的平面波。

二、本发明所采用的技术方案步骤二，具体步骤如下：

1)将0°相与180°相的射频回波信号X₀(t)与X₁₈₀(t)在整个时域加和运算得到射频回波信号X(t)，X(t)＝X₀(t)+X₁₈₀(t)；

2)选定第i条(1≤i≤128～256)待合成扫描线上目标点并根据目标点与中心阵元垂直距离计算波束合成的有效阵元数大小N，并选定阵元；

3)根据第i条待合成扫描线上第d个目标点与N个阵元相对位置，对N个阵元通道接收信号计算每一个通道的信号延时时间τ_id，τ_id＝[τ₁τ₂…τ_N]^T，其中，为第j个阵元的延时时间，D_i为第j个阵元到目标点距离，D_c为目标点中心阵元到目标点距离，c为组织中声速，j＝1,2,3,...,N；

4)根据第d个待合成目标点与N个阵元相对位置，对N个阵元对应射频回波信号X(t)按τ_id依次延时，得到延时后列向量信号X_id(t)＝X(t-τ_id)，X_id(t)＝[x₁(t)x₂(t)…x_N(t)]^T。

三、本发明所采用的技术方案步骤三，具体步骤如下：

1)对延时后时域信号X_id(t)进行傅里叶变换，得到频域信号X_id(ω)，X_id(ω)＝[x₁(ω)x₂(ω)…x_N(ω)]^T；

2)对造影回波的频域信号X_id(ω)进行2f₀谐波滤波，对谐波分量ω_k内的频点构成维数为k×N的谐波信号X_id(ω_k)，其有效带宽为1.4f₀～2.6f₀，k为有效带宽内的频点数；

3)基于该频点数k，构造k阶的主对角线元素为1、其余为0的压缩随机抽取矩阵T，抽取比例为30％～40％；

4)对X_id(ω_k)按行进行随机抽取：抽取比例为30％～40％，得到压缩采样后的造影谐波频域信号

四、本发明所采用的技术方案步骤四，具体步骤如下：

1)通过构造所选频点的稳健自相关矩阵，计算该频点所对应的最小方差自适应波束合成的最优加权系数w_id(ω_k)；

2)对于中的频点，把有效阵元数为N的阵列分为L×k、总数为N-L+1的子阵

3)将所有子阵与其对应的最优加权系数w_id(ω_k)依次点乘求和并平均，得到该频点的最优波束合成的频域输出

${\hat{y}}_{id}\left({\mathrm{\ω}}_k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{w}_{id}\left({\mathrm{\ω}}_k\right){\hat{X}}_{id}^l\left({\mathrm{\ω}}_k\right)---\left(1\right)$

其中，为的第l个子阵。

五、本发明所采用的技术方案中，基于最小方差自适应波束合成的最优加权系数w_id(ω_k)的计算方法，具体步骤如下：

1)基于互相关原理，构造所选频点的稳健自相关矩阵

${\tilde{R}}_{id}\left({\mathrm{\ω}}_k\right)=R_{id}\left({\mathrm{\ω}}_k\right)-{\mathrm{AA}}^H{R}_{id}\left({\mathrm{\ω}}_k\right){\mathrm{AA}}^H/{\left(A^{H}A\right)}^2---\left(2\right)$

其中，[]^H表示矩阵的共轭转置，R_id(ω_k)为的自相关矩阵：

$R_{id}\left({\mathrm{\ω}}_k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{\hat{X}}_{id}^l\left({\mathrm{\ω}}_k\right){\hat{X}}_{id}^l{\left({\mathrm{\ω}}_k\right)}^H---\left(3\right)$

其中，A为长度为L的1维单位向量；AA^HR_id(ω_k)/(A^HA)²为造影回波所选频点的主信号能量项；

2)利用自相关矩阵和单位向量A计算最优加权系数w_id(ω_k)：

$w_{id}\left({\mathrm{\ω}}_k\right)={{\tilde{R}}_{id}}^{-1}\left({\mathrm{\ω}}_k\right)A/A^H{{\tilde{R}}_{id}}^{-1}\left({\mathrm{\ω}}_k\right)A---\left(4\right)$

其中，[]^-1表示矩阵的逆，为避免求逆运算出现不稳定性，将奇异值分解出现的负特征值替换为分解所得的最小正值特征解；同时，从R_id(ω_k)中减掉的主信号能量项，以增强的稳健性。

六、本发明所采用的技术方案步骤五，具体步骤如下：

1)重复步骤四，计算第i条待合成扫描线上所有抽取频点的最优波束合成的频域信号输出

2)构建k×k的傅里叶基F作为稀疏基，利用下式：

${{\hat{y}}_i}^H\left({\mathrm{\ω}}_k\right)={{\mathrm{TFb}}_i}^H\left({\mathrm{\ω}}_k\right)---\left(5\right)$

因抽取矩阵T与傅里叶基F相乘矩阵不满秩，故采用RM-FOCUSS算法和构建的抽取矩阵T与傅里叶基F，从阵元矩阵中估计出稀疏投影系数向量b_i(ω_k)；

3)由y_i(ω_k)＝Fb_i(ω_k)计算阵元输出所有有效带宽内的最优波束合成的频域信号输出y_i(ω_k)；

4)对y_i(ω_k)中造影谐波有效带宽外的频点替换为0点，得到第i条扫描线上整个频域的最优波束合成的频域信号输出y_i(ω)。

七、本发明所采用的技术方案步骤六，具体步骤如下：

1)按扫描线波束合成次序，得到所有扫描线在整个频域的最优波束合成的频域信号输出Y(ω)，Y(ω)＝[y₁(ω)y₂(ω)…y_line(ω)]^T，其中，line为波束合成扫描线数，其取值一般为128～256；

2)将0°相与180°相的射频回波信号脉冲逆转加和后并2f₀谐波滤波后的最优波束合成的频域信号Y(ω)进行逆傅里叶变换，得到整个成像区域射频信号Y(t)；

3)对造影信号Y(t)进行包络检波、对数压缩、坐标变换等图像后处理得到高质量的平面波造影图像。

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

本发明采用压缩感知自适应波束合成(CSBF)算法，其各阵元接收数据的有用频点信息都集中在造影回波信号的谐波分量的有效带宽内，因而对这部分有效带宽内的数据进行压缩感知处理，而非对整个频域进行压缩感知处理，从而在频域自适应波束合成的基础上，进一步降低频域采样点数，在保证图像分辨率不降低的前提下，极大地减少波束合成时所需的计算量。同时，联合脉冲逆转谐波(PIH)造影成像方式，超声阵列换能器交替发射一对相位相反的平面波，在有效抑制组织线性回波、高灵敏地检测微泡非线性谐波分量、大幅度改善图像CTR的同时，使得快速合成一幅平面波造影图像只需2次平面波发射，避免了微泡多次辐射造成的破坏，准确地扑捉了微泡的瞬态分布。

本发明的优点如下：

1、通过在频域仅对微泡谐波分量进行压缩随机采样，大幅度减少波束合成所需计算量，由此，可使平面波造影实时成像成为可能。

2、通过压缩感知理论对自适应波束合成算法性能的优化，完美重建平面波造影回波的时域数据，可得到高横向分辨率的造影图像。

3、联合对数字超声影像平台来说实现较为简单的PIH影成像方法，进一步提高平面波造影图像的CTR。

4、发射2次平面波便可获得一幅高CTR与分辨率平面波造影图像，避免了平面波多次发射造成的微泡破坏及运动伪迹的干扰，准确扑捉微泡的瞬态分布。

附图说明

图1是本发明的算法流程图。

图2是基于压缩感知自适应波束合成的脉冲逆转谐波(PIH-CSBF)狭窄血管流动仿体平面波造影图像，从左至右分别为第100帧、3900帧、4000帧平面波造影图像。

图3是第7000帧在不同波束合成方法下的PIH狭窄血管流动仿体平面波造影图像，(a)为DAS波束合成方法所得图像；(b)为传统自适应波束合成(Regular)方法所得图像；(c)为压缩感知自适应波束合成(CSBF)方法所得图像。

图4是第7000帧不同波束合成方法下的PIH狭窄血管流动仿体平面波造影图像横向(a)与纵向(b)分辨率的对比结果。

图5是第7000帧不同波束合成方法下的PIH狭窄血管流动仿体平面波造影图像的CTR与运行时间的比较；(a)DAS所得结果，(b)传统自适应波束合成(Regular)方法所得结果，(c)CSBF所得结果。

具体实施方式

下面结合附图和实施例(一狭窄血管流动仿体平面波造影成像实例)对本发明作进一步说明。

本发明提出一种基于压缩感知自适应波束合成的脉冲逆转谐波(PIH-CSBF)平面波快速造影成像方法，以克服现有技术在高帧率、超快速成像，尤其是以瞬态信息为特征的造影成像研究中所存在的问题和局限性。

1)超声阵列换能器按脉冲逆转模式交替发射一对中心频率为f₀的相位相反的平面波，而计算机采集超声线阵换能器接收0°相与180°相射频回波信号，仅用2次平面波发射得到一副平面波造影图像，以减少对微泡的辐照，减低微泡的破坏率，同时最大程度避免运动等杂波信号的干扰，准确扑捉微泡的瞬态分布。

2)将0°相与180°相射频信号加和处理，并从时域信号转换为频域信号，按其频域信息分布，选取频域上造影微泡谐波分量有效带宽内的频点，并按比率的压缩随机采样(抽取)；这个过程，通过谐波滤波、压缩采样，大幅度剔除冗余数据、减少自适应波束合成计算量(减少50％～70％)，为平面波造影快速成像奠定数据样本基础。

3)对所抽取的频点构造稳健自相关矩阵，并计算该频点所对应的最小方差自适应波束合成的最优加权系数，从而得到该频点的最优输出；此过程中，考虑阵元信号的差异性，通过压缩感知方法提高自适应波束合成方法性能，在提高运算速度的同时获得高质量的CSBF频域数据。

4)按照步骤2)～3)遍历第一条扫描线上所有抽取的频点，并由压缩感知理论重建谐波分量的该相位频域信息，然后通过填0重建整个频域信息；将第一条扫描线上整个频域信息变换到时域射频信号。

5)按扫描线波束合成次序，重复上述2)～4)步骤，得到所有扫描线相应谐波滤波后的CSBF射频数据，并进行包络检波、对数压缩、坐标变换等射频后处理得到平面波造影图像。

参照图1所示流程图，所述平面波快速造影成像方法，具体步骤如下：

(1)超声阵列换能器按脉冲逆转模式交替发射中心频率为f₀的一对相位相反的平面波，其相位分别为0°相与180°相，探头工作频率f₀范围：2.0MHz～5.0MHz；

(2)计算机采集超声阵列换能器阵元的0°相与180°相射频回波信号X₀(t)与X₁₈₀(t)；

(3)将0°相与180°相的射频回波信号X₀(t)与X₁₈₀(t)在整个时域加和运算得到射频回波信号X(t)：

X(t)＝X₀(t)+X₁₈₀(t)(1)

(4)对造影回波数据X(t)，计算第i条(1≤i≤128～256)待合成扫描线上目标点并根据目标点与中心阵元垂直距离计算波束合成的有效阵元数(即有效孔径)大小，记为N；

(5)根据第i条待合成扫描线上第d个目标点与N个阵元相对位置，对N个阵元通道接收信号计算每一个通道的信号延时时间τ_id；

τ_id＝[τ₁τ₂…τ_N]^T(2)

其中，为第j个阵元的延时时间，D_i为第j个阵元到目标点距离，D_c为为目标点中心阵元到目标点距离，c为组织中声速；

(6)根据第d个待合成目标点与N个阵元相对位置，对N个阵元对应射频回波信号X(t)按τ_id依次延时，得到延时后列向量信号X_id(t)＝X(t-τ_id)：

X_id(t)＝[x₁(t)x₂(t)…x_N(t)]^T(3)

(7)对X_id(t)进行傅里叶变换，得到X_id(ω)：

X_id(ω)＝[x₁(ω)x₂(ω)…x_N(ω)]^T(4)

(8)根据超声发射频率f₀，选取2f₀有效带宽ω_k(1.4f₀～2.6f₀有效频带)的频点构成维数为k×N的信号X_id(ω_k)，k为有效带宽内频点数；

(9)在有效带宽ω_k内，构造维数为k×k、主对角线元素为1、其余为0的抽取矩阵T；

(10)对谐波滤波后的频域信号X_id(ω_k)按行进行随机抽取： ${\hat{X}}_{id}\left({\mathrm{\ω}}_k\right)=T\×X_{id}\left({\mathrm{\ω}}_k\right),$ 抽取比例为30％～40％；

(11)对于抽取的频点把有效阵元数为N的阵列分为长度为L的子阵，子阵的总数为N-L+1，采用空间平滑技术计算的空间平滑相关矩阵R_id(ω_k)：

$R_{id}\left({\mathrm{\ω}}_k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{\hat{X}}_{id}^l\left({\mathrm{\ω}}_k\right){\hat{X}}_{id}^l{\left({\mathrm{\ω}}_k\right)}^H---\left(5\right)$

其中，为的第l个子阵；

(12)基于互相关原理，构造所选频点的稳健自相关矩阵

${\tilde{R}}_{id}\left({\mathrm{\ω}}_k\right)=R_{id}\left({\mathrm{\ω}}_k\right)-{\mathrm{AA}}^H{R}_{id}\left({\mathrm{\ω}}_k\right){\mathrm{AA}}^H/{\left(A^{H}A\right)}^2---\left(6\right)$

其中，[]^Η表示矩阵的共轭转置，A为长度为L的1维单位向量；AA^HR_id(ω_k)/(A^HA)²为造影回波所选频点的主信号能量项；

(13)利用自相关矩阵和单位向量A计算最优加权系数：

$w_{id}\left({\mathrm{\ω}}_k\right)={{\tilde{R}}_{id}}^{-1}\left({\mathrm{\ω}}_k\right)A/A^H{{\tilde{R}}_{id}}^{-1}\left({\mathrm{\ω}}_k\right)A---\left(7\right)$

其中，[]^-1表示矩阵的逆，为避免求逆运算出现不稳定性，将奇异值分解出现的负特征值替换为分解所得的最小正值特征解；同时，从R_id(ω_k)中减掉的主信号能量项，以增强的稳健性。

(14)通过下式计算该频点的最优波束合成的频域输出

${\hat{y}}_{id}\left({\mathrm{\ω}}_k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{w}_{id}\left({\mathrm{\ω}}_k\right){\hat{X}}_{id}^l\left({\mathrm{\ω}}_k\right)---\left(8\right)$

(15)重复步骤(11)至步骤(14)，计算第i条扫描线上所有抽取频点的最优波束合成频域输出

(16)构建k×k的傅里叶基F作为稀疏基，利用下式：

${{\hat{y}}_i}^H\left({\mathrm{\ω}}_k\right)={{\mathrm{TFb}}_i}^H\left({\mathrm{\ω}}_k\right)---\left(9\right)$

因抽取矩阵T与傅里叶基F相乘矩阵不满秩，故采用RM-FOCUSS算法从阵元矩阵中估计出稀疏投影系数向量b_i(ω_k)；

(17)由y_i(ω_k)＝Fb_i(ω_k)计算阵元输出所有有效带宽内的最优波束合成的频域信号输出y_i(ω_k)；

(18)对有效带宽外频点以0填充，得到第i条扫描线上整个频域的最优波束合成频域输出y_i(ω)；

(19)按扫描线波束合成次序，重复步骤(4)至步骤(18)，得到所有扫描线在整个频域的最优波束合成的频域信号输出Y(ω)，Y(ω)＝[y₁(ω)y₂(ω)…y_line(ω)]^T，其中，line为波束合成扫描线数，其取值一般为128～256；

(20)对Y(ω)进行逆傅里叶变换，得到整个造影微泡成像区域的射频信号Y(t)；

(21)对造影信号Y(t)进行包络检波、对数压缩、坐标变换等图像后处理得到高质量的平面波造影图像。

实施平台为Sonix-Touch，选用128阵元线阵换能器为成像探头，按脉冲逆转模式，发射中心频率为3MHz的平面波，其负峰值声压约为178kPa，并在外触发源的控制下，按触发频率200Hz工作。参见图2，实施对象一狭窄血管流动仿体，其正常直径为5mm，按0.2mL的剂量对浓度为2.0×10⁸MB/mL团注。阵元数据被Sonix-DAQ采集约40s的数据。

图3结果显示，CSBF与Regular方法所得PIH平面波造影图像质量相当。

图4结果显示，相比DAS与Regular方法，CSBF在横向与轴向分辨率均取得最佳结果。

图5结果显示，相比基于DAS的PIH成像方法(PIH-DAS)与基于Regular的PIH成像方法(PIH-Regular)，PIH-CSBF所得造影图像的CTR最高可提高约2dB。而在取得与PIH-Regular相当图像质量(CTR)的同时，其分辨率更好，且运算时间仅为PIH-Regular的22.5％。同时，在前期实验中，PIH-CSBF相比基于DAS波束合成的B模式成像方法，CTR大幅度提高。

本发明通过压缩感知自适应波束合成方法(CSBF)大幅度剔除冗余数据，减少波束合成计算量，并实现快速自适应波束合成；同时，联合脉冲逆转谐波(PIH)造影成像技术，使得微泡非线性谐波分量更为灵敏地检测出来，得到高CTR与分辨率的造影图像，也使得波束合成数据仅为压缩采样的阵元接收的造影回波的谐波频带，再一次减少波束合成数据量。其中，方法基于近年兴起的稀疏采样理论，在远低于奈奎斯特采样率的情况下，通过非线性重建算法高概率的完美重建信号，可有效改善波束合成算法的计算性能；其压缩采样比例仅为有效频带内的回波信号的30％～40％，并通过最优稀疏基与阵元加权系数计算，快速频域加和、时域重建造影回波信号。因此，本发明一种基于压缩感知自适应波束合成的脉冲逆转谐波(PIH-CSBF)平面波快速造影成像方法的提出，使得高CTR与分辨率的平面波造影成像的实时实现成为可能。

Claims (7)

1.一种基于压缩感知自适应波束合成的脉冲逆转谐波平面波快速造影成像方法，其特征在于，包括以下步骤：

1)超声阵列换能器交替发射相位分别为0°相与180°相的平面波，计算机采集超声阵列换能器阵元的0°相与180°相射频造影回波信号X₀(t)与X₁₈₀(t)；

2)在0°相与180°相射频造影回波信号加和运算得到射频回波信号X(t)之后，依据第i条待合成扫描线上的第d个目标点空间位置与对应的动态孔径N，对N个有效阵元通道数据延时后得到加和延时时域信号X_id(t)；

3)对加和延时时域信号X_id(t)通过傅立叶变换转换为频域信号X_id(ω)，并选取谐波分量内的频点X_id(ω_k)，进行按比率的压缩随机采样抽取，得到压缩采样后的造影谐波频域频点

4)通过构造所选频点的稳健自相关矩阵，计算该频点所对应的最小方差自适应波束合成的最优加权系数w_id(ω_k)，从而得到该频点的最优波束合成的频域输出

5)按照步骤4)遍历第i条待合成扫描线上所有抽取的造影谐波频点并由多点欠定系统聚焦求解算法重建谐波分量的频域信息，然后重建第i条待合成扫描线的频域信息y_i(ω)；

6)按扫描线波束合成次序，重复步骤2)～步骤5)得到所有待合成扫描线的频域信息Y(ω)，将Y(ω)通过逆傅立叶变换变换到时域，得到Y(t)，对Y(t)进行图像后处理得到平面波造影图像。

2.根据权利要求1所述的方法，其特征在于：步骤1)具体包括以下步骤：

1.1)超声阵列换能器按脉冲逆转模式交替发射中心频率为f₀的一对相位相反的平面波，其相位分别为0°相与180°相，f₀范围为2.0MHz～5.0MHz；

1.2)计算机采集超声阵列换能器阵元的0°相与180°相的射频造影回波信号X₀(t)与X₁₈₀(t)。

3.根据权利要求1所述的方法，其特征在于：步骤2)具体包括以下步骤：

2.1)将0°相与180°相的射频造影回波信号X₀(t)与X₁₈₀(t)在整个时域加和运算得到射频回波信号X(t)；

X(t)＝X₀(t)+X₁₈₀(t)

2.2)选定第i条待合成扫描线上目标点并根据目标点与中心阵元垂直距离计算波束合成的有效阵元数大小N，并选定有效阵元，待合成扫描线数量为128～256；

2.3)根据第i条待合成扫描线上第d个目标点与N个有效阵元相对位置，对N个阵元通道接收信号计算每一个通道的信号延时时间τ_id；

τ_id＝[τ₁τ₂…τ_N]^T

其中，为第j个有效阵元的延时时间，D_j为第j个有效阵元到第d个目标点的距离，D_c为目标点中心阵元到第d个目标点的距离，c为组织中声速，j＝1,2,3,...,N；

2.4)对N个有效阵元对应射频回波信号按τ_id依次延时，得到加和延时时域信号X_id(t)＝X(t-τ_id)

X_id(t)＝[x₁(t)x₂(t)…x_N(t)]^T

其中，x_j(t)表示有效阵元中第j个射频回波信号，j＝1,2,3,...,N。

4.根据权利要求1所述的方法，其特征在于：步骤3)具体包括以下步骤：

3.1)对加和延时时域信号X_id(t)进行傅里叶变换，得到频域信号X_id(ω)；

X_id(ω)＝[x₁(ω)x₂(ω)…x_N(ω)]^T

其中，x_j(ω)是对应于x_j(t)的傅立叶变换后的频域信号，j＝1,2,3,...,N；

3.2)对频域信号X_id(ω)进行2f₀谐波滤波，将谐波分量ω_k内的频点构成维数为k×N的谐波信号X_id(ω_k)，谐波分量有效带宽为1.4f₀～2.6f₀，k为有效带宽内的频点数；

3.3)基于该频点数k，构造k阶的主对角线元素为1，其余为0的压缩随机抽取矩阵T，抽取比例为30％～40％；

3.4)对X_id(ω_k)按行进行随机抽取：得到压缩采样后的造影谐波频域频点

5.根据权利要求1所述的方法，其特征在于：步骤4)具体包括以下步骤：

4.1)通过构造所选频点的稳健自相关矩阵，计算该频点所对应的最小方差自适应波束合成的最优加权系数w_id(ω_k)；

4.2)对于中的频点，把有效阵元数为N的阵列分为L×k、总数为N-L+1的子阵

4.3)将所有子阵与其对应的最优加权系数w_id(ω_k)依次点乘求和并平均，得到中频点的最优波束合成的频域输出

${\hat{y}}_{id}\left({\mathrm{\ω}}_k\right)=\frac{1}{N-L+1}\underset{l=1}{\overset{N-L+1}{\Σ}}{w}_{id}\left({\mathrm{\ω}}_k\right){\hat{X}}_{id}^l\left({\mathrm{\ω}}_k\right)$

其中，为的第l个子阵。

6.根据权利要求1所述的方法，其特征在于：步骤5)具体包括以下步骤：

5.1)计算第i条待合成扫描线上所有抽取频点的最优波束合成的频域信号输出

5.2)构建k×k的傅里叶基F作为稀疏基，k为有效带宽内的频点数，利用下式：

${{\hat{y}}_i}^H\left({\mathrm{\ω}}_k\right)={{\mathrm{TFb}}_i}^H\left({\mathrm{\ω}}_k\right)$

采用多点欠定系统聚焦求解算法和构建的抽取矩阵T与傅里叶基F，从中估计出稀疏投影系数向量b_i(ω_k)；

5.3)由y_i(ω_k)＝Fb_i(ω_k)计算阵元输出所有有效带宽内的最优波束合成的频域信号输出y_i(ω_k)；

5.4)对所有谐波分量ω_k外的频点替换为0点，得到第i条扫描线上整个频域的最优波束合成的频域信息y_i(ω)。

7.根据权利要求1所述的方法，其特征在于：步骤6)具体包括以下步骤：

6.1)按扫描线波束合成次序，得到所有扫描线在整个频域的最优波束合成的频域信息Y(ω)

Y(ω)＝[y₁(ω)y₂(ω)…y_line(ω)]^T

其中，line为波束合成扫描线数；

6.2)将Y(ω)进行逆傅里叶变换，得到整个成像区域射频信号Y(t)；

6.3)对Y(t)进行包络检波、对数压缩以及坐标变换后得到平面波造影图像。