CN111835362B - 一种基于正交基线性表示测量矩阵的压缩感知超声成像方法 - Google Patents

Abstract

本发明涉及一种基于正交基线性表示测量矩阵的压缩感知超声成像方法，属于超声成像技术领域。该方法包括：对超声阵列接收到的回波信号进行处理，得到所需要的超声回波信号x；构造测量矩阵，即正交基线性表示矩阵OBLR对超声回波信号进行压缩采样，得到测量信号y；选取离散余弦变换DCT作为稀疏字典Ψ，对超声回波信号x进行稀疏表示；通过计算得到正交基线性表示矩阵OBLR和稀疏字典Ψ之间的相干系数μ；利用重构算法求解最优化问题，恢复出原始超声信号利用原始超声信号进行波束合成并最终成像；本发明能够显著提高超声信号测量效率，能够对超声信号进行快速、近似最优的重构。

Description

一种基于正交基线性表示测量矩阵的压缩感知超声成像方法

技术领域

本发明属于超声成像技术领域，涉及一种基于正交基线性表示测量矩阵的压缩感知超声成像方法。

背景技术

随着对超声成像质量要求的提高，必然要求更高的采样频率，因此需要存储的回波数据量十分巨大，增加了硬件实现的复杂度。压缩感知理论(CS)是近年来针对高速数据采集与大容量数据存储而提出的一种办法，该理论认为当信号本身或在某个变换域上是稀疏的，就能以远低于奈奎斯特的速率对信号采样的同时进行压缩，再通过重构算法便可以从少量采样数据中以极高的精度重建原始信号，减少需要存储的数据量，降低硬件实现复杂度。尽管压缩感知在理论上取得了巨大的成功，但是CS在感知超声信号方面仍然面临着一些实际的挑战，包括设计简单而有效的测量矩阵。

在压缩感知理论中，测量矩阵性能的优劣直接关系到信号重建精度的高低。测量矩阵性能越好，需要的采样数越少，采样速率越快，重建误差也越小。目前的测量矩阵主要分为随机性测量矩阵和确定性测量矩阵。在压缩感知中常用的随机性测量矩阵包括高斯随机测量矩阵、伯努利随机测量矩阵等。该类矩阵中的每个元素都服从相互独立的同分布，保证了各列向量之间较好的非相关性，重构精度较高。但是存储空间和时间复杂度较大，并且由于其非结构化的本质导致其计算复杂，不利于硬件实现。常用的确定性测量矩阵包括部分哈达玛矩阵，拓普利兹和轮换矩阵等。该类矩阵结构简单，重构效果较好，但是当测量次数较小时，仍要先构造的高维矩阵，再选取行，浪费了存储资源，未能得到广泛应用。特别地，将这些测量矩阵应用到具有重叠性的超声信号时，图像的重构效果很差。

综上所述，现亟需发明一种简单又高效、测量效率高的测量矩阵对超声信号进行压缩采样，提高超声图像重构质量。

发明内容

有鉴于此，本发明的目的在于提供一种基于正交基线性表示测量矩阵的压缩感知超声成像方法。所利用的正交基线性表示测量矩阵与常用稀疏矩阵不相关，并且能够提供快速、近似最优的重构，从而保证重构超声图像质量，在低采样率下的重构超声图像质量优于高斯随机测量矩阵以及以拓普利兹矩阵为代表的确定性测量矩阵。

为达到上述目的，本发明提供如下技术方案：

一种基于正交基线性表示测量矩阵的压缩感知超声成像方法，该方法包括以下步骤：

S1：对超声阵列接收到的回波信号进行处理，得到所需要的超声回波信号x；

S2：构造测量矩阵，即正交基线性表示矩阵OBLR对超声回波信号进行压缩采样，得到测量信号y；

S3：选取离散余弦变换DCT作为稀疏字典Ψ，对超声回波信号x进行稀疏表示；

S4：通过计算得到正交基线性表示矩阵OBLR和稀疏字典Ψ之间的相干系数μ；

S5：利用重构算法求解最优化问题，恢复出原始超声信号

S6：利用原始超声回波信号进行波束合成并最终成像。

可选的，所述步骤S2具体包括以下步骤：

S21：根据信号长度n，以及采样率p，确定测量个数m＝n×p,然后使用随机的方法生成一个全不为零元素的列向量c_i∈±1；

S22：利用列向量C生成m×m维对角阵D＝diag(C),其结构如式表示为：

其中，diag(·)表示对角矩阵，D_m×m表示其大小为m×m维；

S23：采用高斯随机均匀分布在[-1,1]区间上的随机数，生成m×(n-m)维线性表示系数矩阵充当标准正交基系数,结合正交基线性表示方法和正交矩阵D线性表示测量矩阵的剩余n-m个列向量，然后将D和线性表出的n-m列向量拼合构成测量矩阵Φ_m×n,其结构式表示为：

其中：Φ_m×n表示其大小为m×n维；a_1,1,a_1,2,…,a_1,n-m；a_1,1,a_2,2,…,a_2,n-m；a_n,1,a_n,2,…,a_n,n-m为采用高斯随机均匀分布在[-1,1]区间上生成的随机数；

S24：对测量矩阵Φ_m×n中所有列向量进行列归一化处理及相关性优化得到测量矩阵Φ；

S25：用正交基线性表示矩阵OBLR对超声回波信号x进行压缩测量，得到测量信号为：y＝Φx。

可选的，所述步骤S3具体包括以下步骤：

S31：选取的稀疏字典Ψ为离散余弦变换DCT，其表达式为：

其中，k为第k个离散点，x(n)表示超声回波信号x的离散信号，大小为N维，X_c(0)和X_c(k)为变换之后的信号；

S32：对超声回波信号进行稀疏表示为：

x＝Ψα

其中，为n×n维稀疏矩阵，/>是n×1维稀疏系数向量。

可选的，所述步骤S4具体包括：OBLR矩阵和稀疏字典Ψ之间的相干系数μ表达式为：

其中，max(·)表示求取最大元素，Φ_i和Ψ_j分别表示测量矩阵Φ的第i行向量和稀疏矩阵Ψ的第j列向量，其中1≤i≤m,1≤j≤n；|<Φ_i,Ψ_j>|表示计算向量Φ_i和向量Ψ_j的内积的绝对值，||Φ_i||₂表示计算Φ_i的l₂范数，||Ψ_j||₂表示计算Ψ_j的l₂范数。

可选的，所述步骤S5具体包括以下步骤：

S51：通过测量信号y、稀疏字典Ψ以及测量矩阵Φ得到：

y＝Φx＝ΦΨα＝Θα

其中，Θ＝ΦΨ表示为感知矩阵；

S52：计算稀疏系数向量α的逼近值即通过l₁范数最小法求解以下最优化问题：

其中，min(·)表示求取最小元素，表示/>的l₁范数；

S53：通过逼近值恢复出原始超声信号/>

可选的，所述步骤S6具体：利用原始超声信号进行波束合成，计算得到波束信号：

其中，s_DAS表示得到的波束信号，表示第i个阵元上的重建原始回波信号，N₁为超声阵列总数。

本发明的有益效果在于：本发明提供了一种在压缩感知超声成像中测量矩阵的设计方法，该OBLR矩阵由两部分构造而成，第一部分是m×m维的标准正交基矩阵，该矩阵只有对角线上的元素非零，且该元素的值为随机的“1”或“-1”。第二部分是由高斯随机均匀分布在[-1,1]区间上的随机数所生成的m×(n-m)维线性表示系数矩阵，用该矩阵充当标准正交基系数，由于高斯随机数全为非零实数，所以该矩阵的任意一列都只能被标准正交基矩阵的所有列唯一线性表示，保证了测量矩阵最大的近似非相关性。然后对该测量矩阵的所有列进行列归一化及相关性优化得到最终的测量矩阵。与传统的测量矩阵相比，所提出的基于正交基线性表示测量矩阵在一定程度上比它们稀疏，且构造方法简单。本发明采用的OBLR矩阵在低采样率下的重构超声图像质量优于高斯随机测量矩阵以及以拓普利兹矩阵为代表的确定性测量矩阵。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：

图1为本发明所述的基于正交基线性表示测量矩阵的压缩感知超声成像方法的流程图；

图2不同压缩率下稀疏矩阵和三种测量矩阵的相干系数图；

图3为不同测量矩阵和采样数据时的点目标重构图像；

图4为不同测量矩阵采样50％数据时的吸声斑重构图像；

图5为矩阵OBLR采样不同数据量时的吸声斑重构图像；

图6为不同测量矩阵采样50％数据时的geabr_0重构图像；

图7为不同测量矩阵在60mm处横向分辨率曲线对比图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

本发明实施例的附图中相同或相似的标号对应相同或相似的部件；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”、“前”、“后”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本发明的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

请参阅图1～图7，图1为本发明的算法流程图，如图1所示，本发明提供一种在压缩感知超声成像中测量矩阵的设计方法，包括以下步骤：

步骤1：对超声阵列接收到的超声信号进行处理，得到所需要的超声回波信号x。

步骤2：构造测量矩阵，即正交基线性表示矩阵OBLR对超声回波信号进行压缩采样，得到测量信号y，具体包括以下步骤：

S21：根据信号长度n，以及采样率p，确定测量个数m＝n×p,然后使用随机的方法生成一个全不为零元素的列向量c_i∈±1。

S22：利用列向量C生成m×m维对角阵D＝diag(C),其结构如式表示为:

其中，diag(·)表示对角矩阵，D_m×m表示其大小为m×m维。

S23：采用高斯随机均匀分布在[-1,1]区间上的随机数，生成m×(n-m)维线性表示系数矩阵充当标准正交基系数,结合正交基线性表示方法和正交矩阵D线性表示测量矩阵的剩余n-m个列向量,然后将D和线性表出的n-m列向量拼合构成测量矩阵Φ_m×n,其结构式表示为:

其中:Φ_m×n表示其大小为m×n维。a_1,1,a_1,2,…,a_1,n-m；a_1,1,a_2,2,…,a_2,n-m；a_n,1,a_n,2,…,a_n,n-m为采用高斯随机均匀分布在[-1,1]区间上生成的随机数。

S24：对测量矩阵Φ_m×n中所有列向量进行列归一化处理及相关性优化得到测量矩阵Φ。

S25：用测量矩阵OBLR对超声回波信号x进行压缩测量，得到测量信号为：y＝Φx。

步骤3：选取离散余弦变换DCT作为稀疏字典Ψ，对超声回波信号x进行稀疏表示，具体包括以下步骤：

S31：选取的稀疏字典Ψ为离散余弦变换DCT，其表达式为：

其中，k为第k个离散点，x(n)表示超声回波信号x的离散信号，大小为N维，X_c(0)和X_c(k)为变换之后的信号；

S32：对超声回波信号进行稀疏表示为：

x＝Ψα

其中，为n×n维稀疏矩阵，/>是n×1维稀疏系数向量。

步骤4：OBLR矩阵和稀疏字典DCT之间的相干性使用相干系数μ进行数学量化，相干系数μ表达式为：

其中，max(·)表示求取最大元素，Φ_i和Ψ_j分别表示测量矩阵Φ的第i行向量和稀疏矩阵Ψ的第j列向量，其中1≤i≤m,1≤j≤n。|<Φ_i,Ψ_j>|表示计算向量Φ_i和向量Ψ_j的内积的绝对值，||Φ_i||₂表示计算Φ_i的l₂范数，||Ψ_j||₂表示计算Ψ_j的l₂范数。

步骤5：利用重构算法求解最优化问题，恢复出原始超声信号具体包括以下步骤：

S51：通过测量信号y、稀疏字典Ψ以及测量矩阵Φ得到：

y＝Φx＝ΦΨα＝Θα

其中，Θ＝ΦΨ表示为感知矩阵；

S52：计算稀疏系数向量α的逼近值即通过l₁范数最小法求解以下最优化问题：

其中，min(·)表示求取最小元素，表示/>的l₁范数。

S53：通过逼近值恢复出原始超声信号/>

S6：利用原始超声信号进行波束合成并最终成像，计算得到波束信号：

其中，s_DAS表示得到的波束信号，表示第i个阵元上的重建原始回波信号，N₁为超声阵列总数。

FieldII是丹麦理工大学基于声学原理开发的一款超声实验仿真平台，其在理论研究上获得了广泛的认可和使用。为验证所提算法的有效性，利用FieldII对超声成像中常用的点散射目标和吸声斑目标进行成像并利用实际实验数据进行成像对比实验。在点目标仿真实验中，每隔10mm设置了8个散射点目标，均匀分布在30mm～100mm的深度之间，采用发射定点聚焦和接收动态聚焦方式，发射焦点固定在60mm处，设置图像的成像动态范围为50dB。同时，设置了三个半径分别为3mm、2.5mm、2mm的散射暗斑在左边区域，三个半径分别为2mm、2.5mm、3mm的散射亮斑在右边区域，以及3个散射点目标中间位置，外部随机分布100000个散射点。亮斑区域的散射强度是背景区域的10倍，暗斑区域的散射强度定义为零，并设定成像动态范围为50dB。实验所采用的阵元中心频率为3.33MHz，阵元数目为64个，间距为0.2413mm，采样频率为17.76MHz，声速为1500m/s，设成像动态范围为50dB。对上述三个实验目标采用正交基线性表示矩阵OBLR，高斯随机矩阵(Φ_Gaussian)以及拓普利兹矩阵(Φ_Toeplitz)进行对比成像实验。同时，从均方误差和峰值信噪比来评价超声图像重构质量，判断不同测量矩阵的优劣以及重构图像差异。

图2给出了不同采样率下稀疏矩阵Φ_DCT和三种测量矩阵的相关系数μ。由图2可知，矩阵Φ_Toeplitz和稀疏矩阵Ψ_DCT的不相关性最差，这是由于其结构具有一定的确定性且不稀疏。矩阵Φ_Gaussian和稀疏矩阵Ψ_DCT的不相关性较好，这是因为Φ_Gaussian具有随机性，但是存储空间较大。在较低的数据压缩率下μ(Φ_Gaussian,Ψ_DCT)和μ(Φ_OBLR,Ψ_DCT)的值差不多，但在较高数据采样率下μ(Φ_OBLR,Ψ_DCT)要小于μ(Φ_Gaussian,Ψ_DCT)。因此，本文所提矩阵Φ_OBLR不仅更加稀疏，所需存储空间较小，而且和稀疏矩阵Ψ_DCT具有较好的不相关性，能够实现超声信号的重构。

图3给出了不同测量矩阵和采样数据量的点目标重构图像。表1为采集50％数据量时，在三种测量矩阵下重构点目标图像均方误差和峰值信噪比计算值。由图3和表1可知，在采样50％数据量的情况下，三种测量矩阵均能以高质量恢复出原始图像，但是在本文提出的测量矩阵Φ_OBLR下重构出的超声图像均方误差最小，峰值信噪比最大。然而，当采样数据量为30％时，在矩阵Φ_Toeplitz下产生了大量的纵向伪影，在矩阵Φ_Gaussian下也产生较多的纵向伪影。只有矩阵Φ_OBLR能够较准确的恢复出原始超声图像。最后，当采样数据量为20％时，在矩阵Φ_Gaussian和Φ_Toeplitz下产生了大量伪像，以至于分辨不出目标点位置。在矩阵Φ_OBLR下重构出的点目标图像分辨率有所下降，目标点变暗，产生了较少的纵向伪影，但是仍然能够辨别出八个目标点。因此，在不同采样数据量下，通过本文提出的测量矩阵Φ_OBLR重构出的超声图像质量最佳，误差最小。同时，在相同测量矩阵下，采样数据量越大时，图像重构质量也越好，但是重构复杂度也相应增加。

表1三种测量矩阵采样50％数据时的点目标重构图像均方误差和峰值信噪比

测量矩阵	Toeplitz	Gaussian	OBLR
				MSE	0.0095	0.0087	0.0021
PSNR(dB)	68.35	68.74	74.91

图4给出了三种测量矩阵采样50％数据量时的吸声斑重构图像，对应的均方误差和峰值信噪比在表2中列出。由图4可知，在采样50％数据量的情况下，通过矩阵Φ_Toeplitz重构出的吸声斑图像分辨率较低，亮斑处出现大量伪影。在矩阵Φ_Gaussian和Φ_OBLR下的重构图像质量很好，并且图像对比度还有所提升。结合表2可以更加清楚地展现出矩阵Φ_OBLR的优越性。在矩阵Φ_OBLR采样下重构出的超声图像均方误差最小，峰值信噪比最大。

表2三种测量矩阵采样50％数据时的吸声斑重构图像均方误差和峰值信噪比

测量矩阵	Toeplitz	Gaussian	OBLR
				MSE	0.0232	0.0091	0.0024
PSNR(dB)	64.47	68.54	74.33

图5给出了测量矩阵Φ_OBLR采样不同数据量时的吸声斑重构图像。表3列出了吸声斑重构图像均方误差和峰值信噪比。由图5可知，采样20％数据量时重构出的吸声斑图像失真较严重，左下方暗斑很难清楚辨别出。采样30％数据量时能够较准确地重构出原始图像。而采样50％数据量时则能够非常准确地重构出原始图像。同时，由表3可以看出，采样数据量越多，重构出的图像均方误差越小，峰值信噪比越大。

表3矩阵OBLR采样不同数据量时的吸声斑重构图像均方误差和峰值信噪比

数据采样量	20％	30％	50％
				MSE	0.0171	0.0067	0.0024
PSNR(dB)	65.80	69.89	74.33

图6给出了三种测量矩阵采样50％数据量时geabr_0的重构图像，成像动态范围为50dB。同时，表4给出了相应的重构图像均方误差和峰值信噪比。由图6可知，在矩阵Φ_OBLR采样数据下重构出的图像最接近原始图像，图像质量较好，基本能够分辨出所有目标点、亮斑以及暗斑。在矩阵Φ_Gaussian采样数据下重构出的图像出现了失真，暗斑的恢复质量不佳，图像分辨率也下降了。在矩阵Φ_Toeplitz采样数据下重构出的图像质量最差，所有暗斑基本无法分辨出，亮斑处也引入了大量噪声，远场区域的目标点也无法成像。表4清楚地反映出在矩阵Φ_OBLR下重构出的图像均方误差最小，峰值信噪比最大。

表4三种测量矩阵采样50％数据时的geabr_0重构图像均方误差和峰值信噪比

测量矩阵	Toeplitz	Gaussian	OBLR
				MSE	0.0586	0.0271	0.0203
PSNR(dB)	60.45	63.80	65.06

为了更加直观地比较重构图像质量，作出三种测量矩阵在60mm处横向分辨率曲线如图7所示。从图7可以看出，在矩阵Φ_OBLR下的横向分辨率曲线和原始图像分辨率曲线重合度最高，重构图像分辨率也最好。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种基于正交基线性表示测量矩阵的压缩感知超声成像方法，其特征在于：该方法包括以下步骤：

S1：对超声阵列接收到的回波信号进行处理，得到所需要的超声回波信号x；

S2：构造测量矩阵，即正交基线性表示矩阵OBLR对超声回波信号进行压缩采样，得到测量信号y；

S3：选取离散余弦变换DCT作为稀疏字典Ψ，对超声回波信号x进行稀疏表示；

S4：通过计算得到正交基线性表示矩阵OBLR和稀疏字典Ψ之间的相干系数μ；

S5：利用重构算法求解最优化问题，恢复出原始超声信号

S6：利用原始超声回波信号进行波束合成并最终成像；

所述步骤S2具体包括以下步骤：

S21：根据信号长度n，以及采样率p，确定测量个数m＝n×p,然后使用随机的方法生成一个全不为零元素的m×1维列向量c_i∈±1；

S22：利用列向量C生成m×m维对角阵D＝diag(C),其结构如式表示为：

其中，diag(·)表示对角矩阵，D_m×m表示其大小为m×m维；

S23：采用高斯随机均匀分布在[-1,1]区间上的随机数，生成m×(n-m)维线性表示系数矩阵充当标准正交基系数,结合正交基线性表示方法和正交矩阵D线性表示测量矩阵的剩余n-m个列向量，然后将D和线性表出的n-m列向量拼合构成测量矩阵Φ_m×n,其结构式表示为：

其中：Φ_m×n表示其大小为m×n维；a_1,1,a_1,2,…,a_1,n-m；a_2,1,a_2,2,…,a_2,n-m；a_n,1,a_n,2,…,a_n,n-m为采用高斯随机均匀分布在[-1,1]区间上生成的随机数；

S24：对测量矩阵Φ_m×n中所有列向量进行列归一化处理及相关性优化得到测量矩阵Φ；

S25：用正交基线性表示矩阵OBLR对超声回波信号x进行压缩测量，得到测量信号为：y＝Φx。

2.根据权利要求1所述的一种基于正交基线性表示测量矩阵的压缩感知超声成像方法，其特征在于：所述步骤S3具体包括以下步骤：

S31：选取的稀疏字典Ψ为离散余弦变换DCT，其表达式为：

其中，k为第k个离散点，x(n)表示超声回波信号x的离散信号，大小为N维，X_c(0)和X_c(k)为变换之后的信号；

S32：对超声回波信号进行稀疏表示为：

x＝Ψα

其中，为n×n维稀疏矩阵，/>是n×1维稀疏系数向量。

3.根据权利要求1所述的一种基于正交基线性表示测量矩阵的压缩感知超声成像方法，其特征在于：所述步骤S4具体包括：OBLR矩阵和稀疏字典Ψ之间的相干系数μ表达式为：

其中，max(g)表示求取最大元素，Φ_i和Ψ_j分别表示测量矩阵Φ的第i行向量和稀疏矩阵Ψ的第j列向量，其中1≤i≤m,1≤j≤n；|<Φ_i,Ψ_j>|表示计算向量Φ_i和向量Ψ_j的内积的绝对值，||Φ_i||₂表示计算Φ_i的l₂范数，||Ψ_j||₂表示计算Ψ_j的l₂范数。

4.根据权利要求1所述的一种基于正交基线性表示测量矩阵的压缩感知超声成像方法，其特征在于：所述步骤S5具体包括以下步骤：

S51：通过测量信号y、稀疏字典Ψ以及测量矩阵Φ得到：

y＝Φx＝ΦΨα＝Θα

其中，Θ＝ΦΨ表示为感知矩阵；

S52：计算稀疏系数向量α的逼近值即通过l₁范数最小法求解以下最优化问题：

其中，min(·)表示求取最小元素，表示/>的l₁范数；

S53：通过逼近值恢复出原始超声信号/>

5.根据权利要求1所述的一种基于正交基线性表示测量矩阵的压缩感知超声成像方法，其特征在于：所述步骤S6具体：利用原始超声信号进行波束合成，计算得到波束信号：

其中，s_DAS表示得到的波束信号，表示第i个阵元上的重建原始回波信号，N₁为超声阵列总数。