CN106940883A - 基于超声系统点扩散函数仿真和压缩感知的超声成像方法 - Google Patents

Abstract

本发明涉及超声成像领域，更具体地涉及基于超声系统点扩散函数仿真和压缩感知的超声成像方法。本发明的目的是：在保证高帧频和较高成像质量的超声成像的同时，实现只需要较低水平的硬件计算平台就能进行成像。为此本发明提供的所述超声成像方法包括以下步骤：获得所述超声成像系统的点扩散函数；根据所述点扩散函数建立关系矩阵M；建立反映超声回波射频信号s与超声图像的像素点上的散射子散射强度分布I之间关系的方程组：s＝MI；其特征在于，设置一个阈值，将所述关系矩阵M中低于所述阈值的元素都置为0。本发明可以应用于需要高帧频和高分辨率成像且硬件计算平台水平较低的领域，如医学成像等，具有较大的实用价值。

Description

基于超声系统点扩散函数仿真和压缩感知的超声成像方法

技术领域

本发明涉及超声成像领域，更具体地涉及基于超声系统点扩散函数仿真和压缩感知的超声成像方法。

背景技术

传统的超声成像工作原理，以线阵式超声换能器为例，线阵换能器拥有N个可以独立发射/接收的阵元，对应于N个超声发射通道和信号接收通道，超声发射时，利用多个通道的延时发射，使不同通道的超声信号同时到达聚焦位置，形成发射聚焦；接收回波时，对接收到的信号进行类似的延时，将从同一反射物返回的不同通道接收到的信号累加在一起，形成接收聚焦。这样一次发射和一次接收，可以形成一条扫描线。通常超声成像都采用电子扫描的方式，进行M次聚焦发射/接收来获得M条扫描线，再将这些扫描变换成一幅完整的二维图像。因此，传统超声成像的帧频很低，通常在十几帧到几十帧之间。对于高硬度组织弹性成像、大动脉高速血流成像、心脏成像、以及追踪超声造影剂状态变化等成像目标快速运动因而需要极高帧频的应用领域，传统超声成像的帧频已经远远无法满足需要。

超声平面波成像技术包括超声平面波发射和相应的超声回波波束形成技术，是国际上近年来关于提高超声成像帧频的热点研究方向。该技术能将传统的超声成像帧频(一般为十几帧到几十帧)提高几百倍，达到10000-20000帧。该方法一般将线阵换能器的所有阵元都用于发射，采用互相之间没有相对延时的相同的电压脉冲，同时激励线阵换能器各个阵元以产生沿垂直于换能器表面方向向前传播的超声平面波；接收回波信号时，采用基于图像像素点位置的DAS(Delay and Sum，延时叠加法)波束形成技术形成一幅二维图像。这样，只需要一次发射/接收即可完成一次二维成像，极大地提高了成像帧频。但是，由于在使用平面波成像技术时，超声能量均匀地分布在整个二维成像平面，从不同散射物反射的回波会混叠在一起，被各个通道接收，难以区分。因此，由普通的波束形成方法得到的图像会出现很明显的伪影干扰。

为解决这一问题，多角度相干叠加成像法被提出。该方法从2N+1(N为某一正整数)个角度(其中一个角度为通常使用的垂直于超声换能器表面的角度，其他2N个角度围绕这个垂直角度呈对称形态分布，如-2°，-1°，0°，1°，2°)发射超声平面波并同样采用基于图像像素点位置的DAS波束形成技术获得2N+1幅二维图像，将这些图像进行叠加，相当于从多个角度发射的超声平面波之间实现了相干增强，产生类似于聚焦的效果，从而实现了图像分辨率和对比度的增强。N值越大，对分辨率和对比度的提高效果就越显著。利用这一技术，已经实现了高时空分辨率的，能够对全脑微血管响应脑部活动所产生的动态变化进行实时成像的新技术——超声脑功能成像技术(functional ultrasound,fUS)。高达千赫兹数量级的帧频成像效果，是研究动态血流变化情况的关键。此外，该技术还被应用到多个生物医学超声学的前沿研究方向上，如实时三维超声成像、高速多普勒血流流场速度分布成像、二维实时弹性成像、心脏、大动脉应变成像等，具有十分广阔的应用前景。但是，多角度相干叠加成像法相当于再次降低了帧频，例如，采用普通的超声平面波成像方法，可以实现10000帧每秒的帧频，但为了提高图像的分辨率和对比度，改为采用多角度相干叠加成像法，由51个角度的发射/接收结果合成一幅图像，帧频就下降到了低于200帧每秒。因此，多角度相干叠加成像法的应用范围受到严重地制约。

综上所述，如何在保证帧频不下降的同时，尽可能地提高图像的分辨率和对比度，成为超声平面波成像需要解决的重要问题。

近年来，国内外先后有一些关于基于压缩感知的超声平面波成像方法的论文发表。这些方法都分为以下两个步骤：

(1)将图像的每个像素点视为二维平面内的一个网格节点，假设每个网格节点处都存在一个散射子可以造成入射超声的散射，形成回波，则可以认为，我们所要成的超声图像实际反映的是网格节点上散射子的散射强度在二维平面内的分布。首先，需要建立反映超声回波射频信号s与网格节点上散射子散射强度分布I之间关系的数学模型，形成如下形式的方程组：

s＝MI

其中矩阵M为关系矩阵。由于实际成像时超声回波射频信号中噪声的存在，求解该方程组通常情况下是一个不适定问题，无法求得唯一解。

(2)当I稀疏(sparse)时，即其中的非零元素数量远小于零元素数量时，则可以通过压缩感知方法对上述方程组进行求解：

其中β反映了允许多少噪声成分存在。

对于超声成像来说，在第(1)步骤中，如何根据其所遵循的物理原理，建立能够尽可能真实反映s和I之间关系的数学模型，并据此建立便于完成后续迭代计算的矩阵M，是决定超声成像质量和成像方法实用性的关键。而第(2)步骤中，方程组求解的具体计算方法已经有很多成熟的数值迭代方法可供选择，不属于本发明阐述的重点。

对于步骤(1)中所述的反映超声回波射频信号与网格节点上散射子散射强度分布之间关系的数学模型，已发表的主要有以下三种：①基于待成像介质的可压缩性分布情况的比较复杂的模型(德国波鸿鲁尔大学Martin F.Schiffner and Georg Schmitz等)；②基于频域信号延时的比较简单的模型(汕头大学沈民奋等)。

模型①的最终形式为：

其中G是一个N_elN_k×N的矩阵，N_el是超声换能器阵列接收回波信号的通道数，N_k指将宽带的超声回波信号分为N_k个离散的波数k_l，1≤l≤N_k，N＝N_xN_z是图像的总的像素数(或者说是网格节点数)，N_x、N_z分别是x方向(宽度方向)和z方向(深度方向)上图像像素的行数和列数。矩阵G中的每个元素定义为：

其中m表示换能器上的第m个阵元，1≤m≤N_el，i表示图像上的第i个像素，1≤i≤N，代表入射超声的声压，r_el,m表示超声换能器上第m个阵元的空间位置，r_i表示图像上第i个像素的位置，g_l(r_el,m-r_i)是开放空间的格林函数，定义为：

其中j表示虚部，是零阶的第二类Hankel函数。p^sc代表超声回波射频信号，γ_κ代表所要成像的介质可压缩性的分布情况(介质的可压缩性是决定其声散射强度的主要因素)。

模型②的最终形式为：

X(ω)＝A(ω)·S(ω)

由于是基于频域信号进行处理，实际上取ω＝2πf₀，其中f₀为所使用的超声换能器的中心发射频率。X为经过短时傅里叶变换后的超声回波射频信号，S为所要成像的散射子散射强度在频域上对应于f₀的映射，A为K×L的由时间延迟数据组成的关系矩阵，定义为：

[A(ω)_k]_i＝exp[jωτ_k(ρ_i)]

其中K为超声换能器阵列接收回波信号的通道数，L为图像的总的像素数(或者说是网格节点数)，1≤k≤K，1≤i≤L，ρ_i表示图像上的一个像素点(或网格节点)，表示从某个像素点发出的回波信号到达某个超声换能器阵列通道的时间延迟，表示某个超声换能器阵列通道的空间位置，表示某个像素点ρ_i的空间位置。需要指出，X是从全部超声回波射频信号中截取出一小段进行短时傅立叶变换后获得的频域信号，因此若全部超声回波射频信号被分成Q段，则要完成全部成像，需要将后续的求解过程重复Q次。

建立上述两个模型后，通过压缩感知方法进行方程组求解，就可以解出γ_κ(模型①)或者S(模型②)，再将其从向量变换为对应图像像素点数量的矩阵，就可以显示为我们所希望获得的图像。

上述两个反映超声回波射频信号与网格节点上散射子散射强度分布之间关系的数学模型，各有其局限性。

模型①的建立，是从已经被证明比较准确的声传播和散射的数学模型出发，优点是能够比较真实的反映声在介质中的各种物理现象，但缺点也非常明显，就是模型过于复杂。关系矩阵的尺寸过于庞大，需要占用大量内存，同时也造成后续的求解过程的计算量十分巨大。以其论文中进行的成像实验数据为例，当N_x＝400，N_z＝600，N_el＝128，N_k＝1000时，矩阵G占用的内存高达458GB。因此，为实现该算法，只好采用每次调用G时都重新计算其各个元素数值的方法来进行，极大的增加了计算量。而且，实际上上述参数值根本无法满足正常医学超声成像的需要，如果成像深度超过5cm，N_z的取值通常都在3000以上，因此内存占用量还要再增加到5倍，已经完全不是普通计算机能够承担的任务。

模型②由于只考虑超声回波射频信号的时间延迟，同时只考虑超声的中心发射频率f₀而不考虑信号其他频率成分，因此其所使用的关系矩阵A(ω)的规模大为缩小。但是，该模型仍有以下几个问题。首先，矩阵A(ω)的所有元素都是非零的，而且X是从全部超声回波射频信号中截取出一小段进行短时傅立叶变换后获得的频域信号，若全部超声回波射频信号被分成Q段，则要完成全部成像，需要将后续的求解过程重复Q次，因此进行后续的矩阵乘法运算时计算量仍然很大。其次，为了方便的对信号进行时间延迟计算，该模型的所有运算都在频域进行。这就需要首先将时域的超声回波射频信号，通过短时傅里叶变换转换到频域。这个过程不仅增加了计算量，也会引入由于信号长度有限产生的计算误差，进而影响到最终的成像质量。

有鉴于此，需要开发一种新的方法来克服上述缺陷。

发明内容

针对上述现有技术的不足，本发明的目的是，相对于模型①，尽可能的简化关系矩阵，降低运算时的内存存储空间和计算量；相对于模型②，避免使用傅立叶变换和频域计算。

为了达到上述目的，本发明提出了一种基于超声成像系统点扩散函数仿真和压缩感知的超声成像方法，以在保证高帧频和较高成像质量的同时，实现只需要较低水平的硬件计算平台就能使用本方法进行成像，便于实现本发明的产业转化。

本发明提供了一种基于超声成像系统点扩散函数仿真和压缩感知的超声成像方法，其可以包括以下步骤：获得所述超声成像系统的点扩散函数；根据所述点扩散函数建立关系矩阵M；以及建立反映超声回波射频信号s与超声图像的像素点上的散射子散射强度分布I之间关系的方程组：s＝MI；其特征在于，设置一个阈值，将所述关系矩阵M中低于所述阈值的元素都置为0。

在一些实施方式中，通过仿真计算获得所述超声成像系统在某种超声发射模式下在与超声图像上所有像素点对应的空间位置处的所述点扩散函数。

在一些实施方式中，所述仿真计算可以包括以下步骤：(i)针对所述超声成像系统设定各种参数值；(ii)利用成熟的声学仿真方法或软件进行仿真计算，模拟当在与超声图像上某个像素点对应的空间位置上放置尺寸足够小的散射子而其他位置没有任何散射子时形成的被所述超声成像系统接收的超声回波信号，其中被接收的所述超声回波信号即所述系统在所述某种超声发射模式下在所述空间位置的点扩散函数；(iii)对所述超声图像上的所有像素点执行步骤(i)和(ii)，得到所述超声成像系统在某种超声发射模式下在与所述超声图像上所有像素点对应的空间位置处的所述点扩散函数。

在一些实施方式中，所述各种参数值包括但不限于：介质声速/衰减系数、超声探头的类型/阵元数/阵元尺寸和间距/带宽、所使用的超声中心频率、发射波形/焦距、对超声回波信号的采样频率/采样长度、成像范围。

在一些实施方式中，所述成熟的声学仿真方法或软件包括但不限于：时域有限差分法(FDTD)、k空间频谱法、Field II、k-wave。

在一些实施方式中，所述方法还可以包括以下步骤：通过成熟的压缩感知算法对所述方程组s＝MI进行求解，得到向量I，将其从向量变换为对应图像像素点数量的矩阵，再经过调整动态范围和数字扫描变换等步骤，获得希望得到的超声图像。

在一些实施方式中，所述成熟的压缩感知算法包括但不限于匹配追踪法(matching pursuit method)、Bregman算法、operator/variable splitting、FPC(Fixed-point continuation：定点连续)算法、L1-magic算法、牛顿下降法等。

在一些实施方式中，当I不是稀疏的时，可以对I进行稀疏变换Ψ，令θ＝ΨI，其中θ是I在稀疏变换域内的系数，θ是稀疏的，可根据压缩感知的求解公式来求解I，其中β表示允许存在多少噪声成分。

在一些实施方式中，所述稀疏变换Ψ可以包括但不限于离散余弦变换(DCT)、各种小波变换等。

在一些实施方式中，所述超声成像系统可以采用超声平面波发射模式。

在一些实施方式中，所述超声成像系统可以采用超声凸面波发射模式。

在一些实施方式中，所述超声成像系统可以采用超声凹面波发射模式。

在一些实施方式中，所述超声成像系统可以采用任意不规则波形信号激励的发射模式。

本发明在保证高帧频和较高成像质量的超声成像的同时，实现只需要较低水平的硬件计算平台就能进行成像，因而具有较大的实用价值。

本领域技术人员在阅读整个说明书和权利要求书时将理解本发明的这些优点和其它优点。

附图说明

图1示出了采用根据本发明的实施方式的算法对点状的稀疏仿体进行成像的仿真实验结果。

图2示出了采用传统的延时叠加法(DAS)对点状的稀疏仿体进行成像的仿真实验结果。

具体实施方式

下面结合附图对本发明的具体实施例进行说明。在下文所描述的本发明的具体实施例中，为了能更好地理解本发明而描述了一些很具体的技术特征，但显而易见的是，对于本领域的技术人员来说，并不是所有的这些技术特征都是实现本发明的必要技术特征。下文所描述的本发明的一些具体实施例只是本发明的一些示例性的具体实施例，其不应被视为对本发明的限制。另外，为了避免使本发明变得难以理解，对于一些公知的技术没有进行描述。

在一实施方式中，首先使用超声成像系统的计算机控制超声发射电路对超声换能器阵列进行激励，以发射超声信号。当超声换能器阵列的各个通道(每个通道对应于一个阵元)被同时激励时，发射出的超声信号是一组平面波信号，即可以认为其波阵面垂直于超声发射方向，波阵面抵达成像平面内某一深度的时间是一致的。该超声信号在介质中传播，发生散射，形成超声回波信号。超声回波信号被超声换能器阵列接收，再由超声接收电路采样，形成超声回波射频信号。超声回波射频信号被送回到所述计算机中，并在所述计算机中实现本发明所述的基于压缩感知的超声平面波成像。

在该实施方式中，已知超声换能器阵列包括K个阵元，其中第k个阵元的坐标为(x_k,0)。所要成的超声图像的像素点数(即对成像平面划分的网格节点数)为N＝N_x×N_z，其中N_x、N_z分别是x方向(宽度方向)和z方向(深度方向)上图像像素的行数和列数。某一个网格节点处的散射子n的坐标为(x_n,z_n)。超声回波射频信号的采样频率为f_s，每个通道的采样点数为D，则一次超声平面波发射/接收所采集的超声射频回波数据总数为D×K。

首先，通过仿真计算，获得所使用的超声成像系统在某种超声发射模式下的点扩散函数。在一个实施方式中，按照所使用的超声系统设定各种参数值，例如，介质声速/衰减系数、超声探头的类型/阵元数/阵元尺寸和间距/带宽、所使用的超声中心频率、发射波形/焦距、对超声回波信号的采样频率/采样长度、成像范围等(在一个实施方式中，介质声速为1540m/s，衰减系数为0.5dB/[MHz cm]，超声探头的类型为线阵，阵元数为128，阵元间距为0.3mm，带宽为60％，中心频率为5MHz，发射波形为方波，焦距为5cm，采样频率为40MHz，采样长度为5000，成像范围为40×100mm)，利用现有的成熟声学仿真方法(FDTD法，k空间频谱法等)或软件(如Field II，k-wave等)进行仿真计算，模拟当图像上某个像素点所对应的空间位置上放置散射子而其他位置没有任何散射子时所述超声成像系统所接收到的超声回波信号，该超声回波信号即该系统在某种超声发射模式下在该空间位置的点扩散函数，可以看作是一个尺寸为D×K的矩阵，并可进一步变换为一个长度为D×K的向量m。

因此，对每一个网格节点处的散射子n，可以生成一个向量m_n。计算该系统在某种超声发射模式下在图像上所有像素点所对应的空间位置上的点扩散函数，即一共有N个散射子，组合起来，就构造成为一个D×K行，N列的关系矩阵M。设I为所有网格节点上的散射子的散射强度，即I是一个长度为N的向量。则有：

s＝MI

另外，此时的关系矩阵M中，大部分元素的值接近于0。因此可以设置一个阈值，将M中低于该阈值的元素都置为0，则M可以通过稀疏表达的方式来存储和使用，大幅降低其占用的内存存储空间和计算量。

至此，反映超声回波射频信号s与网格节点上散射子散射强度分布I之间关系的方程组建立完毕。最后，通过成熟的压缩感知算法对上述方程组进行求解，得到向量I，将其从向量变换为对应图像像素点数量的矩阵，再经过调整动态范围和数字扫描变换等步骤，即可得到希望获得的超声图像。所述成熟的压缩感知算法包括但不限于匹配追踪法(matching pursuit method)、Bregman算法、operator/variable splitting、FPC(Fixed-point continuation：定点连续)算法、L1-magic算法、牛顿下降法等。

在该实施方式中，所述超声成像系统采用的是超声平面波发射模式。但实际上本方法并不限于只针对于超声平面波发射模式。例如，如果超声系统采用的是凸面波或者凹面波发射模式，甚至是采用任意不规则波形信号激励的发射模式，都可以利用本方法实现成像。

需要注意的是，压缩感知理论要求未知信号I是稀疏的。而在实际进行超声成像时，网格节点上的散射子的散射强度分布本身可能并不满足这一个条件。这种情况下，对I进行稀疏变换Ψ，令θ＝ΨI，其中θ是I在稀疏变换域内的系数。此时，θ是稀疏的，利用压缩感知的求解公式变为：

这样就可以先通过成熟的压缩感知算法对上述方程组进行求解求出I。其中，稀疏变换Ψ包括但不限于离散余弦变换(DCT)、各种小波变换等。

在本发明中所使用的模型，相对于模型①②都大为简化。解决了对信号时间延迟计算的时域表达问题，不再需要进行傅立叶变换和在频域进行计算，避免了因进行傅立叶变换而造成的计算误差。此外，由于本发明模型所建立的关系矩阵M最终被简化为了一个可以进行稀疏表达的矩阵，相比于模型①，大幅降低了其占用的内存存储空间和计算量。相比于模型②，由于需要将全部超声回波射频信号分成Q段(数据的前后两个分段之间需要有一定的重叠以提高深度方向的分辨率，因此4000多个数据点的长度至少需要被分成100段，在第一实施例中为100段，在第二实施例中为200段，在第三实施例中为400段)，对截取出的每一段数据进行短时傅立叶变换以获得频域信号，并重复后续的求解计算。因此模型②在总的计算量上远大于本发明所使用的模型。此外，模型②为了方便的对信号进行时间延迟计算，所有运算都在频域进行。而短时傅里叶变换计算不仅增加了计算量，也会引入由于信号长度有限产生的计算误差，进而影响到最终的成像质量。

现在参考图1和2，其示出了利用超声成像仿真软件Field II对超声成像进行仿真的结果。图1示出了采用根据本发明的实施方式的算法(其中，K＝128，D＝4364，N_x＝256，N_z＝3000)对点状的稀疏仿体进行成像的仿真实验结果。图2示出了采用传统的延时叠加法(DAS)对点状的稀疏仿体进行成像的仿真实验结果。从图1和图2可见，本发明方法所成图像的对比度明显优于传统方法。

因此，本发明所提出的新的超声成像方法，一方面可以实现超高帧频的快速超声成像，另一方面保证了较高的成像质量，同时只需要较低水平的硬件运算平台就能实现，便于实现产业转化。

尽管已经根据优选的实施方案对本发明进行了说明，但是存在落入本发明范围之内的改动、置换以及各种替代等同方案。还应当注意的是，存在多种实现本发明的方法和系统的可选方式。因此，意在将随附的权利要求书解释为包含落在本发明的主旨和范围之内的所有这些改动、置换以及各种替代等同方案。

Claims (10)

1.一种基于超声成像系统点扩散函数仿真和压缩感知的超声成像方法，其包括以下步骤：

获得所述超声成像系统的点扩散函数；

根据所述点扩散函数建立关系矩阵M；以及

建立反映超声回波射频信号s与超声图像的像素点上的散射子散射强度分布I之间关系的方程组：s＝MI；

其特征在于，

设置一个阈值，将所述关系矩阵M中低于所述阈值的元素都置为0。

2.根据权利要求1所述的方法，其特征在于，通过仿真计算获得所述超声成像系统在某种超声发射模式下在与超声图像上所有像素点对应的空间位置处的所述点扩散函数。

3.根据权利要求1或2所述的方法，其特征在于，所述仿真计算包括以下步骤：

(i)针对所述超声成像系统设定各种参数值；

(ii)利用成熟的声学仿真方法或软件进行仿真计算，模拟当在与超声图像上某个像素点对应的空间位置上放置尺寸足够小的散射子而其他位置没有任何散射子时形成的被所述超声成像系统接收的超声回波信号，其中被接收的所述超声回波信号即所述系统在所述某种超声发射模式下在所述空间位置的点扩散函数；

(iii)对所述超声图像上的所有像素点执行步骤(i)和(ii)，得到所述超声成像系统在某种超声发射模式下在与所述超声图像上所有像素点对应的空间位置处的所述点扩散函数。

4.根据权利要求3所述的方法，其特征在于，所述各种参数值包括：介质声速/衰减系数、超声探头的类型/阵元数/阵元尺寸和间距/带宽、所使用的超声中心频率、发射波形/焦距、对超声回波信号的采样频率/采样长度、成像范围。

5.根据权利要求3所述的方法，其特征在于，所述成熟的声学仿真方法或软件包括：时域有限差分法、k空间频谱法、Field II、k-wave。

6.根据权利要求1所述的方法，其特征在于所述方法还包括以下步骤：通过成熟的压缩感知算法对所述方程组s＝MI进行求解，得到向量I，将其从向量变换为对应图像像素点数量的矩阵，再经过调整动态范围和数字扫描变换步骤，获得希望得到的超声图像。

7.根据权利要求6所述的方法，其特征在于，所述成熟的压缩感知算法包括：匹配追踪法、Bregman算法、operator/variable splitting、FPC算法、L1-magic算法、牛顿下降法。

8.根据权利要求5或6所述的方法，其特征在于，当I不是稀疏的时，对I进行稀疏变换Ψ，令θ＝ΨI，其中θ是I在稀疏变换域内的系数，θ是稀疏的，根据压缩感知的求解公式来求解I，其中β表示允许存在多少噪声成分。

9.根据权利要求7所述的方法，其特征在于，所述稀疏变换Ψ包括：离散余弦变换(DCT)、各种小波变换。

10.根据权利要求1或2所述的方法，其特征在于，所述超声成像系统采用超声平面波发射模式、超声凸面波发射模式、超声凹面波发射模式、或任意不规则波形信号激励的发射模式。