技术领域
本发明属于超声造影成像技术领域，具体涉及一种基于压缩感知自适应波束 合成的脉冲逆转谐波(PIH-CSBF)平面波快速造影成像方法。
背景技术
超声造影成像是利用包膜微泡在声场激励下的强非线性响应，以增强超声 对血流的检测能力，有效地从周围组织中鉴别微血管分布的成像方式，可得到高 分辨性能的微血管造影图像。它在实时描绘造影灌注的动态分布过程，同时可从 中提取与生理病理相关的血流动力学参量，使得心血管疾病与恶性肿瘤新生血管 的早期诊断和血流灌注时相分析成为可能，并在临床广泛应用。传统的超声造影 成像技术采取的聚焦波逐条扫描成像方式，以获得对比度、分辨率较佳的造影图 像。但该扫描成像方式存在扫描线间延时导致的伪实时成像、成像帧率难以提高 以及在焦点区域对造影微泡破坏率较高的问题。虽然，可以通过降低成像发射声 压以减少造影微泡的破坏率、减少扫描线次数并降低成像深度以提高图像帧率； 尽管如此，其改善依然是局限性的，未能从技术上突破，且以损失细节分辨率、 信噪比为代价。
因此，为克服上述问题，近年兴起了以平面波为发射方式的造影成像技术。 但因平面波发射声压较低缺乏焦点的问题，其成像造影组织比(Contrast-to-tissue ratio,CTR)和分辨率受限，为此发展了多角度发射的相干复合平面波成像方法、 脉冲逆转微泡小波平面波造影成像等方法，以提高造影图像的CTR和分辨率， 并在血流仿体、兔肾脏造影团注成像中应用。遗憾的是，上述平面波造影成像方 法均为离线处理方法，它们均采用延时叠加(DelayandSum,DAS)波束合成法。 但DAS没有考虑通道数据的差异，加权系数固定，其主瓣宽度过宽，旁瓣高度 过高，所得图像质量较差。为此就有上述脉冲逆转微泡小波变换与多角度发射相 干复合平面造影成像技术的出现，以改善DAS所得造影图像的图像质量。近期 有文献显示已有学者为了提高横向分辨率和对比度，通过自适应波束合成方法以 改善DAS的技术缺陷。
尽管自适应波束合成相比DAS可得到高分辨率的造影图像，但其较复杂的 算法引入更多的计算量。在实际应用中，为获得高CTR与分辨率造影图像而采 用的多角度发射相干复合与脉冲逆转微泡小波变换造影成像方法联合自适应波 束合成技术均极为复杂，在面对成像系统高采样率、平面波快速发射采集的极为 庞大的数据量时，使得超快平面波造影成像难以实时实现。因此，如何在保证造 影图像CTR与分辨率的前提下，实现平面波造影实时成像是本领域一大挑战。
发明内容
本发明的目的在于提供一种基于压缩感知自适应波束合成的脉冲逆转谐波 平面波快速造影成像方法。
为了实现上述目的，本发明采用了以下技术方案：
步骤一超声阵列换能器按脉冲逆转模式发射平面波，其交替发射相位分别 为0°相与180°相，计算机采集超声阵列换能器阵元的0°相与180°相射频回波信 号X0(t)与X180(t)；
步骤二将0°相与180°相射频回波信号加和运算所得射频回波信号X(t)， 依据第i条待合成扫描线上的第d个目标点空间位置与动态孔径N，对N个有 效阵元通道数据延时后得到加和延时时域信号Xid(t)；
步骤三对加和延时时域信号Xid(t)通过傅立叶变换转换为频域信号Xid(ω)，并选取造影回波的谐波分量内的频点Xid(ωk)，进行按比率的压缩随机采样抽取，得到压缩采样后的造影谐波频域信号
步骤四通过构造所选频点的稳健自相关矩阵，计算该频点所对应的最小方差自适应波束合成的最优加权系数wid(ωk)，从而得到该频点的最优波束合成的频域输出
步骤五按照步骤四遍历所有抽取的造影谐波频点并由多点欠定系统聚焦求解(RM-FOCUSS)算法重建谐波分量的频域信息，然后重建第i条扫描线的频域信息yi(ω)；
步骤六按扫描线波束合成次序，重复步骤二～步骤五得到所有扫描线的频 域信息Y(ω)，并将加和射频数据的频域信息Y(ω)通过逆傅立叶变换变换到时 域，得到Y(t)，对Y(t)进行包络检波、对数压缩、坐标变换等图像后处理得到 高质量的平面波造影图像。
具体步骤如下：
一、本发明所采用的技术方案步骤一，具体步骤如下：
1)超声阵列换能器按脉冲逆转模式交替发射中心频率为f0的一对相位相反 的平面波，其相位分别为0°相与180°相，探头工作频率f0范围：2.0MHz～5.0MHz；
2)计算机采集超声阵列换能器阵元的0°相与180°相射频回波信号X0(t)与 X180(t)；
3)其成像特征为，合成一副平面波造影图像，仅需发射2次，即：一对相 位相反的平面波。
二、本发明所采用的技术方案步骤二，具体步骤如下：
1)将0°相与180°相的射频回波信号X0(t)与X180(t)在整个时域加和运算 得到射频回波信号X(t)，X(t)＝X0(t)+X180(t)；
2)选定第i条(1≤i≤128～256)待合成扫描线上目标点并根据目标点 与中心阵元垂直距离计算波束合成的有效阵元数大小N，并选定阵元；
3)根据第i条待合成扫描线上第d个目标点与N个阵元相对位置，对N个阵元通道接收信号计算每一个通道的信号延时时间τid，τid＝[τ1τ2…τN]T，其中，为第j个阵元的延时时间，Di为第j个阵元到目标点距离，Dc为目标点中心阵元到目标点距离，c为组织中声速，j＝1,2,3,...,N；
4)根据第d个待合成目标点与N个阵元相对位置，对N个阵元对应射频 回波信号X(t)按τid依次延时，得到延时后列向量信号Xid(t)＝X(t-τid)， Xid(t)＝[x1(t)x2(t)…xN(t)]T。
三、本发明所采用的技术方案步骤三，具体步骤如下：
1)对延时后时域信号Xid(t)进行傅里叶变换，得到频域信号Xid(ω)， Xid(ω)＝[x1(ω)x2(ω)…xN(ω)]T；
2)对造影回波的频域信号Xid(ω)进行2f0谐波滤波，对谐波分量ωk内的 频点构成维数为k×N的谐波信号Xid(ωk)，其有效带宽为1.4f0～2.6f0，k为 有效带宽内的频点数；
3)基于该频点数k，构造k阶的主对角线元素为1、其余为0的压缩随机 抽取矩阵T，抽取比例为30％～40％；
4)对Xid(ωk)按行进行随机抽取：抽取比例为30％～40％，得到压缩采样后的造影谐波频域信号
四、本发明所采用的技术方案步骤四，具体步骤如下：
1)通过构造所选频点的稳健自相关矩阵，计算该频点所对应的最小方差自适应波束合成的最优加权系数wid(ωk)；
2)对于中的频点，把有效阵元数为N的阵列分为L×k、总数为N-L+1的子阵
3)将所有子阵与其对应的最优加权系数wid(ωk)依次点乘求和并平均，得到该频点的最优波束合成的频域输出
y ^ i d ( ω k ) = 1 N - L + 1 Σ l = 1 N - L + 1 w i d ( ω k ) X ^ i d l ( ω k ) - - - ( 1 ) ]]>
其中，为的第l个子阵。
五、本发明所采用的技术方案中，基于最小方差自适应波束合成的最优加 权系数wid(ωk)的计算方法，具体步骤如下：
1)基于互相关原理，构造所选频点的稳健自相关矩阵
R ~ i d ( ω k ) = R i d ( ω k ) - AA H R i d ( ω k ) AA H / ( A H A ) 2 - - - ( 2 ) ]]>
其中，[]H表示矩阵的共轭转置，Rid(ωk)为的自相关矩阵：
R i d ( ω k ) = 1 N - L + 1 Σ l = 1 N - L + 1 X ^ i d l ( ω k ) X ^ i d l ( ω k ) H - - - ( 3 ) ]]>
其中，A为长度为L的1维单位向量；AAHRid(ωk)/(AHA)2为造影回波所选频点的主信号能量项；
2)利用自相关矩阵和单位向量A计算最优加权系数wid(ωk)：
w i d ( ω k ) = R ~ i d - 1 ( ω k ) A / A H R ~ i d - 1 ( ω k ) A - - - ( 4 ) ]]>
其中，[]-1表示矩阵的逆，为避免求逆运算出现不稳定性，将奇异值分解出现的负特征值替换为分解所得的最小正值特征解；同时，从Rid(ωk)中减掉的主信号能量项，以增强的稳健性。
六、本发明所采用的技术方案步骤五，具体步骤如下：
1)重复步骤四，计算第i条待合成扫描线上所有抽取频点的最优波束合成的频域信号输出
2)构建k×k的傅里叶基F作为稀疏基，利用下式：
y ^ i H ( ω k ) = TFb i H ( ω k ) - - - ( 5 ) ]]>
因抽取矩阵T与傅里叶基F相乘矩阵不满秩，故采用RM-FOCUSS算法和构建的抽取矩阵T与傅里叶基F，从阵元矩阵中估计出稀疏投影系数向量bi(ωk)；
3)由yi(ωk)＝Fbi(ωk)计算阵元输出所有有效带宽内的最优波束合成的 频域信号输出yi(ωk)；
4)对yi(ωk)中造影谐波有效带宽外的频点替换为0点，得到第i条扫描线 上整个频域的最优波束合成的频域信号输出yi(ω)。
七、本发明所采用的技术方案步骤六，具体步骤如下：
1)按扫描线波束合成次序，得到所有扫描线在整个频域的最优波束合成的 频域信号输出Y(ω)，Y(ω)＝[y1(ω)y2(ω)…yline(ω)]T，其中，line为 波束合成扫描线数，其取值一般为128～256；
2)将0°相与180°相的射频回波信号脉冲逆转加和后并2f0谐波滤波后的最 优波束合成的频域信号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)超声阵列换能器按脉冲逆转模式交替发射一对中心频率为f0的相位相反 的平面波，而计算机采集超声线阵换能器接收0°相与180°相射频回波信号，仅 用2次平面波发射得到一副平面波造影图像，以减少对微泡的辐照，减低微泡的 破坏率，同时最大程度避免运动等杂波信号的干扰，准确扑捉微泡的瞬态分布。
2)将0°相与180°相射频信号加和处理，并从时域信号转换为频域信号， 按其频域信息分布，选取频域上造影微泡谐波分量有效带宽内的频点，并按比率 的压缩随机采样(抽取)；这个过程，通过谐波滤波、压缩采样，大幅度剔除冗 余数据、减少自适应波束合成计算量(减少50％～70％)，为平面波造影快速成像 奠定数据样本基础。
3)对所抽取的频点构造稳健自相关矩阵，并计算该频点所对应的最小方差 自适应波束合成的最优加权系数，从而得到该频点的最优输出；此过程中，考虑 阵元信号的差异性，通过压缩感知方法提高自适应波束合成方法性能，在提高运 算速度的同时获得高质量的CSBF频域数据。
4)按照步骤2)～3)遍历第一条扫描线上所有抽取的频点，并由压缩感知 理论重建谐波分量的该相位频域信息，然后通过填0重建整个频域信息；将第一 条扫描线上整个频域信息变换到时域射频信号。
5)按扫描线波束合成次序，重复上述2)～4)步骤，得到所有扫描线相应 谐波滤波后的CSBF射频数据，并进行包络检波、对数压缩、坐标变换等射频后 处理得到平面波造影图像。
参照图1所示流程图，所述平面波快速造影成像方法，具体步骤如下：
(1)超声阵列换能器按脉冲逆转模式交替发射中心频率为f0的一对相位相 反的平面波，其相位分别为0°相与180°相，探头工作频率f0范围： 2.0MHz～5.0MHz；
(2)计算机采集超声阵列换能器阵元的0°相与180°相射频回波信号X0(t) 与X180(t)；
(3)将0°相与180°相的射频回波信号X0(t)与X180(t)在整个时域加和运 算得到射频回波信号X(t)：
X(t)＝X0(t)+X180(t)(1)
(4)对造影回波数据X(t)，计算第i条(1≤i≤128～256)待合成扫描 线上目标点并根据目标点与中心阵元垂直距离计算波束合成的有效阵元数(即有 效孔径)大小，记为N；
(5)根据第i条待合成扫描线上第d个目标点与N个阵元相对位置，对N 个阵元通道接收信号计算每一个通道的信号延时时间τid；
τid＝[τ1τ2…τN]T(2)
其中，为第j个阵元的延时时间，Di为第j个阵元到目标点距离，Dc为为目标点中心阵元到目标点距离，c为组织中声速；
(6)根据第d个待合成目标点与N个阵元相对位置，对N个阵元对应射 频回波信号X(t)按τid依次延时，得到延时后列向量信号Xid(t)＝X(t-τid)：
Xid(t)＝[x1(t)x2(t)…xN(t)]T(3)
(7)对Xid(t)进行傅里叶变换，得到Xid(ω)：
Xid(ω)＝[x1(ω)x2(ω)…xN(ω)]T(4)
(8)根据超声发射频率f0，选取2f0有效带宽ωk(1.4f0～2.6f0有效频 带)的频点构成维数为k×N的信号Xid(ωk)，k为有效带宽内频点数；
(9)在有效带宽ωk内，构造维数为k×k、主对角线元素为1、其余为0 的抽取矩阵T；
(10)对谐波滤波后的频域信号Xid(ωk)按行进行随机抽取： X ^ i d ( ω k ) = T × X i d ( ω k ) , ]]>抽取比例为30％～40％；
(11)对于抽取的频点把有效阵元数为N的阵列分为长度为L的子阵，子阵的总数为N-L+1，采用空间平滑技术计算的空间平滑相关矩阵Rid(ωk)：
R i d ( ω k ) = 1 N - L + 1 Σ l = 1 N - L + 1 X ^ i d l ( ω k ) X ^ i d l ( ω k ) H - - - ( 5 ) ]]>
其中，为的第l个子阵；
(12)基于互相关原理，构造所选频点的稳健自相关矩阵
R ~ i d ( ω k ) = R i d ( ω k ) - AA H R i d ( ω k ) AA H / ( A H A ) 2 - - - ( 6 ) ]]>
其中，[]Η表示矩阵的共轭转置，A为长度为L的1维单位向量；AAHRid(ωk)/(AHA)2为造影回波所选频点的主信号能量项；
(13)利用自相关矩阵和单位向量A计算最优加权系数：
w i d ( ω k ) = R ~ i d - 1 ( ω k ) A / A H R ~ i d - 1 ( ω k ) A - - - ( 7 ) ]]>
其中，[]-1表示矩阵的逆，为避免求逆运算出现不稳定性，将奇异值分解出现的负特征值替换为分解所得的最小正值特征解；同时，从Rid(ωk)中减掉的主信号能量项，以增强的稳健性。
(14)通过下式计算该频点的最优波束合成的频域输出
y ^ i d ( ω k ) = 1 N - L + 1 Σ l = 1 N - L + 1 w i d ( ω k ) X ^ i d l ( ω k ) - - - ( 8 ) ]]>
(15)重复步骤(11)至步骤(14)，计算第i条扫描线上所有抽取频点的最优波束合成频域输出
(16)构建k×k的傅里叶基F作为稀疏基，利用下式：
y ^ i H ( ω k ) = TFb i H ( ω k ) - - - ( 9 ) ]]>
因抽取矩阵T与傅里叶基F相乘矩阵不满秩，故采用RM-FOCUSS算法从阵元矩阵中估计出稀疏投影系数向量bi(ωk)；
(17)由yi(ωk)＝Fbi(ωk)计算阵元输出所有有效带宽内的最优波束合成 的频域信号输出yi(ωk)；
(18)对有效带宽外频点以0填充，得到第i条扫描线上整个频域的最优 波束合成频域输出yi(ω)；
(19)按扫描线波束合成次序，重复步骤(4)至步骤(18)，得到所有扫 描线在整个频域的最优波束合成的频域信号输出Y(ω)， Y(ω)＝[y1(ω)y2(ω)…yline(ω)]T，其中，line为波束合成扫描线数，其 取值一般为128～256；
(20)对Y(ω)进行逆傅里叶变换，得到整个造影微泡成像区域的射频信号 Y(t)；
(21)对造影信号Y(t)进行包络检波、对数压缩、坐标变换等图像后处理 得到高质量的平面波造影图像。
实施平台为Sonix-Touch，选用128阵元线阵换能器为成像探头，按脉冲逆 转模式，发射中心频率为3MHz的平面波，其负峰值声压约为178kPa，并在外 触发源的控制下，按触发频率200Hz工作。参见图2，实施对象一狭窄血管流动 仿体，其正常直径为5mm，按0.2mL的剂量对浓度为2.0×108MB/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 与分辨率的平面波造影成像的实时实现成为可能。