技术领域
本发明涉及医学成像技术领域，特别涉及一种超声成像方法及装置。
背景技术
医学超声成像技术经过半个多世纪的发展，由于其实时性强、对软组织鉴别力高、易于使用和经济性好等优点，已经成为医学临床中应用最广泛的影像技术和临床多种疾病诊断的首选方法。
传统超声成像的帧频很低，通常在十几帧到几十帧之间。对于高硬度组织弹性成像、大动脉高速血流成像、心脏成像、以及追踪超声造影剂状态变化等成像目标快速运动因而需要极高帧频的应用领域，传统超声成像的帧频已经远远无法满足需要。
超声平面波成像技术包括超声平面波发射和相应的超声回波波束形成技术，超声成像技术能将传统的超声成像帧频(一般为十几帧到几十帧)提高几百倍，达到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帧每秒。因此，多角度相干叠加成像法的应用范围受到严重地制约。
综上所述，如何在保证帧频不下降的同时，尽可能地提高图像的分辨率和对比度，成为超声波成像需要解决的重要问题。
发明内容
为解决上述技术的问题，本发明提出一种超声成像方法及装置。
为实现上述目的，本发明提供了一种超声成像方法，该方法包括：
获取N×K个第一时间延迟数据tn,k；其中，所述第一时间延迟数据tn,k为超声信号从超声发射时刻开始，经过第n个散射子的散射，再回到超声换能器阵列的第k个通道位置的总的时间延迟；N表示所要成的超声图像的像素点数，K表示超声换能器阵列中阵元数目；n表示某一个散射子，散射子n＝1，2，……，N；k表示超声换能器阵列中某一个通道，k＝1,2，……，K；
获取D×K个第二时间延迟数据td,k；其中，所述第二时间延迟数据td,k为已知从超声发射时刻开始，经过时间延迟，开始对超声回波射频信号进行采样，得到超声回波射频数字信号，第k个通道接收到所述超声回波射频数字信号的第d个数据点的时刻相对于超声信号发送时刻的时间延迟；D表示每个通道的采样点数，d表示一通道接收到超声回波射频数字信号的某一数据点，d＝1,2，……，D；
将所述超声回波射频数字信号对应地矩阵变换为长度为D×K的向量s；
利用所述向量s中第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k确定矩阵M；其中，所述矩阵M为稀疏矩阵；
根据所述向量s和所述矩阵M建立方程组；其中，所述方程组为：s＝MI；I为表示散射子的散射强度分布的长度为N的向量；
对所述方程组通过压缩感知算法进行求解，得到向量I；
根据所述向量I确定所要成的超声图像。
优选地，所述确定矩阵M的步骤包括：
对第k个通道的数据点d建立长度为N的向量m，获得D×K个向量m；所述向量m中的每个元素确定方法为：如果第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k满足则从第n个散射子散射的超声回波信息包括在第k个通道的数据点d中，则向量m中第n个元素设为1，否则，设为0；
将D×K个向量m构造成一个D×K行、N列的矩阵M。
优选地，所述超声信号为超声平面波信号。
优选地，所述第一时间延迟数据tn,k的表达式为：
其中，第k个通道的空间位置坐标为(xk,0)；所要成的超声图像上的第n个散射子的空间位置坐标为(xn,zn)；c表示声波在介质中的传播速度；x方向是超声图像的宽度方向，z方向是超声图像的深度方向。
优选地，所述超声信号为超声凸面波信号、超声凹面波信号。
优选地，所述第二时间延迟数据td,k的表达式为：
td,k＝t0+(d-1)/fs
其中，t0表示从超声发射时刻到开始对超声回波射频信号进行采样时刻之间的时间延迟值，fs表示采样频率。
优选地，所述对所述方程组通过压缩感知算法进行求解的步骤包括：
如果向量I稀疏时，通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号；
否则，对向量I进行稀疏变换；通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，Ψ表示稀疏变换矩阵，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号。
优选地，所述压缩感知算法为匹配追踪法、Bregman算法、运算符/变量分裂算法、固定点连续算法、L1范数魔术算法、牛顿下降法。
优选地，所述确定所要成的超声图像的步骤包括：
将向量I变换为图像矩阵Φ；其中，所述图像矩阵Φ是一个Nx列Nz行的矩阵；Nx表示三维坐标系中x方向上图像像素的个数，Nz表示三维坐标系中z方向上图像像素的个数，N＝Nx×Nz；
对所述矩阵Φ进行处理，得到所要成的超声图像；其中，处理的方法包括：取信号包络、对数压缩、调整图像显示动态范围和进行数字扫描变换。
为实现上述目的，本发明还提供了一种超声成像装置，该装置包括：
第一时间延迟数据确定单元，用于获取N×K个第一时间延迟数据tn,k；其中，所述第一时间延迟数据tn,k为超声信号从超声发射时刻开始，经过第n个散射子的散射，再回到超声换能器阵列的第k个通道位置的总的时间延迟；N表示所要成的超声图像的像素点数，K表示超声换能器阵列中阵元数目；n表示某一个散射子，散射子n＝1，2，……，N；k表示超声换能器阵列中某一个通道，k＝1,2，……，K；
第二时间延迟数据确定单元，用于获取D×K个第二时间延迟数据td,k；其中，所述第二时间延迟数据td,k为已知从超声发射时刻开始，经过时间延迟，开始对超声回波射频信号进行采样，得到超声回波射频数字信号，第k个通道接收到所述超声回波射频数字信号的第d个数据点相对于超声信号发送时刻的时间延迟；D表示每个通道的采样点数，d表示一通道接收到超声回波射频数字信号的某一数据点，d＝1,2，……，D；
变换单元，用于将所述超声回波射频数字信号对应地矩阵变换为长度为D×K的向量s；
矩阵确定单元，用于利用所述向量s中第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k确定矩阵M；其中，所述矩阵M为稀疏矩阵；
模型建立单元，用于根据所述向量s和所述矩阵M建立方程组；其中，所述方程组为：s＝MI；I表示散射子的散射强度分布、长度为N的向量；
求解单元，用于对所述方程组通过压缩感知算法进行求解，得到向量I；
超声成像单元，用于根据所述向量I确定所要成的超声图像。
优选地，所述矩阵确定单元包括：
向量建立模块，用于对第k个通道的数据点d建立长度为N的向量m，获得D×K个向量m；所述向量m中的每个元素确定方法为：如果第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k满足则从第n个散射子散射的超声回波信息包括在第k个通道的数据点d中，则向量m中第n个元素设为1，否则，设为0；
构造模块，用于将D×K个向量m构造成一个D×K行、N列的矩阵M。
优选地，所述超声成像装置处理的超声信号为超声平面波信号。
优选地，所述第一时间延迟数据确定单元获得的第一时间延迟数据tn,k的表达式为：
其中，第k个通道的坐标为(xk,0)；所要成的超声图像的第n个散射子的坐标为(xn,zn)；c表示声波在介质中的传播速度。
优选地，所述超声成像装置处理的超声信号为超声凸面波信号、超声凹面波信号。
优选地，所述第二时间延迟数据确定单元获得的第二时间延迟数据td,k的表达式为：
td,k＝t0+(d-1)/fs
其中，t0表示时间延迟值，fs表示采样频率。
优选地，所述求解单元包括：
第一求解模块，用于所述向量I稀疏时，通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号；
第二求解模块，用于所述向量I不稀疏时，对向量I进行稀疏变换；通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，Ψ表示稀疏变换矩阵，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号。
优选地，所述求解单元采用的压缩感知算法为匹配追踪法、Bregman算法、运算符/变量分裂算法、固定点连续算法、L1范数魔术算法、牛顿下降法。
优选地，所述超声成像单元包括：
变换模块，用于将向量I变换为图像矩阵Φ；其中，所述图像矩阵Φ是一个Nx行Nz列的矩阵；NX表示三维坐标系中x方向上图像像素的个数，Nz表示三维坐标系中z方向上图像像素的个数，N＝Nx×Nz；
处理模块，用于对所述矩阵Φ进行处理，得到所要成的超声图像；其中，处理的方法包括：取信号包络、对数压缩、调整图像显示动态范围和进行数字扫描变换。
与现有技术方案相比，本发明实现了一次超声发射/接收便可完成一次二维成像，达到了理论上的最快帧频，同时显著地提高了图像的分辨率和对比度。而且，它既可以应用于超声平面波成像，也具有应用于其他非平面波发射方式的超声成像的潜力，为医学超声成像技术的发展提供了一个新的方向，对我国医学超声影像设备打破国际专利垄断，实现掌握自主知识产权的先进技术突破具有重要意义。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为传统地超声成像原理示意图；
图2为本发明提出的一种超声成像方法流程图；
图3为采用本技术方案对点状的稀疏仿体进行成像的仿真实验结果图；
图4为采用传统的延时叠加法对点状的稀疏仿体进行成像的仿真实验结果图；
图5为采用本技术方案对非稀疏仿体进行成像的仿真结果图；
图6为采用传统的延时叠加法对非稀疏仿体进行成像的仿真结果图；
图7为肾脏组织切片原始图像；
图8为采用本技术方案获得的成像结果图；
图9为采用传统的延时叠加法获得的成像结果图；
图10为本发明提出的一种超声成像装置框图。
具体实施方式
下面将参考若干示例性实施方式来描述本发明的原理和精神。应当理解，给出这些实施方式仅仅是为了使本领域技术人员能够更好地理解进而实现本发明，而并非以任何方式限制本发明的范围。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。
本领域技术技术人员知道，本发明的实施方式可以实现为一种系统、装置、设备、方法或计算机程序产品。因此，本公开可以具体实现为以下形式，即：完全的硬件、完全的软件(包括固件、驻留软件、微代码等)，或者硬件和软件结合的形式。
根据本发明的实施方式，提出了一种超声成像方法及装置。
在本文中，需要理解的是，所涉及的术语中：
1、超声成像：以线阵式超声换能器为例，如图1所示，线阵换能器拥有K个可以独立发射/接收的阵元，对应于K个超声发射通道和信号接收通道。超声发射时，利用多个通道的延时发射，使不同通道的超声信号同时到达聚焦位置，形成发射聚焦；接收回波时，对接收到的信号进行类似的延时，将从同一反射物返回的不同通道接收到的信号累加在一起，形成接收聚焦。这样一次发射和一次接收，可以形成一条扫描线。通常超声成像都采用电子扫描的方式，进行P次聚焦发射/接收来获得P条扫描线，再将这些扫描变换成一幅完整的二维图像。
2、压缩感知(Compressive sensing)，也被称为压缩采样(Compressive sampling)，稀疏采样(Sparse sampling)，压缩传感。它作为一个新的采样理论，它通过开发信号的稀疏特性，在远小于Nyquist采样率的条件下，用随机采样获取信号的离散样本，然后通过非线性重建算法完美的重建信号。
此外，附图中的任何元素数量均用于示例而非限制，以及任何命名都仅用于区分，而不具有任何限制含义。
下面参考本发明的若干代表性实施方式，详细阐释本发明的原理和精神。
发明概述
近年来，国内外先后有一些关于基于压缩感知的超声平面波成像方法的论文发表。这些方法都分为两个步骤：
(1)将图像的每个像素点视为二维平面内的一个网格节点，假设每个网格节点处都存在一个散射子可以造成入射超声的散射，形成回波。则可以认为，我们所要成的超声图像实际反映的是网格节点上散射子的散射强度在二维平面内的分布。
首先，需要建立反映超声回波射频信号s与网格节点上散射子散射强度分布I之间关系的数学模型，形成如下形式的方程组：
s＝MI
其中，矩阵M为编码矩阵(encoding matrix)。但求解该方程组通常情况下是一个不适定问题，即方程数量小于未知数数量，无法求得唯一解。
(2)当I稀疏(sparse)时，即其中的非零元素数量远小于零元素数量时，则可以通过压缩感知方法对上述方程组进行求解：
其中β反映了我们允许多少噪声成分存在。
对于超声成像来说，在第(1)步骤中，如何根据其所遵循的物理原理，建立能够尽可能真实反映s和I之间关系的数学模型，并据此建立便于完成后续迭代计算的编码矩阵M，是决定超声成像质量和成像方法实用性的关键。而第(2)步骤中，方程组求解的具体计算方法已经有很多成熟的数值迭代方法可供选择，不属于本申请阐述的重点。
对于步骤(1)中所述的反映超声回波射频信号与网格节点上散射子散射强度分布之间关系的数学模型，已发表的主要有以下两种：
第一种数学模型：基于待成像介质的可压缩性分布情况的比较复杂的模型
第一种数学模型的最终形式为：
其中，G是一个NelNk×N的编码矩阵，Nel是超声换能器阵列接收回波信号的通道数，Nk指将宽带的超声回波信号分为Nk个离散的波数kl，1≤l≤Nk，N＝Nx×Nz是图像的总的像素数(或者说是网格节点数)，Nx，Nz分别是x方向(宽度方向)和z方向(深度方向)上图像像素的行数和列数。矩阵G中的每个元素定义为：
其中，m表示换能器上的第m个阵元，1≤m≤Nel，i表示图像上的第i个像素，代表入射超声的声压，rel,m表示超声换能器上第m个阵元的空间位置，ri表示图像上第i个像素的位置，gl(rel,m-ri)是开放空间的格林函数，定义为：
其中，j表示虚部，是零阶的第二类Hankel函数。psc代表超声回波射频信号，γκ代表所要成像的介质可压缩性的分布情况(介质的可压缩性是决定其声散射强度的主要因素)。
第二种数学模型：基于频域信号延时的比较简单的模型
第二种数学模型的最终形式为：X(ω)＝A(ω)·S(ω)
由于是基于频域信号进行处理，实际上取ω＝2πf0，其中，f0为所使用的超声换能器的中心发射频率。X为经过短时傅里叶变换后的超声回波射频信号，S为所要成像的散射子散射强度在频域上对应于f0的映射，A为K×L的由时间延迟数据组成的编码矩阵，定义为：
[A(ω)k]i＝exp[jωτk(ρi)]
其中，K为超声换能器阵列接收回波信号的通道数，L为图像的总的像素数(或者说是网格节点数)，1≤k≤K，1≤i≤L，ρi表示图像上的一个像素点(或网格节点)，表示从某个像素点发出的回波信号到达某个超声换能器阵列通道的时间延迟，rk表示某个超声换能器阵列通道的空间位置，表示某个像素点ρi的空间位置。需要指出，X是从全部超声回波射频信号中截取出一小段进行短时傅立叶变换后获得的频域信号，因此若全部超声回波射频信号被分成Q段，则要完成全部成像，需要将后续的求解过程重复Q次。
建立上述两个模型后，通过压缩感知方法进行方程组求解，就可以解出γκ(第一种数学模型)或者S(第二种数学模型)，再将其从向量变换为对应图像像素点数量的矩阵，就可以显示为我们所希望获得的图像。
上述两个反映超声回波射频信号与网格节点上散射子散射强度分布之间关系的数学模型，各有其局限性。
第一种数学模型的建立，是从已经被证明比较准确的声传播和散射的数学模型出发，优点是能够比较真实的反映声在介质中的各种物理现象，但缺点也非常明显，就是模型过于复杂。编码矩阵G的尺寸过于庞大，需要占用大量内存，同时也造成后续的求解过程的计算量十分巨大。以其论文中进行的成像实验数据为例，当Nx＝400，Nz＝600，Nel＝128，Nk＝1000时，矩阵G占用的内存高达458GB。因此，为实现该算法，只好采用每次调用G时都重新计算其各个元素数值的方法来进行，极大的增加了计算量。而且，实际上上述参数值根本无法满足正常医学超声成像的需要，如果成像深度超过5cm，Nz的取值通常都在3000以上，因此内存占用量还要再增加到5倍，已经完全不是普通计算机能够承担的任务。
第一种数学模型由于只考虑超声回波射频信号的时间延迟，同时只考虑超声的中心发射频率f0而不考虑信号其他频率成分，因此其所使用的编码矩阵A(ω)的规模大为缩小。但是，该模型仍有以下几个问题。首先，矩阵A(ω)的所有元素都是非零的，而且X是从全部超声回波射频信号中截取出一小段进行短时傅立叶变换后获得的频域信号，若全部超声回波射频信号被分成Q段，则要完成全部成像，需要将后续的求解过程重复Q次，因此进行后续的矩阵乘法运算时计算量仍然很大。其次，为了方便的对信号进行时间延迟计算，该模型的所有运算都在频域进行。这就需要首先将时域的超声回波射频信号，通过短时傅里叶变换转换到频域。这个过程不仅增加了计算量，也会引入由于信号长度有限产生的泄漏误差和旁瓣误差，进而影响到最终的成像质量。
本技术方案的目的是，相对于第一数学模型，尽可能的简化编码矩阵，降低运算时的内存存储空间和计算量；相对于第二数学模型，避免使用傅立叶变换和频域计算；在保证高帧频和较高成像质量的同时，实现只需要较低水平的硬件计算平台就能使用本方法进行成像，便于实现本技术方案的产业转化。
本技术方案通过超声回波射频信号相对于超声波发射时刻的时间延迟信息与散射子分布的空间位置信息之间的关系，建立在时域表达的编码矩阵M，形成如下的反映超声回波射频数字信号s与网格节点上散射子散射强度分布I之间关系的方程组s＝MI；由于此时编码矩阵M是一个稀疏矩阵，因此可以通过稀疏表达的方式大幅度降低其占用的内存存储空间，同时在利用其运算时的计算量也非常低。最后，通过成熟的压缩感知算法对上述方程组进行求解，得到向量I，将其从向量变换为对应图像像素点数量的矩阵，即为所希望获得的超声图像。
在介绍了本发明的基本原理之后，下面具体介绍本发明的各种非限制性实施方式。
应用场景总览
超声成像是利用超声声束扫描人体，通过对反射信号的接收、处理，以获得体内器官的图象。常用的超声仪器有多种：A型(幅度调制型)是以波幅的高低表示反射信号的强弱，显示的是一种“回声图”。M型(光点扫描型)是以垂直方向代表从浅至深的空间位置，水平方向代表时间，显示为光点在不同时间的运动曲线图。以上两型均为一维显示，应用范围有限。B型(辉度调制型)即超声切面成像仪，简称“B超”。是以亮度不同的光点表示接收信号的强弱，在探头沿水平位置移动时，显示屏上的光点也沿水平方向同步移动，将光点轨迹连成超声声束所扫描的切面图，为二维成像。由于B型超声图象清晰、直观，层次感强，故在临床广为应用。至于D型是根据超声多普勒原理制成.C型则用近似电视的扫描方式，显示出垂直于声束的横切面声象图。近年来，超声成像技术不断发展，如灰阶显示和彩色显示、实时成像、超声全息摄影、穿透式超声成像、超声计并机断层圾影、三维成像、体腔内超声成像等。
超声成像方法常用来判断脏器的位置、大小、形态，确定病灶的范围和物理性质，提供一些腺体组织的解剖图，鉴别胎儿的正常与异常，在眼科、妇产科及心血管系统、消化系统、泌尿系统的应用十分广泛。
本技术方案首先由计算机控制超声发射电路对超声换能器阵列进行激励，发射超声信号。当超声换能器阵列的各个通道(每个通道对应于一个阵元)被同时激励时，发射出超声信号，超声信号在介质中传播，发生散射，形成超声回波信号。超声回波信号被超声换能器阵列接收，形成超声回波射频信号，再由超声接收电路采样，形成超声回波射频数字信号。超声回波射频数字信号被送回到计算机中，并在计算机中实现超声波成像。
示例性方法
接下来，参考图2对本发明示例性实施方式进行介绍。
如图2所示，为本发明提出的一种超声成像方法流程图。该方法包括：
步骤201)：获取N×K个第一时间延迟数据tn,k；其中，所述第一时间延迟数据tn,k为超声信号从超声发射时刻开始，经过第n个散射子的散射，再回到超声换能器阵列的第k个通道位置的总的时间延迟；N表示所要成的超声图像的像素点数，K表示超声换能器阵列中阵元数目；n表示某一个散射子，散射子n＝1，2，……，N；k表示超声换能器阵列中某一个通道，k＝1,2，……，K；
步骤202)：获取D×K个第二时间延迟数据td,k；其中，所述第二时间延迟数据td,k为已知从超声发射时刻开始，经过时间延迟，开始对超声回波射频信号进行采样，得到超声回波射频数字信号，第k个通道接收到所述超声回波射频数字信号的第d个数据点的时刻相对于超声信号发送时刻的时间延迟；D表示每个通道的采样点数，d表示一通道接收到超声回波射频数字信号的某一数据点，d＝1,2，……，D；
步骤203)：将所述超声回波射频数字信号对应地矩阵变换为长度为D×K的向量s；
步骤204)：利用所述向量s中第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k确定矩阵M；其中，所述矩阵M为稀疏矩阵；
步骤205)：根据所述向量s和所述矩阵M建立方程组；其中，所述方程组为：s＝MI；I为表示散射子的散射强度分布的长度为N的向量；
步骤206)：对所述方程组通过压缩感知算法进行求解，得到向量I；
步骤207)：根据所述向量I确定所要成的超声图像。
本方法并不限于只针对于超声平面波成像。例如，如果超声换能器阵列发出的是凸面波或者凹面波，只要根据具体情况，修改上述的计算超声信号从超声发射时刻开始，经过某个散射子的散射，再回到超声换能器阵列的某个阵元位置的第一时间延迟数据tn,k的计算公式即可。其他均相同。
以超声信号为超声平面波信号为例，已知超声换能器阵列包括K个阵元，其中第k个阵元的坐标为(xk,0)。所要成的超声图像的像素点数(即对成像平面划分的网格节点数)为N＝Nx×Nz，其中，Nx、Nz分别是x方向(宽度方向)和z方向(深度方向)上图像像素的行数和列数。某一个网格节点处的散射子n的坐标为(xn,zn)。则有，超声信号从超声发射时刻开始，经过这个散射子的散射，再回到超声换能器阵列的第k个阵元位置，总的第一时间延迟为：
其中，c是声波在介质中的传播速度。对应的，可以得到N×K个第一时间延迟数据tn,k。
已知从超声发射时刻开始，经过t0的时间延迟，开始对超声回波射频信号进行采样，采样频率为fs，每个通道的采样点数为D，则第k个通道接收到超声回波射频信号的第d个数据的时间点，相对于超声发射时刻的时间延迟为：
td,k＝t0+(d-1)/fs
对应的，可以得到D×K个第二时间延迟数据td,k。
对于本实施例来说，确定矩阵M的步骤包括：
对第k个通道的数据点d建立长度为N的向量m，获得D×K个向量m；所述向量m中的每个元素确定方法为：如果第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k满足则从第n个散射子散射的超声回波信息包括在第k个通道的数据点d中，则向量m中第n个元素设为1，否则，设为0；
将D×K个向量m构造成一个D×K行、N列的矩阵M。
对于本实施例来说，所述对所述方程组通过压缩感知算法进行求解的步骤包括：
如果向量I稀疏时，通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号；
否则，对向量I进行稀疏变换；通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，Ψ表示稀疏变换矩阵，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号。
对于本技术方案来说，成熟的压缩感知算法包括但不限于匹配追踪法(matching pursuitmethod)、Bregman算法、运算符/变量分裂算法(operator/variable splitting)、固定点连续算法(Fixed-point continuation)、L1范数魔术算法(L1-magic)、牛顿下降法等。
对于本实施例来说，确定所要成的超声图像的步骤包括：
将向量I变换为图像矩阵Φ；其中，所述图像矩阵Φ是一个NX行Nz列的矩阵；NX表示三维坐标系中x方向上图像像素的个数，Nz表示三维坐标系中z方向上图像像素的个数，N＝Nx×Nz；
对所述矩阵Φ进行处理，得到所要成的超声图像；其中，处理的方法包括：取信号包络、对数压缩、调整图像显示动态范围和进行数字扫描变换。
在本实施例中，利用Field II超声成像仿真软件对本技术方案进行了验证。如图3所示，为采用本技术方案对点状的稀疏仿体进行成像的仿真实验结果图。如图4所示，为采用传统的延时叠加法(DAS)对点状的稀疏仿体进行成像的仿真实验结果图。经过对比可以发现，本技术方案可以去除传统方法中出现的横向伪影。
如图5所示，为采用本技术方案对非稀疏仿体进行成像的仿真结果图。如图6所示，为采用传统的延时叠加法(DAS)对非稀疏仿体进行成像的仿真结果图。经过对比可以发现，本技术方案可以获得更好的对比度，并且可以去除传统方法近场中出现的伪影。
如图7所示，为肾脏组织切片原始图像。如图8所示，为采用本技术方案获得的成像结果图。如图9所示，为采用传统的延时叠加法(DAS)获得的成像结果图。经过对比发现，本技术方案所成的图像，显示肾脏的轮廓、结构更清晰，软组织层次更清楚，组织对比度增加，并且可以去除传统方法近场中出现的伪影。
本技术方案考虑超声回波射频信号相对于超声波发射时刻的时间延迟，不区分其中的不同频率成分，相对于第一数学模型、第二数学模型都大为简化。本技术方案将超声回波射频信号相对于超声波发射时刻的时间延迟信息与散射子分布的空间位置信息联系起来，通过建立更大规模的编码矩阵M，解决了对信号时间延迟计算的时域表达问题，不再需要进行傅立叶变换和在频域进行计算，避免了因进行傅立叶变换而造成的计算误差。
此外，由于本技术方案所建立的编码矩阵M是一个稀疏矩阵，可以通过稀疏表达的方式大幅度降低其占用的内存存储空间，同时在利用其运算时的计算量也非常低。以我们进行的一次仿真成像实验为例，超声回波射频信号的数据点数为128×4364，所成图像像素点数为256×3000，则矩阵M的大小为558592×768000，假设矩阵M中的所有元素都为64bit双精度数，则需要内存存储空间为3423GB。实际上，由于矩阵M是稀疏的，其中的非零元素的数量为98×106，因此当采用稀疏表达时，M实际占用的内存存储空间是1.57GB。相比于第一数学模型来说，如果Nk取1000，则其编码矩阵G需要786GB内存存储空间。可见本技术方案所占用的内存储存空间大为缩减。
相比于第二数学模型，虽然其编码矩阵A只需要786MB的存储空间，但由于需要将全部超声回波射频信号分成Q段，对截取出的每一段数据进行短时傅立叶变换以获得频域信号，并重复后续的求解计算。因此，第二数学模型在总的计算量上仍然大于本技术方案(数据的前后两个分段之间需要有一定的重叠以提高深度方向的分辨率，因此4000多个数据点的长度至少需要被分成100段)。此外，第二数学模型为了方便的对信号进行时间延迟计算，所有运算都在频域进行。而短时傅里叶变换计算不仅增加了计算量，也会引入由于信号长度有限产生的泄漏误差和旁瓣误差，进而影响到最终的成像质量。
应当注意，尽管在附图中以特定顺序描述了本发明方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。
示例性设备
下面参考图10对本发明示例性实施方式的装置进行介绍。
如图10所示，本发明提出的一种超声成像装置框图。该装置包括：
第一时间延迟数据确定单元101，用于获取N×K个第一时间延迟数据tn,k；其中，所述第一时间延迟数据tn,k为超声信号从超声发射时刻开始，经过第n个散射子的散射，再回到超声换能器阵列的第k个通道位置的总的时间延迟；N表示所要成的超声图像的像素点数，K表示超声换能器阵列中阵元数目；n表示某一个散射子，散射子n＝1，2，……，N；k表示超声换能器阵列中某一个通道，k＝1,2，……，K；
第二时间延迟数据确定单元102，用于获取D×K个第二时间延迟数据td,k；其中，所述第二时间延迟数据td,k为已知从超声发射时刻开始，经过时间延迟，开始对超声回波射频信号进行采样，得到超声回波射频数字信号，第k个通道接收到采样后的超声回波射频信号的第d个数据点相对于超声信号发送时刻的时间延迟；D表示每个通道的采样点数，d表示一通道接收到超声回波射频数字信号的某一数据点，d＝1,2，……，D；
变换单元103，用于将所述超声回波射频数字信号对应地矩阵变换为长度为D×K的向量s；
矩阵确定单元104，用于利用所述向量s中第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k确定矩阵M；其中，所述矩阵M为稀疏矩阵；
模型建立单元105，用于根据所述向量s和所述矩阵M建立方程组；其中，所述方程组为：s＝MI；I表示散射子的散射强度分布、长度为N的向量；
求解单元106，用于对所述方程组通过压缩感知算法进行求解，得到向量I；
超声成像单元107，用于根据所述向量I确定所要成的超声图像。
对于本实施例来讲，所述矩阵确定单元104包括：
向量建立模块，用于对第k个通道的数据点d建立长度为N的向量m，获得D×K个向量m；所述向量m中的每个元素确定方法为：如果第k个通道的数据点d对应的第二时间延迟数据td,k以及相同通道下对应的第一时间延迟数据tn,k满足则从第n个散射子散射的超声回波信息包括在第k个通道的数据点d中，则向量m中第n个元素设为1，否则，设为0；
构造模块，用于将D×K个向量m构造成一个D×K行、N列的矩阵M。
对于本实施例来讲，所述超声成像装置处理的超声信号为超声平面波信号。这种情况下，所述第一时间延迟数据确定单元获得的第一时间延迟数据tn,k的表达式为：
其中，第k个通道的坐标为(xk,0)；所要成的超声图像的第n个散射子的坐标为(xn,zn)；c表示声波在介质中的传播速度。
本方法并不限于只针对于超声平面波成像。例如，如果超声换能器阵列发出的是凸面波或者凹面波，只要根据具体情况，修改上述的计算超声信号从超声发射时刻开始，经过某个散射子的散射，再回到超声换能器阵列的某个阵元位置的第一时间延迟数据tn,k的计算公式即可。其他均相同。
所述第二时间延迟数据确定单元获得的第二时间延迟数据td,k的表达式为：
td,k＝t0+(d-1)/fs
其中，t0表示时间延迟值，fs表示采样频率。
对于本实施例来讲，所述求解单元106包括：
第一求解模块，用于所述向量I稀疏时，通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号；
第二求解模块，用于所述向量I不稀疏时，对向量I进行稀疏变换；通过压缩感知算法对方程组s＝MI进行求解，即满足条件的I值为方程组的最优解；其中，Ψ表示稀疏变换矩阵，β表示允许多少噪声成分存在；s.t.是表示逻辑关系的公式符号。
对于本实施例来讲，成熟的压缩感知算法包括但不限于匹配追踪法(matching pursuitmethod)、Bregman算法、运算符/变量分裂算法(operator/variable splitting)、固定点连续算法(Fixed-point continuation)、L1范数魔术算法(L1-magic)、牛顿下降法等。
优选地，所述超声成像单元107包括：
变换模块，用于将向量I变换为图像矩阵Φ；其中，所述图像矩阵Φ是一个NX行Nz列的矩阵；NX表示三维坐标系中x方向上图像像素的个数，Nz表示三维坐标系中z方向上图像像素的个数，N＝Nx×Nz；
处理模块，用于对所述矩阵Φ进行处理，得到所要成的超声图像；其中，处理的方法包括：取信号包络、对数压缩、调整图像显示动态范围和进行数字扫描变换。
有上述示例性方法描述可知，本技术方案一方面可以实现超高帧频的快速超声成像，另一方面保证了较高的成像质量，同时只需要较低水平的硬件运算平台就能实现，便于实现产业转化。
应当注意，尽管在上文详细描述中提及了超声成像装置为一种应用计算机软件代码实现的装置，该装置中包括的若干单元或模块，但是这种划分仅仅并非强制性的。实际上，根据本发明的实施方式，上文描述的两个或更多装置的特征和功能可以在一个系统中具体化。反之，上文描述的一个系统的特征和功能可以进一步划分为由多个装置来具体化。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。