Dynamic reconstruction algorithm based on the split Bregman iteration for electrical capacitance tomography

Abstract

Electrical capacitance tomography (ECT) is considered a promising visualization measurement technique, in which image reconstruction algorithms play an important role in real applications. In this paper, a dynamic reconstruction model, which integrates the ECT measurement information and the dynamic evolution information of the reconstruction objects, is presented. A generalized objective function, which simultaneously considers the ECT measurement information, the dynamic evolution information of the objects of interest, the temporal constraints and the spatial constraints, is proposed. An iteration scheme that integrates the beneficial advantages of the split Bregman iteration technique and the homotopy algorithm is developed for solving the proposed objective function. Numerical simulations are implemented to evaluate the feasibility and effectiveness of the proposed algorithm. For the cases simulated in this paper, the accuracy of the images reconstructed by the proposed algorithm is improved and the artifacts in the reconstructed images can be removed effectively, which indicates that the proposed algorithm is successful in solving ECT inverse problems.

Keywords

Dynamic reconstruction electrical capacitance tomography homotopy method inverse problems split Bregman iteration technique

Introduction

Applying process tomography technologies to investigate the complicated dynamic behaviours of the multiphase systems in process equipments to improve system efficiency and reduce pollutant emission has attracted increasing attention. Electrical capacitance tomography (ECT) is a non-invasive image-based measurement technique, which aims at acquiring information about the spatial material distribution of inaccessible objects to achieve the visualization monitoring of industrial processes. Owing to the distinct merits such as the non-intrusive sensing, high speed, low cost, easy implementation and high safety, which are important in monitoring the objects in a fast dynamic process, ECT is considered a promising visualization measurement technique. In recent years, ECT technology has been accepted as a potential tool for visualizing the distribution of the flame in the porous media, identifying the two-phase flow patterns and exploring the complicated dynamic behaviours of the multiphase system or process (Bennett et al., 2002; Du et al., 2006; Hamidipour and Larachi, 2010; Liu et al., 2008; Makkawi and Ocone, 2007; Niedostatkiewicz et al., 2009; Rimpilainen et al., 2011; Wang and Yang, 2010; Xiong et al., 2010; Zhao et al., 2010).

ECT technology attempts to reconstruct the permittivity distribution of the cross-section via an appropriate algorithm based on the capacitance measurement data, in which reconstructing high-quality images plays an important role in real applications. Owing to the ill-posed nature of the inverse problem, the ‘soft-field’ effect and the underdetermined problem in the process of ECT image reconstruction, exactly reconstructing original objects is challenging when the measurement objects are in a complicated dynamic process. In the past few years, improving the quality of the reconstructed images has attracted intensive attention. Consequently, various algorithms have been developed for ECT image reconstruction. Overall, ECT image reconstruction algorithms can be roughly divided into two categories: the static image reconstruction algorithms and the dynamic image reconstruction algorithms. Presently, various static image reconstruction algorithms have been proposed for ECT image reconstruction, such as the linear back-projection (LBP) method (Xie et al., 1992), the Tikhonov regularization method (Tikhonov and Arsenin, 1977), the Landweber iteration algorithm (Jang et al., 2006; Landweber, 1951; Yang et al., 1999), the neural network algorithm (Warsito and Fan, 2001), the truncated singular value decomposition (TSVD) method (Yang and Peng, 2003), the algebraic reconstruction technique (ART) and the simultaneous iterative reconstruction technique (SIRT) (Yang and Peng, 2003), the offline iteration and online reconstruction (OIOR) algorithm (Liu et al., 2004), the simulated annealing (SA) algorithm (Ortiz-Aleman et al., 2004), the generalized Tikhonov regularization methods (Fang, 2004; Lei et al., 2011; Soleimani and Lionheart, 2005; Wang et al., 2007), the genetic algorithm (GA) (Mou et al., 2005), the generalized vector sampled pattern matching method (Takei, 2006) and the level set algorithm (Banasiak and Soleimani, 2010; Kortschak et al., 2007). Owing to the fact that the image reconstruction process only involves a simple matrix–vector multiplication, the advantages of the LBP method include easy implementation, low computational complexity and cost; however, the quality of the images reconstructed by the LBP method is relatively low for the complicated reconstruction tasks. As a result, the LBP method is often applied to qualitative analysis. The standard Tikhonov regularization (STR) method is an efficient approach for solving inverse problems, and has found numerous applications in various fields. In essence, a Tikhonov regularization solution is a result of balancing the accuracy and stability of an inversion solution. When the STR method is applied to ECT image reconstruction, however, the quality of the images reconstructed by the algorithm is far from perfect, owing to the excessive smoothness effect. At present, different methods had been developed to improve the STR method (Fang, 2004; Soleimani and Lionheart, 2005; Wang et al., 2007). In the past years, the Landweber iteration algorithm has enjoyed popularity in the field of ECT image reconstruction. The advantages of the algorithm include the easy implementation and low computational complexity and cost, owing to the fact that only the gradient information of the objective function is utilized in the process of implementing image reconstruction. From the viewpoint of the numerical optimization, however, it is only the steepest descent algorithm and the rate of convergence is relatively low. In essence, the neural network method belongs to the optimization-based reconstruction technique, which had been applied for imaging two- and three-phase flow systems (Warsito and Fan, 2001). The TSVD method achieves the numerical stability of the solution by truncating the small singular values of the coefficient matrix. The advantages of the method involve the easy implementation and low computational complexity; in practice, however, determining the truncated singular values is challenging, especially when the singular values of the coefficient matrix are continuously decreasing. Additionally, when a large-scale problem is considered, the computational burden in implementing the SVD is heavy. Therefore, the method is seldom adopted in the field of ECT image reconstruction. The GA and the SA algorithms belong to the stochastic global optimal algorithms. A distinct advantage of this kind of algorithms is to obtain a possible global optimal solution. The OIOR algorithm is executed in two stages. The first stage is an offline process, in which the generalized inverse of the coefficient matrix is iteratively solved. In the second stage, the generalized inverse matrix is used for online image reconstruction. A distinct advantage of the OIOR method is that the algorithm can achieve online reconstruction; however, for the relatively complicated reconstruction objects, the quality of the images reconstructed by the OIOR method is not satisfactory. The ART method and the SIRT algorithm are commonly used for image reconstruction in X-ray computerized tomography. They are simple and effective, especially when the system matrix is very large (Yang and Peng, 2003). The ART method and the SIRT algorithm have many variations; for the complicated reconstruction of objects in ECT applications, however, the quality of the images reconstructed by both algorithms is far from satisfactory. The generalized vector sampled pattern matching method is an iterative technique; it has one characteristic that other iterative methods do not have: the inclusion of an objective function in the general solution (Takei, 2006). The level set method presents a novel insight for ECT image reconstruction and shows distinct advantages in shape reconstruction (Banasiak and Soleimani, 2010). More discussions on numerical performances and efficiencies of other algorithms can be found in Yang and Peng (2003) and Lei (2008).

In general, the above-mentioned algorithms have played an important role in promoting the development of ECT technology and found numerous applications. It is worth mentioning that static reconstruction algorithms are often used to image a dynamic object (Liu et al., 2008; Xiong et al., 2010). Unfortunately, these algorithms fail to consider the information about the temporal dynamics when reconstruction objects are in a dynamic process such as the multiphase flow and the visualization of combustion flame. Applications indicate that ECT measurement tasks often involve the time-varying objects. It would be more appropriate to image a dynamic object using a dynamic reconstruction algorithm that can take into account the dynamic behaviours of the objects of interest. In the field of ECT image reconstruction, however, the dynamic reconstruction algorithms do not yet received enough attention. Fortunately, several pioneers have investigated this subject, such as the particle filter (PF) method (Watzenig et al., 2007), the Kalman filter (KF) method (Soleimani et al., 2007) and the four-dimensional imaging method (Soleimani et al., 2009). In particular, in the fields of the electrical resistance tomography and the electrical impedance tomography, the dynamic reconstruction algorithms, such as the KF method (Kaipio et al., 1999; Seppänen et al., 2001; Vauhkonen et al., 2000), the extended Kalman filter (EKF) method (Kim et al., 2001; Tossavainen et al., 2006) and the unscented Kalman filter (UKF) method (Khambampati et al., 2009), have been proposed to improve the reconstruction quality. Overall, the investigations on the dynamic reconstruction algorithms in the field of ECT are far from perfect, and finding an effective dynamic reconstruction algorithm is highly desirable.

Studies reveal that one of drawbacks of ECT image reconstruction is the lack of enough information. Static reconstruction algorithms often consider ECT measurement information, while failing to pay attention to the dynamic evolution information when the measurement objects of interest are in a dynamic process. Therefore, introducing other useful information, such as the dynamic evolution information of a dynamic object, in the process of ECT image reconstruction may improve the reconstruction quality. Additionally, applications indicate that there is a close correlation between the images at different time instants when the measurement object is in a dynamic process (Xiong et al., 2010). So, considering such a temporal correlation will be necessary to improve reconstruction quality. According to the above discussions, simultaneously considering the temporal constraints, the spatial constraints, the ECT measurement information and the dynamic evolution information of a dynamic object may be essential to improve the quality of the reconstructed images. In practice, additionally, the ECT image reconstruction process is often reformulated into an optimization problem and developing an efficient algorithm to solve the optimization problem is essential. This paper presents a dynamic reconstruction model, which integrates the ECT measurement information and the dynamic evolution information of the objects of interest. A generalized objective function, which simultaneously takes into account the dynamic evolution information of the objects of interest, the ECT measurement information, the temporal constraints and the spatial constraints, is proposed. An iterative scheme that integrates the merits of the split Bregman iteration technique and the homotopy method is designed to solve the proposed objective function. Numerical simulations are implemented to validate the feasibility and effectiveness of the proposed algorithm.

The rest of this paper is organized as follows. The next section presents the static and dynamic models for ECT image reconstruction and a concise comparison between the both models is given. Then, a generalized objective function that simultaneously considers the dynamic evolution information of the objects of interest, the ECT measurement information, the temporal constraints and the spatial constraints is proposed. An iteration scheme, which integrates the advantages of the split Bregman iteration technique and the homotopy method, is developed for solving the proposed objective function. Numerical simulations are implemented to evaluate the feasibility and effectiveness of the proposed algorithm, and a concise discussion on the numerical results is given and conclusions are presented.

Model representation

Static reconstruction model

Within an ECT sensor, the relationship between the capacitance and the permittivity distribution is governed by (Yang and Peng, 2003):

C = \frac{Q}{V} = - \frac{1}{V} \int \int_{Γ} ε (x, y) \nabla ϕ (x, y) d Γ

(1)

where $Q$ is the electric charge; $ε (x, y)$ and $ϕ (x, y)$ represent the permittivity distribution and the electrical potential distribution; $V$ and $Γ$ stand for the potential difference between two electrodes forming the capacitance and the electrode surface.

In practice, the linearization approximation reconstruction model is often used to achieve fast reconstruction. In particular, when the inaccurate properties on the measurement data are considered, the static linearization reconstruction model can be formulated by (Yang and Peng, 2003):

S G = C + r

(2)

where $C$ represents an $m \times 1$ dimensional vector indicating the normalized capacitance values; $G$ is an $n \times 1$ dimensional vector standing for the normalized permittivity distribution, which indicates the grey level values in the reconstructed image; $S$ stands for a matrix of dimension $m \times n$ , which is called the sensitivity matrix in the field of ECT image reconstruction; $r$ is an $m \times 1$ dimensional vector indicating the noise in the capacitance data.

Dynamic reconstruction model

Equation (2) only considers ECT measurement information and fails to take into account the dynamic evolution information of the objects of interest. ECT measurement objects are often in a dynamic process such as the two-phase flow and the multiphase flow systems that can be depicted by the fluid dynamics equations; therefore, developing a reconstruction model that considers not only the ECT measurement information but also the dynamic evolution information of the measurement objects will be highly desirable for real applications. So, a dynamic imaging model can be formulated by:

G_{k + 1} = f (G_{k}, v_{k})

(3)

y_{k} = h (G_{k}, ϖ_{k})

(4)

where $f (\cdot)$ describes the process evolution information of the objects of interest, which can be formulated by a set of partial differential equations in the multiphase flow measurement; $h (\cdot)$ is called the measurement equation, which can be specified by Equation (1); $y_{k}$ represents the capacitance measurement data at time instant $k$ ; $v_{k}$ and $ϖ_{k}$ depict the uncertainties on the process evolution equation and the measurement equation; the subscript $k$ stands for the discrete time index.

In real applications, to achieve a fast reconstruction, Equations (3) and (4) can be approximated by the following linearization formulas:

G_{k} = F_{k} G_{k - 1} + v_{k}

(5)

C_{k} = S_{k} G_{k} + ϖ_{k}

(6)

where $F_{k}$ is the state-transition operator at time instant $k$ ; $S_{k}$ represents the measurement operator at time instant $k$ , which can be called the sensitivity matrix. If set $F_{k} = I$ , where $I$ is an identity matrix, Equation (5) can be named a purely random-walk evolution model, which is often used in practice when a better temporal dynamic model is not known (Butala et al., 2010). Ideally, the non-linear dynamic model can be incorporated into the evolution operator $F_{k}$ . In Equation (6), only the inaccuracy of the capacitance measurement data is considered. In fact, the model approximation distortions derived from the approximation and discretization of a real problem and physically implementing an ECT sensor in real applications may produce errors. As a result, developing a model that takes into account the inaccurate properties of the capacitance measurement data, the model approximation and the sensitivity matrix may be highly desirable for real applications. In this paper, a generalized measurement model considering the above-mentioned inaccuracies is proposed, which can be formulated by:

C_{k} = (S_{k} + E_{k}) G_{k} + B_{k} - n_{k}

(7)

where $B_{k}$ is an $m \times 1$ dimensional vector representing the linearization approximation errors or the model errors at time instant $k$ ; $E_{k}$ represents a matrix of dimension $m \times n$ indicating the inaccurate property on the sensitivity matrix at time instant $k$ ; $n_{k}$ is an $m \times 1$ dimensional vector indicating the measurement noises at time instant $k$ .

Comparisons of the both reconstruction models

Compared with the static reconstruction model (Equation 2), the dynamic reconstruction model (Equations 5 and 7) has the following properties:

One of the main drawbacks of ECT image reconstruction is the lack of adequate information. The static model only considers ECT measurement information; however, the dynamic model simultaneously focuses on the dynamic evolution information of a dynamic object and the ECT measurement information, which increases the quantity of information in the process of image reconstruction; consequently, the reconstruction accuracy may be improved in theory. A distinct advantage of the dynamic model is that it integrates the dynamic evolution information of a dynamic object and the ECT measurement information.

In the case of the static reconstruction model, the solution only reflects an instantaneous measurement $C$ , whereas in the case of the dynamic reconstruction model it reflects a sequence of temporally successive measurement $C_{k}$ , such that the temporal correlations of the dynamic objects of interest should be imposed. In other words, the dependence of the measurement $C_{k}$ on the evolution of $G_{k}$ is explicitly considered in the dynamical reconstruction model. If the evolution of $G_{k}$ does not follow any dynamics, the dynamic reconstruction model is reduced to the static reconstruction model. Obviously, the static reconstruction model is a special case of the dynamic reconstruction model.

The dynamic reconstruction model simultaneously takes into account the inaccurate properties of the process evolution equation of the objects of interest and the ECT measurement information, which is especially suitable for real applications owing to the fact that it is difficult to depict the complicated dynamic behaviours of a dynamic object in practice and to achieve an exact capacitance measurement exactly.

The dynamic reconstruction model considers the model approximation errors derived mainly from the approximation and discretization of a real problem, the inaccuracy on the sensitivity matrix derived from physically implementing an ECT sensor in real applications, and the measurement noise, which distinctly differs from traditional liner and non-linear image reconstruction algorithms.

Design of the objective function

Directly solving dynamic reconstruction model Equations (5) and (7) is challenging. One of the popular approaches is to formulate the solution of Equations (5) and (7) into an optimization problem. Similar to the derivation process in Beck and Ben-Tal (2006), according to the Tikhonov regularization method and the optimization theory, Equations (5) and (7) can be reformulated as:

{\begin{matrix} min_{G_{k}, E_{k}, B_{k}, n_{k}} {‖ E_{k} ‖_{F}^{2} + ‖ n_{k} ‖_{2}^{2} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2}} \\ s . t . C_{k} = (S_{k} + E_{k}) G_{k} + B_{k} - n_{k} \end{matrix}

(8)

where operator $‖ \cdot ‖_{F}$ defines the Frobenius norm for a matrix and $α > 0$ stands for the regularization parameter. According to the optimization theory, Equation (8) can be rewritten as a double minimization problem (Beck and Ben-Tal, 2006):

{\begin{matrix} min_{G_{k}, B_{k}} min_{E_{k}, n_{k}} {‖ E_{k} ‖_{F}^{2} + ‖ n_{k} ‖_{2}^{2} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2}} \\ s . t . C_{k} = (S_{k} + E_{k}) G_{k} + B_{k} - n_{k} \end{matrix}

(9)

Taking the inner minimization problem into account, the following expression can be obtained:

{\begin{matrix} min_{E_{k}, n_{k}} {‖ E_{k} ‖_{F}^{2} + ‖ n_{k} ‖_{2}^{2} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2}} \\ s . t . C_{k} = (S_{k} + E_{k}) G_{k} + B_{k} - n_{k} \end{matrix}

(10)

The Lagrangian function of Equation (10) can be formulated by:

min L (E_{k}, n_{k}, λ_{k}) = {\begin{matrix} ‖ E_{k} ‖_{F}^{2} + ‖ n_{k} ‖_{2}^{2} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2} + \\ 2 λ_{k}^{T} ((S_{k} + E_{k}) G_{k} + B_{k} - (n_{k} + C_{k})) \end{matrix}}

(11)

where $λ_{k}$ is a Lagrangian multiplier.

According to the optimization theory, we can obtain the following expressions:

2 E_{k} + 2 λ_{k} G_{k}^{T} = 0, (\nabla_{E_{k}} L = 0)

(12)

2 n_{k} - 2 λ_{k} = 0, (\nabla_{n_{k}} L = 0)

(13)

C_{k} = (S_{k} + E_{k}) G_{k} + B_{k} - n_{k}, (feasibility)

(14)

The equality $λ_{k} = n_{k}$ can be obtained from Equation (13). Submitting this equality into Equation (12) yields:

E_{k} = - n_{k} G_{k}^{T}

(15)

Submitting Equation (15) to Equation (14) yields:

(S_{k} - n_{k} G_{k}^{T}) G_{k} + B_{k} = C_{k} + n_{k}

(16)

Thereby,

n_{k} = \frac{S_{k} G_{k} - C_{k} + B_{k}}{‖ G_{k} ‖_{2}^{2} + 1}

(17)

Submitting Equation (17) to Equation (15), the following result can be obtained:

E_{k} = - \frac{(S_{k} G_{k} - C_{k} + B_{k}) G_{k}^{T}}{‖ G_{k} ‖_{2}^{2} + 1}

(18)

Similarly, submitting Equations (17) and (18) into Equation (10), we can show that its value is equal to $\frac{‖ S_{k} G_{k} - C_{k} + B_{k} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2}$ . As a consequence, Equation (8) can be reduced to:

min_{G_{k}, B_{k}} {\frac{‖ S_{k} G_{k} - C_{k} + B_{k} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2}}

(19)

For concise depiction, Equation (19) can be reformulated as:

min_{G_{k}, B_{k}} {‖ D_{k} ‖_{2}^{2} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2}}

(20)

where

{D_{k,}}_{j} = \frac{{S_{k,}}_{j} G_{k} - {C_{k,}}_{j} + {B_{k,}}_{j}}{\sqrt{‖ G_{k} ‖_{2}^{2} + 1}}

(21)

where $S_{k, j}$ represents the $j th$ row of matrix $S_{k}$ at time instant $k$ ; $C_{k, j}$ is the $j th$ element of vector $C_{k}$ at time instant $k$ ; $B_{k, j}$ stands for the $j th$ element of vector $B_{k}$ at time instant $k$ . Equation (19) is the total least squares (TLS) estimation when $B_{k} = 0$ ; obviously, the TLS estimation is a special case of Equation (19).

ECT image reconstruction process is a typical ill-posed problem, methods that ensure a stable numerical solution while improving the reconstruction quality should be employed. According to the Tikhonov regularization method, Equation (19) can be reformulated as:

min_{G_{k}, B_{k}} {\frac{‖ S_{k} G_{k} - C_{k} + B_{k} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2} + α_{1} Ω_{1} (B_{k}) + α_{2} Ω_{2} (G_{k})}

(22)

where $α_{1} > 0$ and $α_{2} > 0$ are the regularization parameters; $Ω_{1} (B_{k})$ and $Ω_{2} (G_{k})$ represent the stabilizing functions (or the constrained functions).

In Equation (22), merely the spatial constraints are considered; however, the temporal constraints are not considered. Applications indicate that ECT measurement objects are often in a dynamic process and there is a close correlation in the reconstructed images at different time instants. As a result, considering such a temporal correlation may be essential to improve the reconstruction quality. Schmitt and Louis (2002) proposed the temporal constraints. Therefore, the temporal constraint is introduced to Equation (22) in order to emphasize the temporal correlation of a dynamic object, which can be formulated as:

min_{G_{k}, B_{k}} {\begin{matrix} \frac{‖ S_{k} G_{k} - C_{k} + B_{k} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1 |} ‖_{2}^{2} + \\ α_{1} Ω_{1} (B_{k}) + α_{2} Ω_{2} (G_{k}) + α_{3} Ω_{3} (G_{k} - G_{k - 1}) \end{matrix}}

(23)

where $α_{3} > 0$ is a regularization parameter and $Ω_{3} (G_{k} - G_{k - 1})$ stands for a temporal constraint function.

The design of functions $Ω_{1}$ and $Ω_{2}$ is vital to the applications of Equation (23). For easy computation, functions $Ω_{1}$ and $Ω_{2}$ are defined as:

Ω_{1} (B_{k}) = ‖ B_{k} ‖_{1}

(24)

Ω_{2} (G_{k}) = ‖ W_{2} G_{k} ‖_{2}^{2}

(25)

where operator $‖ \cdot ‖_{1}$ stands for the 1-norm and $W_{2}$ represents a weighted matrix.

Furthermore, according to the idea in Schmitt and Louis (2002), the temporal constraint function is defined as:

Ω_{3} (G_{k} - G_{k - 1}) = ‖ W_{3} (G_{k} - G_{k - 1}) ‖_{1}

(26)

where $W_{3}$ is a weighted matrix.

Following the above discussions, a generalized objective function for ECT image reconstruction can be formulated as:

min_{G_{k}, B_{k}} {\begin{matrix} \frac{‖ S_{k} G_{k} - C_{k} + B_{k} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2} + α_{1} ‖ B_{k} ‖_{1} + \\ α_{2} ‖ W_{2} G_{k} ‖_{2}^{2} + α_{3} ‖ W_{3} (G_{k} - G_{k - 1}) ‖_{1} \end{matrix}}

(27)

In essence, Equation (27) is a generalized Tikhonov regularization function, and the selection of the regularization parameter plays an important role in real applications. At present, various methods, such as the discrepancy principle method (Morozov, 1984), the generalized cross-validation method (Golub et al., 1979), the generalized Arcangeli criteria method (Engl, 1987), the L-curve method (Hansen, 1992) and the Engl criteria method (Engl et al., 1996), are available for determining the regularization parameter, and more discussions on the selection of the regularization parameter can be found in Wang (2007). Owing to the complexity and particularity of the problem, determining an optimal regularization parameter is challenging at present and in the field of ECT image reconstruction the regularization parameter is often chosen empirically (Yang and Peng, 2003). In this paper, we determined the regularization parameter empirically.

Several desirable properties can be found from Equation (27):

Equation (27) takes into account the dynamic evolution information of the objects of interest and the ECT measurement information, which increases the quantity of information in the process of ECT image reconstruction.

Equation (27) considers the spatial constraints and the temporal constraints, which is a distinct feature that differs from other image reconstruction algorithms.

Equation (27) takes into account the inaccurate properties on the dynamic evolution equation of the measurement objects and the measurement equation owing to the fact that it is hard to describe the dynamic evolution information of a dynamic object and to achieve the exact capacitance measurements exactly.

ECT image reconstruction process is a typical ill-posed problem, and its solution is unstable. Methods that ensure the numerical stability of a solution while improving the quality of reconstruction should be used. In Equation (27), the Tikhonov regularization technique is introduced to stabilize a numerical solution, which is especially important to real applications.

Solving of the objective function

Previously, the original dynamic reconstruction model is reformulated into an optimization problem and developing an efficient algorithm is important for the applications of Equation (27). In this section, the split Bregman iteration method and the homotopy algorithm are introduced and an iterative scheme that integrates the advantages of the both methods is developed to solve Equation (27).

The split Bregman iteration technique

The split Bregman iteration technique is an efficient algorithm for solving the $ℓ_{1}$ -regularization problems, which effectively separates the $ℓ_{2}$ and $ℓ_{1}$ portions of the objective function (Cai et al., 2009a, 2009b; Goldstein and Osher, 2009; Goldstein et al., 2010; Osher et al., 2010; Yin et al., 2008; Zhang et al., 2011). Owing to the high computational efficiency and the satisfactory numerical performance, the split Bregman iteration technique has found wide applications in the fields such as signal and image processing, and compressed sensing. In essence, the split Bregman iteration algorithm has been developed for solving the following optimization problem (Cai et al., 2009a; Goldstein and Osher, 2009):

min_{u} {{‖ Φ (u) ‖}_{1} + J (u)}

(28)

where $J$ is a convex function. Introducing an additional variable $d$ and an equality constraint $d = Φ (u)$ , Equation (28) can be reformulated into:

{\begin{matrix} min_{u, d} {{‖ d ‖}_{1} + J (u)} \\ s . t . d = Φ (u) \end{matrix}

(29)

Finally, minimizing Equation (29) is equivalent to solving the following sub-problems (Goldstein and Osher, 2009):

(u^{κ + 1}, d^{κ + 1}) = min_{u, d} {{‖ d ‖}_{1} + J (u) + \frac{λ}{2} ‖ d - Φ (u) - b^{κ} ‖_{2}^{2}}

(30)

b^{κ + 1} = b^{κ} + (Φ (u^{κ + 1}) - d^{κ + 1})

(31)

where $κ$ is the index of iterations.

Equation (30) can be minimized efficiently by iteratively minimizing with respect to $u$ and $d$ , in two steps (Goldstein and Osher, 2009):

u^{κ + 1} = min_{u} {J (u) + \frac{λ}{2} ‖ d^{κ} - Φ (u) - b^{κ} ‖_{2}^{2}}

(32)

d^{κ + 1} = min_{d} {{‖ d ‖}_{1} + \frac{λ}{2} ‖ d - Φ (u^{κ + 1}) - b^{κ} ‖_{2}^{2}}

(33)

Clearly, the speed of the split Bregman iteration technique is dependent mainly on how fast we can solve each of the two sub-problems in Equations (32) and (33). Equation (32) can be solved by the traditional optimization algorithms; however, Equation (33) can be solved quickly by the shrinkage algorithm (Yin et al., 2008; Osher et al., 2010), which can be formulated by:

d_{i}^{κ + 1} = shrinkage ({(Φ (u^{κ + 1}) + b^{κ})}_{i}, 1 / λ), i = 1, 2, \dots, n

(34)

where for $y$ and $a$ , we define:

shrinkage (y, a) = {\begin{matrix} y - a, & y \in (a, \infty) \\ 0, & y \in [- a, a] \\ y + a, & y \in (- \infty, - a) \end{matrix}

(35)

The proposed iterative scheme

According to the alternating direction iteration optimization method, Equation (27) can be reformulated into:

B_{k}^{κ + 1} = min_{B_{k}} {\frac{‖ S_{k} G_{k}^{κ} - C_{k} + B_{k} ‖_{2}^{2}}{‖ G_{k}^{κ} ‖_{2}^{2} + 1} + α_{1} {‖ B_{k} ‖}_{1}}

(36)

G_{k}^{κ + 1} = min_{G_{k}} {\begin{matrix} \frac{‖ S_{k} G_{k} - C_{k} + B_{k}^{κ + 1} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2} + \\ α_{2} ‖ W_{2} G_{k} ‖_{2}^{2} + α_{3} ‖ W_{3} (G_{k} - G_{k - 1}) ‖_{1} \end{matrix}}

(37)

Equation (37) can be reformulated into a constrained optimization problem, which can be expressed by:

{\begin{matrix} G_{k}^{κ + 1} = min_{G_{k}} {\frac{‖ S_{k} G_{k} - C_{k} + B_{k}^{κ + 1} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2} + α_{2} ‖ W_{2} G_{k} ‖_{2}^{2} + {‖ d_{k, 1} ‖}_{1}} \\ s . t . d_{k, 1} = α_{3} W_{3} (G_{k} - G_{k - 1}) \end{matrix}

(38)

Following the discussions presented in previous sections, a sequential iteration scheme can be obtained for solving Equation (27), which can be summarized as follows:

Step 1. Specify the algorithmic parameters and the initial values.

Step 2. Solve Equation (36) using the shrinkage algorithm and obtain the estimation of variable $B$ .

Step 3. Update variable $G$ by solving Equation (39).

G_{k}^{κ + 1} = min_{G_{k}} {\begin{matrix} \frac{‖ S_{k} G_{k} - C_{k} + B_{k}^{κ + 1} ‖_{2}^{2}}{‖ G_{k} ‖_{2}^{2} + 1} + α ‖ G_{k} - F_{k} G_{k - 1} ‖_{2}^{2} + \\ α_{2} ‖ W_{2} G_{k} ‖_{2}^{2} + \frac{λ_{1}}{2} ‖ d_{k, 1} - α_{3} W_{3} (G_{k} - G_{k - 1}) - b_{k, 1}^{κ} ‖_{2}^{2} \end{matrix}}

(39)

Step 4. Update variable $d$ by solving Equation (40).

d_{k, 1}^{κ + 1} = min_{d_{k}} {{‖ d_{k, 1} ‖}_{1} + \frac{λ_{1}}{2} ‖ d_{k, 1} - α_{3} W_{3} ({G_{k}}^{κ + 1} - G_{k - 1}) - b_{k, 1}^{κ} ‖_{2}^{2}}

(40)

Step 5. Update variable $b$ using Equation (41).

b_{k, 1}^{κ + 1} = b_{k, 1}^{κ} + α_{3} W_{3} (G_{k}^{κ + 1} - G_{k - 1}) - d_{k, 1}^{κ + 1}

(41)

Step 6. Loop to Step 2 until a predetermined stopping criterion is satisfied.

Several desirable properties can be found from the above iterative scheme:

It can be found from the above iterative scheme that a complicated optimization problem is separated into several simple sub-problems, which effectively reduces the computational complexity and cost.

Equation (39) can be easily solved by traditional optimization algorithms such as the conjugate gradient method, the Newton algorithm and the homotopy method.

Equations (36) and (40) can be fast solved by the shrinkage algorithm, which further reduces the computational complexity and cost. This feature is highly suitable for real applications.

Homotopy method

In the previous subsection, several unconstrained optimization problems, Equations (36), (39) and (40), need to be sequentially solved. Equations (36) and (40) can be solved quickly by the shrinkage algorithm, and developing an efficient optimization algorithm to solve Equation (39) plays a vital role in the applications of the iteration algorithm. Owing to advantages such as the easy implementation and the low computational complexity and cost, the homotopy method is used to solve the unconstrained optimization problem. For concise depiction, the homotopy method is introduced according to variable $G$ , and the superscript and subscript are omitted.

According to the optimization theory, minimizing Equation (39) is equivalent to solving the following system of non-linear Equations (Guo et al., 2008):

ℜ (G) = 0

(42)

where $ℜ (G)$ is the gradient vector of Equation (39).

The general homotopy paradigm involves embedding the equation to be solved, $ℜ (G) = 0$ , in a system of equations of one higher dimension, $ℵ (G, ƛ) = 0$ , with the introduction of one more variable $ƛ$ , called the homotopy parameter. Typically, $ƛ$ is restricted to the range [0, 1], and the embedding is done so that the augmented system $ℵ (G, ƛ) = 0$ is easy to solve and reduces to the original system when $ƛ = 1$ (Li et al., 1987; Huang et al., 2004; Ma, 2005). To solve the equation $ℵ (G, ƛ) = 0$ , homotopy parameter $ƛ$ can be divided into (Huang et al., 2004; Li et al., 1987; Ma, 2005):

0 = ƛ_{0} < ƛ_{1} < \dots < ƛ_{ν} = 1

(43)

At each point $ƛ_{i}$ , solve the following discrete equations:

ℵ (G, ƛ_{i}) = 0, i = 1, 2, \dots, ν

(44)

and set the solution $G (ƛ_{i})$ as the initial value of equation $ℵ (G, ƛ_{i + 1}) = 0$ ; repeat this process until $ƛ_{i} = 1$ .

For simplicity, the fixed-point homotopy is employed in this paper, which can be formulated by (Li et al., 1987; Ma, 2005):

ℵ (G, ƛ_{i}) = (1 - ƛ_{i}) (G - G_{0}) + ƛ_{i} ℜ (G) = 0 r

(45)

For concise notation, Equation (45) can be reformulated by:

G = Φ (G, ƛ_{i})

(46)

where $Φ (\cdot)$ is a function of variables $G$ and $ƛ_{i}$ .

Equation (46) is a system of non-linear equations, and the fixed-point iterative technique is used to solve Equation (46) in this work, which can be described by (Huang et al., 2004):

G^{κ + 1} = Φ (G^{κ}, ƛ_{i})

(47)

It can be known in advance that the inversion solution belongs in the range $[Θ_{1}, Θ_{2}]$ ; therefore, the iterative scheme is slightly modified according to prior information. Consequently, a projected operator is introduced to the iterative scheme:

G^{κ + 1} = Project {Φ (G^{κ}, ƛ_{i})}

(48)

where

Project [f (x)] = {\begin{matrix} Θ_{1}, & if f (x) < Θ_{1} \\ f (x), & if Θ_{1} \leq f (x) \leq Θ_{2} \\ Θ_{2} & if f (x) > Θ_{2} \end{matrix}

(49)

Numerical simulations and discussion

For concise depiction, the proposed iterative scheme is called the dynamic image reconstruction (DIR) algorithm. In this section, the static and dynamic reconstruction cases are implemented to evaluate the feasibility and effectiveness of the DIR algorithm, and the quality of the images reconstructed by the DIR algorithm are compared with the projected Landweber iteration (PLI) method, the ART algorithm and the KF method. A 12-electrode square ECT sensor is selected for simulations and an image is presented using 32×32 pixels. For easy computation, $S_{k} = S$ and $F_{k}$ is defined as an identity matrix (Soleimani et al., 2007). The initial values for the compared algorithms are computed by the STR method. All algorithms are implemented using the MATLAB 7.0 software on a PC with a Pentium IV 2.4-GHz CPU with 4 Gbytes of memory.

Case 1

Obviously, the DIR algorithm can achieve the single-frame reconstruction (or the static reconstruction) and the multi-frame reconstruction (or the dynamic reconstruction). In this section, the single-frame reconstruction cases are first implemented to evaluate the feasible and effectiveness of the DIR algorithm and the reconstruction quality is compared with the PLI algorithm and the ART algorithm.

Figure 1 shows three static permittivity distributions chosen for simulations, where the black colour stands for the high permittivity materials with a value of 2.6 and the white colour represents the low permittivity materials with a value of 1.0. The diameters of the cylinders in Figure 1(a) and (b) are 20 mm, the permittivity of the cylinders is 2.6, and the permittivity of the rest of the reconstruction region is 1.0. The diameter of the cylinder in Figure 1(c) is 30 mm, the permittivity of the cylinder is 2.6 and the permittivity of the rest of the reconstruction region is set as 1.0.

Tables 1 and 2 list the algorithmic parameters for the PLI algorithm and the ART algorithm, respectively. In the DIR algorithm, $ƛ = 1$ , $α = α_{3}$ ; $W_{2} = diag (1 / β_{1}, \dots, 1 / β_{n})$ and $W_{3} = diag (1 / τ_{1}, \dots, 1 / τ_{n})$ , in which $τ_{j} = 1 + (G_{j}^{2} + ξ)^{1 / 2}$ , $β_{j} = ((1 + (G_{j}^{2} + ξ)^{ϑ / 2}) \times (G_{j}^{2} + ξ)^{(2 - ϑ) / 2} + ξ)$ , $0 < ϑ < 2$ , and $G$ is a latest previous solution; the values of $ϑ$ , $α$ and $α_{2}$ are selected according to a specific reconstruction task. The rest of the algorithmic parameters for the DIR algorithm are listed in Table 3. Figures 2 –4 present the images reconstructed by the PLI algorithm, the ART method and the DIR algorithm. The image errors are shown in Table 4.

Table 1

Algorithmic parameters for the projected Landweber iteration (PLI) algorithm

Algorithmic parameters	Figure 2(a)	Figure 2(b)	Figure 2(c)
Relaxation factor	1	1	1
Number of iteration	480	500	111

Table 2

Algorithmic parameters for the algebraic reconstruction technique (ART) algorithm

Algorithmic parameters	Figure 3(a)	Figure 3(b)	Figure 3(c)
Relaxation factor	1	1	1
Number of iteration	530	440	114

Table 3

Algorithmic parameters for the dynamic image reconstruction (DIR) algorithm

Algorithmic parameters	Figure 4(a)	Figure 4(b)	Figure 4(c)
$ϑ$	0.1	0.1	1.0
$α_{2}$	0.0005	0.00035	0.0006
$α$	0.00003	0.00004	0.00004
Number of iteration	8	8	8

Table 4

Image error (%)

Algorithms	Figure 1(a)	Figure 1(b)	Figure 1(c)
PLI	12.61	17.71	15.22
ART	12.44	18.47	15.07
DIR	Near zero	Near zero	Near zero

PLI, projected Landweber iteration; ART, algebraic reconstruction technique; DIR, dynamic image reconstruction.

Figure 2 shows the images reconstructed by the PLI algorithm. Numerical simulation results indicate that one of the distinct advantages of the PLI algorithm is that the numerical implementation is simple, and the computational complexity and cost are relatively low. For the cases simulated in this paper, unfortunately, the quality of the images reconstructed by the PLI method is far from perfect and the distortions are relatively serious, which indicate that it is difficult for the PLI algorithm to reconstruct the original objects exactly.

The images reconstructed by the ART algorithm are illustrated in Figure 3. Similarly, numerical simulation results reveal that the implementation of the ART algorithm is relatively easy, and the computational complexity and cost are low. However, it can be found from Figure 3 that for the cases simulated in this section, the quality of the images reconstructed by the ART algorithm is not satisfactory and the distortions of the reconstructed images are relatively large, which indicates that it is difficult for the ART algorithm to reconstruct detailed information on the original objects.

Figure 4 shows the images reconstructed by the DIR algorithm. Numerical simulation results indicate that the DIR algorithm shows a satisfactory numerical performance and can ensure the numerical stability of the solution owing to the fact that the Tikhonov regularization technique is used to stabilize the numerical solution. It can be seen from Figure 4 that the quality of the images reconstructed by the DIR algorithm is improved compared with the PLI algorithm and the ART method, the artifacts in the reconstructed images can be effectively eliminated and the distortion of the reconstructed images is very slight, which indicates that the DIR algorithm is successful in solving ECT inverse problems.

The image errors are listed in Table 4. The second and third rows in Table 4 are the image errors for the PLI algorithm and the ART algorithm. It can be observed that the quality of the images reconstructed by the PLI algorithm and the ART algorithm is not satisfactory and the distortions are relatively large. The fourth row in Table 4 shows the image errors for the DIR algorithm. As can be expected, the quality of the images reconstructed by the DIR algorithm is higher than the PLI method and the ART algorithm. For the cases simulated in this section, the DIR algorithm gives the smallest image errors, which indicates that the DIR algorithm is a promising candidate for solving ECT image reconstruction problems.

The CPU times for reconstructing subfigures a, b and c in Figure 2 using the PLI algorithm are about 0.68, 0.71 and 0.18 s. The CPU times for reconstructing subfigures a, b and c in Figure 3 using the ART algorithm are about 6.86, 5.71 and 1.48 s. The CPU times for reconstructing subfigures a, b and c in Figure 4 using the DIR algorithm are about 3.28 s. It can be found that for the cases simulated in this section, the reconstruction speed of the PLI algorithm is faster than the ART method and the DIR algorithm; however, the quality of the image reconstructed by the DIR algorithm is higher than the PLI method and the ART algorithm. In the future, more investigations on accelerating the DIR algorithm should be further implemented.

Case 2

Studies indicate that the ECT image reconstruction process is a typically ill-posed problem, and a small perturbation in the input data may bring a large fluctuation in the final results, which may make the inversion solution meaningless. Applications reveal that the noises in the capacitance measurement data are inevitable and complicated. A successful algorithm should have the advantage of dealing with the inaccurate properties on the capacitance measurement data. In this section, the noise-contaminated capacitance data are used to evaluate the numerical performance and robustness of the DIR algorithm, and the definition of the noise level can be found in Yang and Peng (2003). Algorithmic parameters for the DIR algorithm are the same as Table 3. Figures 5 –7 illustrate the images reconstructed by the DIR algorithm under the noise levels of 1%, 6% and 15%. The image errors under different noise levels for the DIR algorithm are listed in Table 5.

Table 5

Image errors for the dynamic image reconstruction (DIR) algorithm under different noise levels (%)

Noise levels	Figure 1(a)	Figure 1(b)	Figure 1(c)
1%	Near zero	Near zero	Near zero
6%	Near zero	Near zero	Near zero
15%	7.97	13.36	5.78

Figures 5 –7 present the images reconstructed by the DIR algorithm under the noise levels of 1%, 6% and 15%. As can be expected, the DIR algorithm shows satisfactory numerical performance and the favourable robustness, and the artifacts in the reconstructed images can be effectively removed, which indicates that the DIR algorithm is successful in treating with the inaccuracies on the measurement data. Meanwhile, it can be observed from Table 5 that the quality of the images reconstructed by the DIR algorithm is acceptable under the noise levels simulated in the section. In particular, when the noise level is 15%, the image errors for the original reconstruction objects Figures 1(a–c) are 7.97%, 13.36% and 5.78%, which indicates that the DIR algorithm is robust. Additionally, it can be seen from Figures 5 –7 and Table 5 that the quality of the reconstructed images decreases with the augmentation of the noise levels, which indicates that the inaccurate properties on the capacitance measurement data should be treated seriously, and this issue should be further investigated in the future.

Case 3

Applications indicate that ECT reconstruction objects are often in a dynamic process; a successful algorithm should be able to treat with such reconstruction tasks. In order to evaluate further the feasibility and efficiency of the DIR algorithm, a complicated dynamic reconstruction case, where the shapes and locations of the permittivity distributions change with the index of time, is simulated in this section. Subfigures (a), (b) and (c) in Figure 8 show the different permittivity distributions at discrete time indexes $k$ , $k + 1$ and $k + 2$ . In Figure 8(a) and (c), the diameters of the cylinders are 20 mm, the permittivity is defined as 2.6 and the permittivity of the rest of reconstruction region is 1.0. In the case of Figure 8(b), the diameter of the cylinder is 30 mm, the permittivity is defined as 2.6 and the permittivity of the rest of reconstruction region is 1.0.

The specific procedures on simulating the dynamic image reconstruction process can be summarized as: (1) the capacitance data (such as $C_{k}$ , $C_{k + 1}$ and $C_{k + 2}$ ) at different time instants (such as $k$ , $k + 1$ and $k + 2$ ) is firstly computed according to the presupposed permittivity distributions; (2) the capacitance data at time instant $k$ is used to implemented image reconstruction by solving Equation (27), and the permittivity distribution $G_{k}$ at time instant $k$ can be obtained; (3) based on the permittivity distribution $G_{k}$ at the time instant $k$ and the capacitance data $C_{k + 1}$ at time instant $k + 1$ , the permittivity distribution $G_{k + 1}$ at the time instant $k + 1$ can be obtained by solving Equation (27) under the condition of replacing the discrete time index $k$ with $k + 1$ ; (4) increase the index of time, and solve Equations (27) until the maximum discrete time index is reached. More details on the implementation of the sequential dynamic reconstruction methods can be found in Soleimani et al. (2007) and Khambampati et al. (2009).

Tables 6 and 7 list the algorithmic parameters for the PLI algorithm and the ART method. Table 8 shows the algorithmic parameters for the DIR algorithm and the rest of the algorithmic parameters are the same as previously given. The images reconstructed by the PLI algorithm, the ART method and the DIR algorithm are illustrated in Figures 9 –11. The image errors for the compared algorithms are presented in Table 9.

Table 6

Algorithmic parameters for the projected Landweber iteration (PLI) algorithm

Algorithmic parameters	Figure 9(a)	Figure 9(b)	Figure 9(c)
Relaxation factor	1	1	1
Number of iteration	530	111	500

Table 7

Algorithmic parameters for the algebraic reconstruction technique (ART) algorithm

Algorithmic parameters	Figure 10(a)	Figure 10(b)	Figure 10(c)
Relaxation factor	1	1	1
Number of iteration	405	114	320

Table 8

Algorithmic parameters for the dynamic image reconstruction (DIR) algorithm

Algorithmic parameters	Figure 11(a)	Figure 11(b)	Figure 11(c)
$ϑ$	0.1	1.0	0.1
$α_{2}$	0.00025	0.0006	0.00045
$α$	0.000035	0.00004	0.00004
Number of iteration	8	8	8

Table 9

Image error (%)

Algorithms	Figure 8(a)	Figure 8(b)	Figure 8(c)
PLI	18.84	15.22	15.09
ART	20.45	15.07	16.37
DIR	Near zero	Near zero	Near zero

PLI, projected Landweber iteration; ART, algebraic reconstruction technique; DIR, dynamic image reconstruction.

Figures 9 and 10 present the images reconstructed by the PLI algorithm and the ART method under different time indexes. It can be observed that the distinct advantages of the PLI algorithm and the ART method include the easy implementation, the low computational complexity and cost. For the dynamic reconstruction cases simulated in the section, however, the quality of the images reconstructed by the PLI algorithm and the ART method is far from satisfactory and the distortions are relatively serious. Additionally, it is worth mentioning that ECT reconstruction objects are often in a complicated dynamic evolution process; the PLI algorithm and the ART method merely consider the ECT measurement information and fail to consider the dynamic evolution information of the objects of interest. As a result, it is hard for the PLI algorithm and the ART method to reconstruct the complicated dynamic objects exactly.

The images reconstructed by the DIR algorithm under different time indexes are shown in Figure 11. As can be expected, when the dynamic reconstruction objects are considered, the DIR algorithm shows satisfactory numerical performance, the quality of the images reconstructed by the DIR algorithm is improved compared with the PLI algorithm and the ART method, and the artifacts in the reconstructed images can be effectively eliminated. In addition, it can be found from Figures 9 –11 that the quality of the images reconstructed by the DIR algorithm is higher than the PLI algorithm and the ART method, which indicates that the DIR algorithm is successful in solving ECT inverse problems.

Table 9 lists the image errors for the PLI algorithm, the ART method and the DIR algorithm. It can be observed from Table 9 that for the dynamic reconstruction objects, the quality of the images reconstructed by the PLI algorithm and the ART method is not satisfactory and the image errors are relatively large. It can be seen from Table 9 that for the dynamic reconstruction cases simulated in this section, the DIR algorithm gives the smallest image errors, which indicates that the DIR algorithm is successful in solving ECT image reconstruction problems.

Case 4

In order to evaluate further the feasibility and effectiveness of the DIR algorithm, another dynamic reconstruction case is simulated in this section. Especially, the quality of the images reconstructed by the DIR algorithm is compared with the KF method (one of the typical dynamic reconstruction algorithms). Subfigures (a), (b) and (c) in Figure 12 represents the different locations of the dynamic reconstruction objects at discrete time indexes $k$ , $k + 1$ and $k + 2$ . The diameters of the cylinders in Figure 12 are 20 mm, the permittivity is defined as 2.6 and the permittivity of the rest of reconstruction region is 1.0. Table 10 lists the algorithmic parameters for the GIR algorithm, and the rest of the algorithmic parameters are the same as previously shown. In the KF method, the covariance of the process evolution is defined as $Q_{k} = 1 \times 10^{- 4} I_{n}$ , where $I_{n}$ is a $n \times n$ identity matrix. In order to simulate a real environment, the variance of the measurement noise is defined as $R_{k} = 1 \times 10^{- 6} I_{m}$ . Figures 13 and 14 illustrate the images reconstructed by the KF method and the DIR algorithm. The image errors for the KF method and the DIR algorithm are presented in Table 11.

Table 10

Algorithmic parameters for the dynamic image reconstruction (DIR) algorithm

Algorithmic parameters	Figure 14(a)	Figure 14(b)	Figure 14(c)
$ϑ$	0.1	0.1	0.1
$α_{2}$	0.0008	0.0004	0.0005
$α$	0.0001	0.00006	0.00006
Number of iteration	8	8	8

Table 11

Image error (%)

Algorithms	Figure 12(a)	Figure 12(b)	Figure 12(c)
KF	15.57	23.48	21.17
DIR	Near zero	Near zero	Near zero

KF, Kalman filter; DIR, dynamic image reconstruction.

The images reconstructed by the KF method under different time indexes are illustrated in Figure 13. It can be observed that the distinct advantages of the KF method include the easy numerical implementation and the low computational complexity and cost. For the dynamic reconstruction cases simulated in the section, however, the quality of the images reconstructed by the KF method is far from satisfactory, the distortions are relatively serious and the artifacts are comparatively more. In addition, the results in Figure 13 indicate that it is hard for the KF method to reconstruct the original dynamic objects exactly.

The images reconstructed by the DIR algorithm under different time indexes are presented in Figure 14. As can be expected, when the dynamic reconstruction objects are considered, the DIR algorithm shows satisfactory numerical performance, the quality of the images reconstructed by the DIR algorithm is improved compared with the KF method, and the artifacts in the reconstructed images can be eliminated effectively. Meanwhile, it can be observed from Figures 13 and 14 that the quality of the images reconstructed by the DIR algorithm is higher than the KF method, which confirms that the DIR algorithm is successful in imaging the dynamic objects.

Table 11 lists the image errors for the KF method and the DIR algorithm. The second row in Table 11 is the image errors for the KF method. It can be observed that for the dynamic reconstruction objects, the quality of the images reconstructed by the KF method is far from perfect. The third row in Table 11 lists the image errors for the DIR algorithm. As can be expected, in the cases of the dynamic reconstruction objects, the quality of the images reconstructed by the DIR algorithm is higher than the KF method. Meanwhile, it can be observed from Table 11 that for the dynamic reconstruction cases simulated in this section, the DIR algorithm gives the smallest image errors, which indicates that the DIR algorithm is a promising candidate for solving the dynamic reconstruction problems.

Case 5

In this section, to evaluate further the feasibility and robustness of the DIR algorithm, another dynamic reconstruction case is simulated. In order to simulate real measurement environment, the noise-contaminated capacitance data is used. Figure 15 shows the dynamic process of the objects of interest under different time indexes, and subfigures a, b and c in Figure 15 represents the different locations of the dynamic reconstruction objects at discrete time indexes $k$ , $k + 1$ and $k + 2$ . The diameters of the cylinders in Figure 15 are 20 mm, the permittivity is defined as 2.6 and the permittivity of the rest of reconstruction region is 1.0. Table 12 lists the algorithmic parameters for the GIR algorithm, and the rest of the algorithmic parameters are the same as given above. Figures 16 and 17 illustrate the images reconstructed by the DIR algorithm under the noise levels of 5% and 13%. The image errors for the DIR algorithm under the different noise levels are presented in Table 13.

Table 12

Algorithmic parameters for the dynamic image reconstruction (DIR) algorithm

Algorithmic parameters	Figure 15(a)	Figure 15(b)	Figure 15(c)
$ϑ$	0.5	1.1	0.5
$α_{2}$	0.00036	0.0008	0.0009
$α$	0.000032	0.00004	0.00004
Number of iteration	8	8	8

Table 13

Image error (%)

Noise levels	Figure 15(a)	Figure 15(b)	Figure 15(c)
5%	7.63	Near zero	Near zero
13%	12.06	10.50	4.60

Figures 16 and 17 illustrate the images reconstructed by the GIR algorithm under the noise levels of 5% and 13%. As can be expected, when the measurement noises in the dynamic reconstruction objects are considered, the DIR algorithm shows satisfactory robustness and the quality of the images reconstructed by the DIR algorithm is acceptable. It can be seen from Figures 16 and 17 and Table 13 that the quality of the reconstructed images decreases with the increase of the noise levels, which indicates that the inaccurate properties on the measurement data should be treated seriously, and the issue should be further investigated in the future.

Conclusions

ECT is considered a promising visualization measurement technology, in which reconstructing high-quality images is highly desirable for real applications. Applications indicate that ECT measurement objects are often in a complicated dynamic process, and developing a dynamic reconstruction algorithm to image a dynamic object is important for real applications. In this paper, a dynamic reconstruction model, which integrates the ECT measurement information and the dynamic evolution information of the reconstruction objects, is presented. A generalized objective function, which simultaneously considers the ECT measurement information, the dynamic evolution information of the dynamic objects of interest, the temporal constraints and the spatial constraints, is proposed. An iteration scheme that integrates the advantages of the split Bregman iteration technique and the homotopy method is developed for solving the proposed objective function. Numerical simulations are implemented to evaluate the feasibility and effectiveness of the proposed algorithm. For the cases simulated in this paper, the accuracy of the images reconstructed by the proposed algorithm is improved and the artifacts in the reconstructed images can be removed effectively, which indicates that the proposed algorithm is successful in solving ECT inverse problems. As a result, a promising algorithm is introduced for solving ECT image reconstruction problems.

Applications indicate that each algorithm may illustrate different numerical performances to different reconstruction tasks, and the selection of a specific image reconstruction algorithm depends mainly on the measurement requirements and the prior information of a specific reconstruction task. Our work provides an alternative approach for solving ECT inverse problems, which needs to be further validated by more cases in the future.

Footnotes

Funding

This work was supported by the Fundamental Research Funds for the Central Universities (No. 10MG20); the China Postdoctoral Science Foundation (Nos. 20090460263 and 201003088); the National Natural Science Foundation of China (Nos. 50736002, 51006106 and 50906086); and the Program for the Changjiang Scholars and Innovative Research Team in University (No. IRT0952).

References

Banasiak

Soleimani

(2010) Shape based reconstruction of experimental data in 3D electrical capacitance tomography. NDT & E International 43: 241–249.

Beck

Ben-Tal

(2006) On the solution of the Tikhonov regularization of the total least squares problems. SIAM Journal on Optimization 17: 98–118.

Bennett

West

Luke

Williams

(2002) The investigation of bubble column and foam processes using electrical capacitance tomography. Minerals Engineering 15: 225–234.

Butala

Hewett

Frazin

Kamalabadi

(2010) Dynamic three-dimensional tomography of the solar corona. Solar Physics 262: 495–509.

Cai

Osher

Shen

(2009a) Split Bregman methods and frame based image restoration. Multiscale Modeling & Simulation 8: 337–369.

Cai

Osher

Shen

(2009b) Linearized Bregman iterations for compressed sensing. Mathematics of Computation 78: 1515–1536.

Warsito

Fan

(2006) Imaging the choking transition in gas–solid risers using electrical capacitance tomography. Industrial & Engineering Chemistry Research 45: 5384–5395.

Engl

(1987) Discrepancy principle for Tikhonov regularization of ill-posed problems leading to optimal convergence rates. Journal of Optimization Theory and Applications 52: 209–215.

Engl

Hanke

Neubauer

(1996) Regularization of Inverse Problems. Dordrecht: Kluwer Academic Publishing.

10.

Fang

(2004) A nonlinear image reconstruction algorithm for electrical capacitance tomography. Measurement Science and Technology 15: 2124–2132.

11.

Goldstein

Osher

(2009) The split Bregman method for L1-regularized problems SIAM Journal on Imaging Sciences 2: 323–343.

12.

Goldstein

Bresson

Osher

(2010) Geometric applications of the split Bregman method: segmentation and surface reconstruction. Journal of Scientific Computing 45: 272–293.

13.

Golub

Heath

Wahba

(1979) Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 21: 215–223.

14.

Guo

Wang

Cai

(2008) Engineering Optimization: Principle, Algorithm and Implementation. Beijing: China Machine Press.

15.

Hamidipour

Larachi

(2010) Dynamics of filtration in monolith reactors using electrical capacitance tomography. Chemical Engineering Science 65: 504–510.

16.

Hansen

(1992) Analysis of discrete ill-posed problems by means of the L-curve. SIAM Review 34: 561–580.

17.

Huang

Zeng

(2004) Theory and Methods for Nonlinear Numerical Analysis. Wuhan: Wuhan University Press.

18.

Jang

Lee

Kim

Choi

(2006) Modified iterative Landweber method in electrical capacitance tomography. Measurement Science and Technology 17: 1909–1917.

19.

Kaipio

Karjalainen

Somersalo

Vauhkonen

(1999) State estimation in time-varying electrical impedance tomography. Annals of the New York Academy Sciences 873: 430–439.

20.

Khambampati

Rashid

Ijaz

Kim

Soleimani

Kim

(2009) Unscented Kalman filter approach to tracking a moving interfacial boundary in sedimentation processes using three-dimensional electrical impedance tomography. Philosophical Transactions of the Royal Society 367: 3095–3120.

21.

Kim

Lee

Vauhkonen

(2001) Image reconstruction in time-varying electrical impedance tomography based on the extended Kalman filter. Measurement Science and Technology 12: 1032–1039.

22.

Kortschak

Wegleiter

Brandstatter

(2007) Formulation of cost functionals for different measurement principles in nonlinear capacitance tomography. Measurement Science and Technology 18: 71–78.

23.

Landweber

(1951) An iteration formula for Fredholm integral equations of the first kind. American Journal of Mathematics 73: 615–624.

24.

Lei

(2008) Research on image reconstruction algorithms of electrical capacitance tomography for the multiphase flow. PhD thesis, Institute of Engineering Thermophysics, Chinese Academy of Sciences.

25.

Lei

Liu

Guo

Han

(2011) An image reconstruction algorithm based on the semiparametric model for electrical capacitance tomography. Computers and Mathematics with Applications 61: 2843–2853.

26.

(1987) Numerical Methods for System of Nonlinear Equations. Beijing: Science Press.

27.

Liu

Chen

Xiong

Zhang

Lei

(2008) Preliminary study on ECT imaging of flames in porous media. Measurement Science and Technology 19: 1–7.

28.

Liu

Yang

Wang

Jiang

(2004) Prior-online iteration for image reconstruction with electrical capacitance tomography. IEE Proceedings – Science, Measurement and Technology 151: 195–200.

29.

(2005) Handbook of Modern Applied Mathematics. Beijing: Tsinghua University Press.

30.

Makkawi

Ocone

(2007) Integration of ECT measurement with hydrodynamic modelling of conventional gas–solid bubbling bed. Chemical Engineering Science 62: 4304–4315.

31.

Morozov

(1984) Methods for Solving Incorrectly Posed Problems. New York: Springer-Verlag.

32.

Mou

Peng

Yao

Xiao

(2005) Image reconstruction using a genetic algorithm for electrical capacitance tomography. Tsinghua Science and Technology 10: 587–592.

33.

Niedostatkiewicz

Tejchman

Chaniecki

Grudzien

(2009) Determination of bulk solid concentration changes during granular flow in a model silo with ECT sensors. Chemical Engineering Science 64: 20–30.

34.

Ortiz-Aleman

Martin

Gamio

(2004) Reconstruction of permittivity images from capacitance tomography data by using very fast simulated annealing. Measurement Science and Technology 15: 1382–1390.

35.

Osher

Mao

Dong

Yin

(2010) Fast linearized Bregman iteration for compressive sensing and sparse denoising. Communications in Mathematical Sciences 8: 93–111.

36.

Rimpilainen

Poutiainen

Heikkinen

Savolainen

Vauhkonen

Ketolainen

(2011) Electrical capacitance tomography as a monitoring tool for high-shear mixing and granulation. Chemical Engineering Science 66: 4090–4100.

37.

Schmitt

Louis

(2002) Efficient algorithms for the regularization of dynamic inverse problems: I. Theory Inverse Problems 18: 645–658.

38.

Seppänen

Vauhkonen

Somersalo

Kaipio

(2001) State estimation with fluid dynamical evolution models in process tomography – an application to impedance tomography. Inverse problems 17: 467–484.

39.

Soleimani

Lionheart

WRB

(2005) Nonlinear image reconstruction for electrical capacitance tomography using experimental data. Measurement Science and Technology 16: 1987–1996.

40.

Soleimani

Mitchell

Banasiak

Wajman

Adler

(2009) Four-dimensional electrical capacitance tomography imaging using experimental data. Progress in Electromagnetics Research-PIER 90: 171–186.

41.

Soleimani

Vauhkonen

Yang

Peyton

Kim

(2007) Dynamic imaging in electrical capacitance tomography and electromagnetic induction tomography using a Kalman filter. Measurement Science and Technology 18: 3287–3294.

42.

Takei

(2006) GVSPM image reconstruction for capacitance CT images of particles in a vertical pipe and comparison with the conventional method. Measurement Science and Technology 17: 2104–2112.

43.

Tikhonov

Arsenin

(1977) Solution of Ill-posed Problems. New York: V.H. Winston & Sons.

44.

Tossavainen

Vauhkonen

Kolehmainen

Kim

(2006) Tracking of moving interfaces in sedimentation processes using electrical impedance tomography. Chemical Engineering Science 61: 7717–7729.

45.

Vauhkonen

Mäkinen

Karjalainen

Kaipio

(2000) Dynamic electrical impedance tomography – phantom studies. Inverse Problem in Engineering 8: 495–510.

46.

Wang

Yang

(2010) Measurement of fluidised bed dryer by different frequency and different normalisation methods with electrical capacitance tomography. Powder Technology 199: 60–69.

47.

Wang

Tang

Cao

(2007) An image reconstruction algorithm based on total variation with adaptive mesh refinement for ECT. Flow Measurement and Instrumentation 18: 262–267.

48.

Wang

(2007) Computational Methods for Inverse Problems and Their Applications. Beijing: Higher Education Press.

49.

Warsito

Fan

(2001) Neural network based multi-criterion optimization image reconstruction technique for imaging two-and three-phase flow systems using electrical capacitance tomography. Measurement Science and Technology 12: 2198–2210.

50.

Watzenig

Brandner

Steiner

(2007) A particle filter approach for tomographic imaging based on different state–space representations. Measurement Science and Technology 18: 30–40.

51.

Xie

Huang

Hoyle

Thorn

Lenn

Snowden

Beck

(1992) Electrical capacitance for flow imaging: system model for development of image reconstruction algorithms and design of primary sensors. IEE Proceedings – G 139: 89–98.

52.

Xiong

Zhang

Liu

Lei

(2010) Wavelet enhanced visualization of solids distribution in the top of a CFB. Chemical Engineering Journal 158: 61–68.

53.

Yang

Peng

(2003) Image reconstruction algorithms for electrical capacitance tomography. Measurement Science and Technology 14: L1–L13.

54.

Yang

Spink

York

McCann

(1999) An image reconstruction algorithm based on Landweber’s iteration method for electrical capacitance tomography. Measurement Science and Technology 10: 1065–1069.

55.

Yin

Osher

Goldfarb

Darbon

(2008) Bregman iterative algorithms for ℓ₁ minimization with applications to compressed sensing. SIAM Journal on Imaging Sciences 1: 143–168.

56.

Zhang

Burger

Osher

(2011) A unified primal-dual algorithm framework based on Bregman iteration. Journal of Scientific Computing 46: 20–46.

57.

Zhao

Takei

Doh

(2010) ECT measurement and CFD-DEM simulation of particle distribution in a down-flow fluidized bed. Flow Measurement and Instrumentation 21: 212–218.