CN104614411A - Electrical tomography Lp regularized reconstruction method based on p vector equal difference decline - Google Patents

Abstract

The invention relates to an electrical tomography Lp regularized reconstruction method based on p vector equal difference decline and is applicable to bubble flow tomography. Lp regularized inverse problem solution is performed in each step of iteration by utilizing Gauss-Newton iteration, a p vector consisting of p values on each pixel point in the image is updated according to the obtained solution, p distribution with field domain object space distribution characteristics is obtained, and finally, calculation is finished to acquire a reconstruction image. The method comprises the following steps: acquiring a relative boundary measurement value vector b and a sensitivity matrix A needed by reconstruction; establishing an Lp regularized objective function; calculating an equal difference decline factor; solving by utilizing a Gauss-Newton iteration formula; updating the p vector in each step of iteration by utilizing the solution; and imaging. According to the method disclosed by the invention, accurate solution of the electrical tomography inverse problems is promoted, and the image reconstruction quality is improved.

Description

Based on the electricity tomography Lp regularization reconstruction method that p vector equal difference declines

Technical field

The invention belongs to electricity chromatography technical field of imaging, relate to and utilize Lp regularization method to realize image reconstruction and Gauss-Newton alternative manner.

Background technology

Polyphasic flow refers to comprise obvious interfacial fluid system, as the liquid (gas) of bubbles (drop), immiscible liquid, gas or liquid etc. containing solid particle, they often appear in the processes such as power, chemical industry, oil, nuclear energy, metallurgical engineering, have very important effect to commercial production and scientific research.The flow pattern of polyphasic flow refers to the geometry that presents in its pipeline and the different nowed forming of dynamic characteristic, and it carrys out qualitative description by the form of component or phase, and flow pattern common in two-phase flow comprises bubble flow, slug flow, annular flow etc.

Electricity chromatography imaging technique (Electrical Tomography, ET) be a kind of process tomographic imaging technology based on electrical characteristics sensitive mechanism newly occurred from the later stage eighties in last century, its physical basis is that different mediums has different electrical characteristics (conductivity/dielectric coefficient/complex admittance/magnetic permeability), by judging that the electrical characteristics distribution of object in sensitivity field just can know the distribution situation of this middle medium by inference.Electricity chromatography imaging technique mainly comprises Electrical Resistance Tomography (Electrical Resistance Tomography, ERT), capacitance chromatography imaging (Electrical Capacitance Tomography, ECT), electrical impedance tomography (Electrical Impedance Tomography, and electromagnetic chromatographic (Electrical Magnetic Tomography, EMT) EIT).Electricity chromatography is imaged on polyphasic flow and biomedical sector has wide practical use, and can realize for a long time, continue to monitor.

Electricity tomography inverse problem (i.e. image reconstruction problem) solve have non-linear.By linearization process, problem can be converted into linear inverse problem solving.For the ill-posedness of reverse temperature intensity, usually choose regularization method process inverse problem.The thought of regularization method is that the stable disaggregation that searching one is retrained by prior imformation carrys out approaching to reality solution.The difference that is different and Regularization function form of choosing of prior imformation makes regularization method have different application forms, such as with 2 norms of separating for Regularization function realizes the stable L2 regularization method solved of inverse problem: the people such as VauhkonenM were published in " IEEE medical imaging " (Medical Imaging in 1998, IEEE Transactions) the 17th volume, 285-293 page, be entitled as the article of " Tikhonov regularization and prior imformation selection based on electrical impedance tomography " (Tikhonov regularization and prior information in electrical impedance tomography), stablize for Regularization function realizes inverse problem the L1 regularization method solved: Jin with 1 norm of separating, the people such as Bangti were published in " numerical evaluation in engineering " (International Journal For Numerical Methods In Engineering) the 89th volume in 2012,337-353 page, is entitled as the article of " the electrical impedance tomography reconstruction algorithm based on sparse regularization " (A reconstruction algorithm for electrical impedance tomography based on sparsity regularization).

But adopt L2 regularization to solve inverse problem gained solution and there will be smooth phenomenon, become image to have larger tail shadow; And L1 regularization to have smooth object distribution field domain solve and there will be sparse problem, the actual size of field domain object can not be fully demonstrated, therefore the Lp regularization method (generally getting p ∈ [1,2] in electricity tomography) that to introduce with p norm be Regularization function.The people such as Daubechies I were published in " Mathematics and Applied Mathematics " (Communications on Pure and Applied Mathematics) the 57th volume in 2004,1413-1457 page, the article being entitled as " the iteration threshold algorithm of linear inverse problem sparse constraint " (An iterative thresholding algorithm for linear inverse problems with a sparsity constraint) provides the iterative algorithm solving Lp regularization.The people such as Zhang Lingling were published in 2013 " polyphasic flow detect and instrument and meter " (Flow Measurement and Instrumentation) the 33rd volume, 244-250 page, the article being entitled as " Electrical Resistance Tomography inverse problem iteration threshold algorithm " (An iterative thresholding algorithm for the inverse problem of electrical resistance tomography) by iteration threshold algorithm application in Electrical Resistance Tomography, and imaging results during p=1.5 is discussed, demonstrate the validity of method.

But in existing research, utilize Lp regularization to carry out electricity tomography reverse temperature intensity and usually a fixing p value is selected to whole field domain, and the field domain of different objects distribution needs given different p value, to obtain better stable solution.This method have ignored the spatial character of the field domain of different objects distribution, can not give prominence to the feature of field domain self, and the adjustment of p value brings extra work amount to solving of problem, is unfavorable for the popularization of method.

Summary of the invention

The object of the invention is to the above-mentioned deficiency overcoming prior art, a kind of electricity tomography Lp regularization reconstruction method is proposed, the present invention is based on Lp regularization, in conjunction with Gauss-Newton iteration, solve L2 Regularization Solution and cross smooth and L1 Regularization Solution is excessively sparse problem, improve solving precision and the image reconstruction quality of electricity tomography inverse problem.Technical scheme of the present invention is as follows:

A kind of electricity tomography Lp regularization reconstruction method declined based on p vector equal difference, be applicable to bubble flow tomography, the method regards electricity chromatography imaging problem as a linear ill-posed problem Ax=b, wherein A is sensitivity matrix, b is retive boundary measured value vector, x is the imaging gray value vectors distributing corresponding with field domain material electrical properties, is called solution vector, adopts the method for Lp regularization inverse problem solving to carry out image reconstruction.It is characterized in that,

Utilize Gauss-Newton iteration to carry out often walking in iteration of Lp regularization inverse problem solving and upgrade according to gained solution the p vector be made up of the p value on each pixel in image, obtain the p distribution with field domain object space distribution character, finally complete to calculate to obtain and rebuild image, step is as follows:

(1) according to the measurement to tested field domain, the retive boundary measured value vector b and sensitivity matrix A needed for rebuilding is obtained;

(2) objective function of Lp regularization is set up;

(3) initiation parameter is set: the initial value x of solution vector x ₀=0, p vector initial value p ₀=2, p vector stop value p _stop=1; Setting iterations N;

(4) equal difference descending factors α=(p is calculated ₀-p _stop)/N;

(5) Gauss-Newton iterative formula is utilized to solve;

(6) utilize institute to solve and upgrade p vector: judge in solution vector, whether each element is zero, if then the p value of corresponding pixel points keeps back p value constant; If not, then the p value of corresponding pixel points is updated to the difference of back p value and equal difference descending factors;

(7) judge whether iteration completes, if then iteration ends, carry out next step operation, if not, then rebound (5) step continues to solve;

(8) according to finally solving gained gray-scale value, imaging is carried out.

Preferably, the objective function of described Lp regularization is: wherein λ > 0 is regularization coefficient, || || be Euclid norm, p vector meets arbitrary element and belongs to [1,2]; Objective function non-differentiability when considering p=1, utilizes approximate above-mentioned objective function, wherein n is the dimension of solution vector x, and j is the counting integer from 1 to n, x _jfor the element of jth in solution vector x, β is small adjustable parameter, meets β > 0.

Gauss-Newton iterative formula is utilized to be:

$x_k=x_{k-1}-{\left[{\▿}^{2}F\left(x_{k-1}\right)\right]}^{-1}\▿F\left(x_{k-1}\right)$

Wherein k is current iterations, meets 1≤k≤N; x _kthe solution that kth time iteration obtains, x _k-1it is the solution that (k-1) secondary iteration obtains; for working as x=x _k-1time objective function first differential, for x=x _k-1time objective function second-order differential, obtain respectively by two formula below:

$\▿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)$

Wherein p _k-1that the secondary iteration of back i.e. (k-1) upgrades the p vector obtained; Diag () is the function by vector structure diagonal matrix, and each element of vector forms the diagonal element of diagonal matrix, and the off-diagonal element of diagonal matrix is zero.

The invention has the beneficial effects as follows the process based on utilizing Gauss-Newton iterative Lp regularization, propose a kind ofly to utilize the equal difference of p in iterative process to decline to realize Lp adaptive electricity tomographic image method for reconstructing, for rebuild each pixel of smooth domain of image provide value be 2 or be close to 2 p value; For rebuild image have each pixel of object area provide value be 1 or be close to 1 p value, the p obtained with field domain object space distribution character distributes, and completes and finally solves.The present invention effectively overcomes the shortcoming of L2 regularization and L1 regularization, improves reverse temperature intensity precision and image reconstruction quality, has stronger robustness and larger applicability.

Accompanying drawing explanation

Fig. 1 is the FB(flow block) of a kind of electricity tomography Lp self-adaptation method for reconstructing based on the decline of p equal difference of the present invention;

Fig. 2 is the circular tested field domain in single cross section of Electrical Resistance Tomography system of the present invention and distribution of electrodes;

Fig. 3 is the true distribution of three models that example of the present invention is chosen: (a) is two roundlets model (b) be three round models (c) is four roundlet models

Fig. 4 is the imaging results schematic diagram of three models under L2 regularization solves in example of the present invention: the model (a-c) wherein in (a-c) difference corresponding diagram 3;

Fig. 5 is the imaging results schematic diagram of three models under L1 regularization solves in example of the present invention: the model (a-c) wherein in (a-c) difference corresponding diagram 3;

Fig. 6 is the imaging results schematic diagram of three models under the self-adaptation Lp regularization that this method proposes solves in example of the present invention: the model (a-c) wherein in (a-c) difference corresponding diagram 3;

In figure:

1, tested field domain 2, electrode

Embodiment

In conjunction with the accompanying drawings and embodiments a kind of electricity tomography Lp self-adaptation method for reconstructing declined based on p equal difference of the present invention is illustrated.

A kind of electricity tomography Lp self-adaptation method for reconstructing declined based on p equal difference of the present invention, based on Lp regularization, for the solution that L2 regularization is tried to achieve the smooth and problem that L1 regularization is tried to achieve solution is excessively sparse, in conjunction with Gauss-Newton iterative formula, the p that proposition utilizes back result of calculation to upgrade in an iterative process and is made up of the p value on each pixel in image is vectorial, and utilize the p vector after upgrading to complete current calculating, until the method for reconstructing of iteration ends.

As shown in Figure 1, be a kind of electricity tomography Lp self-adaptation method for reconstructing process flow diagram declined based on p equal difference of the present invention.Be illustrated in figure 2 the circular tested field domain in single cross section of Electrical Resistance Tomography system and the distribution of electrodes of one of electricity tomography, adopt 16 electrodes to be evenly distributed on field domain outer wall.Choosing three typical bubble flow models is embodiment, and in field domain, object truly distributes as shown in Fig. 3 (a-c).In order to embody the different of self-adaptation Lp regularization and L2 regularization and L1 regularization in the present invention better, provide the solving result of three models under these three kinds of regularization methods respectively.Embodiment comprises following concrete steps:

A kind of electricity tomography Lp self-adaptation method for reconstructing declined based on p equal difference, be applicable to bubble flow tomography, the method regards electricity chromatography imaging problem as a linear ill-posed problem Ax=b, wherein A is sensitivity matrix, b is retive boundary measured value vector, and x is the imaging gray value vectors distributing corresponding with field domain material electrical properties.

The objective function of Lp regularization is:

$\min F\left(x\right)={\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\mathrm{\λ}{\left|\right|x\left|\right|}_p^p$

Wherein λ > 0 is regularization coefficient, || || be Euclid norm, p is the constant vector be made up of the p value of rebuilding on each pixel of image.Objective function non-differentiability during consideration p=1, utilizes new objective function:

$\min F\left(x\right)={\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\mathrm{\λ}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(\sqrt{x_j^2+\mathrm{\β}}\right)}^p$

Approximate former objective function.Method for reconstructing includes following steps:

(1) for three typical bubble flow models, the boundary survey value vector sum sensitivity matrix needed for rebuilding separately is obtained respectively:

Boundary survey value vector measurand is placed in electricity chromatography imaging measurement system, 16 electrodes (as shown in Figure 2) are uniformly distributed outside tested field domain, adopt the pattern that current excitation voltage measurement and exciting electrode are not measured, gather the boundary voltage under cycle motivation circulation measurement, obtain the vector that 208 measured values are formed altogether; Inverse problem right-hand vector b is not containing the barnyard boundary voltage vector b1 of inclusions and the difference (i.e. right-hand vector retive boundary measured value vector b=b1-b2) having the boundary survey voltage vector b2 of thing field containing inclusions;

Sensitivity matrix is according to not containing the boundary survey voltage of the barnyard of inclusions, in conjunction with sensitivity theory, and meter sensitivity matrix, computing formula is:

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

Wherein A _ijbe a jth electrode pair to the sensitivity coefficient of i-th electrode pair, φ _i, φ _jbe respectively i-th electrode pair and a jth electrode pair is I at exciting current _i, I _jtime field domain Potential Distributing, x, y are distributed as the transverse and longitudinal coordinate of field domain;

(2) initiation parameter is set:

The initiation parameter arranged comprises: preset parameter, solution vector initial value x ₀=0, p vector initial value p ₀=2, p vector stop value p _stop=1; The iterations N=5 rule of thumb chosen, regularization coefficient λ=1 × 10 ^-4, fine setting parameter beta=1 × 10 ^-16;

(3) calculate equal difference descending factors α, computing formula is:

α＝(p ₀-p _stop)/N；

(4) utilize Gauss-Newton iterative formula to solve, its Iteration is:

$x_k=x_{k-1}-{\left[{\▿}^{2}F\left(x_{k-1}\right)\right]}^{-1}\▿F\left(x_{k-1}\right)$

Wherein k is the iterations of current step, meets 1≤k≤N; x _kthe solution that kth time iteration obtains, x _k-1it is the solution that (k-1) secondary iteration obtains; for working as x=x _k-1time objective function first differential, for x=x _k-1time objective function second-order differential, obtain respectively by two formula below:

$\▿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)$

Wherein p _k-1that the secondary iteration of back i.e. (k-1) upgrades the p vector obtained; Diag () is the function by vector structure diagonal matrix, and each element of vector forms the diagonal element of diagonal matrix, and the off-diagonal element of diagonal matrix is zero.And calculating second-order differential time eliminate higher order term;

(5) the p vector in iterative process is determined:

In order to compare the different of this method and L2 regularization and L1 regularization better, the p vector needed for L2 regularization, L1 regularization and self-adaptation Lp regularization is provided respectively to the model of three in embodiment:

To L2 regularization, in iterative process, p vector keeps all elements to meet p=2 constant;

To L1 regularization, in iterative process, p vector keeps all elements to meet p=1 constant;

To self-adaptation Lp regularization, utilize institute to solve in iterative process and upgrade p vector, update strategy is:

$p_k\left(l\right)=\left\{\begin{array}{cc}{p}_{k-1}\left(l\right)& x_{k-1}\left(l\right)=0\\ p_{k-1}\left(l\right)-\mathrm{\α}& x_{k-1}\left(l\right)>0\end{array}\right.$

Wherein p _kit is the p vector that kth time iteration upgrades; L is the position of each pixel of field domain;

(6) judge whether iteration completes, if then iteration ends, carry out next step operation, if not, then rebound step (4) continues to solve;

(7) according to finally solving gained gray-scale value, imaging is carried out:

Gained Regularization Solution is corresponded on the pixel of Electrical Resistance Tomography, carries out Grey imaging.Figure 4 shows that L2 regularization result of calculation image, Figure 5 shows that L1 regularization result of calculation image, Fig. 6 is self-adaptation Lp regularization result of calculation image, and in figure, (a-c) is corresponding with model (a-c) respectively.

Can find out, under identical optimum configurations, the solution of L2 regularization is excessively smooth, become image tail shadow large; And L1 regularization is excessively sparse, the size of object can not be described well; The result of self-adaptation Lp regularization, in the middle of L2 regularization and L1 regularization acquired results, improves reverse temperature intensity precision and image reconstruction quality, describes the distribution of object in field domain more accurately.

The above embodiment is several better model of the present invention, and the present invention is not limited to the content disclosed in this embodiment and accompanying drawing.Every do not depart from spirit disclosed in this invention under the equivalence that completes or amendment, all in the scope of protection of the invention.

Claims (3)

1. the electricity tomography Lp regularization reconstruction method declined based on p vector equal difference, be applicable to bubble flow tomography, the method regards electricity chromatography imaging problem as a linear ill-posed problem Ax=b, wherein A is sensitivity matrix, b is retive boundary measured value vector, x is the imaging gray value vectors distributing corresponding with field domain material electrical properties, is called solution vector, adopts the method for Lp regularization inverse problem solving to carry out image reconstruction.It is characterized in that,

Utilize Gauss-Newton iteration to carry out often walking in iteration of Lp regularization inverse problem solving and upgrade according to gained solution the p vector be made up of the p value on each pixel in image, obtain the p distribution with field domain object space distribution character, finally complete to calculate to obtain and rebuild image, step is as follows:

(1) according to the measurement to tested field domain, the retive boundary measured value vector b and sensitivity matrix A needed for rebuilding is obtained;

(2) objective function of Lp regularization is set up;

(3) initiation parameter is set: the initial value x of solution vector x ₀=0, p vector initial value p ₀=2, p vector stop value p _stop=1; Setting iterations N;

(4) equal difference descending factors α=(p is calculated ₀-p _stop)/N;

(5) Gauss-Newton iterative formula is utilized to solve;

(6) utilize institute to solve and upgrade p vector: judge in solution vector, whether each element is zero, if then the p value of corresponding pixel points keeps back p value constant; If not, then the p value of corresponding pixel points is updated to the difference of back p value and equal difference descending factors;

(7) judge whether iteration completes, if then iteration ends, carry out next step operation, if not, then rebound (5) step continues to solve;

(8) according to finally solving gained gray-scale value, imaging is carried out.

2. the electricity tomography Lp regularization reconstruction method declined based on p vector equal difference according to claim 1, it is characterized in that, the objective function of described Lp regularization is: wherein λ > 0 is regularization coefficient, || || be Euclid norm, p vector meets arbitrary element and belongs to [1,2]; Objective function non-differentiability when considering p=1, utilizes approximate above-mentioned objective function, wherein n is the dimension of solution vector x, and j is the counting integer from 1 to n, x _jfor the element of jth in solution vector x, β is small adjustable parameter, meets β > 0.

3. the electricity tomography Lp regularization reconstruction method declined based on p vector equal difference according to claim 2, is characterized in that: utilize Gauss-Newton iterative formula to be:

x _k＝x _k-1-[▽ ²F(x _k-1)] ^-1▽F(x _k-1)

Wherein k is current iterations, meets 1≤k≤N; x _kthe solution that kth time iteration obtains, x _k-1it is the solution that (k-1) secondary iteration obtains; ▽ F (x _k-1) for work as x=x _k-1time objective function first differential, ▽ ²f (x _k-1) be x=x _k-1time objective function second-order differential, obtain respectively by two formula below:

$\▿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)$

Wherein p _k-1that the secondary iteration of back i.e. (k-1) upgrades the p vector obtained; Diag () is the function by vector structure diagonal matrix, and each element of vector forms the diagonal element of diagonal matrix, and the off-diagonal element of diagonal matrix is zero.