CN111047663A - 电学层析成像伪影抑制图像重建方法 - Google Patents

Abstract

本发明公开了电学层析成像伪影抑制图像重建方法，包含如下步骤：(1)根据被测场域获取重建所需的相对边界测量值向量b_meas和灵敏度矩阵A；(2)设置初始化参数；(3)更新自适应正则化参数λ^k；(4)更新辅助变量p^k+1；(5)更新权重因子ω^k；(6)更新辅助变量q^k+1；(7)更新辅助变量r^k+1；(8)更新电导率分布g^k+1；(9)更新增广拉格朗日乘子γ₁ ^T,γ₂ ^T,和γ₃ ^T；(10)判断迭代是否符合迭代终止条件

或者k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步，继续迭代求解；(11)根据最终求解所得灰度值进行成像。本发明可以有效地降低阶梯伪影及提高抗噪性，同时简化了算法结构，提高了算法的适用性。

Description

电学层析成像伪影抑制图像重建方法

技术领域

本发明属于电学层析成像技术领域，涉及利用电学层析成像伪影抑制图像重建方法实现图像重建，具体涉及电学层析成像伪影抑制图像重建方法。

背景技术

电学层析成像技术(Electrical Tomography，ET)出现于20世纪80年代后期，是一种基于电特性敏感机理的过程层析成像技术，它基于边界测量值对测量区域介质的电特性(电导率/介电系数/复导纳/磁导率)分布信息进行成像，进而得到介质的分布信息。电学层析成像在多相流、地质勘探以及医学成像领域有广泛的应用前景，可以实现长期、持续监测，进而实现功能性成像。

电学层析成像逆问题(即图像重建问题)求解具有非线性，通过线性化处理可以将问题转化为线性逆问题求解。针对逆问题求解的不适定性，通常选取正则化方法处理逆问题，例如Tikhonov正则化方法因其简单、稳定、直接的重建而不需要迭代而引起了人们的广泛关注。如Y B Xu等人2016年发表于《流量测量和仪器》(Flow Measurement andInstrumentation)第50卷，第1-12页，题为电阻层析成像自适应Tikhonov正则化参数选择方法(An adaptive Tikhonov regularization parameter choice method forelectrical resistance tomography)，这种方法能很好地恢复边缘光滑的物体，而当恢复边缘锐利的物体时则会出现过度的平滑。为了缓解Tikhonov正则化方法带来的过度光滑性，研究了用全变分(Total variation regularization method,TV)正则化方法重建尖锐不连续物体的方法。如JLXu等人2016年发表于《AEU-国际电子和通讯杂志》(AEU-International Journal of Electronics and Communications)第70卷，第1128-1133页，题为通过改进的变分模型对图像进行去模糊和去噪(Image deblurring and denoisingby an improved variational model)，然而，无论原始目标是否连续，都会产生分段不变的图像，称为“阶梯伪影”。为了抑制阶梯伪影，对TV 正则化方法的改进进行了大量的工作，其中，混合正则化方法在克服阶梯伪影和保留锐利边缘方面提供了一种有效的解决方案。如S Wang等人2018年发表于《数值算法》 (Numerical Algorithms)第78卷，第513-533页，题为利用一阶和二阶全变差去除超声图像中的散斑噪声(Speckle noise removal inultrasound images by first and second-order total variation)和TAdam等人2019年发表于《多维系统与信号处理》(Multidimensional Systems and Signal Processing)第30卷，第503-527页，题为基于重叠群稀疏的高阶非凸总变分的图像去噪(Image denoisingusing combined higher order non-convex total variation with overlapping groupsparsity)等文章针对混合方法在去除阶梯伪影方面给出了相应的改进策略。

在现有的研究中，这些混合正则化算法都局限在对噪声敏感的L₂范数保真度项和在一定程度上产生阶梯伪影的各向同性TV惩罚项，这使得成像结果容易受到噪声的影响，并且无法很好地去除阶梯伪影。因此，这种混合方式仍旧无法得到很好的成像效果，而且还复杂了算法的结构，使得很难在实际运用中推广使用。

发明内容

本发明的目的在于提出一种电学层析成像伪影抑制图像重建方法，该方法可以有效地降低阶梯伪影及提高抗噪性，同时简化了算法结构，提高了算法的适用性。

本发明为实现上述目的采用如下技术方案，电学层析成像伪影抑制图像重建方法，其特征在于具体过程为：将电学层析成像看作一个线性不适定问题Ag＝b_meas，其中，A为灵敏度矩阵，b_meas为相对边界测量值向量，g为所求成像灰度值；设计最小化的目标函数为：

其中，λ为正则化参数，ω为权重因子，||·||₁为向量L₁范数，

为TV各向异性惩罚项，▽为梯度算子，l_x和l_y分别为x和y方向上等距分布的离散网格数，

和

为一阶离散微分，

为高阶惩罚项，

为二阶离散微分；

在分裂变量法的基础上引入3个辅助变量p、q和r，将最小化的目标函数中的无约束问题转化为约束极小化形式：

s.t. p＝Ag-b_meas,q＝▽g,r＝▽²g,

上述约束极小化形式可通过交替方向乘子法表示为增广拉格朗日函数：

其中，γ₁，γ₂和γ₃为增广拉格朗日乘子，上标T是转置算子，||·||₂为向量L₂范数；

图像重建算法即对上述增广拉格朗日函数的求解包含以下步骤：

(1)根据被测场域获取重建所需的相对边界测量值向量b_meas和灵敏度矩阵A；边界测量值的获取是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布L个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压，相对边界测量值向量b_meas为不含内含物的空场边界测量电压向量b_meas1和含有内含物的有物场的边界测量电压向量b_meas2之差；(2)设置初始化参数，调节最优权重因子的正参数ξ，最小迭代阈值θ_min，最大迭代次数k_max，惩罚参数β₁,β₂,β₃，增广拉格朗日乘子γ₁,γ₂,γ₃和调节最优正则化参数的范围因子H，0＜H＜100；(3)更新自适应正则化参数λ^k；(4)更新辅助变量p^k+1；(5)更新权重因子ω^k；(6)更新辅助变量q^k+1；(7)更新辅助变量r^k+1；(8)更新电导率分布g^k+1；(9)更新增广拉格朗日乘子γ₁ ^T,γ₂ ^T和γ₃ ^T；(10) 判断迭代是否符合迭代终止条件

或者k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步，继续迭代求解；(11)根据最终求解所得灰度值进行成像。

进一步优选，所述图像重建算法的具体步骤为：

(1)针对三种典型模型分别获取各自重建所需的边界测量值b_meas和灵敏度矩阵A：边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值；逆问题右端项b_meas为不含内含物的空场边界电压b_meas1和含有内含物的有物场的边界测量电压b_meas2之差即右端项相对边界测量值b_meas＝b_meas2-b_meas1；

灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理论，计算灵敏度矩阵，计算公式为：

其中，A_ij是第j个电极对对第i个电极对的灵敏度系数，

分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)设置初始化参数：调节最优权重因子的正参数ζ＝0.3，最小迭代阈值θ_min＝10^-6，初始迭代次数k＝0，最大迭代次数k_max＝200，惩罚参数β₁＝300，β₂＝β₃＝50，增广拉格朗日乘子

和调节最优正则化参数的范围因子H＝0.03，初始灰度值g⁰＝A^Tb_meas；

(3)更新自适应正则化参数λ^k：

其中，

和

分别代表第 k次迭代时g的最小值和最大值，d代表求解逆问题时有限元剖分的平均网格长度，H为调节最优正则化参数的因子0＜H＜100；

(4)更新辅助变量p^k+1：

因此，其解可表示为

其中，|·|and

分别表示分量绝对值和点态乘积，“sgn”表示为正负号函数；

(5)更新权重因子ω^k：

其中，ξ是一个正参数，它被调整以确保ω最优值；

(6)更新辅助变量q^k+1：

因此，其解可以表示为

(7)更新辅助变量r^k+1：

因此，其解可以表示为

(8)更新电导率分布g^k+1：

为了提高计算效率和提高解的稳定性，采用快速傅里叶变换求解g子问题：

其中，上标T是转置算子，▽^T▽,(▽²)^T▽²和A^TA是具有循环块结构的块循环矩阵；

(9)更新增广拉格朗日乘子γ₁ ^T,γ₂ ^T,和γ₃ ^T：在每次迭代中，均更新增广拉格朗日乘子，

和

(10)判断迭代是否符合迭代终止条件

或者k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步，继续迭代求解；

(11)根据最终求解所得灰度值进行成像。

本发明具有以下有益效果：本发明的电学层析成像伪影抑制图像重建方法通过自适应加权因子来控制一阶和高阶项之间的权重和自适应正则化参数平衡了保真度项和混合TV 项之间的权衡。另外，采用L₁范数作为保真度项来增强算法的抗噪性能。为了解决所提出的极小化问题，研究了一种高效、快速的交替方向乘子迭代算法。相比于各向同性TV和各向异性TV，本发明提出的一种电学层析成像伪影抑制图像重建方法在降低阶梯伪影和提高抗噪性能等方面有着很好的效果。同时自适应参数的选择增强了参数选择的客观性和简易性和交替方向乘子迭代算法的使用简化了算法的复杂度，提高了算法的运行速度，增强了算法的适用性。

附图说明

图1为本发明的电学层析成像伪影抑制图像重建方法的流程框图；

图2为本发明的电阻层析成像系统圆形单截面被测场域，激励电流和测量电压的模式以及电极分布；

图3为本发明的实施例在选取三种模型的真实分布时，各向同性TV算法，各向异性TV算法和电学层析成像伪影抑制图像重建方法的图像重建结果示意图；

图4为权重因子与梯度之间的关系；

图5为三种方法在不同噪声水平下对模型B的重建图像；

图6为三种方法在不同噪声水平下模型B的图像相对误差和相关系数。

图中：1-被测场域；2-电极；3-激励电流；4-测量电压。

具体实施方式

结合附图和实施例对本发明的电学层析成像伪影抑制图像重建方法加以具体说明。

一种电学层析成像伪影抑制图像重建方法，针对传统全变分正则化算法产生的阶梯伪影问题，以各向异性TV惩罚项和高阶惩罚项为基础，结合采用L₁范数作为保真度项，通过自适应选择正则化参数和权重因子来选取最优值，以一种高效、快速的交替方向乘子迭代算法来求解本发明提出的目标函数，完成最终的逆问题求解。

如图1所示，为本发明的电学层析成像伪影抑制图像重建方法流程图。

如图2所示，为电学层析成像之一的电阻抗层析成像系统圆形单截面被测场域1，激励电流3和测量电压4的模式以及电极分布，采用16电极2均匀分布在被测场域1外壁。

选取三种典型的介质模型为实施例，场域内物体真实分布如图3左侧一竖列所示，其他三列从左到右分别表示为各向同性TV算法，各向异性TV算法和电学层析成像伪影抑制图像重建方法。为了较好地体现本发明中一种电学层析成像伪影抑制图像重建方法与其它两种方法的不同，分别给出了三种典型模型在这三种正则化算法下的求解结果。

实施例包括如下具体步骤：电学层析成像伪影抑制图像重建方法，该方法将电学层析成像看作一个线性不适定问题Ag＝b_meas，其中，A为灵敏度矩阵，b_meas为相对边界测量值向量，g为所求成像灰度值，逆问题可以用最小二乘优化形式的目标函数表示：

其中，f(g)为目标函数，正则化方法为提高解的稳定性提供了一种新的选择,其形式可以描述为：

其中，λ是正则化参数，它控制最小二乘项(通常称为数据保真度项)

和正则化项R(g)之间的权重。

TV正则化算法在保持锐利边缘方面表现出良好的效果。其一般形式可以表示为：

其中，||▽g||₁表示TV正则化算法中的惩罚项，可以定义为：

其中，

在这里我们假设灰度值g被定义在等距网格上， {(x_i,y_j)|x_i＝iΔx,y_j＝jΔy,i＝1,2,...,l_x,j＝1,2,...,l_y}，其中，g_i,j表示的是g在点(x_i,y_j)处的离散值。

在边缘保持方面，TV正则化方法有着明显的优势。然而，以L₂范数为保真度项的正则化方法对测量噪声敏感，鲁棒性差。此外，由于各向同性罚项及其分段逼近，在恢复图像中观察到明显的阶梯伪影。为了克服传统TV正则化方法的不足，提高重建图像的质量，

为了克服传统TV正则化方法的不足，提高重建图像的质量，本发明针对传统TV正则化算法产生的阶梯伪影问题，以各向异性TV惩罚项和高阶惩罚项为基础，结合采用L₁范数作为保真度项，通过自适应选择正则化参数和权重因子来选取最优值，以一种高效、快速的交替方向乘子迭代算法来求解。其最小化的目标函数为：

其中，

这里，λ是自适应正则化参数，ω∈[0,1]是自适应权重因子。

在分裂变量法的基础上，引入3个辅助变量p、q和r，将目标函数中的无约束问题转化为约束极小化形式：

s.t.p＝Ag-b_meas,q＝▽g,r＝▽²g,

为了得到约束问题的最优解，研究了一种有效的迭代算法，以解决所提出的约束极小化模型。基于交替方向乘子法，增广拉格朗日函数可以表示为：

图像重建算法包含以下步骤：

(1)针对三种典型模型分别获取各自重建所需的边界测量值b_meas和灵敏度矩阵A：边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极(如图2所示)，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值；逆问题右端项b_meas为不含内含物的空场边界电压b_meas1和含有内含物的有物场的边界测量电压b_meas2之差(即右端项相对边界测量值b_meas＝ b_meas2-b_meas1)；

灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理论，计算灵敏度矩阵，计算公式为：

其中A_ij是第j个电极对对第i个电极对的灵敏度系数，

分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)设置初始化参数：调节最优权重因子的正参数ζ＝0.3,最小迭代阈值θ_min＝10^-6,初始迭代次数k＝0，最大迭代次数k_max＝200,惩罚参数β₁＝300,β₂＝β₃＝50,增广拉格朗日乘子

和调节最优正则化参数的范围因子H＝0.03，初始灰度值g⁰＝A^Tb_meas；

(3)更新自适应正则化参数λ^k：

其中，

和

分别代表第k次迭代时g的最小值和最大值，d代表求解逆问题时有限元剖分的平均网格长度，H为调节最优正则化参数的因子0＜H＜100；

(4)更新辅助变量p^k+1：

因此，其解可以表示为

其中，|·|and

分别表示分量绝对值和点态乘积，“sgn”表示为正负号函数；

(5)更新权重因子ω^k：

其中，ξ是一个正参数，它被调整以确保ω最优值。图4为权重因子与梯度之间的关系，可以看出ω选择为1，当||▽g^k||₁相对较小或较大的，如红色方框所示。ω选择为0，当||▽g^k||₁取中间值时，如蓝色方框所示。当阶梯伪影出现时，该方法为如何确定ω提供了一种有效的解决方案；

(6)更新辅助变量q^k+1:

因此，其解可以表示为

(7)更新辅助变量r^k+1：

因此，其解可以表示为

(8)更新电导率分布g^k+1：

为了提高计算效率和提高解的稳定性，采用快速傅里叶变换求解了g子问题：

其中，上标T是转置算子，▽^T▽,(▽²)^T▽²和A^TA是具有循环块结构的块循环矩阵；

(9)更新增广拉格朗日乘子γ₁ ^T,γ₂ ^T,和γ₃ ^T：在每次迭代中，我们都会更新增广拉格朗日乘子：

和

(10)判断迭代是否符合迭代终止条件

或者k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步，继续迭代求解；

(11)根据最终求解所得灰度值进行成像。

图3为本发明的实施例选取三种典型的介质模型，各向同性TV算法，各向异性TV算法和电学层析成像伪影抑制图像重建方法的图像重建结果示意图。可以看出，三种典型模型中，采用各向同性TV正则化方法时，背景中存在明显的阶梯伪影，严重影响了图像恢复的质量。此外，它还不能准确地识别正方形包含物(例如，模型B)和具有方形和圆的混合包含(例如，模型C)的形状。与各向同性TV相比，各向异性TV有效地降低了阶梯伪影。然而，由于存在较大的几何畸变，具有圆形的夹杂物(如模型A和模型C)被扭曲为方形。相比于前面两种正则化算法，一种电学层析成像伪影抑制图像重建方法的重建图像的阶梯伪影得到了很好的抑制，并保持了夹杂物的形状，并且有着更高的图像分辨率，明显的提高了图像重建的质量。在实际应用中，电学层析成像系统容易受到噪声的影响。为了验证本发明提出的一种电学层析成像伪影抑制图像重建方法的抗噪声性能，图5 比较了三种方法在不同噪声水平下对模型B的重建图像。可以看出，当噪声水平逐渐增大时，本发明提出的方法在抑制阶梯伪影和保持锐利边缘方面优于其它两种方法。同时，在电学层析成像中，通常采用图像相对误差(Relative Error，RE)和相关系数(Correlation Coefficient，CC)评价算法来定量图像重建质量，表达式如①、②所示，图像相对误差越小，相关系数越大，表明图像重建质量越好。

其中，σ是重建区域的计算电导率，σ^*是实际电导率，t表示像素数，

和

表示σ和σ^*的平均值，σ_i和σ_i ^*表示的是σ和σ^*的第i个三角形单元。

图6给出了这三种方法不同噪声水平下模型B的图像相对误差和相关系数，可以看出，本发明提出的一种电学层析成像伪影抑制图像重建方法与各向同性TV算法，各向异性TV 算法相比，具有最低的相对误差和最高的相关系数，能够准确的描绘出被测场域内部分布，明显地提高了电学层析成像逆问题求解精度和抗噪声能力。

以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.电学层析成像伪影抑制图像重建方法，其特征在于具体过程为：将电学层析成像看作一个线性不适定问题Ag＝b_meas，其中，A为灵敏度矩阵，b_meas为相对边界测量值向量，g为所求成像灰度值；设计最小化的目标函数为：

其中，λ为正则化参数，ω为权重因子，||·||₁为向量L₁范数，

为TV各向异性惩罚项，

为梯度算子，l_x和l_y分别为x和y方向上等距分布的离散网格数，

和

为一阶离散微分，

为高阶惩罚项，

为二阶离散微分；

在分裂变量法的基础上引入3个辅助变量p、q和r，将最小化的目标函数中的无约束问题转化为约束极小化形式：

上述约束极小化形式可通过交替方向乘子法表示为增广拉格朗日函数：

其中，γ₁，γ₂和γ₃为增广拉格朗日乘子，上标T是转置算子，||·||₂为向量L₂范数；

图像重建算法即对上述增广拉格朗日函数的求解包含以下步骤：(1)根据被测场域获取重建所需的相对边界测量值向量b_meas和灵敏度矩阵A；边界测量值的获取是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布L个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压，相对边界测量值向量b_meas为不含内含物的空场边界测量电压向量b_meas1和含有内含物的有物场的边界测量电压向量b_meas2之差；(2)设置初始化参数，调节最优权重因子的正参数ξ，最小迭代阈值θ_min，最大迭代次数k_max，惩罚参数β₁,β₂,β₃，增广拉格朗日乘子γ₁,γ₂,γ₃和调节最优正则化参数的范围因子H，0＜H＜100；(3)更新自适应正则化参数λ^k；(4)更新辅助变量p^k+1；(5)更新权重因子ω^k；(6)更新辅助变量q^k+1；(7)更新辅助变量r^k+1；(8)更新电导率分布g^k+1；(9)更新增广拉格朗日乘子γ₁ ^T,γ₂ ^T和γ₃ ^T；(10)判断迭代是否符合迭代终止条件

或者k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步，继续迭代求解；(11)根据最终求解所得灰度值进行成像。

2.根据权利要求1所述的电学层析成像伪影抑制图像重建方法，其特征在于所述图像重建算法的具体步骤为：

(1)针对三种典型模型分别获取各自重建所需的边界测量值b_meas和灵敏度矩阵A：边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值；逆问题右端项b_meas为不含内含物的空场边界电压b_meas1和含有内含物的有物场的边界测量电压b_meas2之差即右端项相对边界测量值b_meas＝b_meas2-b_meas1；

灵敏度矩阵是根据不含内含物的空场的边界测量电压，结合灵敏度理论，计算灵敏度矩阵，计算公式为：

其中，A_ij是第j个电极对对第i个电极对的灵敏度系数，

分别为第i个电极对及第j个电极对在激励电流为I_i，I_j时场域电势分布；

(2)设置初始化参数：调节最优权重因子的正参数ζ＝0.3，最小迭代阈值θ_min＝10^-6，初始迭代次数k＝0，最大迭代次数k_max＝200，惩罚参数β₁＝300，β₂＝β₃＝50，增广拉格朗日乘子

和调节最优正则化参数的范围因子H＝0.03，初始灰度值g⁰＝A^Tb_meas；

(3)更新自适应正则化参数λ^k：

其中，

和

分别代表第k次迭代时g的最小值和最大值，d代表求解逆问题时有限元剖分的平均网格长度，H为调节最优正则化参数的因子0＜H＜100；

(4)更新辅助变量p^k+1：

因此，其解可表示为

其中，|·|and

分别表示分量绝对值和点态乘积，“sgn”表示为正负号函数；

(5)更新权重因子ω^k：

其中，ξ是一个正参数，它被调整以确保ω最优值；

(6)更新辅助变量q^k+1：

因此，其解可以表示为

(7)更新辅助变量r^k+1：

因此，其解可以表示为

(8)更新电导率分布g^k+1：

为了提高计算效率和提高解的稳定性，采用快速傅里叶变换求解g子问题：

其中，上标T是转置算子，

和A^TA是具有循环块结构的块循环矩阵；

(9)更新增广拉格朗日乘子γ₁ ^T,γ₂ ^T,和γ₃ ^T：在每次迭代中，均更新增广拉格朗日乘子，

和

(10)判断迭代是否符合迭代终止条件

或者k≤k_max，若是则迭代终止，进行下一步操作；若否，设置k＝k+1并跳回第(3)步，继续迭代求解；

(11)根据最终求解所得灰度值进行成像。

Images