Image Reconstruction Based on Homotopy Perturbation Inversion Method for Electrical Impedance Tomography

Abstract

The image reconstruction for electrical impedance tomography (EIT) mathematically is a typed nonlinear ill-posed inverse problem. In this paper, a novel iteration regularization scheme based on the homotopy perturbation technique, namely, homotopy perturbation inversion method, is applied to investigate the EIT image reconstruction problem. To verify the feasibility and effectiveness, simulations of image reconstruction have been performed in terms of considering different locations, sizes, and numbers of the inclusions, as well as robustness to data noise. Numerical results indicate that this method can overcome the numerical instability and is robust to data noise in the EIT image reconstruction. Moreover, compared with the classical Landweber iteration method, our approach improves the convergence rate. The results are promising.

1. Introduction

Electrical impedance tomography (EIT) is an imaging modality which seeks to reconstruct the electrical conductivity distribution inside the interested objection from the surface electrical measurements. The process of generating an image from the measurements is called image reconstruction. Being a noninvasive, portable, and inexpensive methodology, the EIT technique has been extensively applied in medical imaging, geophysical exploration, nondestructive testing of materials, and so forth [1–3].

However, the image reconstruction for EIT mathematically is by nature a nonlinear and ill-posed inverse problem, which makes the image reconstruction a challenging task. Nonlinearity increases the requirement for using iterative image reconstruction techniques rather than noniterative ones, while ill-posedness makes any numerical reconstruction method require some form of regularization technique to obtain stable solution. Image quality greatly depends on the performance of the forward solver and the measurement errors but also depends on the regularization technique and the converging nature of iterative algorithm used in inverse problem. At present, various reconstruction algorithms have been developed to improve the image quality and increase the speed of inversion algorithms [4–10]. In the past few years, iterative image reconstruction with an updated forward solution has been investigated. In particular, the Landweber iteration method is one of the most straightforward and popular iterative algorithms for image reconstruction problem [11–14]. The advantages of this algorithm include the simple implementation and low computational cost, due to the fact that only the gradient information of the data fitting is used in the process of implementing the image reconstruction. Nevertheless, the ill-conditioning of the sensitivity matrix makes the convergence of Landweber iteration algorithm relatively slow and need a lot of iterations to get a good-quality image. Landweber iteration has a semiconvergence characteristic; that is, reconstructed image error often starts to increase gradually after getting to local minimum. Some preconditioning methods have been proposed to speed up the convergence rate [15–17].

Homotopy perturbation method is a novel and effective iteration method, which has been successfully applied to solve various of nonlinear equations [18–20] in various branches of science and engineering. Recently, it has also been extended for dealing with nonlinear ill-posed inverse problems [21, 22]. The essential idea of this method is to introduce a homotopy parameter and combine the traditional perturbation method with the homotopy technique. One of the most notable features of homotopy perturbation method is that usually only a few perturbation terms can be sufficient to obtain a reasonably accurate solution. In this paper, a novel iteration regularization scheme based on the homotopy perturbation technique, namely, homotopy perturbation inversion method, is applied to investigate the EIT inverse problem for the first time. Numerical simulations are carried out to illustrate the feasibility and effectiveness of this method. Simulation results indicate that this method can overcome the numerical instability and is robust to data noise in the EIT image reconstruction. Moreover, in contrast to the classical Landweber iteration method, our approach improves the convergence rate. As a result an efficient method for EIT image reconstruction is introduced.

The outline of this paper is as follows. In Section 2, we briefly describe the mathematical model for the EIT forward and inverse problem. In Section 3, the homotopy perturbation inversion scheme is presented, which is an iteration regularization scheme. A stable reconstruction is guaranteed provided the step length parameter is chosen appropriately and a stopping rule is formulated. To illustrate the well performance of this method, numerical simulations are provided in Section 4. Finally, a short conclusion is drawn in Section 5.

2. Mathematical Formulation of the EIT Problem

In this section we give a brief account on the mathematical model for EIT, which is composed of forward problem and inverse problem.

We first consider a typical acquisition experiment for EIT. Assume that 𝐿 electrodes with perfect conductor have been fixed around the surface 𝜕Ω of the object. Current is injected through some subset of these electrodes to excite the object under investigation, and the induced voltages are measured at all 𝐿 electrodes. This procedure is repeated many times with different electrodes until a sufficient amount of data has been obtained. Mathematical models for the measurement are described in [23, 24]. The complete electrode model, nowadays the standard model for medical applications, is presently the most accurate model for EIT since it takes into account the presence of the discrete electrodes and the effect of the contact impedance and can match the measurement precision of the experiment well [24]. The model is composed of the following equations: ()𝜎∇⋅𝜎∇𝑢=0,inΩ,𝜕𝑢𝜕𝑛=0,on𝜕Ω⋃^𝐿_𝑙=1𝑒_𝑙𝑒𝑙,𝜎𝜕𝑢𝜕𝑛𝑑𝑆=𝐼_𝑙,for𝑙=1,2,…,𝐿,𝑢+𝑧_𝑙𝜎𝜕𝑢𝜕𝑛=𝑈_𝑙,on𝑒_𝑙,𝑙=1,2,…,𝐿, (1) where 𝜎 is the electrical conductivity, 𝑢 is the electrical potential, 𝑒_𝑙 denotes the 𝑙th electrode, 𝑧_𝑙 is the effective contact impedance between the 𝑙th electrode and tissue, and 𝑛 is the unit outward normal to 𝜕Ω and 𝜕𝑢/𝜕𝑛=∇𝑢⋅𝑛. Furthermore, 𝑈_𝑙 and 𝐼_𝑙 are the electrical potential and electrical current on the 𝑙th electrode, respectively. In addition, in order to ensure the solvability of the elliptic boundary value problem (1) and the uniqueness of the solution, the following conditions are required to be satisfied: ^𝐿_𝑙=1𝐼_𝑙=0,^𝐿_𝑙=1𝑈_𝑙=0. (2) We denote vectors consisting of the electrical potentials and currents on the electrodes by 𝑈 and 𝐼; that is, 𝑈=(𝑈₁,𝑈₂,…,𝑈_𝐿)∈ℝ^𝐿 and 𝐼=(𝐼₁,𝐼₂,…,𝐼_𝐿)∈ℝ^𝐿, respectively.

2.1. The Forward Problem and Finite Element Approximation

For known conductivity distribution 𝜎, injecting current 𝐼, and contact impedance 𝑧_𝑙, the forward problem is to find the internal electrical potentials 𝑢 in Ω and the voltages 𝑈 on the electrodes that satisfy the complete electrode model (1). In practice the system (1) cannot be solved analytically for any arbitrary distribution, and therefore numerical approximations must be used to find an approximate solution. Finite element method [25, 26] is particularly well suited to the discretization of this problem.

As we all know, the procedure in the finite element method starts with the so-called variational (Galerkin) formulation of the problem (1), also called the weak form. We first need define some notation. By 𝐻¹(Ω) we denote the 𝐿²(Ω)-based Sobolev space with smoothness index 1, and 𝐻¹(Ω)={𝑢∈𝐻¹(Ω)∶∫_𝜕Ω𝑢𝑑𝑆=0}. For convenience, let us denote the space of length 𝐿 by ℝ^𝐿_⋄={𝑉∈ℝ^𝐿∶∑^𝐿_𝑙=1𝑉_𝑙=0}. The weak formulation of the forward problem can be stated as follows: given a current vector 𝐼∈ℝ^𝐿_⋄, a conductivity distribution 𝜎, and positive contact impedance 𝑧_𝑙, to find the solution (𝑢,𝑈)∈𝐻¹(Ω)⊕ℝ^𝐿_⋄ such that 𝑏((𝑢,𝑈),(𝑣,𝑉))=^𝐿_𝑙=1𝐼_𝑙𝑉_𝑙,∀(𝑣,𝑉)∈𝐻¹(Ω)⊕ℝ^𝐿_⋄, (3) where the strictly elliptic bilinear form 𝑏 is defined by 𝑏),(=((𝑢,𝑈𝑣,𝑉))_Ω𝜎∇𝑢⋅∇𝑣𝑑𝑥+^𝐿_𝑙=11𝑧_𝑙𝑒𝑙𝑢−𝑈_𝑙𝑣−𝑉_𝑙𝑑𝑆. (4) Indeed, the existence and uniqueness of a solution (𝑢,𝑈)∈𝐻¹(Ω)⊕ℝ^𝐿_⋄ have been shown using the Lax-Milgram lemma in [24].

Now we consider the piecewise linear function {𝜙_𝑖} as the finite element basis function, and the solution domain Ω will be divided into small triangles elements. We approximate the potential distribution 𝑢 within Ω as 𝑢^ℎ=^𝑁_𝑖=1𝛼_𝑖𝜙_𝑖 (5) and potentials 𝑈 on the electrodes as 𝑈^ℎ=^𝐿−1_𝑙=1𝛽_𝑗𝑛_𝑗, (6) where 𝑁 is the number of nodes in finite element mesh, 𝑛₁=[1,−1,0,…,0]^𝑇, 𝑛₂=[1,0,−1,0,…,0]^𝑇∈ℝ^𝐿, and so forth and 𝛼_𝑖 and 𝛽_𝑗 are the coefficients to be determined. For notational reasons we identify 𝑢^ℎ and 𝑈^ℎ with their coordinates and write 𝑢^ℎ=(𝛼₁,…,𝛼_𝑁), 𝑈^ℎ=(𝛽₁,…,𝛽_𝐿−1). By substituting these two approximations into variational equation (3) and by choosing 𝑣=𝜙_𝑖 and 𝑉=𝑛_𝑗, we obtain a discrete linear system 𝐴(𝜎)𝑏=𝑓, (7) where the data vector 𝑓=(𝟎,𝐼)̂^𝑇 with 𝟎∈ℝ^𝑁 and 𝐼̂=(𝐼₁−𝐼₂,𝐼₁−𝐼₃,…,𝐼₁−𝐼_𝐿), the determined coefficient vector 𝑏=(𝑢^ℎ,𝑈^ℎ)^𝑇, and the stiffness matrix 𝐴(𝜎)=𝐶𝐵𝐶_𝑇𝐷, (8) with 𝐵_𝑖𝑗=_Ω𝜎∇𝜙_𝑖⋅∇𝜙_𝑗+𝑑𝑥_𝐿𝑙=11𝑧_𝑙𝑒𝑙𝜙_𝑖𝜙_𝑗𝐶𝑑𝑆,𝑖,𝑗=1,2,…,𝑁,_𝑖𝑙1=−𝑧_𝑙𝑒𝑙𝜙_𝑖𝐷𝑑𝑆,𝑖=1,2,…,𝑁,𝑙=1,2,…,𝐿−1,_𝑙𝑙=1𝑧_𝑙𝑒𝑙𝑑𝑆,𝑙=1,2,…,𝐿−1. (9) After solving the determined coefficients 𝑏=𝐴(𝜎)⁻¹𝑓, we can use (6) to find the approximate potentials 𝑈^ℎ on the electrodes.

2.2. The Inverse Problem

In practice, the conductivity distribution 𝜎 is not known. What we know is merely all pairs of injected current data 𝐼 and resulting voltage data 𝑈 on the electrodes. The inverse problem under the complete electrode model (1) is to estimate the conductivity distribution 𝜎 from these data.

By Ohm’s law, the output is indeed a measurement (discretization) of the linear current-to-voltage operator that is, 𝑅_𝜎∶𝐼∈ℝ^𝐿_⋄→𝑈∈ℝ^𝐿_⋄, which maps the applied electrode current data 𝐼 to the computed electrode voltage data 𝑈 or denoted by 𝑈(𝜎) to emphasize the dependence on the conductivity parameter as follows: 𝑈(𝜎)=𝑅_𝜎𝐼 (10) and 𝑅_𝜎𝐼=Υ_𝐷𝑢∈ℝ^𝐿_⋄, where Υ_𝐷∶𝐻¹(Ω)→ℝ^𝐿_⋄ denotes the Dirichlet trace operator, which projects the solution 𝑢 of the forward problem to the electrode voltages 𝑈.

Now we denote the solution operator by 𝐹(𝜎;𝐼), such that 𝑈(𝜎)=𝐹(𝜎;𝐼), (11) where 𝐹 relates the conductivity distribution 𝜎 with the corresponding virtual measurements 𝑈, often referred to as the parameter-to-observation map, also called the forward operator. Notably, 𝐹 depends linearly on 𝐼 but nonlinearly on 𝜎. For a fixed current vector 𝐼, we can view 𝐹(𝜎;𝐼) as a function of 𝜎 only. It is important to distinguish between the PDE operator 𝐴(𝜎) and the forward operator 𝐹(𝜎). Though both are associated with the model problem (1), their meaning is not the same. Indeed, 𝐴(𝜎) describes the forward problem, while 𝐹(𝜎) actually represents its solution.

However, in practical applications, actual observations for EIT are typically potential differences between certain electrodes, denoted by 𝑈^obs, which may not be known exactly and unavoidably noisy due to the measurement or model error. Thus the observation model response of EIT can be formulated as the following nonlinear operator equation: 𝑈^obs=𝐹(𝜎)+𝑒, (12) where 𝑒 is the noise in the measurements.

Thus the inverse problem for EIT becomes that of finding the approximation conductivity distribution 𝜎 from the observed partial noisy knowledge 𝑈^obs, also called image reconstruction.

3. Homotopy Perturbation Inversion Method

One popular approach to deal with the inverse problem (12) is the nonlinear least-square method min_𝜎𝐹(𝜎)−𝑈^obs2. (13)

The corresponding Euler equation of (13) is given by 𝐹^′(𝜎)^𝑇𝐹(𝜎)−𝑈^obs=0. (14)

By the homotopy perturbation technique, a fixed point homotopy function 𝐻(𝑣,𝑝)∶𝑅×[0,1]→𝑅 is constructed for (14), which satisfies 𝐻𝑣,𝑝=𝑝𝐹^′(𝑣)^𝑇𝐹(𝑣)−𝑈^obs+1−𝑝𝑣−𝜎₀=0, (15) where 𝑝∈[0,1] is an embedding homotopy parameter and 𝜎₀ is an initial guess which incorporate a prior estimate of the true conductivity distribution 𝜎^∗. Obviously, from these definitions we will have 𝐻(𝑣,0)=𝑣−𝜎₀𝐻()=0,𝑣,1=𝐹^′(𝑣)^𝑇𝐹(𝑣)−𝑈^obs=0. (16) As the parameter 𝑝 changed continuously from zero to one, the trivial problem 𝐻(𝑣,0)=0 is continuously deformed to the original problem 𝐻(𝑣,1)=0 that is, 𝑣 changed from 𝜎₀ to the approximation solution 𝜎. In topology, this is called deformation, and 𝜎−𝜎₀ and 𝐹^′(𝜎)^𝑇(𝐹(𝜎)−𝑈^obs) are called homotopic.

Next, according to the homotopy perturbation technique, we can use the parameter 𝑝 as a small parameter and assume that the solution of (15) can be expressed as a power series in 𝑝, 𝑣=𝜎₀+𝑝𝜎₁+𝑝²𝜎₂+𝑝³𝜎₃+⋯. (17) As 𝑝→1, the approximate solution of (14) is obtained with 𝜎=lim_𝑝→1𝑣=𝜎₀+𝜎₁+𝜎₂+𝜎₃+⋯. (18)

We rewrite (15) with the Taylor expansion of 𝐹(𝜎) at 𝜎₀ and ignore the higher order term, which yields 𝐻𝑣,𝑝=𝑝𝐹^′𝜎₀^𝑇𝐹𝜎₀+𝐹^′𝜎₀𝑣−𝜎₀−𝑈^obs+1−𝑝𝑣−𝜎₀=0. (19) Substituting (17) into (19), we will have 𝑝𝐹^′𝜎₀^𝑇𝐹𝜎₀+𝐹^′𝜎₀×𝑝𝜎₁+𝑝²𝜎₂+𝑝³𝜎₃+⋯−𝑈^obs+1−𝑝𝑝𝜎₁+𝑝²𝜎₂+𝑝³𝜎₃+⋯=0. (20) Equating the coefficients for the different power of 𝑝, we get the following formulations: 𝑝¹∶𝐹^′𝜎₀^𝑇𝐹𝜎₀−𝑈^obs+𝜎₁𝑝=0,_𝑘+1∶𝐹^′𝜎₀^𝑇𝐹^′𝜎₀𝜎_𝑘+𝜎_𝑘+1−𝜎_𝑘=0,𝑘=1,2,3,…. (21) Then starting with an initial prior estimate 𝜎₀ and in turn solving (21), we can obtain 𝜎₁=−𝐹^′𝜎₀^𝑇𝐹𝜎₀−𝑈^obs,𝜎_𝑘+1=𝐼−𝐹^′𝜎₀^𝑇𝐹^′𝜎₀^𝑘−𝐹^′𝜎₀^𝑇𝐹𝜎₀−𝑈^obs,𝑘=1,2,3,…. (22)

Thus by (18) we can get the approximate solution of (14); that is, 𝜎=lim_𝑝→1𝑣=𝜎₀+𝜎₁+𝜎₂+𝜎₃+⋯=𝜎₀+−𝐹^′𝜎₀^𝑇𝐹𝜎₀−𝑈^obs+𝐼−𝐹^′𝜎₀^𝑇𝐹^′𝜎₀−𝐹^′𝜎₀^𝑇𝐹𝜎₀−𝑈^obs+𝐼−𝐹^′𝜎₀^𝑇𝐹^′𝜎₀²−𝐹^′𝜎₀^𝑇𝐹𝜎₀−𝑈^obs+⋯. (23)

The above formulation will allow us to construct two types of iteration methods for solving the nonlinear inverse problem (12). The first iteration scheme is to use the first-order approximation, as follows, for a given prior 𝜎₀, 𝜎_𝑘+1=𝜎_𝑘−𝛼_𝑘𝐹^′𝜎_𝑘^𝑇𝐹𝜎_𝑘−𝑈^obs, (24) which is the well-known classical Landweber iteration method (CLIM) [27], viewed as one of the variations of the steepest gradient descent method, and the second iteration scheme is to use the second-order approximation, as follows, for a given prior 𝜎₀, 𝜎_𝑘+1=𝜎_𝑘−𝛼_𝑘2𝐼−𝐹^′𝜎_𝑘^𝑇𝐹^′𝜎_𝑘×𝐹^′𝜎_𝑘^𝑇𝐹𝜎_𝑘−𝑈^obs, (25) where 𝛼_𝑘 is an appropriately chosen step length. Since the iteration-varying step length usually shows better performance than the constant value, the step length 𝛼_𝑘 for iterations (24) and (25) was chosen with the same criterion; that is, 𝛼_𝑘=1/𝜆_max(𝐽_𝑘^𝑇𝐽_𝑘), where 𝐽_𝑘 is the Jacobian matrix of 𝐹(𝜎) calculated at the 𝑘th iteration 𝜎_𝑘 and 𝜆_max(𝐽_𝑘^𝑇𝐽_𝑘) denotes the maximum eigenvalue of 𝐽_𝑘^𝑇𝐽_𝑘. We here call iteration scheme (25) the homotopy perturbation inversion method (HPIM). It has faster convergence than the classical Landweber iteration (24), which has been analyzed and proved in [28].

Considering the measurement errors of the observation data 𝑈^obs we employ the discrepancy principle as a stopping rule for HPIM, which determines the stopping index 𝑘_∗=𝑘_∗(𝛿,𝑈^obs) by 𝑈^obs−𝐹𝜎_𝑘∗≤𝜏𝛿<𝑈^obs−𝐹𝜎_𝑘,0≤𝑘<𝑘_∗, (26) for some sufficiently large 𝜏>0. In fact, this stopping rule renders this method as a regularization method.

4. Numerical Simulations

In this section, some numerical results for simulated data were presented to demonstrate the feasibility and effectiveness of the homotopy perturbation inversion method (HPIM). The quality of the images reconstructed by HPIM is compared with that of the classical Landweber iteration method (CLIM). All the simulations are 2D, as this allows an easier visualization of the results and faster simulation. A 16-electrode system was selected for simulations. We used constant injection current between adjacent electrodes and adjacent voltage measurement between all other electrodes. To avoid committing what is referred to as an inverse crime, two different finite element meshes were used for the forward inverse solvers, as shown in Figure 1, and the small green strips on the perimeter of the meshes show the position of the ideal electrodes. A mesh with 1968 elements and 1049 nodes was used for the forward simulations. For the inverse computations another mesh with 492 elements and 279 nodes was used to reduce the computational burden. Considering the simplicity’s sake of mesh refinement, we divided one element used in the inverse calculation into four elements to generate the forward mesh. In the reconstructions, contact impedances under the electrodes were assumed to be known and are set to be 0.05 for all electrodes.

(a)

(b)

Figure 1 

The forward and inverse models. (a) the forward model with 1968 triangulation elements; (b) the inverse model with 492 elements.

Figure 2 shows three original phantom models used for our simulations. The phantoms from left to right have various inclusions at different locations, occupying 36 cells, 16 cells, and 28 cells on the total size (492 cells) of the unknown discrete conductivity vector, respectively. The higher (background) and lower (inclusions) conductivity distribution values are set to be 1.0 and 6.0, respectively.

(a)

(b)

(c)

Figure 2 

Three simulation phantoms.

In the following, images for these three phantom models were reconstructed by the homotopy perturbation inversion method (HPIM) and classical Landweber iteration method (CLIM), respectively. For all simulations, the synthetic ideal data were generated using the finite element forward model based on the EIDORS 2D EIT software developed by Vauhkonen et al. [29]. There are a total of 256 data points being generated for the entire survey. The initial conductivity distribution values for both methods are set to the background values.

In addition to the reconstructed images, in order to compare the reconstruction performance more quantitatively, the following criteria are used: DE=𝐹𝜎^†−𝑈^obs𝜎,RE=^†−𝜎^∗‖𝜎^∗‖, (27) where 𝜎^† represents the reconstructed conductivity distribution, 𝜎^∗ is the true one, and DE and RE represent the discrepancy error and relative image error, respectively.

4.1. Noise-Free Case

The numerical performances of HPIM and CLIM were first tested by the noise-free data. Figures 3(a), 3(b), 4(a), 4(b), 5(a), and 5(b) show the reconstructed results obtained from HPIM and CLIM for three different test phantoms, respectively. Columns from left to right in each figure represent the reconstructions at iterations 10, 20, and 50, respectively. It can be seen that HPIM is superior to CLIM and can reconstruct the conductivity distribution even after 10-step iteration. After 20 iterations, the reconstruction performance with HPIM shows great improvement over CLIM. To compare the performance of the above two methods more quantitatively, the relative errors of reconstructed images in Figures 3(a), 3(b), 4(a), 4(b), 5(a), and 5(b), using (27), are calculated and listed in Table 1. It can be observed that HPIM has less relative error than CLIM at the same iterations. To visibly illustrate the convergence behavior of both methods, we also draw the curves from discrepancy errors and relative errors versus iteration time for the reconstructions of these three phantoms, as shown in Figures 3(c)–5(c). It is clear that the discrepancy errors and relative errors of HPIM are consistently lower than those of CLIM. As can be expected, the convergence rate of the proposed HPIM is accelerated significantly compared to that of CLIM.

Table 1 

Comparison of relative errors for the reconstructions displayed in Figures 3(a), 3(b), 4(a), 4(b), 5(a), and 5(b).

(a)

(b)

(c)

Figure 3 

Reconstruction results with HPIM (see (a)) and CLIM (see (b)) for the first phantom at different iteration number, as well as convergence behavior (see (c)). (c) Left: discrepancy error; right: relative error. In both plots, the horizontal axes denote the time of the iteration process.

(a)

(b)

(c)

Figure 4 

Reconstruction results with HPIM (see (a)) and CLIM (see (b)) for the second phantom at different iteration number, as well as convergence behavior (see (c)). (c) Left: discrepancy error; right: relative error. In both plots, the horizontal axes denote the time of the iteration process.

(a)

(b)

(c)

Figure 5 

Reconstruction results with HPIM (see (a)) and CLIM (see (b)) for the third phantom at different iteration number, as well as convergence behavior (see (c)). (c) Left: discrepancy error; right: relative error. In both plots, the horizontal axes denote the time of the iteration process.

4.2. Noise-Contaminated Case

Studies indicate that EIT image reconstruction process is a typical ill-posed problem. A small perturbation in the input data may yield a large fluctuation of the solution, which may make the inversion results meaningless. In the data acquisition processes, the measurement data unavoidably contain noise. The robustness of reconstruction technique with respect to noise is crucial. In order to simulate real measurement environment, the noise-contaminated input data are used to evaluate the numerical performance and robustness to noise of HPIM. Usually, the synthetic Gaussian noise was generated to simulate systematic and random errors existing in the data acquisition processes. The resulting noisy data was indicated by adding noise to the synthetic ideal data, defined as 𝑈^noise=𝑈^syn+𝜀⋅𝑛, (28) where 𝜀 is the noise level and 𝑛 is a noise vector obeying a Gaussian distribution with zero mean and standard deviation of 1.

In the following, we illustrate the performance of HPIM for the noise-contaminated case by comparing against the common CLIM. According to the stopping rule (26) with 𝜏=20, for each phantom, the corresponding reconstructed images by HPIM under four different noise levels of 5%, 3%, 1%, and 0.5% are listed in Figures 6(a)–6(d), respectively. As can be expected, HPIM shows the favorable robustness. Since the images reconstructed by CLIM have the similar qualities of those by HPIM, we do not list them in Figure 6. We summarized more detailed computational results in Table 2, which is against the second phantom with two inclusions, and the corresponding reconstructed images are partially displayed in the middle column of Figure 6. Here 𝑘_∗ and 𝑇 denote the stooping steps and the CPU time (in seconds), respectively. We find from Table 2 that images reconstructed by HPIM have the similar error qualities to those by CLIM, but in less steps and computational time. Additionally, it also can be seen from Figure 6 and Table 2 that with the increase of the noise levels the reconstruction quality decreases, which indicates that the accuracy of the measurement data has an effect on the imaging quality.

Table 2 

Comparison of reconstructions for the second phantom with 𝜏=20 at four different noise levels.

(a)

(b)

(c)

(d)

Figure 6 

Reconstruction results by HPIM for three phantoms with different noise levels. (a) 𝜀=5%; (b) 𝜀=3%; (c) 𝜀=1%; (d) 𝜀=0.5%.

5. Conclusion

This paper presents an application of homotopy perturbation inversion method to the image reconstruction for EIT. Simulation results indicate the feasibility and effectiveness of this method. Compared with the classical Landweber iteration method, our approach reduces the computational time and produces more satisfactory effect.

Acknowledgment

This work was supported by the National Natural Science Foundation of China under Grant no. 91230119.