CN104306023A - 基于压缩感知的超声成像快速实现方法 - Google Patents

Abstract

本发明公开了一种基于压缩感知理论的超声成像快速实现方法，主要解决压缩感知框架下快速迭代收缩阈值算法计算迭代步长时间较慢的问题。其实现步骤为：1.根据设定的分辨率，对探测区域离散化；2.对离散化探测区域进行宽带脉冲扫描，得到回波向量和观测矩阵；3.根据回波向量和观测矩阵建立超声成像数学模型；4.通过求观测矩阵每一列元素模值和的二范数确定迭代步长；5.将迭代步长带入快速迭代收缩阈值算法解出重建观测场景散射强度；6.将该散射强度取模值排列成一个二维矩阵，得到重建的超声图像。本发明相对传统的快速迭代收缩阈值算法，计算迭代步长的时间得到了大幅度减少，提高了超声成像的实时性，可用于超声实时处理领域。

Description

基于压缩感知的超声成像快速实现方法

技术领域

本发明属于图像处理技术领域，特别涉及一种快速成像的方法，可用于B超成像。

背景技术

医学超声成像经过60余年的发展，它具有相对安全、实时性好、无创、便携、价格低廉等优点，其与X射线诊断技术、电子计算机断层扫描CT成像技术、核磁共振成像技术一起称为现代医学四大影像技术，已令亿万患者受益。

但是超声成像仍然存在一些不足，如分辨率不高，多为毫米级；受噪声干扰严重，图像质量较差；实时性一般。

近年来，在信号处理领域兴起的压缩感知理论吸引了诸多学者的关注，该理论指出，只要信号在某一个空间Ψ上具有稀疏性，就可以利用观测矩阵对其以远低于奈奎斯特采样速率进行观测，进而利用优化手段高概率地从混叠观测中重建原信号，这使得传感器的采样成本大大降低。而且通过恰当选择空间Ψ，信号的稀疏性越强，精确恢复原信号的可能性就越大，这样就会在提高图像分辨率、抑制噪声方面有出色的表现。从近几年国内外发表的文献来看，对压缩感知理论的研究已经涉及众多领域如压缩感知CS雷达成像、医学图像处理、光谱分析、遥感图像处理等，具有非常广阔的应用前景。

由于病灶区域与正常组织的密度特征有明显差别，可认为超声图像在空间域是稀疏的，将压缩感知理论应用到超声成像可以较好的解决超声成像分辨率不高和噪声干扰严重的问题，但是压缩感知理论的问题在于重建过程中观测矩阵维度巨大、导致计算复杂度高，图像的重建时间过长。

针对这个问题，以色列学者A.Beck等人在论文“A fast iterativeshrinkage-thresholding algorithm for linear inverse problems”SIAMJ.IMAGE SCIENCES,Vol.2.No.1,pp.183-202,2009中提出了快速迭代收缩阈值算法，简称FISTA。利用该算法，超声成像的基本框架可以表示为：

$X^{*}=\underset{x}{\arg \min }\{{\left|\right|\mathrm{\Ψx}-b\left|\right|}_2^2+\mathrm{\λ}{\left|\right|x\left|\right|}_1\}$

其中，X^*为重建观测场景散射强度，x为目标向量，b是采样后超声阵元接收的回波数据，λ为正则化参数，表示向量Ψx-b二范数的平方，||x||₁表示目标向量x的一范数。

FISTA算法复杂度低，适合求解大规模矩阵的重建问题，同时具备全局收敛性，收敛速度快。但是，这个方法的不足是在采用固定步长确定梯度下降序列时，用观测矩阵Ψ最大奇异值的倒数作为迭代步长，当观测矩阵维度较大时，求最大奇异值的时间将比较久，不利于满足超声成像实时性的要求。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于压缩感知的超声成像快速实现方法，以缩短压缩感知理论应用于超声成像时的重建时间。

本发明的技术方案是这样实现的：

一.技术原理

在经典的梯度下降算法中，迭代步长μ的取值范围为μ∈(0,1/||Ψ^TΨ||]，其中||Ψ^TΨ||表示观测矩阵Ψ的最大奇异值，μ的值越接近||Ψ^TΨ||，算法收敛速度越快，因此设计简单易求的迭代步长μ，使μ尽可能的接近||Ψ^TΨ||，将有助于提高超声成像重建的速度。

二.技术方案

根据上述原理，本发明的实现步骤如下：

(1)将超声探测区域二维离散化，得到N个离散化的像素点，其中N＝T×S，T表示轴向像素的个数，S表示侧向像素的个数；

(2)将超声宽带脉冲发射信号在频域均匀采样得到W个频点，按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描，每次扫描得到一个长度为A的局部观测回波向量b_t，并将这W个局部观测回波向量按从上到下顺序合成长度为M＝A×W的观测回波向量b＝{b₁,...,b_t,...,b_W}，同时保存由这W个频点产生的回波声场强度矩阵Ψ₁,...,Ψ_t,...,Ψ_W，其中，矩阵Ψ_t的宽度为A，长度为N，A表示超声线阵的阵元个数，1≤t≤W；

(3)将回波声场强度矩阵Ψ₁,...,Ψ_t,...,Ψ_W，按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ，Ψ的第i行第j列元素表示为Ψ[i，j],1≤i≤M,1≤j≤N；将离散化的二维探测区域按照行优先的顺序排列成一个目标向量x；

(4)根据回波向量b和观测矩阵Ψ定义基于压缩感知的超声成像数学模型：

$X^{*}=\underset{x}{\arg \min }\{{\left|\right|\mathrm{\Ψx}-b\left|\right|}_2^2+\mathrm{\λ}{\left|\right|x\left|\right|}_1\}---1)$

其中X^*为重建观测场景散射强度，λ为正则化参数，表示向量Ψx-b二范数的平方，||x||₁表示目标向量x的一范数；

(5)对上述数学模型进行求解，得到重建观测场景散射强度X^*：

(5a)初始化：n＝0，ε＝10^-3，n表示第n次迭代，ε表示迭代终止条件；

(5b)根据观测矩阵Ψ，构造基于梯度下降算法的迭代步长μ：

μ＝1/norm(sum_x,2)   2)

其中norm()表示向量的二范数；长度为N，表示实数域，它的每一个元素sum_x[j]，1≤j≤N表示观测矩阵Ψ的每一列所有元素模值的和，计算如下：

$\mathrm{sum}\_x\left[j\right]=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{abs}\left(\mathrm{\Ψ}\right[i,j\left]\right)---3)$

其中abs()表示取元素的模值,1≤i≤M,1≤j≤N；

(5c)将回波b、观测矩阵Ψ和迭代步长μ带入快速迭代收缩阈值算法中，经过多次梯度下降和快速阈值收缩过程，直到目标向量满足迭代终止条件，得到重建观测场景散射强度X^*；

(6)对重建场景散射强度X^*取模值，并按照先行后列的顺序排列成一个二维矩阵，即得到重建的超声图像。

本发明由于在计算重建观测场景散射强度的过程中采用简单易求的方法计算迭代步长，相比求最大奇异值的倒数作为迭代步长的方案，在观测矩阵维度超过1GB的时候可以减少约96％的计算时间，提高了超声成像重建的速度。

附图说明

图1是本发明的主流程图；

图2是本发明对探测区域离散化的示意图；

图3是本发明使用的超声阵列与探测区域位置示意图；

图4是本发明将回波数据排列成列向量的示意图；

图5是本发明中进行快速迭代收缩阈值算法的子流程图；

图6是本发明中仿真使用的20个理想点目标的二维场景图。

具体实施方式

以下结合附图对本发明作进一步详细描述。

参照图1，本发明实现步骤如下：

步骤一：对探测区域进行离散化。

本发明对探测区域离散化采用等间隔采样，得到离散化的二维探测区域N＝T×S，如图2所示，其中N表示离散化总像素点数，T表示轴向像素点数，S表示侧向像素点数。

步骤二：对二维离散化的探测区域进行宽带脉冲平面波扫描，得到观测回波向量b和观测矩阵Ψ。

具体实现如下：

(2.1)根据图3所示的超声阵列与探测区域的相对位置建立直角坐标系(x,y)，其中x表示侧向，y表示轴向；

(2.2)在直角坐标系中，将超声阵列固定在轴向坐标为零的位置，即y＝0，且将阵列中心与探测区域中心对齐，超声阵列的长度为(A-1)×d，则第l个阵元的横坐标x_l为：

$x_l=(1-\frac{A+1}{2})\×d,$

其中A为超声线阵的阵元个数，d为相邻两个阵元之间的间隔，1≤l≤A；

(2.3)将超声宽带脉冲发射信号在频域均匀采样得到W个频点，按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描，每次扫描得到一个长度为A的局部观测回波向量b_t，1≤t≤W，并将这W个局部观测回波向量按图4所示从上到下顺序合成长度为M＝A×W的观测回波向量b＝{b₁,...,b_t,...,b_W}；

(2.4)计算探测区域中各个像素点在第m个超声线阵阵元处产生的声场强度

其中A_in(k_t)表示超声宽带脉冲发射信号在频点取值为k_t时的幅度；表示超声宽带脉冲发射信号在频点取值为k_t时离散二维探测区域中各个像素点返回的相位，表示超声宽带脉冲发射信号方位角的单位向量，指定为轴向方向；表示超声宽带脉冲发射信号从超声线阵轴心位置到离散二维探测区域各个像素点距离的向量；表示格林函数，表示超声线阵轴心位置到到超声线阵各个阵元距离的向量，1≤m≤A；

(2.5)由所述声场强度计算出当频点为k_t时，宽度为A，长度为N的回波声场强度矩阵Ψ_t：

其中1≤t≤W,1≤j≤N,1≤m≤A；

(2.6)将回波声场强度矩阵Ψ_t按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ，Ψ的第i行第j列元素表示为Ψ[i，j],1≤i≤M,1≤j≤N；

(2.7)将离散化的二维探测区域按照行优先顺序排列成一个长度为N的目标向量x；

步骤三：定义基于压缩感知的超声成像数学模型为：

其中X^*为重建观测场景散射强度，λ为正则化参数，表示向量Ψx-b二范数的平方，||x||₁表示目标向量x的一范数。

步骤四：对上述数学模型进行求解，得到重建观测场景散射强度X^*。

具体实现如下：

(4.1)初始化：n＝0，ε＝10^-3，n表示第n次迭代，ε表示迭代终止条件；

(4.2)根据观测矩阵Ψ，构造基于梯度下降算法的迭代步长μ：

μ＝1/norm(sum_x,2)，

其中norm()表示向量的二范数；长度为N，它的每一个元素sum_x[j]，1≤j≤N表示观测矩阵Ψ的每一列所有元素模值的和，计算如下：

$\mathrm{sum}\_x\left[j\right]=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\mathrm{abs}\left(\mathrm{\Ψ}\right[i,j\left]\right),$

其中abs()表示取元素的模值，1≤i≤M,1≤j≤N；

(4.3)将回波b、观测矩阵Ψ和迭代步长μ带入快速迭代收缩阈值算法中，进行梯度下降和快速阈值收缩，得到重建观测场景散射强度X^*。

参见图5，本发明步骤具体步骤如下：

(4.3.a)更新梯度下降序列u_n：

u_n＝y_n-μΨ^H(Ψy_n-b)，

其中u_n的长度为N；y_n是快速收缩向量，初始值为0,长度为N；Ψ^H表示观测矩阵Ψ的共轭转置矩阵；

(4.3.b)将梯度下降序列u_n带入下式，得到当前迭代的目标向量x_n：

x_n＝S_Γ(u_n)，

其中S_Γ为阈值函数：

其中Γ为阈值，Γ＝λμ，λ取值在2e⁴～5e⁴之间，e表示科学计数，取值为10；

sign()表示取符号函数；

(4.3.c)判断迭代终止条件||x_n-x_n-1||₂<ε是否成立：

若成立，则停止计算，重建观测场景散射强度X^*＝x_n；

若不成立，令n＝n+1，更新快速收缩向量y_n为：

y_n＝x_n-1+(x_n-2-x_n-1)×(1-t₁)/t₂,

其中当n＝1时，x₀＝0；j表示系数向量x_n和x_n-1中的第j个元素，x_n[j]表示向量x_n的第j个元素的值，x_n-1[j]表示向量x_n-1的第j个元素值；t₁、t₂是两个数值不同的加速因子，t₁初始值为1；将t₁更新为t₁＝t₂，同时返回步骤(4.3.a)。

步骤五：对重建场景散射强度X^*取模值，并按照先行后列的顺序排列成一个二维矩阵，即得到重建的超声图像。

本发明的效果可以通过以下仿真实验说明：

1.仿真条件

仿真测试平台为CPU：Intel Core i7 4770，内存：16GB，windows7 64位操作系统，仿真程序采用VS2010下的C语言实现。仿真一个有128个阵元的超声阵列，阵元间距为300um，二维离散化区域的像素点轴向与侧向间距均为100um，超声宽带脉冲发射信号的带宽为3MHz，中心频率为4.5MHz，采样频点为64个。

2.仿真内容

构造大小为128×128的二维探测场景，如图6所示，它有20个理想点目标，由此产生的观测矩阵Ψ的总数据量为1GB。理想点目标位于探测区域中不同位置，用传统的FISTA算法和本发明中的方法计算迭代步长的时间如表1所示.

表1

从表1可见,本发明计算迭代步长μ的时间，要比传统FISTA算法采用乘幂法来计算最大奇异值，再对最大奇异值取倒数得到迭代步长μ节省96％的时间。

Claims (3)

1.一种基于压缩感知的超声成像快速实现方法，包括以下步骤：

(1)将超声探测区域二维离散化，得到N个离散化的像素点，其中N＝T×S，T表示轴向像素的个数，S表示侧向像素的个数；

(2)将超声宽带脉冲发射信号在频域均匀采样得到W个频点，按频点顺序依次对已离散化的二维探测区域进行一次平面波扫描，每次扫描得到一个长度为A的局部观测回波向量b_t，并将这W个局部观测回波向量按从上到下顺序合成长度为M＝A×W的观测回波向量b＝{b₁,...,b_t,...,b_W}，同时保存由这W个频点产生的回波声场强度矩阵Ψ₁,...,Ψ_t,...,Ψ_W，其中，矩阵Ψ_t的宽度为A，长度为N，A表示超声线阵的阵元个数，1≤t≤W；

(3)将回波声场强度矩阵Ψ₁,...,Ψ_t,...,Ψ_W，按照从上到下顺序排列成一个大小为M×N的观测矩阵Ψ，Ψ的第i行第j列元素表示为Ψ[i，j],1≤i≤M,1≤j≤N；将离散化的二维探测区域按照行优先的顺序排列成一个目标向量x；

(4)根据回波向量b和观测矩阵Ψ定义基于压缩感知的超声成像数学模型：

$X^{*}=\underset{x}{\arg \min }\{{\left|\right|\mathrm{\&Psi;x}-b\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|x\left|\right|}_1\}---1)$

其中X^*为重建观测场景散射强度，λ为正则化参数，表示向量Ψx-b二范数的平方，||x||₁表示目标向量x的一范数；

(5)对上述数学模型进行求解，得到重建观测场景散射强度X^*：

(5a)初始化：n＝0，ε＝10^-3，n表示第n次迭代，ε表示迭代终止条件；

(5b)根据观测矩阵Ψ，构造基于梯度下降算法的迭代步长μ：

μ＝1/norm(sum_x,2)   2)

其中norm()表示向量的二范数；长度为N，表示实数域，它的每一个元素sum_x[j]，1≤j≤N表示观测矩阵Ψ的每一列所有元素模值的和，计算如下：

$\mathrm{sum}\_x\left[j\right]=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\mathrm{abs}\left(\mathrm{\&Psi;}\right[i,j\left]\right)---3)$

其中abs()表示取元素的模值,1≤i≤M,1≤j≤N；

(5c)将回波b、观测矩阵Ψ和迭代步长μ带入快速迭代收缩阈值算法中，经过多次梯度下降和快速阈值收缩过程，直到目标向量满足迭代终止条件，得到重建观测场景散射强度X^*；

(6)对重建场景散射强度X^*取模值，并按照先行后列的顺序排列成一个二维矩阵，即得到重建的超声图像。

2.如权利要求1所述的压缩感知的超声成像快速实现方法，其特征在于，所述步骤2中，频点取k_t时的回波声场强度矩阵Ψ_t，通过下式计算：

${\mathrm{\&Psi;}}_t=\left\{\begin{array}{c}{k}_t^2{A}_{\mathrm{in}}\left(k_t\right)[\left(e^{{-\mathrm{jk}}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_1}\right)\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,l}-{\stackrel{\&RightArrow;}{r}}_1),...,e^{-{\mathrm{jk}}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_j}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,l}-{\stackrel{\&RightArrow;}{r}}_j),...,e^{-j{k}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_N}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,l}-{\stackrel{\&RightArrow;}{r}}_N)]\\ ...\\ k_t^2{A}_{\mathrm{in}}\left(k_t\right)[e^{-{\mathrm{jk}}_i{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_1}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,m}-{\stackrel{\&RightArrow;}{r}}_1),...,e^{-j{k}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_j}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,m}-{\stackrel{\&RightArrow;}{r}}_j),...,e^{-j{k}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_N}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,m}-{\stackrel{\&RightArrow;}{r}}_N)]\\ ...\\ k_t^2{A}_{\mathrm{in}}\left(k_t\right)[e^{-j{k}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_1}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,A}-{\stackrel{\&RightArrow;}{r}}_1),...,e^{-{\mathrm{jk}}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_1}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,A}-{\stackrel{\&RightArrow;}{r}}_j),...e^{-{\mathrm{jk}}_t{\stackrel{\&RightArrow;}{e}}_{\mathrm{\&theta;}}\&CenterDot;{\stackrel{\&RightArrow;}{r}}_N}\&CenterDot;G({\stackrel{\&RightArrow;}{r}}_{A,A}-{\stackrel{\&RightArrow;}{r}}_N)]\end{array}\right.---4)$

其中，A_in(k_t)表示超声宽带脉冲发射信号在频点取值为k_t时的幅度；表示超声宽带脉冲发射信号在频点取值为kt时离散二维探测区域中各个像素点返回的相位，表示超声宽带脉冲发射信号方位角的单位向量，指定为轴向方向；表示超声宽带脉冲发射信号从超声线阵轴心位置到离散二维探测区域各个像素点距离的向量；

表示格林函数，表示超声线阵轴心位置到到超声线阵各个阵元距离的向量,1≤t≤W,1≤j≤N,1≤m≤A。

3.如权利要求1所述的压缩感知的超声成像快速实现方法，其特征在于，所述步骤5中，将回波b、观测矩阵Ψ和迭代步长μ带入快速阈值迭代收缩算法中，经过多次梯度下降和快速阈值收缩过程，其步骤如下：

(5.1)更新梯度下降序列u_n：

u_n＝y_n-μΨ^H(Ψy_n-b)   5)

其中u_n的长度为N；y_n是快速收缩向量，初始值为0,长度为N；Ψ^H表示观测矩阵Ψ的共轭转置矩阵；

(5.2)将梯度下降序列u_n带入下式，得到当前迭代的目标向量x_n：

x_n＝S_Γ(u_n)   6)

其中S_Γ为阈值函数：

其中Γ为阈值，Γ＝λμ，λ取值在2e⁴～5e⁴之间，e表示科学计数，取值为10；sign()表示取符号函数；

(5.3)判断迭代终止条件||x_n-x_n-1||₂<ε是否成立：

若成立，则停止计算，重建观测场景散射强度X^*＝x_n；

若不成立，令n＝n+1，更新快速收缩向量y_n为：

y_n＝x_n-1+(x_n-2-x_n-1)×(1-t₁)/t₂   7)

其中当n＝1时，x₀＝0；j表示系数向量x_n和x_n-1中的第j个元素，x_n[j]表示向量x_n的第j个元素的值，x_n-1[j]表示向量x_n-1的第j个元素值；t₁、t₂是两个数值不同的加速因子，t₁初始值为1，将t₁更新为t₁＝t₂，同时返回步骤(5.1)。