CN104634829B - 基于p向量等比收缩的电学层析成像Lp正则化重建方法 - Google Patents

Abstract

本发明涉及一种基于p向量等比收缩的电学层析成像Lp正则化重建方法，适用于泡状流层析成像，利用Gauss‑Newton迭代进行Lp正则化逆问题求解的每步迭代中根据所得解更新由图像中各个像素点上的p值构成的p向量，得到具有场域物体空间分布特性的p分布，最终完成计算获取重建图像，步骤如下：获取重建所需的相对边界测量值向量b和灵敏度矩阵A；建立Lp正则化的目标函数；计算等比收缩因子；利用Gauss‑Newton迭代公式进行求解；在每次迭代中，利用所求解更新p向量；成像。本发明有利于电学层析成像逆问题的精确求解，提高图像重建质量。

Description

基于p向量等比收缩的电学层析成像Lp正则化重建方法

技术领域

本发明属于电学层析成像技术领域，涉及利用Lp正则化方法实现图像重建的方法。

背景技术

多相流指包含明显分界面的流体系统，如含气泡(液滴)的液体(气体)、不混溶的液体、含固体颗粒的气体或液体等，它们经常出现在动力、化工、石油、核能、冶金工程等过程中，对工业生产与科学研究有着十分重要的作用。多相流的流型指的是其管道中呈现出的几何与动力特征各异的流动形态，它可通过组分或相的形态来定性描述，两相流中常见的流型包括泡状流、弹状流、环状流等。

电学层析成像技术(Electrical Tomography，ET)是自上世纪80年代后期出现的一种新的基于电特性敏感机理的过程层析成像技术，它的物理基础是不同的媒质具有不同的电特性(电导率/介电系数/复导纳/磁导率)，通过判断敏感场内物体的电特性分布便可推知该场中媒质的分布情况。电学层析成像技术主要包括电阻层析成像(ElectricalResistance Tomography,ERT)、电容层析成像(Electrical Capacitance Tomography,ECT)、电阻抗层析成像(Electrical Impedance Tomography,EIT)和电磁层析成像(Electrical Magnetic Tomography,EMT)。电学层析成像在多相流及生物医学领域有广泛的应用前景，可以实现长期、持续监测。

电学层析成像逆问题(即图像重建问题)求解具有非线性。通过线性化处理，可以将问题转化为线性逆问题求解。针对逆问题求解的不适定性，通常选取正则化方法处理逆问题。正则化方法的思想是寻找一个由先验信息约束的稳定解集来逼近真实解。先验信息的选取不同和正则化函数形式的不同使得正则化方法具有不同的应用形式，例如以解的2范数为正则化函数实现逆问题的稳定求解的L2正则化方法：Vauhkonen·M等人在1998年发表于《IEEE医学成像》(Medical Imaging,IEEE Transactions)第17卷，第285-293页，题为《基于电阻抗层析成像的Tikhonov正则化及先验信息选择》(Tikhonov regularizationand prior information in electrical impedance tomography)的文章；以解的1范数为正则化函数实现逆问题稳定求解的L1正则化方法：Jin,Bangti等人在2012年发表于《工程中的数值计算》(International Journal For Numerical Methods In Engineering)第89卷，第337-353页，题为《基于稀疏正则化的电阻抗层析成像重建算法》(A reconstructionalgorithm for electrical impedance tomography based on sparsityregularization)的文章。

但是采用L2正则化求解逆问题所得解会出现过光滑现象，所成图像具有较大的尾影；而L1正则化对具有光滑物体分布的场域所求解会出现过稀疏的问题，不能充分体现场域物体的真实大小，因此引入以p范数为正则化函数的Lp正则化方法(一般在电学层析成像中取p∈[1,2])。Daubechies I等人于2004年发表于《数学与应用数学》(Communicationson Pure and Applied Mathematics)第57卷，第1413-1457页，题为《针对线性逆问题稀疏约束的迭代阈值算法》(An iterative thresholding algorithm for linear inverseproblems with a sparsity constraint)的文章提供了求解Lp正则化的迭代算法。张玲玲等人于2013年发表于《多相流检测与仪器仪表》(Flow Measurement andInstrumentation)第33卷，第244-250页，题为《电阻层析成像逆问题迭代阈值算法》(Aniterative thresholding algorithm for the inverse problem of electricalresistance tomography)的文章将迭代阈值算法应用于电阻层析成像中，并对p＝1.5时的成像结果进行讨论，验证了方法的有效性。

然而现有研究中，利用Lp正则化进行电学层析成像逆问题求解通常对整个场域选择一个固定的p值，且不同物体分布的场域需要给定不同的p值，以获得更好的稳定解。这种方法忽略了不同物体分布的场域的空间特性，不能突出场域自身的特点，且p值的调节给问题的求解带来额外的工作量，不利于方法的推广。

发明内容

本发明的目的在于克服现有技术的上述不足，提出一种电学层析成像Lp重建方法。本发明以Lp正则化为基础，结合Gauss-Newton迭代，解决了L2正则化解过光滑而L1正则化解过稀疏的问题，提高电学层析成像逆问题的求解精度和图像重建质量。本发明的技术方案如下：

一种基于p向量等比收缩的电学层析成像Lp正则化重建方法，适用于泡状流层析成像，该方法将电学层析成像问题看作一个线性不适定问题Ax＝b，其中A为灵敏度矩阵，b为相对边界测量值向量，x为与场域物质电特性分布对应的成像灰度值向量，称其为解向量，采用Lp正则化逆问题求解的方法进行图像重建，其特征在于，

利用Gauss-Newton迭代进行Lp正则化逆问题求解的每步迭代中根据所得解更新由图像中各个像素点上的p值构成的p向量，得到具有场域物体空间分布特性的p分布，最终完成计算获取重建图像，步骤如下：

(1)根据对被测场域的测量，获取重建所需的相对边界测量值向量b和灵敏度矩阵A；

(2)建立Lp正则化的目标函数；

(3)设置初始化参数：解向量x的初始值x₀＝0，p向量初始值p₀＝2，p向量终止值p_stop＝1；设定迭代次数N；

(4)计算等比收缩因子

(5)利用Gauss-Newton迭代公式进行求解；

(6)利用所求解更新p向量：判断解向量中各个元素是否为零，若是则对应像素点的p值保持前一步p值不变；若否，则对应像素点的p值更新为前一步p值与收缩因子的乘积；

(7)判断迭代是否完成，若是则迭代终止，进行下一步操作，若否，则跳回第(5)步继续求解；

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

作为优选实施方式，所述的Lp正则化的目标函数为：其中λ＞0是正则化系数，||·||为欧几里得范数，p向量满足任一元素属于[1,2]；考虑到p＝1时目标函数不可微，利用近似上述的目标函数，其中n为解向量x的维数，j为从1到n的计数整数，x_j为解向量x中第j个元素，β是微小的可调参数，满足β＞0。

利用Gauss-Newton迭代公式为：

x_k＝x_k-1-[▽²F(x_k-1)]^-1▽F(x_k-1)

其中k是当前的迭代次数，满足1≤k≤N；x_k是第k次迭代得到的解，x_k-1是第(k-1)次迭代得到的解；▽F(x_k-1)为当x＝x_k-1时目标函数的一阶微分，▽²F(x_k-1)为x＝x_k-1时目标函数的二阶微分，分别通过下面两个公式获得：

其中p_k-1是前一步即第(k-1)次迭代更新获得的p向量；diag(·)是通过向量构造对角阵的函数，向量的各个元素构成对角阵的对角元素，对角阵的非对角元素均为零。

本发明的有益效果是基于利用Gauss-Newton迭代求解Lp正则化的过程，提出一种利用迭代过程中p的等比收缩实现Lp自适应的电学层析成像图像重建方法，为重建图像的光滑区域每个像素点提供值为2或接近为2的p值；为重建图像有物体区域每个像素点提供值为1或接近为1的p值，获得了带有场域物体空间分布特性的p分布，完成最终求解。本发明有效地克服了L2正则化和L1正则化的缺点，提高了逆问题求解精度和图像重建质量，具有更强的鲁棒性和更大的适用性。

附图说明

图1为本发明的一种基于p等比收缩的电学层析成像Lp自适应重建方法的流程框图；

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

图3为本发明的实例选取的三个模型的真实分布：(a)为一个小圆模型(b)为两个圆模型(c)为三个小圆模型

图4为本发明的实例中三个模型在L2正则化求解下的成像结果示意图：其中(a-c)分别对应图3中的模型(a-c)；

图5为本发明的实例中三个模型在L1正则化求解下的成像结果示意图：其中(a-c)分别对应图3中的模型(a-c)；

图6为本发明的实例中三个模型在本方法提出的自适应Lp正则化求解下的成像结果示意图：其中(a-c)分别对应图3中的模型(a-c)；

图中：

1、被测场域 2、电极

具体实施方式

结合附图和实施例对本发明的一种基于p等比收缩的电学层析成像Lp自适应重建方法加以说明。

本发明的基于p等比收缩的电学层析成像Lp自适应重建方法，以Lp正则化为基础，针对L2正则化求得的解过光滑而L1正则化求得的解过稀疏的问题，结合Gauss-Newton迭代公式，提出在迭代过程中利用前一步计算结果由图像中各个像素点上的p值构成的p向量，并利用更新后的p向量完成当前计算，直至迭代终止的重建方法。

如图1所示，为本发明的一种基于p等比收缩的电学层析成像Lp自适应重建方法流程图。如图2所示为电学层析成像之一的电阻层析成像系统圆形单截面被测场域及电极分布，采用16电极均匀分布在场域外壁。选取三个典型的泡状流模型为实施例，场域内物体真实分布如图3(a-c)所示。为了更好地体现本发明中自适应Lp正则化与L2正则化和L1正则化的不同，分别给出三个模型在这三种正则化方法下的求解结果。实施例包括如下具体步骤：

一种基于p等比收缩的电学层析成像Lp自适应重建方法，适用于泡状流层析成像，该方法将电学层析成像问题看作一个线性不适定问题Ax＝b，其中A为灵敏度矩阵，b为相对边界测量值向量，x为与场域物质电特性分布对应的成像灰度值向量。

Lp正则化的目标函数为：

其中λ＞0为正则化系数，||·||为欧几里得范数，p为由重建图像各个像素点上的p值构成的常数向量。考虑p＝1时目标函数不可微，利用新的目标函数：

近似原目标函数。重建方法包含有以下步骤：

(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₀＝0，p向量初始值p₀＝2，p向量终止值p_stop＝1；根据经验选取的迭代次数N＝5，正则化系数λ＝1×10^-4，微调参数β＝1×10^-16；

(3)计算等比收缩因子γ，计算公式为：

(4)利用Gauss-Newton迭代公式进行求解，其迭代格式为：

x_k＝x_k-1-[▽²F(x_k-1)]^-1▽F(x_k-1)

其中k是当前步的迭代次数，满足1≤k≤N；x_k是第k次迭代得到的解，x_k-1是第(k-1)次迭代得到的解；▽F(x_k-1)为当x＝x_k-1时目标函数的一阶微分，▽²F(x_k-1)为x＝x_k-1时目标函数的二阶微分，分别通过下面两个公式获得：

其中p_k-1是前一步即第(k-1)次迭代更新获得的p向量；diag(·)是通过向量构造对角阵的函数，向量的各个元素构成对角阵的对角元素，对角阵的非对角元素均为零。且计算二阶微分▽²F(x_k-1)时省略了高阶项；

(5)确定迭代过程中的p向量：

为了更好地比较本方法与L2正则化和L1正则化的不同，对实施例中的三个模型分别给出L2正则化、L1正则化和自适应Lp正则化所需的p值：

对L2正则化，迭代过程中p向量保持所有元素满足p＝2不变；

对L1正则化，迭代过程中p向量保持所有元素满足p＝1不变；

对自适应Lp正则化，迭代过程中利用所求解更新p向量，更新策略为：

其中p_k是第k次迭代更新的p向量；l是场域各个像素点的位置；

(6)判断迭代是否完成，若是则迭代终止，进行下一步操作，若否，则跳回步骤(4)继续求解；

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

将所得正则化解对应到电阻层析成像的像素点上，进行灰度成像。图4所示为L2正则化计算结果成像图，图5所示为L1正则化计算结果成像图，图6为自适应Lp正则化计算结果成像图，图中(a-c)分别与模型(a-c)对应。

可以看出，在相同的参数设置下，L2正则化的解过光滑，所成图像尾影大；而L1正则化过稀疏，不能很好地描述物体的大小；自适应Lp正则化的结果介于L2正则化和L1正则化所得结果中间，提高了逆问题求解精度和图像重建质量，更加准确的描述了场域内物体的分布。

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

Claims (2)

1.一种基于p向量等比收缩的电学层析成像Lp正则化重建方法，适用于泡状流层析成像，该方法将电学层析成像问题看作一个线性不适定问题Ax＝b，其中A为灵敏度矩阵，b为相对边界测量值向量，x为与场域物质电特性分布对应的成像灰度值向量，称其为解向量，采用Lp正则化逆问题求解的方法进行图像重建.其特征在于，

利用Gauss-Newton迭代进行Lp正则化逆问题求解的每步迭代中根据所得解更新由图像中各个像素点上的p值构成的p向量，得到具有场域物体空间分布特性的p分布，最终完成计算获取重建图像，步骤如下：

(1)根据对被测场域的测量，获取重建所需的相对边界测量值向量b和灵敏度矩阵A；

(2)建立如下的Lp正则化的目标函数：

其中λ＞0是正则化系数，||·||为欧几里得范数，p向量满足任一元素属于[1,2]；考虑到p＝1时目标函数不可微，利用近似上述的目标函数，其中n为解向量x的维数，j为从1到n的计数整数，x_j为解向量x中第j个元素，β是微小的可调参数，满足β＞0；

(3)设置初始化参数：解向量x的初始值x₀＝0，p向量初始值p₀＝2，p向量终止值p_stop＝1；设定迭代次数N；

(4)计算等比收缩因子

(5)利用Gauss-Newton迭代公式进行求解；

(6)利用所求解更新p向量：判断解向量中各个元素是否为零，若是则对应像素点的p值保持前一步p值不变；若否，则对应像素点的p值更新为前一步p值与收缩因子的乘积；

(7)判断迭代是否完成，若是则迭代终止，进行下一步操作，若否，则跳回第(5)步继续求解；

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

2.根据权利要求1所述的基于p向量等比收缩的电学层析成像Lp正则化重建方法，其特征在于：利用Gauss-Newton迭代公式为：

x_k＝x_k-1-[▽²F(x_k-1)]^-1▽F(x_k-1)

其中k是当前的迭代次数，满足1≤k≤N；x_k是第k次迭代得到的解，x_k-1是第(k-1)次迭代得到的解；▽F(x_k-1)为当x＝x_k-1时目标函数的一阶微分，▽²F(x_k-1)为x＝x_k-1时目标函数的二阶微分，分别通过下面两个公式获得：

$\▿F\left(x_{k-1}\right)=A^T({\mathrm{Ax}}_{k-1}-b)+\mathrm{\λ}diag\left(p_{k-1}{\left(\sqrt{x_{k-1}^2+\mathrm{\β}}\right)}^{p_{k-1}-2}{x}_{k-1}\right)$

${\▿}^{2}F\left(x_{k-1}\right)=A^{T}A+\mathrm{\λ}diag\left(p_{k-1}{\left(\sqrt{x_{k-1}^2+\mathrm{\β}}\right)}^{p_{k-1}-2}\right)$

其中p_k-1是前一步即第(k-1)次迭代更新获得的p向量；diag(·)是通过向量构造对角阵的函数，向量的各个元素构成对角阵的对角元素，对角阵的非对角元素均为零。