CN104535294A - 一种基于二阶微分的修正l曲线电学层析成像重建方法 - Google Patents

Abstract

本发明提供一种基于二阶微分的修正L曲线电学层析成像重建方法，适用于泡状流层析成像，包括：根据被测场域，获取重建所需的相对边界测量值向量b和灵敏度矩阵A；利用Tikhonov正则化，计算并绘制L-曲线；判断L-曲线是否存在局部拐点；若不存在局部拐点，则通过L-曲线法确定优化选取的正则化系数；若存在局部拐点，通过修正的L-曲线法确定优化选取的正则化系数，方法为：通过计算L-曲线二阶微分的第二大峰值确定曲线的局部拐点，并将此局部拐点对应的正则化系数作为优化选取的系数；将正则化系数代入Tikhonov正则化中进行图像重建逆问题求解；成像。本发明有利于电学层析成像逆问题的精确求解，提高了图像重建质量。

Description

一种基于二阶微分的修正L曲线电学层析成像重建方法

技术领域

本发明属于电学层析成像技术领域，涉及利用Tikhonov正则化方法实现图像重建，特别是基于L-曲线实现Tikhonov正则化系数的优化选取。

背景技术

多相流在工业生产与科学研究中有着十分重要的作用，它们经常出现在动力、化工、石油、核能、冶金工程等工业过程中，涉及安全与经济问题，对其精确测量也给工程和科研人员提出了挑战。管道中的多相流呈现出的几何与动力特征各异的流动形态称为流型，它可通过组分或相的形态来定性描述，两相流中常见的流型包括泡状流、弹状流、环状流等。

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

电学层析成像逆问题(即图像重建问题)求解具有非线性。通过线性化处理，可以将问题转化为线性逆问题求解。针对逆问题求解的不适定性，通常选取正则化方法处理逆问题，例如W Q Yang等人2002年发表于《测量科学与技术》(Measurement Science&Technology)第14卷，第R1-R13页，题为《电容层析成像图像重建算法》(Image reconstruction algorithms for electrical capacitancetomography)的综述文章介绍的部分利用正则化实现图像重建的方法。在利用正则化方法求解过程中，存在一个不可忽略的问题就是正则化系数的选取，一个合理的正则化系数直接关系到方法是否能够收敛到问题的真实解以及收敛速度。

目前针对正则化系数选取方法的研究已经十分广泛，学者们提出了各种方法用于正则化系数的选取，例如：偏差原则、广义交叉验证、残差法、L-曲线法等。其中L-曲线法作为最早提出的正则化系数选取方法之一，获得了广泛的应用。

C.L.LAWSON和R.J.HANSON在1974年出版于普伦蒂斯·霍尔出版社的《最小二乘问题的求解》(Solving Least Squares Problems)一书中第25和26章首次介绍了L-曲线方法；P.C.Hansen在1992年发表于《工业与应用数学》(Society for Industrial and Applied Mathematics)第34卷，第561-580页，题为《基于L-曲线法的离散不适定问题的分析》(Analysis of discrete ill–posed problems bymeans of the L–curve)将L-曲线法应用到逆问题求解中。自此展开了基于L-曲线的正则化系数选取方法的研究。

在传统L-曲线法中，选取曲线的全局拐点对应的正则化系数作为最优的正则化系数，研究全局拐点确定方法的文献有：P.C.Hansen等人1993年发表于《工业与应用数学》(Society for Industrial andApplied Mathematics)第14卷，第1487-1503页，题为《L-曲线在离散不适定问题正则化方法中的应用》(The use of the L-curve in the regularization of discrete ill-posed problems)的文章；J.Lian等人1998年发表于《医学和生物工程协会第20届国际会议》(Engineering in Medicine andBiology Society,Proceedings of the 20th Annual International Conference of the IEEE)第4卷，第2155-2158页，题为《一种新的实现皮层电位成像的正则化的方法》(A new method forimplementation of regularization in cortical potential imaging)的文章；J.L.Castellanos等人2002年发表于《应用数值数学》(Applied Numerical Mathematics)第43卷，第359-373页，题为《寻找L-曲线拐点的三角形法》(The triangle method for finding the corner of the L-curve)的文章；P.C.Hansen等人2007年发表于《计算和应用数学杂志》(Journal of Computational and AppliedMathematics)第198卷，第483-492页，题为《离散L-曲线准则自适应修剪算法》(An adaptive pruningalgorithm for the discrete L-curve criterion)的文章，等等。

但是传统的L-曲线法存在一定的局限性，例如M.Hanke在1996年发表于《BIT计算数学》(BITNumerical Mathematics)第36卷，第287-301页，题为《不适定问题L-曲线法的局限性》(Limitationsof the L-curve Method in Ill-Posed Problems)的文章给出了L-曲线法局限性的证明。S.Morigi等人2001年发表于《应用数学与计算》(Applied Mathematics and Computation)第121卷，第55-73页，题为《欠定线性系统的L曲线正则化分块计算法》(A regularizing L-curve Lanczos method forunderdetermined linear system)；L.Reichel等人2008年发表于《计算和应用数学杂志》(Journal ofComputational and Applied Mathematics)第219卷，第493-508页，题为《解决不适定问题的新L-曲线法》(A new L-curve for ill-posed problems)；M.Rezghi等人2009年发表于《计算和应用数学杂志》(Journal of Computational and Applied Mathematics)第231卷，第914-924页，题为《基于Tikhonov正则化的L-曲线法的新的变形方法》(A new variant of L-curve for Tikhonovregularization)等文章针对L-曲线计算量和L-曲线的变形问题给出了相应的改进策略。

现有研究中，以L-曲线为基础进行正则化系数选取的方法及其改进方法都是将曲线的全局拐点对应的系数值作为优化选取的正则化系数，但是这种选取方法应用到Tikhonov正则化参数选取中，再进行电学层析成像逆问题重建时，并不能保证对所有被测介质分布均适用，也就是说，在某些被测介质分布模型下，选取曲线的全局拐点对应的系数作为最优的正则化系数的L-曲线法在图像重建求解时会失效。

发明内容

本发明的目的在于提出一种基于二阶微分的修正L曲线电学层析成像重建方法，方法以L-曲线为基础，解决传统L-曲线法选取的系数在Tikhonov正则化进行泡状流电学层析成像逆问题求解时的不适用问题，提高电学层析成像逆问题的求解精度和图像重建质量。本发明的技术方案如下：

一种基于二阶微分的修正L曲线电学层析成像重建方法，适用于泡状流层析成像，该方法将电学层析成像问题看作一个线性不适定问题Ax＝b。其中，A为灵敏度矩阵，b为相对边界测量值向量，x为所求成像灰度值，包含有以下步骤：

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

(2)利用Tikhonov正则化，计算并绘制L-曲线；

(3)判断L-曲线是否存在局部拐点；

(4)若不存在局部拐点，则通过传统L-曲线法确定优化选取的正则化系数；若存在局部拐点，通过修正的L-曲线法确定优化选取的正则化系数，方法为：通过计算L-曲线二阶微分的第二大峰值确定曲线的局部拐点，并将此局部拐点对应的正则化系数作为优化选取的系数；

(5)将正则化系数代入Tikhonov正则化中进行图像重建逆问题求解；

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

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

本发明的有益效果是基于L-曲线提出了一种基于二阶微分的修正L曲线电学层析成像重建方法，对传统L-曲线进行了补充修正，解决了传统L-曲线法选取的正则化系数在泡状流电学层析成像逆问题求解时的不适用问题，丰富了基于L-曲线选取正则化系数的策略，拓宽了L-曲线在正则化系数选取问题上的模型适用性，提高了泡状流电学层析成像逆问题的求解精度和图像重建质量。

附图说明

图1为本发明的一种基于二阶微分的修正L曲线电学层析成像重建方法的流程框图；

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

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

图4为本发明的实例中三个典型模型的L-曲线及曲线上的全局拐点和局部拐点：其中(a-c)分别对应图3中的模型(a-c)；

图5为本发明的实例中三个模型在传统的L-曲线法选取的正则化系数下Tikhonov正则化成像结果示意图：其中(a-c)分别对应图3中的模型(a-c)；

图6为本发明的实例中三个模型在修正的L-曲线法选取的正则化系数下Tikhonov正则化成像结果示意图：其中(a-c)分别对应图3中的模型(a-c)；

图中：

1、被测场域    2、电极

具体实施方式

结合附图和实施例对本发明的一种基于二阶微分的修正L曲线电学层析成像重建方法加以说明。

本发明的一种基于二阶微分的修正L曲线电学层析成像重建方法，针对传统L-曲线法不适用的模型对应的L-曲线存在局部拐点的特点，利用曲线拐点处二阶微分存在突变的性质，通过计算曲线二阶微分的第二大峰值确定曲线的局部拐点，并将此局部拐点对应的系数作为优化选取的Tikhonov正则化系数，完成最终的逆问题求解。

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

(1)针对泡状流的三个典型模型，分别获取各自重建所需的边界测量值和灵敏度矩阵：

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

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

$A_{\mathrm{ij}}=-\∫\frac{{\▿\mathrm{\φ}}_i}{I_i}\·\frac{{\▿\mathrm{\φ}}_j}{I_j}\mathrm{dxdy}$

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

(2)利用Tikhonov正则化，分别计算并绘制相应的L-曲线：

Tikhonov正则化求解逆问题的正则化解为：

x_λ＝(A^TA+λ²I)^-1A^Tb

其中A为灵敏度矩阵，b为相对边界测量值，λ为正则化系数，I为单位矩阵。L-曲线为正则化解对应的残差范数与解范数的关系曲线，一般绘制在对数坐标下，即如图4所示，其中L-曲线(a-c)分别对应模型(a-c)；

(3)正则化系数的优化选取：

为了更好地体现传统L-曲线法与修正的L-曲线法的不同，对实例中的三个模型分别按照传统L-曲线法和修正的L-曲线法进行参数的优化选取：

传统的L-曲线法：计算曲线的曲率

$\hat{c}=\frac{{\widehat{\mathrm{\ρ}}}^{\′}{\widehat{\mathrm{\ξ}}}^{\′\′}-{\widehat{\mathrm{\ρ}}}^{\′\′}{\widehat{\mathrm{\ξ}}}^{\′}}{{({\left({\widehat{\mathrm{\ρ}}}^{\′}\right)}^2+{\left({\widehat{\mathrm{\ξ}}}^{\′}\right)}^2)}^{3/2}}$

其中 $\widehat{\mathrm{\ρ}}=\log \left({\left|\right|{\mathrm{Ax}}_{\mathrm{\λ}}-b\left|\right|}_2^2\right),\widehat{\mathrm{\ξ}}=\log \left({\left|\right|x_{\mathrm{\λ}}\left|\right|}_2^2\right),$ 符号'和”表示符号下变量对正则化系数λ的一阶微分和二阶微分，选择最大曲率点为曲线的全局拐点，并将此全局拐点对应的系数作为传统L-曲线法优化选取的正则化系数；

修正的L-曲线法：计算L-曲线的二阶微分

$\hat{k}=\frac{d^2\widehat{\mathrm{\ξ}}}{d{\widehat{\mathrm{\ρ}}}^2}=\frac{{\widehat{\mathrm{\ρ}}}^{\′}{\widehat{\mathrm{\ξ}}}^{\′\′}-{\widehat{\mathrm{\ρ}}}^{\′\′}{\widehat{\mathrm{\ξ}}}^{\′}}{{\left({\widehat{\mathrm{\ρ}}}^{\′}\right)}^3}$

其 $\widehat{\mathrm{\ρ}}=\log \left({\left|\right|{\mathrm{Ax}}_{\mathrm{\λ}}-b\left|\right|}_2^2\right),\widehat{\mathrm{\ξ}}=\log \left({\left|\right|x_{\mathrm{\λ}}\left|\right|}_2^2\right),$ '和”表示变量对正则化系数λ的一阶微分和二阶微分，选择二阶微分的第二大峰值点为曲线的局部拐点，并将此拐点对应的系数作为修正的L-曲线法优化选取的正则化系数；

(4)将优化选取的正则化系数代入Tikhonov正则化中进行图像重建逆问题求解：

将步骤(3)中通过传统L-曲线法和修正的L-曲线法获得的正则化系数分别代入Tikhonov正则化求解公式中，计算逆问题的正则化解；

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

将步骤(4)计算所得正则化解对应到电阻层析成像的像素点上，进行灰度成像。图5所示为传统L-曲线法选取参数的Tikhonov正则化计算结果成像图，图6所示为修正的L-曲线法选取参数的Tikhonov正则化计算结果成像图，图中(a-c)分别与模型(a-c)对应。

可以看出，传统L-曲线法在这些模型中不适用，所求解偏离真实值，不能成像；而本发明提出的修正的L-曲线法在这些模型中具有较好的适用性，能够基本描绘出被测场域内部分布，提高了电学层析成像逆问题求解精度和图像重建质量。

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

Claims (2)

1.一种基于二阶微分的修正L曲线电学层析成像重建方法，适用于泡状流层析成像，该方法将电学层析成像问题看作一个线性不适定问题Ax＝b。其中，A为灵敏度矩阵，b为相对边界测量值向量，x为所求成像灰度值，包含有以下步骤：

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

(2)利用Tikhonov正则化，计算并绘制L-曲线；

(3)判断L-曲线是否存在局部拐点；

(4)若不存在局部拐点，则通过L-曲线法确定优化选取的正则化系数；若存在局部拐点，通过修正的L-曲线法确定优化选取的正则化系数，方法为：通过计算L-曲线二阶微分的第二大峰值确定曲线的局部拐点，并将此局部拐点对应的正则化系数作为优化选取的系数；

(5)将正则化系数代入Tikhonov正则化中进行图像重建逆问题求解；

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

2.根据权利要求1所述的电学层析成像重建方法，其特征在于：所述步骤(1)边界测量值的获取，是指将被测对象置于电学层析成像测量系统中，被测场域外均匀分布n个电极，采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压，相对边界测量值向量b为不含内含物的空场边界测量电压向量b1和含有内含物的有物场的边界测量电压向量b2之差。