# A Superresolution Image Reconstruction Algorithm Based on Landweber in Electrical Capacitance Tomography

Authors: Chen Deyun; Li Zhiqiang; Gao Ming; Wang Lili; Yu Xiaoyang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

According to the image reconstruction accuracy influenced by the “soft field” nature and ill-conditioned problems in electrical capacitance tomography, a superresolution image reconstruction algorithm based on Landweber is proposed in the paper, which is based on the working principle of the electrical capacitance tomography system. The method uses the algorithm which is derived by regularization of solutions derived and derives closed solution by fast Fourier transform of the convolution kernel. So, it ensures the certainty of the solution and improves the stability and quality of image reconstruction results. Simulation results show that the imaging precision and real-time imaging of the algorithm are better than Landweber algorithm, and this algorithm proposes a new method for the electrical capacitance tomography image reconstruction algorithm.


---

## Body

## 1. Introduction

Electrical capacitance tomography (ECT) is a high-tech cross-section detection technology in recent years, and it can achieve the visualization measurement of the multiphase media component parameter in the closed pipe and process container equipment. The ECT technology could be able to obtain the permittivity distribution inside the pipeline through measuring the capacitance values between the different plates of the capacitive sensing electrode array around the pipeline. It has the advantages of noninvasive, simple structure, low cost, fast response, good safety performance, wide application range, and so forth. So, it can be used for measuring the two-phase flow in the petroleum, chemical and electric power industries, and so forth.There are three principal factors with image reconstruction for the electrical capacitance tomography system. Firstly, it is the “soft field” effect; the sensitivity performance is in homogeneity affected by the medium distribution. Secondly, the measurements are less than the pixels of the image, and the underdetermined equation has numerous solutions. Thirdly, the solution stability is poor; even small perturbation in the measured capacitance data can cause great error. Due to the nonlinear and the “soft field” effect, the stability of the ECT system is poor and seriously pathological. So, it is difficult for image reconstruction [1]. At present, there are many methods about the ECT image reconstruction, such as the linear back projection (LBP), regularization method, Landweber iterative method [2], genetic algorithms [3], projection Landweber iterative method, conjugate gradient method [4], improved Gauss-Newton method [5], and acceleration algorithm based on polynomial [6]. The characteristic of the LBP is simple and fast reconstruction. But because its imaging quality is relatively poor, strictly speaking, the algorithm is only a qualitative algorithm. Regular methods commonly use empirical values, because the regular selection of parameters influences the imaging quality. It is specially worth mentioning that the online preiterative method (OIOR) based on Landweber iterative method has the same imaging quality with several Landweber iterations under the same speed with LBP, but its spatial resolution has yet to exceed Landweber method.


## 2. Basic Principle of ECT System

Electrical capacitance tomography system mainly consists of three parts: capacitive sensor, data acquisition and signal processing, and image reconstruction, as shown in Figure1. Pairs of the metal electrode plates are mounted uniformly on the outside of the insulating pipe, and exterior uses shielding cover. This constitutes the capacitive sensor. Each phase separation medium of the multiphase fluid in the pipeline has a different dielectric constant. When the phase concentration and distribution of the fluid change, the equivalent permittivity of the multiphase flow mixture will change. So, the measuring capacitance values between the plates will be changing too. The size of the capacitor values reflects the size and distribution of the multiphase flow medium concentration. The data acquisition and signal-processing section samples the capacitance values between each pair of the plates and obtains the data in the different angular positions after superposition. These data will be sent to the computer after processing. These measurement values reflect the distribution of the dielectric constant in the pipeline. The distribution information of the fluid medium in the pipeline is obtained through some image reconstruction algorithm by the computer.Figure 1
Composed of ECT system.n is the total number of independent electrode pairs, which is available for the ECT system with the N-electrode plate. Consider
(1)n=CN2=N·(N-1)2.
In this paper, the typical 8-electrode capacitance sensor is the research object, we can get 28 independent measured values as C1,C2,C3,…,C28. The capacitance value of degree of influence on the same plate pair has different effects on the medium of the different locations in the pipeline. That is the soft field characteristic. On the other hand, the same point of the medium within the pipeline has different effects on the capacitance values of different plate pairs. The capacitance values between a pair of capacitor plates can actually be approximately regarded as the varying degrees superposition by the capacitance of all the points within the pipeline. The capacitance measurement values can be seen as the projection data of the multiphase fluid on a direction within the pipeline. The electrical capacitance tomography obtains each phase parameter of the multiphase flow from the projection data; then reconstruct the image. According to the electrical principles, ignoring the impact of the distribution changes of each subphase of the pipe radial multiphase fluid and the shield, the measure value Cj which is for the capacitance between any two electrodes in the array type of the 8-electrode capacitance sensor is
(2)Cj=∬Dε(x,y)·Sj(x,y,ε(x,y))dxdy,
where j=1,2,…,28, D is the pipeline cross-sections, ε(x,y) represents the distribution function of the dielectric in the pipeline cross-section, and Sj(x,y,ε(x,y)) is the sensitivity distribution function of the measuring capacitance Cj, namely, the sensitivity of the capacitance Cj  to the media on point (x,y).The sensitivity distribution data is the priori data required for the image reconstruction. The sensitivity distribution which is also known as the sensitive field can reduce the soft field error of the sensitive field and improve the quality of the image reconstruction. It is extremely necessary to carry out in-depth analysis on the sensitivity distribution.The calculation formula of the sensitivity function is(3)Sj(k)≈Cj(k)-Cj(oil)(ε2-ε1)·δk,
where Cj(oil) is empty-pipe capacitance, that is, the capacitance when the pipeline is filled with the discrete phase (oil), and it is a constant, where j=1,2,…,28. When the dielectric constant of the kth microelement is ε2 in the pipeline and the dielectric constant of other microelements  is ε1, Cj(k)  is the capacitance value.


## 3. The Principle of the Landweber Iterative Algorithms

The ECT image reconstruction is a typical ill-posed problem, and its solution is unstable. In order to obtain meaningful results of reconstruction, we should use the method that can improve the quality and stability of the solution. Regularization technique is an effective method to deal with pathological problems, and Landweber iteration method is a typical regularization techniques. The Landweber iterative format can be expressed as [7]
(4)G(0)=0,G(m)=(I-aSTS)G(m-1)+aSTC.The iterative method actually uses the steepest decent method with the step size asa to solve the minimum of the quadratic functional ∥SG-C∥2. Actually, if
(5)minJ(G)=12∥SG-C∥2,G
(
m
)
=
(
I
-
a
S
T
S
)
G
(
m
-
1
)
+
a
S
T
C is  the steepest descent iterative sequence with the step size as a. On the other hand, constant a in the iterative scheme can also be seen as a relaxation factor, and it is a modified form from G(m)=G(m-1)-S(SG(m-1)-C).TheG(m) is expressed as G(m)=RmC; the operator Rm:Y→X is
(6)Rm:=a∑k=0m-1(I-aSTS)kST,m=1,2,….Use the singular system(μj,xj,yj) of compact operator S to expand. The expression of Rm is as follows:
(7)RmC=a∑j=1∞μj∑k=0m-1(1-aμj2)k(C,yj)μjxj=∑j=1∞1μj[1-(1-aμj2)m](C,yj)xj=∑j=1∞q(m,μj)μj(C,yj)xj,
where the filtering function may be expressed as follows:
(8)q(m,μ)=1-(1-aμ2)m.When constructing the regularization solution sequenceGm by Landweber iteration, the number of iterations m=m(δ) should be given, and it is as the iteration stopping criteria. One common method is that iteration stops when the first m which meets formula (9) is disappeared:
(9)∥SGm-C∥≤γδ.The so-called semiconvergence phenomenon will disappear when solves the inverse problem with Landweber iterative method. That is, in the early stages of the iteration approximate solution can be stably improved, and the effect of self-regularization is shown, but the iterations will tend to diverge when iterations exceed a certain threshold [8]. So, the key question is to find a suitable termination principle in order to match the number of iterations and the error level of original data. However, in practical applications, the error level of the original data is difficult to obtain. So, there is a common method to determine the number of iterations by empirical methods.Landweber iterative method has various modifications and most improvements concentrate on accelerating the problem. From the numerical optimization point, Landweber iterative method belongs to the most descent methods, and its convergence rate is relatively slow.v-method is an accelerated method of Landweber algorithm. The method by Brakhage is a polynomial acceleration method whose applications are broader than the method by Chebyshev. v-method can be expressed as the following iterative scheme [9]:
(10)Gk=Gk-1+μk(Gk-1-Gk-2)-ωkST(C-SCk-1),
where μ1=0,ω1=(4v+2)/(4v+1); when k≥1, iterative scheme is as follows:
(11)μk=(k-1)(2k-3)(2k+2v-1)(k+2v-1)(2k+4v-1)(2k+2v-3),ωk=4(2k+2v-1)(k+v-1)(k+2v-1)(2k+4v-1).In fact, standard Landweber iterative method can be seen as a special kind of steepest descent method. Standard landweber iterative method steepest descent method have the same iteration format, and the only difference is the selected aspect of the iteration step. The steepest descent method requires searching from each of iteration by the 1-D to obtain the optimal step size. Landweber iteration method makes the step size in the interval[0,∥S∥-2].


## 4. Superresolution Image Reconstruction Algorithm

### 4.1. Superresolution Reconstruction Model

In the ECT image reconstruction, small changes of capacitance value which is measured can cause the fluctuation of the solution. This can cause the deterioration of the imaging quality. Sequence images about the same object with different angles can make up such issues to a certain extent. The destination of image restoration is to recover image that has been degraded and make sure that the processed image is close to the original as possible Figure2 simulated the degradation processing of sequence images.Figure 2
Simulated degradation processing of sequence images.In order to establish mathematical model of super-resolution image reconstruction, suppose thatyi is ith low-resolution image, f is high-resolution image which needs to be reconstruct, Ei is the geometry transformation matrix, Ci is the fuzzy matrix, D is the downsampling matrix, and ni is the additional noise. We can obtain superresolution reconstruction model which is as follows [10]:
(12)yj=DCiEif+ni,1≤i≤p,
where f is the MN×1 vector, which is super-resolution image vectors permutation that we need to solve. p shows a low-resolution image with a total of p frames. yi is the MLNL×1 vector, which is said to get the low-resolution image whose size is ML×NL from the i-frame. M/ML=N/NL=L is the multiple of the corresponding super-resolution reconstruction. Ei is the motion compensation matrix of MN×MN which does pixel registration. And it is obtained by general multilevel block-matching registration. Ci is blurring function matrix MN×MN. D is downsampling matrix MLNL×MN. ni is the Gaussian noise obtained by the low-resolution image of the ith frame of the image. Array formula (12) by a column vector of the image and as follows:
(13)[y1y2⋮yp]=[DC1E1DC2E2⋮DCpEp]f+[n1n2⋮np],
abbreviated as
(14)Y=Hf+n.


### 4.2. Regularization of Solutions

Due to the arbitrariness of the displacement among the low-resolution images and the number of the reconstruction low-resolution images with the effect of the drop clarity matrix, the previous reconstruction problem is usually a pathological problem. That is, the above reconstruction requires introducing the regularization techniques for solving the ill-posed problem. The minimized formula solved by formula (14) is shown as follows:
(15)minfβ∥Y-Hf∥2+∥Lf∥2,
where L is the smooth operator for solving ill-posed problem and it usually uses the Laplacian in super-resolution reconstruction. β is the regularization parameter and it can be obtained by experiences in engineering approximation or other methods. Transformation is shown as follows:
(16)x=Lf,(17)A=HL-1.Formula (15) becomes
(18)minxβ∥Y-Ax∥2+∥x∥2.Minimizing the previous equation is equivalent to solving the following equations:(19)(βATA+I)x=βATY.Since the matrixA cannot be diagonalized, the solution of the previously mentioned equation is computationally intensive part in the super-resolution image reconstruction. Approximate the registration results among the images by iterative solution algorithm. And solve the previous equation by preconditioned conjugate gradient method. The computational complexity of the algorithm is high, and the computing complexity with the number of low-resolution images increases exponentially. Then, analyze the previous equation, and derive a direct calculation algorithm, so that the complexity of computation increases linearly with the number of low-resolution image.


### 4.3. The Export of Closed Solution

The solutions of (19) can be directly expressed as follows:
(20)x=β(βATA+I)-1ATY.According to the matrix inversion lemma,(21)(F+BCD)-1=F-1-F-1B(DF-1B+C-1)-1DF-1.
The equation obtained by formula (20) is shown as follows:
(22)x=β(I-AT(AAT+β-1I)-1A)ATY.The formula obtained by solution is shown as follows:(23)x=AT(AAT+β-1I)-1Y.According to formula, (14) and (17), AATis a p×p block matrix, and the (i,j) block is DCiEiL-1L-TEjTCjTDT whose size is MLNL×MLNL. The calculation uses the calculation model of cyclic matrix, so Ci,Ei, and L are all the block cyclic matrixes. According to the nature of the block cyclic matrix, CiEiL-1L-TEjTCjT still is a block cyclic matrix whose size is MN×MN. Since D is the next sampling matrix, the super resolution multiples on the horizontal and vertical directions are the same. DCiEiL-1L-TEjTCjTDT is still a block cyclic matrix whose size is MLNL×MLNL. Since block cyclic matrix plus the unit matrix is still a faster cycle matrix, matrix AAT+β-1I is a p×p partitioned matrix; each block of the matrix is still the block cyclic matrix whose size is MLNL×MLNL. According to the nature that the block cyclic matrix can be diagonalized, let Qi,j be the (i,j) block cyclic matrix of AAT+β-1I:
(24)Qi,j=WΛi,jW-1,
where Λi,j is the diagonal matrix which is diagonalized by the corresponding block circulant matrix. W is the corresponding orthogonal matrix, so the following formula can be obtained:
(25)AAT+β-1I=[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ1,1Λ1,2⋯Λ1,pΛ2,1Λ2,2⋯Λ2,p⋮⋮⋮Λp,1Λp,2⋯Λp,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1],(26)(AAT+β-1I)-1=[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ1,1Λ1,2⋯Λ1,pΛ2,1Λ2,2⋯Λ2,p⋮⋮⋮Λp,1Λp,2⋯Λp,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1].The corresponding calculation mainly converts to seek the inverse of the corresponding block diagonalization matrix:(27)Λ=[Λ1,1Λ1,2⋯Λ1,pΛ2,1Λ2,2⋯Λ2,p⋮⋮⋮Λp,1Λp,2⋯Λp,p].LetP be an element permutation matrix of block matrix [11]. PTP=I,[PTΛP]i,j;k,l=[Λ]k,l;i,j. The function is to substitute the (i,j) element of the first (k,l) block matrix for the (k,l) element of the (i,j) block matrix.So,(28)PTΛP=[A~10⋯00A~2⋯0⋮⋮⋮00⋯A~ML×NL],
where every A~i is a p×p matrix. The inverse of formula (27) is shown as follows:
(29)Λ-1=P[A~10⋯00A~2⋯0⋮⋮⋮00⋯A~ML×NL]-1PT.According to the nature and inverse of the block diagonal matrix in formula (29), its inverse is also the block diagonalization of the corresponding matrix inverse. That is,
(30)Λ-1=P[A~1-10⋯00A~2-1⋯0⋮⋮⋮00⋯A~ML×NL-1]PT.Thus, the inversion of the matrix whose size is(p×ML×NL)×(p×ML×NL) transforms into the inversion of ML×NLp×p matrix. According to the nature of the permutation matrix, formula (31) is shown as follows:
(31)Λ-1=[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p],
where every Λ~i,j is (ML×NL)×(ML×NL) diagonal matrix and substituted into formula (26):
(32)(ATA+β-1I)-1=[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1].According to formula (14), (16), (17), (23), and (32), formula (33) can be obtained. 								(33)f=L-1L-T[(C1E1)T(C2E2)T⋯(CpEp)T]×[DT0⋯00DT⋯0⋮⋮⋮00⋯DT]×[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1]Y.Diagonalize the corresponding block circulant matrix:(34)L=WΛLW-1,Ci=WΛCiW-1,Ei=WΛEiW-1.Substitute formulae (14) and (35) into formula (34):
(35)f=WΛL-1ΛL-*[(ΛC1ΛE1)*(ΛC2ΛE2)*⋯(ΛCpΛEp)*]×[W-1DT0⋯00W-1DT⋯0⋮⋮⋮00⋯W-1DT]×[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1]Y.The simplified formula is shown as follows:(36)f=WΛL-1ΛL-*∑i=1p[(ΛCiΛEp)*W-1DT∑j=1p(WΛ~i,pW-1yj)].So, the superresolution image reconstruction computable closed formula is derived. In the derivation process, there are many methods to be used, such as the diagonalization of the block circulation diagonal matrix, fast Fourier transform, downsampling matrix transpose multiplying, small matrix inversion, and vector product. The block circulant matrix can be diagonalized by fast Fourier transform of the convolution kernel. Down sampling matrix transpose multiplication is equivalent to add zero in the corresponding position. For the small size of the matrix inversion, it can be calculated directly.


## 4.1. Superresolution Reconstruction Model

In the ECT image reconstruction, small changes of capacitance value which is measured can cause the fluctuation of the solution. This can cause the deterioration of the imaging quality. Sequence images about the same object with different angles can make up such issues to a certain extent. The destination of image restoration is to recover image that has been degraded and make sure that the processed image is close to the original as possible Figure2 simulated the degradation processing of sequence images.Figure 2
Simulated degradation processing of sequence images.In order to establish mathematical model of super-resolution image reconstruction, suppose thatyi is ith low-resolution image, f is high-resolution image which needs to be reconstruct, Ei is the geometry transformation matrix, Ci is the fuzzy matrix, D is the downsampling matrix, and ni is the additional noise. We can obtain superresolution reconstruction model which is as follows [10]:
(12)yj=DCiEif+ni,1≤i≤p,
where f is the MN×1 vector, which is super-resolution image vectors permutation that we need to solve. p shows a low-resolution image with a total of p frames. yi is the MLNL×1 vector, which is said to get the low-resolution image whose size is ML×NL from the i-frame. M/ML=N/NL=L is the multiple of the corresponding super-resolution reconstruction. Ei is the motion compensation matrix of MN×MN which does pixel registration. And it is obtained by general multilevel block-matching registration. Ci is blurring function matrix MN×MN. D is downsampling matrix MLNL×MN. ni is the Gaussian noise obtained by the low-resolution image of the ith frame of the image. Array formula (12) by a column vector of the image and as follows:
(13)[y1y2⋮yp]=[DC1E1DC2E2⋮DCpEp]f+[n1n2⋮np],
abbreviated as
(14)Y=Hf+n.


## 4.2. Regularization of Solutions

Due to the arbitrariness of the displacement among the low-resolution images and the number of the reconstruction low-resolution images with the effect of the drop clarity matrix, the previous reconstruction problem is usually a pathological problem. That is, the above reconstruction requires introducing the regularization techniques for solving the ill-posed problem. The minimized formula solved by formula (14) is shown as follows:
(15)minfβ∥Y-Hf∥2+∥Lf∥2,
where L is the smooth operator for solving ill-posed problem and it usually uses the Laplacian in super-resolution reconstruction. β is the regularization parameter and it can be obtained by experiences in engineering approximation or other methods. Transformation is shown as follows:
(16)x=Lf,(17)A=HL-1.Formula (15) becomes
(18)minxβ∥Y-Ax∥2+∥x∥2.Minimizing the previous equation is equivalent to solving the following equations:(19)(βATA+I)x=βATY.Since the matrixA cannot be diagonalized, the solution of the previously mentioned equation is computationally intensive part in the super-resolution image reconstruction. Approximate the registration results among the images by iterative solution algorithm. And solve the previous equation by preconditioned conjugate gradient method. The computational complexity of the algorithm is high, and the computing complexity with the number of low-resolution images increases exponentially. Then, analyze the previous equation, and derive a direct calculation algorithm, so that the complexity of computation increases linearly with the number of low-resolution image.


## 4.3. The Export of Closed Solution

The solutions of (19) can be directly expressed as follows:
(20)x=β(βATA+I)-1ATY.According to the matrix inversion lemma,(21)(F+BCD)-1=F-1-F-1B(DF-1B+C-1)-1DF-1.
The equation obtained by formula (20) is shown as follows:
(22)x=β(I-AT(AAT+β-1I)-1A)ATY.The formula obtained by solution is shown as follows:(23)x=AT(AAT+β-1I)-1Y.According to formula, (14) and (17), AATis a p×p block matrix, and the (i,j) block is DCiEiL-1L-TEjTCjTDT whose size is MLNL×MLNL. The calculation uses the calculation model of cyclic matrix, so Ci,Ei, and L are all the block cyclic matrixes. According to the nature of the block cyclic matrix, CiEiL-1L-TEjTCjT still is a block cyclic matrix whose size is MN×MN. Since D is the next sampling matrix, the super resolution multiples on the horizontal and vertical directions are the same. DCiEiL-1L-TEjTCjTDT is still a block cyclic matrix whose size is MLNL×MLNL. Since block cyclic matrix plus the unit matrix is still a faster cycle matrix, matrix AAT+β-1I is a p×p partitioned matrix; each block of the matrix is still the block cyclic matrix whose size is MLNL×MLNL. According to the nature that the block cyclic matrix can be diagonalized, let Qi,j be the (i,j) block cyclic matrix of AAT+β-1I:
(24)Qi,j=WΛi,jW-1,
where Λi,j is the diagonal matrix which is diagonalized by the corresponding block circulant matrix. W is the corresponding orthogonal matrix, so the following formula can be obtained:
(25)AAT+β-1I=[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ1,1Λ1,2⋯Λ1,pΛ2,1Λ2,2⋯Λ2,p⋮⋮⋮Λp,1Λp,2⋯Λp,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1],(26)(AAT+β-1I)-1=[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ1,1Λ1,2⋯Λ1,pΛ2,1Λ2,2⋯Λ2,p⋮⋮⋮Λp,1Λp,2⋯Λp,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1].The corresponding calculation mainly converts to seek the inverse of the corresponding block diagonalization matrix:(27)Λ=[Λ1,1Λ1,2⋯Λ1,pΛ2,1Λ2,2⋯Λ2,p⋮⋮⋮Λp,1Λp,2⋯Λp,p].LetP be an element permutation matrix of block matrix [11]. PTP=I,[PTΛP]i,j;k,l=[Λ]k,l;i,j. The function is to substitute the (i,j) element of the first (k,l) block matrix for the (k,l) element of the (i,j) block matrix.So,(28)PTΛP=[A~10⋯00A~2⋯0⋮⋮⋮00⋯A~ML×NL],
where every A~i is a p×p matrix. The inverse of formula (27) is shown as follows:
(29)Λ-1=P[A~10⋯00A~2⋯0⋮⋮⋮00⋯A~ML×NL]-1PT.According to the nature and inverse of the block diagonal matrix in formula (29), its inverse is also the block diagonalization of the corresponding matrix inverse. That is,
(30)Λ-1=P[A~1-10⋯00A~2-1⋯0⋮⋮⋮00⋯A~ML×NL-1]PT.Thus, the inversion of the matrix whose size is(p×ML×NL)×(p×ML×NL) transforms into the inversion of ML×NLp×p matrix. According to the nature of the permutation matrix, formula (31) is shown as follows:
(31)Λ-1=[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p],
where every Λ~i,j is (ML×NL)×(ML×NL) diagonal matrix and substituted into formula (26):
(32)(ATA+β-1I)-1=[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1].According to formula (14), (16), (17), (23), and (32), formula (33) can be obtained. 								(33)f=L-1L-T[(C1E1)T(C2E2)T⋯(CpEp)T]×[DT0⋯00DT⋯0⋮⋮⋮00⋯DT]×[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1]Y.Diagonalize the corresponding block circulant matrix:(34)L=WΛLW-1,Ci=WΛCiW-1,Ei=WΛEiW-1.Substitute formulae (14) and (35) into formula (34):
(35)f=WΛL-1ΛL-*[(ΛC1ΛE1)*(ΛC2ΛE2)*⋯(ΛCpΛEp)*]×[W-1DT0⋯00W-1DT⋯0⋮⋮⋮00⋯W-1DT]×[W0⋯00W⋯0⋮⋮⋮00⋯W]×[Λ~1,1Λ~1,2⋯Λ~1,pΛ~2,1Λ~2,2⋯Λ~2,p⋮⋮⋮Λ~p,1Λ~p,2⋯Λ~p,p]×[W-10⋯00W-1⋯0⋮⋮⋮00⋯W-1]Y.The simplified formula is shown as follows:(36)f=WΛL-1ΛL-*∑i=1p[(ΛCiΛEp)*W-1DT∑j=1p(WΛ~i,pW-1yj)].So, the superresolution image reconstruction computable closed formula is derived. In the derivation process, there are many methods to be used, such as the diagonalization of the block circulation diagonal matrix, fast Fourier transform, downsampling matrix transpose multiplying, small matrix inversion, and vector product. The block circulant matrix can be diagonalized by fast Fourier transform of the convolution kernel. Down sampling matrix transpose multiplication is equivalent to add zero in the corresponding position. For the small size of the matrix inversion, it can be calculated directly.


## 5. Experimental Results and Analysis

In order to verify the effectiveness of the algorithm, do numerical simulation test by the aforementioned algorithm, and compare the reconstruction results with the Landweber iterative algorithm. Matlab simulation is used on a computer with Core (TM) 2 Duo CPU, 2 GHz, and 1 G memory. The ECT system uses the 8-electrode; display uses32×32 subdivision, and the imaging region is divided into 1024 pixels; there are 856 imaging units on the effective area of the pipeline cross-section. The research object, respectively, uses the typical stratified flow, core flow, and trickle flow. The image reconstruction results are as shown in Table 1 (Black area is water, and white area is oil).Table 1
Three streams experimental results contrast.
Stratified flow (a)
Core flow (b)
Trickle flow (c)
Landweber algorithm
The algorithmIn order to evaluate the imaging performance of the super resolution image reconstruction algorithm based on Landweber in Electrical Capacitance Tomography, Landweber algorithm is used for comparative evaluation. This paper uses three objective evaluation methods: mean square error, signal to noise ratio, and peak signal to noise ratio. The details are as follows.(1)
Mean Square Error.Mean square error indicates the proximity of two images of the same size. Usually, in the super resolution reconstruction experiment, the reconstructed images are compared with the known images after obtaining the low-resolution image sequences using the original images:(37)MSE=1MN∑i=1M∑j=1N[f(i,j)-g(i,j)]2,
where M,N are the length and width of the image.  f(i,j) represents the pixels of the original image. g(i,j) represents the pixels of the reconstructed image.The exact figures of the evaluations simulated by Matlab are shown in Table2.(2)
Signal to Noise Ratio.Table 2
Mean square error.
Flow pattern
(a)
(b)
(c)
Landweber algorithm
27.86
42.12
81.37
The algorithm
16.81
20.16
41.72Consider(38)SNR=10log10{(1/MN)∑i=1M∑j=1Nf(i,j)2MSE},
where M,N are the length and width of the image.  f(i,j) represents the pixels of the original image.Signal to noise ratio is used to measure the ratio between the reconstructed image and the original image. The larger the ratio is, the more information the image contains. On the contrary, the smaller the ratio is, the less information the image contains.The exact figures of the evaluations simulated by Matlab are shown in Table3.(3)
Peak Signal to Noise Ratio.Table 3
Signal to noise ratio.
Flow pattern
(a)
(b)
(c)
Landweber algorithm
14.23
12.28
14.79
The algorithm
23.02
20.76
21.60Peak signal to noise ratio indicates the ratio between the maximum possible power signal and the destructive noise power which affects the representation accuracy. The larger the ratio is, the less the distortion is:(39)PSNR=10log10(2552MSE).The exact figures of the evaluations simulated by Matlab are shown in Table4.Table 4
Peak signal to noise ratio.
Flow pattern
(a)
(b)
(c)
Landweber algorithm
24.83
19.21
24.61
The algorithm
33.54
28.79
32.39In Table1, the quality of the reconstructed images for the flow pattern (a), (b), (c) has improved by using the algorithm proposed in this paper. It can be seen that the resolution in the periphery and center of the imaging area is better than Landweber reconstruction algorithm. The data listed in Tables 2, 3, and 4 are consistent with the results of the corresponding images listed in Table 1, from the MSE, SNR, and PSNR.


## 6. Conclusions

This paper presents the super-resolution image reconstruction algorithm based on Landweber iteration and gives the mathematical theory calculation steps of the algorithm based on the basic principle of the electrical capacitance tomography. Through the experimental simulation analysis, the algorithm can overcome the “soft field” nature of the electrical capacitance tomography and has a higher imaging accuracy. Experimental results show that the super resolution image reconstruction algorithm based on Landweber iterative can significantly improve imaging quality, stability, and so on. So, it proposes a new method for ECT image reconstruction algorithm.



---

*Source: 128172-2013-10-21.xml*