CN104318619A - 面向无损检测的自适应压缩感知的重建方法 - Google Patents

Abstract

本发明提出一种面向无损检测的自适应压缩感知的重建方法，将高维几何关系矩阵M进行降维，降维后的矩阵大大减少了内存存储空间，结合压缩感知理论，在利用l₁范数最优化方法求解方程的迭代过程中加快计算速度。和现有重建方法比较具有精度高，内存占用少，对硬件配置要求不高的特点。

Description

面向无损检测的自适应压缩感知的重建方法

技术领域

本发明涉及无损检测成像技术领域，具体为一种面向无损检测的自适应压缩感知的重建方法，用基于模型的重建算法结合压缩感知理论和自适应理论进行图像的重建。

背景技术

无损检测技术是指：利用物质的声、光、磁和电等特性，在不损害或不影响被检测对象使用性能的前提下，检测被检对象中是否存在缺陷或不均匀性，给出缺陷大小，位置，性质和数量等信息。它与破坏性检测想比，无损检测具有非破坏性、全面性、全程性等优点。近年来，无损检测技术快速发展，基于不同物理特性的检测技术相继被提出和应用。本专利主要针对无损检测技术中的一种基于光电效应的光声层析成像技术展开研究。

光声层析成像技术中涉及的重建算法的研究在最近几年获得快速发展。最简单的重建方法是延迟迭代算法(DAS)，这是一种在超声波成像中常用的波束成像技术，该技术中引入相关因子和最小方差方法可以进一步提高重建质量；滤波反投影重建算法(FBP)因为其快速高效的重建速度成为在光声成像技术中使用较多的一种重建算法，该算法基于逆向Radon变换并且以目标成像物体位于检测区域的中心和目标物体远离超声波探测器为两个重要前提；上述两种方法方便容易计算，当探测器收集到足够的超声波信号时能够恢复出高质量、高信噪比的图像。但是，如果收集获得的超声波信号不足，重建的图像很容易出现假象。除了上述这两种近似重建方法外，一些解析重建方法相继被提出：如通用反投影算法(UBP)和傅里叶域重建算法，这两种算法为了能够重建出精确的图像也需要超声波信号的充分获取，在某些特殊情况下同样存在着数据扫描不足的问题。为了克服这样的问题基于模型的重建算法(迭代算法)被提出，该算法将重建原始图像的过程转化为一个求解最优化问题，并且能够解决非理想的物理条件和测量环境下的图像重建问题，因此，由超声波的不均匀性和衰减所导致的重建问题也可以得到很好的解决。

可以看到在数据不足，测量环境不理想的情况下，基于模型的重建算法具有其他传统重建算法无可比拟的高质量图像重建效果。但是该方法也有其自身缺点：在重建过程产生几何观测矩阵时，需要计算机相当庞大的内存消耗(特别是对于3D图像重建)，以及在求解目标函数时几十甚至上百次的迭代运算过程中需要相当长的计算时间。有些学者从计算机数据存储结构以及利用计算机图形加速处理器(GPU)加快运算，但是这需要昂贵的硬件支持。

发明内容

为了克服现有技术的不足，本发明提供一种基于压缩感知理论和自适应建模的迭代重建算法，不需要昂贵的硬件支持，基于重建算法本身，在原有的高精度重建基础上，快速的进行图像重建。

本发明解决其技术问题所采用的技术方案为：

所述一种面向无损检测的自适应压缩感知的重建方法，其特征在于：包括以下步骤：

步骤1：根据所获得的超声波信号利用解析法—通用反投影重建算法进行图像的初步重建；

步骤2：将步骤1得到的初步重建结果作为先验模型，依次采用图像特征映射提取算法、自适应网格节点生成算法以及图像的Delaunay三角化方法，产生自适应的三角网格；

步骤3：根据步骤2得到的三角形网格结合三角线性插值方法建立正向模型P＝MH，其中P指换能器检测获得的超声波能量信号矩阵，矩阵M指初始声波压力与探测信号之间的几何关系矩阵，H指自适应三角网格节点的像素值；

步骤4：利用步骤3建立的正向模型，根据换能器检测获得的超声波能量信号矩阵P和初始声波压力与探测信号之间的几何关系矩阵M求解自适应三角网格节点的像素值H：

将正向模型P＝MH求解问题转换为优化问题，其中P_d为换能器获得的超声波信号，H^*＝wH，w为小波变换基，λ＝0.05；采用l₁范数最小化方法求解优化问题，得到自适应三角网格节点的像素值H，再利用三角线性插值方法得到每个像素点的值，完成原始图像的重建。

有益效果

本发明实现了一种光声层析技术中图像重建的方法，该方法可以将高维几何关系矩阵M进行降维，降维后的矩阵大大减少了内存存储空间，结合压缩感知理论，在利用l1范数最优化方法求解方程的迭代过程中加快计算速度。和现有重建方法比较具有精度高，内存占用少，对硬件配置要求不高的特点。

本发明之所以具有上述的有益效果其原因在于：在使用基于模型的重建方法前通过一系列算法构建内存占用较小的矩阵，这在3D图像重建过程中优势尤为明显。先采用通用反投影重建算法(UBP)快速初步重建出原始图像的大概模型，该模型包含了原始图像的基本内容。接着，根据这初步图像利用图像特征映射提取算法、自适应网格节点生成算法以及图像的Delaunay三角化方法，产生自适应的三角网格。然后，根据形成的自适应网格采用三角插值的方法构建几何关系矩阵M，这里M相比较基于像素以及基于均匀网格构建的几何关系矩阵内存占有小很多。接着利用压缩感知理论(CS)将重建信号H利用小波变换基投射到频域空间使得信号尽可能的稀疏，这相当于在求解方程组时降低了未知数的个数，因此为求解原始图像提供了一个数学解决思路。最终达到在测量环境以及硬件配置并不理想的情况下能够较快的高精度的重建出原始图像。

附图说明

图1是本发明的方法流程图。

图2是本发明中Delaunay三角化算法示意图：(a)点集合；(b)外面边界三角形带有尾巴的点集的delaunay三角剖分；(c)演示外接圆效应的特征圆。

图3是本发明中生成自适应三角化网格的示意图。

图4是本发明中列举的相关算法的对比实验：(a)：点模型的图像重建实验；(b)：圆模型的图像重建实验；(c)线模型的图像重建实验；左上角为反投影重建算法(UBP)的重建图像，右上角为迭代算法重建图像，左下角为自适应算法重建图像，右下角为原始图像。

具体实施方式

下面结合具体实施例描述本发明：

附图1展示了本发明所实现基于模型的自适应压缩感知重建算法的光声层析重建方法的总流程图，该总流程图包含了实现最终重建所需的各个主要步骤。本发明的目的是通过对基于模型的重建算法中涉及的几何关系矩阵Γ不断降维，达到减少存储空间，使重建信号尽可能稀疏，减少计算时间消耗的效果。下面以2D图像重建为例，扫描方式围绕目标物体作圆周扫描。

具体的实现步骤为：

步骤1：初步图像重建：根据所获得的超声波信号利用解析法—通用反投影重建算法(UBP)进行图像的初步重建。该步骤可以在很短时间内快速重建原始图像的大概模型。

通用反投影算法：

从物理学角度，光声层析成像技术表示的是一种逆向求源问题。即根据超声波探测器检测到的超声波信息复原出原始图像的过程。在宽束短脉冲激光激励下，即H(r,t)＝H_s(r)δ(t)，满足热封闭和声压封闭条件，光声波动方程可写成：

${\▿}^{2}p(r,t)-\frac{1}{v_s^2}\frac{{\∂}^{2}p(r,t)}{{\∂t}^2}=-\frac{p_0\left(r\right)}{v_s^r}\frac{\∂\mathrm{\δ}\left(t\right)}{\∂t}---\left(1\right)$

上式中，初始声压p₀(r)＝Γ(r)H_s(r)即为光声层析技术所要重建的图像。在大多数情况下，我们将超声波信号表达式简化为如下公式：

$p(r,t)=\frac{\mathrm{\Γ}}{4\mathrm{\πc}}\frac{\∂}{\∂t}{\∫}_{{\mathrm{\θ}}_1}^{{\mathrm{\θ}}_2}{H}_r(x_0+R\cos \mathrm{\θ},y_0+R\sin \mathrm{\θ})\mathrm{d\θ}---\left(2\right)$

式中，R为探测器与重建像素点之间的距离。而在该算法中我们需要对式(1)两边同时进行傅里叶变换，得：

${\▿}^{2}p(r,w)+{\left(\frac{w}{v_s}\right)}^{2}p(r,w)=j\frac{w}{v_s^2}{p}_0\left(r\right)---\left(3\right)$

上式中，p(r,w)为p(r,t)的傅里叶变换，即：

$p(r,w)=\underset{-\∞}{\overset{+\∞}{\∫}}p(r,t)e^{-\mathrm{jwt}}\mathrm{dt}---\left(4\right)$

式(4)中，令k＝w/v_s，则有：

$p(r,w)=\underset{-\∞}{\overset{+\∞}{\∫}}\frac{1}{v_s}p(r,\stackrel{\‾}{t}/v_s)e^{-\mathrm{jk}\stackrel{\‾}{t}}d\stackrel{\‾}{t}---\left(5\right)$

令则有p(r,k)＝v_sp(r,w)，于是式(5)改写为：

${\▿}^{2}p(r,k)+k^{2}p(r,k)={\mathrm{jkp}}_0\left(r\right)---\left(6\right)$

由格林定理，测量声压谱p(r₀,k)可表示为：

$p(r_0,k)=\mathrm{jk}\underset{V}{\∫}G(r_0,r,k)p_0\left(r\right)\mathrm{dr}---\left(7\right)$

式中，V为包含组织域的封闭检测面所含体积，G(r₀,r,k)为对应式的格林函数，且有：

$G(r_0,r,k)=\frac{e^{\mathrm{jk}|r_0-r|}}{4\mathrm{\π}|r_0-r|}---\left(8\right)$

成像积分方程式(8)可由解析或者数值方法求解，反投影解析解为：

$p_0\left(r\right)=\frac{-1}{{2\mathrm{\πc}}^2\mathrm{\Γ}}\underset{t=|r-r_0|/c}{\∫}\frac{1}{t}[\frac{\∂p(r,t)}{\∂t}-\frac{p(r,t)}{t}]\mathrm{dr}---\left(9\right)$

通过式(9)，我们可以通过获得的超声波信号方便快速的求出原始图像的光吸收系数分布从而达到初步图像(先验模型)重建的效果。

步骤2：自适应三角化网格生成：将步骤1得到的初步重建结果作为先验模型，依次采用图像特征映射提取算法、自适应网格节点生成算法以及图像的Delaunay三角化方法(参考文献A fast Approach for Accurate Content-Adaptive Mesh Generation)，产生自适应的三角网格；该网格依据先验模型产生，能够体现原始图像的基本结构。这样在图像内容较复杂的区域三角形网格划分的更加精细，在图像内容简单单一的区域三角形网格划分较大。由此可以集中硬件资源来重建我们需要重建的图像区域，有效避免传统均匀网格化所带来的大量不必要的内存存储和计算时间的消耗。

图像特征映射提取：

根据反投影算法初步得到原始图像的先验模型，我们将图像函数用f(x,y)表示，指在坐标x,y处图像的函数值。我们定义G(x,y)为点x,y处最大的二阶导数值。

$G(x,y)=\underset{\mathrm{\θ}\∈\left[\mathrm{0,2}\mathrm{\π}\right]}{\max }\left|f_{\mathrm{\θ}}^{\′\′}(x,y)\right|---\left(10\right)$

G(x,y)的计算方法如下：

G(x,y)＝2max{|f″_xx(x,y)|,|f″_xy(x,y)|,f″_yy(x,y)}   (11)

我们定义特征函数：

$\mathrm{\σ}(x,y)={\left(\frac{G(x,y)}{A}\right)}^{\mathrm{\γ}}---\left(12\right)$

式(12)中，A指图中所有像素点中二阶导数最大的值，γ为一个正常数，通常将其设置为γ＝1。可以理解为每个像素点归一化的过程。这样直到遍历完整个图像上所有的像素点，每个点的特征函数值即可求得。

自适应网格节点的产生：

弗洛伊德-斯坦伯格算法是在图像数字半色调处理过程中一种经典的基于偏差融合的算法。我们使用修改后的弗洛伊德-斯坦伯格算法依据步生成的特征映射函数σ(x,y)自适应生成网格节点，具体算法流程如下。

1)、定义指示函数b(i,j)来表示点i,j是否有节点，若该位置有网格节点我们使其值为1，否则使其为0。

$b(i,j)=\left\{\begin{array}{c}1,\mathrm{if\σ}(i,j)\≥q\\ 0,\mathrm{else}\end{array}\right.---\left(13\right)$

式中q为阀值。

2)、计算点(i,j)处的数值偏差。定义偏差函数e(i,j)：

e(i,j)＝σ(i,j)-(2q)b(i,j)   (14)

3)、点(i,j)处四个相邻的像素点(i,j+1)、(i+1,j-1)、(i+1,j)、(i+1,j+1)的值会根据点(i,j)的偏差函数e(i,j)进行调整。

σ(i,j+1)＝σ(i,j+1)+w₁e(i,j),

σ(i+1,j-1)＝σ(i+1,j-1)+w₂e(i,j),

σ(i+1,j)＝σ(i+1,j)+w₃e(i,j),

σ(i+1,j+1)＝σ(i+1,j+1)+w₄e(i,j).   (15)

式中w_i(i＝1,2,3,4)，分别为7/16,3/16,5/16,1/16，一种更加合理的确定wi的方法是依据相邻像素之间的几何距离。

4)、重复1)、2)、3)步，直到遍历所有的图像像素点。即可求得初始图像上的自适应节点。

三角化：

根据生成的自适应节点将图像进行三角网格化。这里我们使用的是德劳内(Delaunay)三角网格剖分方法。其目的是使剖分得到的每个三角形中最小的角最大化。这就意味着Delaunay三角剖分技术力图避免出现瘦长三角形。如图2我们可以看到三角剖分的技术要点，那就是任何三角形的外接圆都不包任何括其他点这就叫做外接圆性质。从计算效率上来讲，Delaunay三角化从一个远离三角形边界的外轮廓开始运行，图2(b)中代表虚构的外三角由虚线引至最高点。图2(c)中是表明一些外接圆性质的例子，包括一个连接两个界外点的实时数据的一个顶点虚拟外部三角形。现在又很多种算法计算Delaunay三角剖分，其中一些很有效，但内部描述却很难描述。其中一个较简单的算法要点如下：

1)、添加外部三角形，并从它的一个顶点(在这里会产生一个确定的外部起点)处开始。

2)、加入一个内部点；在三角形的外接圆内搜索该点，去掉包含该点的三角剖分。

3)、重新构造三角图，在刚刚去掉的三角形外接圆内包括新的点。

4)、返回2)直到没有新的点加入。

这样，即可根据先验模型得到自适应三角形网格，网格示意图如图2所示。

步骤3：根据步骤2得到的三角形网格结合三角线性插值方法建立正向模型P＝MH，其中P指换能器检测获得的超声波能量信号矩阵，矩阵M指初始声波压力与探测信号之间的几何关系矩阵，H指自适应三角网格节点的像素值；利用这种方式我们得到的矩阵相比传统基于像素构建的几何关系矩阵以及基于均匀三角网格构建的几何关系矩阵大大降低了存储空间。

根据图3所示的产生的自适应三角网格的示意图，我们引入三角插值方法。对于图3中任意一个三角形来说，已知(x₁,y₁,z₁)，(x₂,y₂,z₂)，(x₃,y₃,z₃)，则三角形的插值函数H_r(x,y)可以方便求得：

$H_r(x,y)=-\frac{\mathrm{Ax}+\mathrm{By}+D}{C}---\left(16\right)$

式中

$A=\left|\begin{array}{ccc}1& y_2& z_2\\ 1& y_3& z_3\\ 1& y_1& z_2\end{array}\right|B=\left|\begin{array}{ccc}{x}_2& 1& z_2\\ x_3& 1& z_3\\ x_1& 1& z_1\end{array}\right|$

$C=\left|\begin{array}{ccc}{x}_2& y_2& 1\\ x_3& y_3& 1\\ x_1& y_1& 1\end{array}\right|D=\left|\begin{array}{ccc}{x}_2& y_2& z_2\\ x_3& y_3& z_3\\ x_1& y_1& z_1\end{array}\right|---\left(17\right)$

|.|为矩阵行列式的值，A、B、C、D为三角插值函数系数，这样我们可以通过插值公式(16)求解三角形内任意一点的像素值。因为本例扫描方式为圆周扫描如图3所示，x₀,y₀为探测器所在位置的横坐标和纵坐标，R为像素点与探测器S之间的距离，θ表示圆心角。根据几何关系可以知道图像上任何一点(x,y)的像素值，x＝x₀+Rcosθ，y＝y₀+Rsinθ。所以我们可以将公式2用如下公式来表示：

$p(r,t)=\frac{\mathrm{\Γ}}{4\mathrm{\πc}}\frac{\∂}{\∂t}{\∫}_{{\mathrm{\θ}}_1}^{{\mathrm{\θ}}_2}{H}_r(x_0+R\cos \mathrm{\θ},y_0+R\sin \mathrm{\θ})\mathrm{d\θ}---\left(18\right)$

根据上述所求得的自适应三角形网格即可对每一个三角形进行上述插值，模拟超声波信号p(r,t)进一步变为：

$p(r,t)=-\frac{\mathrm{\Γ}}{4\mathrm{\πc}}\frac{\∂}{\∂t}\underset{{\mathrm{\θ}}_1}{\overset{{\mathrm{\θ}}_2}{\∫}}\frac{A(x_0+R\cos \mathrm{\θ})+B(y_0+R\sin \mathrm{\θ})+D}{C}\mathrm{d\θ}---\left(19\right)$

则正向模型用矩阵的形式来表示：

P＝MH   (20)

式中，P指换能器检测获得的超声波能量信号矩阵，矩阵M指初始声波压力与探测信号之间的几何关系矩阵，H指自适应三角网格节点的像素值。

步骤4：利用步骤3建立的正向模型，根据换能器检测获得的超声波能量信号矩阵P和初始声波压力与探测信号之间的几何关系矩阵M求解自适应三角网格节点的像素值H，再结合三角线性插值公式(16)就可以恢复出图像上所有的像素点的值。

现在我们已经将图像重建问题转化为式(20)，需要根据P和M去求解H，但是，我们发现矩阵M是稀疏的，我们不能够简单的利用M^-1P求解H，因此我们将式(20)转化为优化如下的l₂范数问题。

$H^{\′}=\arg \underset{H}{\min }{\left|\right|P_d-\mathrm{MH}\left|\right|}_2---\left(21\right)$

式中，H'为估计的最优原始信号值，P_d为换能器获得的超声波信号。一般情况下，相比较原始图像像素的个数N，换能器获得的探测信号数量S要远小的多。因此式(21)是一个欠定问题，简单的最小二乘法难以求解。

考虑到几何关系矩阵A的稀疏性，引入压缩感知理论即可以进一步减小矩阵的大小，又为求解提供很好的解决方法。因此将式(21)用下式表达：

$\min {\left|\right|P_d-{\mathrm{Mw}}^{-1}{H}^{*}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|H^{*}\left|\right|}_1---\left(22\right)$

其中P_d为换能器获得的超声波信号，H^*＝wH，w为小波变换基，λ＝0.05；这样我们将问题转化为求解凸优化方程的问题，这里我们用l₁范数最小化方法进行求解。当探测器获得信号与模拟获得的信号之间的误差达到一定阀值时或者循环迭代次数达到一定次数时，即跳出循环，输出自适应节点的像素值H，此时的H即为所要求解的目标物体的光能量吸收分布即原始图像的像素点。根据自适应三角网格节点的像素值利用三角线性插值的算法，利用式(22)求解每个像素点的值，最终达到在测量环境不理想的情况下快速高质量的重建图像。

本发明算法归纳如下：

(1)：首先根据探测器获得的超声波信号结合通用反投影算法(UBP)初步重建原始图像的大概模型。

(2)：用得到的初始重建图像相继用图像特征映射提取算法、自适应网格节点生成算法以及图像的Delaunay三角化方法，产生自适应的三角网格，自适应三角网格示意图如图2所示。

(3)：用基于模型的重建方法对图像进行进一步高精度重建，在建立正向模型P＝MH时，根据产生的自适应三角网格引入三角线性插值方法，构建几何关系矩阵M，再结合压缩感知理论使得重建信号在另一个空间域里尽可能稀疏，由此将问题转化为求解式(22)的凸优化问题，最后用l1范数最小化方法进行求解，当迭代误差小于一定阀值或者达到一定迭代次数即跳出循环。

为了评价本发明的预测性能和比较仿真实验的结果，我们使用以下指标来说明：均方根误差(RMSE)：

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{F}{\mathrm{\Σ}}}{(y_i-{\hat{y}}_i)}^2}{F}}---\left(23\right)$

F指数据点的总数量，y_i和分别表示i个实际值和预测值，RMSE的值越小，表明预测精度越高。

以下给出仿真试验：

仿真条件：仿真实验基于linux平台，使用C++计算机语言编程。本实验是以128乘128的2D图像进行实验。探测器在图像周围做圆周扫描，扫描360个角度，每个角度采样700次数据。本次试验共使用三种模型—点、圆、线模型，每种模型使用反投影重建(UBP)、迭代算法、自适应算法进行实验，以下分别从时间、内存、信噪比三个性能项目对迭代算法和自适应算法进行分析。

表1.不同模型在两种不同算法下关于时间、内存、信噪比的性能对比

从上表我们可以看出在可接受的时间范围内，自适应算法保持高质量的重建效果的同时，大大降低了内存存储空间的消耗。从图4所示的结果，我们可以发现自适应算法相比较反投影重建算法(UBP)具有更好的图像重建质量。

Claims (1)

1.一种面向无损检测的自适应压缩感知的重建方法，其特征在于：包括以下步骤：

步骤1：根据所获得的超声波信号利用解析法—通用反投影重建算法进行图像的初步重建；

步骤2：将步骤1得到的初步重建结果作为先验模型，依次采用图像特征映射提取算法、自适应网格节点生成算法以及图像的Delaunay三角化方法，产生自适应的三角网格；

步骤3：根据步骤2得到的三角形网格结合三角线性插值方法建立正向模型P＝MH，其中P指换能器检测获得的超声波能量信号矩阵，矩阵M指初始声波压力与探测信号之间的几何关系矩阵，H指自适应三角网格节点的像素值；

步骤4：利用步骤3建立的正向模型，根据换能器检测获得的超声波能量信号矩阵P和初始声波压力与探测信号之间的几何关系矩阵M求解自适应三角网格节点的像素值H：

将正向模型P＝MH求解问题转换为优化问题，其中P_d为换能器获得的超声波信号，H^*＝wH，w为小波变换基，λ＝0.05；采用l₁范数最小化方法求解优化问题，得到自适应三角网格节点的像素值H，再利用三角线性插值方法得到每个像素点的值，完成原始图像的重建。