CN106530367B - 一种基于Firm阈值迭代的电学层析成像稀疏重建方法 - Google Patents

Abstract

本发明涉及一种基于Firm阈值迭代的电学层析成像稀疏重建方法，用于非连续介质的电学层析成像，该方法将电学层析成像问题看作求解一个线性不适定问题，确定目标函数；重建方法为：根据被测场域，获取重建所需的相对边界测量值向量和灵敏度矩阵；设置初始化参数；进行加速预处理；根据上一循环中计算得到的正则化参数，即本次阈值迭代计算中所选择的阈值大小，对经过加速处理后的上一步迭代结果进行阈值迭代计算；计算本次迭代得到的解中非零元素的个数，更新稀疏度；判断迭代是否结束；根据最终求解所得灰度值，进行成像。本发明提高了算法的适用性和易用性，扩展了稀疏重建方法的应用范围。

Description

一种基于Firm阈值迭代的电学层析成像稀疏重建方法

技术领域

本发明属于电学层析成像技术领域，涉及利用稀疏正则化方法实现图像重建及阈值迭代求解方法。

背景技术

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

通过线性化处理，可以将电学层析成像中逆问题转化为线性逆问题并基于灵敏度矩阵来求解。针对逆问题求解的病态性，可选取正则化方法处理。正则化算法通用目标函数为：

其中，λ为正则化参数，用于调节误差项与正则化项在目标函数中的比重，G(x)为正则化项，根据具体应用的需要可以设计为不同的形式。在电学层析成像领域中，经典的Tikhonov算法，即以解的2范数作为正则化项的L₂正则化算法，目标函数为是由Vauhkonen等人在1998年发表于《IEEE医学成像》(IEEETransactions on Medical Imaging)第17卷，第285-293页，题为《基于电阻抗层析成像的Tikhonov正则化及先验信息选择》(Tikhonov regularization and prior informationin electrical impedance tomography)的文章中；近几年，随着压缩传感理论的发展，稀疏正则化算法重新受到多个领域的重视，成为了研究的热点，其中以解的1范数作为正则化项的L₁正则化算法应用最为广泛，其目标函数为：由Jin Bangti等人在2012 年发表于《工程中的数值计算》(International Journal For NumericalMethods In Engineering) 第89卷，第337-353页，题为《基于稀疏正则化的电阻抗层析成像重建算法》(A reconstruction algorithm for electrical impedance tomographybased on sparsity regularization)的文章中。采用 L₂正则化方法求解逆问题所得的解会出现过光滑现象，重建图像具有较大的伪影；而L₁正则化方法在图像的边界锐化和去伪影方面相比L₂正则化优势明显，但由于其目标函数的不可微性，求解起来相对困难，并且其对具有光滑物体分布的场域所求解时会出现过稀疏的问题，导致重建图像小于场域物体的真实大小。但目前L₁正则化算法已成为稀疏化重建算法中最受欢迎的方法之一，不断有新的求解算法出现，其中凸优化类算法在求解L₁正则化问题中应用最为广泛。

凸优化算法是将目标函数的求解问题转化为凸优化问题来进行求解的方法。常见的求解方法包括迭代阈值法、内点法、映射梯度法等。其中迭代阈值算法由于其计算简便、原理清晰而成为研究热点，但同时也以迭代速度慢著称。常用的迭代阈值算法包括迭代软阈值算法、迭代硬阈值算法以及一些基于迭代阈值算法的改进算法等。Amir Beck等人于2009年发表于《工业和应用数学学会杂志：成像科学》(SIAM Journal on ImagingSciences)第二卷，第 183-202页，题为《一种线性逆问题求解的快速迭代收缩阈值方法》(Afast iterative shrinkage-thresholding algorithm for linear inverse problems)提供了一种加速的软阈值迭代方法来求解L₁正则化问题，提高了求解速度。

张玲玲等人于2011年在发表于《世界科学与工程学术学会：电路和系统处理》(WSEAS Transactions on Circuits and Systems)第10卷，第393-402页，题为《一种用于电阻层析成像的快速迭代收缩阈值方法》(A Fast Iterative Shrinkage-ThresholdingAlgorithm for Electrical Resistance Tomography)的文章将该快速迭代阈值算法应用于电阻层析成像中，并对成像结果进行讨论，验证了方法的有效性。

然而在现有研究中，利用基于L₁正则化的软阈值迭代方法仍然如下缺点：(1)迭代速度慢；(2)重建图像的效果虽然在去伪影和边界锐化放面好于L₂正则化方法，但是依然存在所求解过于稀疏、重建图像小于场域物体的真实分布大小的问题；(3)正则化参数的选取与场域分布密切相关，当物体分布发生变化时，通常需要重新计算选择更合适的正则化参数才能保证所求解逼近真实解，而且正则化参数的选取对图像重建的质量有很大的影响。但在现有研究中，基于迭代软阈值的L₁正则化方法在进行电学层析成像逆问题求解时，通常在整个迭代过程始终采用一个固定的正则化系数(λ值)，在实际应用中，对于一些未知量经常变化的应用场景，重复计算正则化参数将使逆问题求解过程更加复杂，难以得到一个合适的正则化参数值。虽然目前很多学者已经针对正则化系数选取问题做了很多研究，但是并没有提出一种具有普遍适用性的最佳选取策略，多数正则化系数选取方法都只能在某些领域或某些情况下才能够达到较好地应用效果，导致算法的应用范围较小。所以研究的重点在于实现针对不同场域分布的正则化系数自适应选取，在每一步迭代中自适应地调节正则化参数，在一定范围内保证图像重建的质量不受初始正则化参数的影响，并且使重建图像更接近于实物分布。

发明内容

本发明的目的在于，基于一种新的阈值迭代方法，提出一个新的稀疏正则化目标函数，并且自适应地更新正则化参数，来解决电学层析成像的稀疏重建问题。本发明以Firm迭代阈值方法为基础，求解新的稀疏正则化目标函数，并根据场域分布，在每一步迭代计算中自适应地更新正则化参数值，解决稀疏正则化参数难以选取、求解计算速度较慢的问题，在不改变正则化参数的前提下，仍然可以在不同场域分布的图像重建中得到高质量的图像；同时也解决了L₁正则化图像重建结果过于稀疏、重建图像小于真实分布的问题，提高电学层析成像逆问题的求解精度和图像质量。本发明的技术方案如下：

一种基于Firm阈值迭代的电学层析成像稀疏重建方法，用于非连续介质的电学层析成像，该方法将电学层析成像问题看作求解一个线性不适定问题Ax＝b，其中A为灵敏度矩阵，b为相对边界测量值向量，x为需要求解的成像灰度值。

所建立的目标函数为：其中τ为初始人为选择的正则化系数，||·||为欧几里得范数，ρ是每一步迭代中实际采用的正则化参数的大小；是基于Firm阈值迭代函数的新的稀疏正则化项。重建方法包含有以下步骤：

(1)根据被测场域，获取重建所需的相对边界测量值向量b和灵敏度矩阵A；边界测量值的获取，通常是指将被测对象置于电学层析成像测量系统中，被测场域外均匀分布n个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压，相对边界测量值向量b为不含内含物的空场边界测量电压向量b1和含有内含物的有物场的边界测量电压向量b2之差。

(2)设置初始化参数：解初始值x₀＝A^Tb，稀疏度即非零元素个数K，初始正则化参数τ，常数t₁＝1，初始解x₀，迭代终止条件tol，解向量的元素总数N。

(3)对x_k进行加速预处理，方式为其中预处理常数 t₁＝1，t按照在每一步中进行迭代更新，在接下来的运算中用y_k代替x_k进行迭代，以实现加速迭代求解的作用。

(4)根据上一循环中计算得到的正则化参数ρ，即本次阈值迭代计算中所选择的阈值大小，对经过加速处理后的上一步迭代结果进行阈值迭代计算，代入迭代格式 x_k+1＝S_ρ,τ(y_k-A^T(Ay_k-b))中，其中S_ρ,τ的阈值迭代函数算子为：将ρ代入进行计算。

(5)计算本次迭代得到的解中非零元素的个数，更新稀疏度K，根据迭代解的稀疏度更新计算下一步迭代计算中使用的正则化参数ρ；正则化参数ρ的更新策略为：

(6)判断迭代是否符合迭代终止条件，若是则迭代终止，进行下一步操作；若否，则跳回第(3)步继续迭代求解；

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

本发明的有益效果是：基于一种新的稀疏正则化目标函数，在计算过程中选择与新目标函数相对应的Firm阈值函数进行迭代求解，同时在每一步迭代过程中自适应地根据解稀疏度来更新正则化参数，是一种新的图像稀疏重建方法。方法有效地克服了传统L₁正则化阈值迭代方法成像结果过于稀疏的缺点，提高了图像重建的精度；也解决了正则化参数难以选取的问题，提高了算法的适用性和易用性，扩展了稀疏重建方法的应用范围。

附图说明

图1为本发明的一种新的稀疏正则化目标函数的电学层析成像重建方法的流程框图；

图2为本发明的电阻层析成像系统圆形单截面被测场域及电极分布；

图3为本发明的实例，在选取一个小圆模型的真实分布时，传统的基于L₁正则化的软阈值迭代求解(FISTA)方法和基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法，分别在一组初始正则化参数条件下，图像重建结果的示意图，其中基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法即本专利发明的方法；

图4为本发明的实例，在选取两个小圆模型的真实分布时，传统的基于L₁正则化的软阈值迭代求解(FISTA)方法和基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法，分别在一组初始正则化参数条件下，图像重建结果的示意图，其中基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法即本专利发明的方法；

图5为本发明的实例中，在两个小圆模型中，传统的基于L₁正则化的软阈值迭代求解 (FISTA)方法和基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法，分别在一组初始正则化参数(τ₁～τ₆)条件下，分别所对应的迭代次数示意图，其中初始正则化参数τ₁～τ₆是从大到小的排序；

图中：

1、被测场域 2、电极

以上成像结果均在1％的随机噪声条件下测得。

具体实施方式

结合附图和实施例对本发明的一种新的稀疏正则化目标函数的电学层析成像重建方法加以说明。

本发明的一种新的稀疏正则化目标函数的电学层析成像重建方法，以稀疏正则化为基础，针对基于L₁正则化的软阈值迭代方法所得解过稀疏以及参数难以选取的问题，结合一种基于基于新的稀疏正则化目标函数的Firm阈值迭代方法，提出在迭代过程中利用前一步计算结果的稀疏度更新下一步迭代中的阈值，直至迭代终止的重建方法。

如图1所示，为本发明的一种基于新的稀疏正则化目标函数的电学层析成像重建方法流程图。如图2所示为电学层析成像之一的电阻层析成像系统圆形单截面被测场域及电极分布，采用16电极均匀分布在场域外壁。选取两种典型的非连续介质模型为实施例，场域内物体真实分布如图3左侧一竖列所示。为了较好地体现本发明中基于新的稀疏正则化目标函数(采用Firm阈值函数进行迭代计算并在迭代中自适应更新阈值大小)的方法与基于传统 L₁正则化(采用软阈值函数进行迭代计算并在迭代中不更新阈值大小)的不同，分别给出这两种算法在一组不同初始正则化参数条件下的求解结果。实施例包括如下具体步骤：

一种基于Firm阈值迭代的电学层析成像稀疏重建方法，适用于非连续介质的电学层析成像，该方法将电学层析成像问题看作求解一个线性不适定问题Ax＝b，其中A为灵敏度矩阵，b为相对边界测量值向量，x为需要求解的成像灰度值。

阈值迭代方法是梯度下降方法的一种拓展，其内在遵循的原理是一致的，均沿梯度的方向进行迭代计算。和传统Landweber迭代方法沿目标函数的梯度方向求解相比，由于目标函数中增加了一项正则化项，变成的形式，此时阈值迭代方法可以视为一步Landweber迭代之后再加上一个“阈值衰减”的步骤，沿着新的目标函数的梯度方向进行迭代求解。传统的阈值迭代方法的迭代格式为 x_k+1＝S_λ(x_k-tA^T(Ax_k-b))。

这种基于Firm阈值迭代方法的电学层析成像稀疏重建方法的目标函数为：其中τ为初始人为选择的正则化系数，||·||为欧几里得范数，ρ是每一步迭代中实际采用的正则化参数的大小。是基于Firm阈值迭代函数得到的新的稀疏正则化项。

重建方法包含有以下步骤：

(1)针对非连续介质的两个典型模型，分别获取各自重建所需的边界测量值和灵敏度矩阵：

边界测量值是将被测对象置于电学层析成像测量系统中，被测场域外均匀分布16个电极(如图2所示)，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下的边界电压，共获得208个测量值；逆问题右端项b为不含内含物的空场边界电压b1和含有内含物的有物场的边界测量电压b2之差(即右端项相对边界测量值b＝b1-b2)；

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

其中A_ij是第j个电极对对第i个电极对的灵敏度系数，φ_i,φ_j分别为第i个电极对及第 j个电极对在激励电流为I_i,I_j时场域电势分布，x,y分布为场域的横纵坐标；

(2)设置初始化参数：

设置的初始化参数包括：解初始值x₀＝A^Tb，初始预估稀疏度(非零元素个数) K＝N/5，正则化系数τ通常从开始从大到小选取以方便两种算法进行比较， t₁＝1，迭代终止条件tol＝2e^-2；

(3)对上一步得到的解x_k进行预处理，目的是为了加速收敛，公式为：

其中，t_k的更新策略是：

(4)对预处理之后的解进行阈值迭代计算，x_k+1＝S_ρ,τ(y_k-A^T(Ay_k-b))；其中

迭代公式的迭代格式为：

(5)计算非零元素的个数，更新稀疏度K，计算下一步迭代计算中使用的阈值ρ；阈值ρ的更新策略为：

(6)根据两步迭代所求解之间残差(||x_k+1-x_k||≤tol是否成立)来判断迭代是否完成，若是则迭代终止，进行下一步操作，若否，则跳回第(3)步继续求解；

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

图3为本发明的实例，在选取一个小圆模型的真实分布时，传统的基于L₁正则化的软阈值迭代求解(FISTA)方法和基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法，分别在一组初始正则化参数条件下，图像重建结果的示意图；图4为本发明的实例，在选取两个小圆模型的真实分布时，传统的基于L₁正则化的软阈值迭代求解(FISTA)方法和基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法，分别在一组初始正则化参数条件下，图像重建结果的示意图；图5为本发明的实例中，在两个小圆模型中，传统的基于L₁正则化的软阈值迭代求解(FISTA)方法和基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法，分别在一组初始正则化参数(τ₁～τ₆)条件下，分别所对应的迭代次数示意图，其中初始正则化参数τ₁～τ₆是从大到小的排序；

可以看出，两个模型中，在相同的正则化参数条件下，基于L₁正则化的软阈值迭代求解(FISTA)方法的解过于稀疏，不能很好地描述物体的大小，相比实物偏小；基于新的稀疏正则化目标函数的Firm阈值迭代(FIVTA)方法更加准确的描述了场域内物体的分布，提高了逆问题求解精度和图像重建质量。同时，在两个模型中还可以看出，随着正则化参数选择的不同，FIVTA(本算法)成像的结果稳定性较高，基本不受初始正则化参数选择的影响，在实际应用中可以较好地解决正则化参数难以选取的问题。并且，本算法的迭代次数和计算时间相比传统基于L₁正则化的软阈值迭代求解(FISTA)方法都有较大的提高。

以上所述实施例为本发明的几个较佳模型，本发明不局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都在本发明保护的范围。

Claims (1)

1.一种基于Firm阈值迭代的电学层析成像稀疏重建方法，用于非连续介质的电学层析成像，该方法将电学层析成像问题看作求解一个线性不适定问题Ax＝b，其中A为灵敏度矩阵，b为相对边界测量值向量，x为需要求解的成像灰度值；所建立的目标函数为：其中τ为初始人为选择的正则化系数，||·||为欧几里得范数，ρ是每一步迭代中实际采用的正则化参数的大小；是基于Firm阈值迭代函数的稀疏正则化项，重建方法包含有以下步骤：

(1)根据被测场域，获取重建所需的相对边界测量值向量b和灵敏度矩阵A；边界测量值的获取，通常是指将被测对象置于电学层析成像测量系统中，被测场域外均匀分布n个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压，相对边界测量值向量b为不含内含物的空场边界测量电压向量b1和含有内含物的有物场的边界测量电压向量b2之差；

(2)设置初始化参数：解初始值x₀＝A^Tb，稀疏度即非零元素个数K，初始正则化参数τ，常数t₁＝1，初始解x₀，迭代终止条件tol，解向量的元素总数N；

(3)对x_k进行加速预处理，方式为其中预处理常数t₁＝1，t按照在每一步中进行迭代更新，在接下来的运算中用y_k代替x_k进行迭代，以实现加速迭代求解的作用；

(4)根据上一循环中计算得到的正则化参数ρ，即本次阈值迭代计算中所选择的阈值大小，对经过加速处理后的上一步迭代结果进行阈值迭代计算，代入迭代格式x_k+1＝S_ρ,τ(y_k-A^T(Ay_k-b))中，其中S_ρ,τ为阈值迭代函数算子，表示为：将ρ代入进行计算；

(5)计算本次迭代得到的解中非零元素的个数，更新稀疏度K，根据迭代解的稀疏度更新计算下一步迭代计算中使用的正则化参数ρ；正则化参数ρ的更新策略为：

(6)判断迭代是否符合迭代终止条件，若是则迭代终止，进行下一步操作；若否，则跳回第(3)步继续迭代求解；

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