CN107204021A

CN107204021A - 基于高斯函数探头响应模型和压缩感知的超声成像方法

Info

Publication number: CN107204021A
Application number: CN201710274228.3A
Authority: CN
Inventors: 王丛知; 杨新新; 刘佳妹; 肖杨; 郑海荣
Original assignee: Shenzhen Institute of Advanced Technology of CAS
Current assignee: Shenzhen Institute of Advanced Technology of CAS
Priority date: 2017-04-25
Filing date: 2017-04-25
Publication date: 2017-09-26
Anticipated expiration: 2037-04-25
Also published as: CN107204021B

Abstract

本发明涉及超声成像领域，更具体地涉及基于高斯函数探头响应模型和压缩感知的超声成像方法。本发明的目的是：在保证高帧频和较高成像质量的超声成像的同时，实现只需要较低水平的硬件计算平台就能进行成像。具体地，本发明提供了通过超声回波射频信号相对于超声平面波发射时刻的时间延迟信息与散射子分布的空间位置信息之间的关系，同时利用高斯函数近似描述超声信号经探头发射和接收时的响应模型，并考虑超声信号在介质中传播时的幅度衰减，建立在时域表达的关系矩阵M。本发明可以应用于需要高帧频和高分辨率成像且硬件计算平台水平较低的领域，如医学成像等，具有较大的实用价值。

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)中所述的反映超声回波射频信号与网格节点上散射子散射强度分布之间关系的数学模型，已发表的主要有以下两种：①基于待成像介质的可压缩性分布情况的比较复杂的模型(德国波鸿鲁尔大学MartinF.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。

在一些实施方式中，计算所述超声成像系统的超声回波信号的时间延迟可以包括计算所述超声回波信号相对于散射子分布的空间位置的时间延迟t_n,k，以及计算所述超声回波信号相对于超声发射时刻的时间延迟t_d,k，其中，所述时间延迟t_n,k的计算的具体过程可以为：所述超声成像系统包括k个阵元，第k个阵元的坐标为(x_k,0)，所要生成的超声图像的像素点数为N＝N_x×N_z，其中N_x，N_z分别是x方向和z方向上图像像素的行数和列数，某个网格节点处的散射子n的坐标为(x_n,z_n)，则超声信号从超声发射时刻开始，经过散射子n的散射，再回到所述超声成像系统的第k个阵元位置的所述时间延迟t_n,k表示为：

其中c是声波在介质中的传播速度；

其中，所述时间延迟t_d,k的计算的具体过程可以为：从超声发射时刻开始，经过t₀的时间延迟，开始对超声回波信号进行采样，采样频率为f_s，每个通道的采样点数为D，则第k个通道接收到超声回波信号的第d个数据的时间点，相对于超声发射时刻的时间延迟t_d,k表示为：

t_d,k＝t₀+(d-1)/f_s。

在一些实施方式中，利用高斯函数近似描述超声信号的响应函数模型的具体过程为：采用高斯函数近似描述所述超声信号经探头发射或接收时的响应函数h_trans：

其中：

f₀是超声信号的中心频率，B是超声探头的响应带宽比值，综合考虑发射和接收的全过程，得到总的响应函数h_TR为：

在一些实施方式中，计算所述超声回波信号在介质中传播时的幅度衰减的具体过程可以为：设衰减系数为α，单位为dB/Hz.m，则对于某一个网格节点处的散射子n，得到幅度衰减的倍数A_n为：

在一些实施方式中，建立关系矩阵M的具体过程可以为：对每一个网格节点处的散射子n，生成长度为D×K的向量m：

针对N个散射子，组合向量m以建立为D×K行、N列的矩阵M。

在一些实施方式中，所述方法还可以包括以下步骤：设置一个阈值，将所述关系矩阵M中低于所述阈值的元素都置为0；通过成熟的压缩感知算法对所述方程组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个阵元的坐标为(xk,0)。所要成的超声图像的像素点数(即对成像平面划分的网格节点数)为N＝N_x×N_z，其中N_x，N_z分别是x方向(宽度方向)和z方向(深度方向)上图像像素的行数和列数。某一个网格节点处的散射子n的坐标为(x_n,z_n)。则有，超声信号从超声发射时刻开始，经过这个散射子的散射，再回到超声换能器阵列的第k个阵元位置，总的时间延迟为：

其中c是声波在介质中的传播速度。对应的，可以得到N×K个时间延迟数据。

在该实施方式中，已知从超声发射时刻开始，经过t₀的时间延迟，开始对超声回波射频信号进行采样，采样频率为f_s，每个通道的采样点数为D，则第k个通道接收到超声回波射频信号的第d个数据的时间点，相对于超声发射时刻的时间延迟为：

t_d,k＝t₀+(d-1)/f_s，

相应地，可以得到D×K个时间延迟数据。

超声信号经探头发射或接收时的响应函数h_trans，通常采用高斯函数来进行近似：

其中：

其中：f₀是超声信号的中心频率，B是超声探头的响应带宽比值(小于1)。则综合考虑发射和接收的全过程，总的响应函数h_TR为：

进一步地，考虑超声信号在介质中传播时的幅度衰减，设其衰减系数为α，单位为dB/Hz.m，则对于某一个网格节点处的散射子n，幅度衰减的倍数为：

如前所述，某一个网格节点处的散射子n对应于超声换能器阵列上通道k的时间延迟为t_n,k，则经该散射子反射并被通道k接收的信号可以通过时间延迟t_n,k，t_d,k，A_n和响应函数h_TR计算得到：

因此，对每一个网格节点处的散射子n，可以生成一个向量m，长度为D×K。一共有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段)，对截取出的每一段数据进行短时傅立叶变换以获得频域信号，并重复后续的求解计算。因此模型②在总的计算量上远大于本发明所使用的模型。此外，模型②为了方便的对信号进行时间延迟计算，所有运算都在频域进行。而短时傅里叶变换计算不仅增加了计算量，也会引入由于信号长度有限产生的计算误差，进而影响到最终的成像质量。

现在参考图1和2，其示出了利用超声成像仿真软件Field II对超声成像进行仿真的结果。图1示出了采用根据本发明的实施方式的算法对点状的稀疏仿体进行成像的仿真实验结果。图2示出了采用传统的延时叠加法(DAS)对点状的稀疏仿体进行成像的仿真实验结果。从图1和图2可见，本发明方法可以去除传统方法中出现的横向伪影。

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

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

Claims

1.一种基于高斯函数探头响应模型和压缩感知的超声成像方法，其特征在于，包括以下步骤：

计算超声成像系统的超声回波信号的时间延迟；

利用高斯函数近似描述超声信号经过探头发射和接收时的响应函数模型；

计算所述超声回波信号在介质中传播时的幅度衰减；

基于所述超声回波信号的所述时间延迟、所述高斯函数探头响应模型、和所述超声回波信号的幅度衰减建立关系矩阵M；以及

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

2.根据权利要求1所述的方法，其特征在于，计算所述超声成像系统的超声回波信号的时间延迟包括计算所述超声回波信号相对于散射子分布的空间位置的时间延迟t_n,k，以及计算所述超声回波信号相对于超声发射时刻的时间延迟t_d,k，

其中，所述时间延迟t_n,k的计算具体为：所述超声成像系统包括K个阵元，第k个阵元的坐标为(x_k,0)，所要生成的超声图像的像素点数为N＝N_x×N_z，其中N_x，N_z分别是x方向和z方向上图像像素的行数和列数，某个网格节点处的散射子n的坐标为(x_n,z_n)，则超声信号从超声发射时刻开始，经过散射子n的散射，再回到所述超声成像系统的第k个阵元位置的所述时间延迟t_n,k表示为：

其中c是声波在介质中的传播速度；

其中，所述时间延迟t_d,k的计算具体为：从超声发射时刻开始，经过t₀的时间延迟，开始对超声回波信号进行采样，采样频率为f_s，每个通道的采样点数为D，则第k个通道接收到超声回波信号的第d个数据的时间点，相对于超声发射时刻的时间延迟t_d,k表示为：

t_d,k＝t₀+(d-1)/f_s。

3.根据权利要求1或2所述的方法，其特征在于，利用高斯函数近似描述超声信号的响应函数模型的具体过程为：采用高斯函数近似描述所述超声信号经探头发射或接收时的响应函数h_trans：

其中：

4.根据权利要求1或2所述的方法，其特征在于，计算所述超声回波信号在介质中传播时的幅度衰减的具体过程为：设衰减系数为α，单位为dB/Hz.m，则对于某一个网格节点处的散射子n，得到幅度衰减的倍数A_n为：

<mrow> <msub> <mi>A</mi> <mi>n</mi> </msub> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msup> <mn>10</mn> <mrow> <msub> <mi>&alpha;f</mi> <mn>0</mn> </msub> <msub> <mi>ct</mi> <mrow> <mi>n</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mo>/</mo> <mn>20</mn> </mrow> </msup> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mo>.</mo> </mrow>

5.根据权利要求1或2所述的方法，其特征在于，建立关系矩阵M的具体过程为对每一个网格节点处的散射子n，生成长度为D×K的向量m：

针对N个散射子，组合向量m以建立D×K行、N列的矩阵M。

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

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

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

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

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

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

10.根据权利要求1或2所述的方法，其特征在于，所述超声信号为超声平面波信号、超声凸面波信号或超声凹面波信号。