Dimensionality reduction-based dynamic reconstruction algorithm for electrical capacitance tomography

Abstract

The image reconstruction task in electrical capacitance tomography (ECT) is an ill-posed problem, in which image reconstruction methods play a vital role in real applications. In this paper, a multiple measurement vector-based dimensionality reduction dynamic reconstruction model that simultaneously utilizes the ECT measurement information and the dynamic evolution information of a dynamic object is proposed. A robust sparse orthogonal projective non-negative matrix factorization (RSOPNMF) method is proposed for extracting the basis vectors from a set of snapshots, and the split Bregman iteration (SBI) algorithm is used to solve the RSOPNMF model. The original unknown variables are projected onto the subspaces spanned by a set of the basis vectors extracted by the RSOPNMF method from a set of snapshots to obtain a low-dimensional model, where the images are indirectly reconstructed by solving the corresponding low-dimensional coefficient vector, and thus the dimensionality of the unknown variables is reduced and the computational cost are decreased. Based on the multiple measurement vectors and the dimensionality reduction model, an objective functional that incorporates the ECT measurement information, the dynamic evolution information of a dynamic object, the spatial constraint and the temporal constraint is proposed, in which the unknown variables are solved in a batching pattern. An iterative scheme that integrates the beneficial advantages of the SBI method and the forward–backward splitting algorithm is developed for solving the proposed objective functional. Numerical simulation results validate the feasibility of the proposed algorithm.

Keywords

Dimensionality reduction dynamic image reconstruction electrical capacitance tomography multiple measurement vectors projective non-negative matrix factorization method

Introduction

The electrical capacitance tomography (ECT) measurement method attempts to acquire the permittivity distribution information of the cross-section from the given capacitance measurement data via an appropriate algorithm, in which reconstruction algorithms play a crucial role in real applications. The ECT image reconstruction task consists of two major procedures: the forward problem and the inverse problem. The former is a well-posed problem with a purpose of computing the capacitance values from the given permittivity distributions, which can be easily solved by different numerical methods such as the finite element method and the finite difference method. According to the inverse problem theory (Engl et al., 2000; Han and Li, 2011; Huang and Zhao, 2005; Kirsch, 1996; Liu, 2005; Tarantola, 2005; Tikhonov and Arsenin, 1977; Vogel, 2002; Wang, 2007), the latter is an ill-posed problem with a motivation of estimating the permittivity distributions from the finite capacitance data.

In a mathematical notation, 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 electric charge; $ε (x, y)$ and $ϕ (x, y)$ represent the permittivity and electrical potential distributions, respectively; $V$ indicates the potential difference between two electrodes forming the capacitance and $Γ$ stands for the electrode surface. Under the considerations of the measurement noises and the linearization approximation of Equation (1), finally, the common static ECT reconstruction model can be approximated as (Yang and Peng, 2003):

Sx = C + r

(2)

where $S$ stands for a matrix of dimension $m \times n$ , and it is called as the sensitivity matrix in ECT image reconstruction; $x$ is an $n \times 1$ dimensional vector standing for the permittivity distribution, which indicates the grey level values in the visualization images; $C$ represents an $m \times 1$ dimensional vector indicating the capacitance values and $r$ is an $m \times 1$ dimensional vector standing for the measurement noises.

Reconstructing high-quality images is desired for real applications. In the past few years, various algorithms, including static reconstruction algorithms and dynamic reconstruction methods, have been developed for ECT image reconstruction. Common static reconstruction algorithms include the linear back-projection algorithm (Xie et al., 1992), the standard Tikhonov regularization (STR) method (Tikhonov and Arsenin, 1977), the Landweber iteration technique (Jang et al., 2006; Landweber, 1951; Yang et al., 1999), the neural network algorithm (Warsito and Fan, 2001), the truncated singular value decomposition method (Yang and Peng, 2003), the algebraic reconstruction technique (ART; Yang and Peng, 2003), the simultaneous iterative reconstruction technique (SIRT; Yang and Peng, 2003), the offline iteration and online reconstruction technique (Liu et al., 2004), the generalized Tikhonov regularization algorithms (Fang, 2004; Soleimani and Lionheart, 2005; Wang et al., 2007), the simulated annealing method (Ortiz-Aleman et al., 2004), the genetic algorithm (Mou et al., 2005), the generalized vector sampled pattern matching (GVSPM) method (Takei, 2006), the non-linear Landweber iteration algorithm (Li and Yang, 2008), the level set method (Banasiak and Soleimani, 2010; Kortschak et al., 2007), the Calderon method (Cao et al., 2011) and the image fusion-based method (Olmos et al., 2012). Popular dynamic reconstruction algorithms include the particle filter algorithm (Watzenig et al., 2007), the Kalman filter method (Soleimani et al., 2007), the four-dimensional imaging technique (Soleimani et al., 2009), the Ensemble Kalman filter method (Lei et al., 2012) and the split Bregman iteration (SBI)-based reconstruction algorithm (Lei et al., 2013); more details on the filter-based methods can be found in Zhu (2010). In the past few years, static reconstruction methods have been intensively studied. Owing to the fact that the considerations of the dynamic evolution information of a dynamic object and the temporal correlations are absent, however, static reconstruction methods may not be optimal for dynamic remonstration tasks. Due to the complexities and particularities of the dynamic reconstruction tasks, currently the investigations on dynamic reconstruction algorithms are far from perfect, and seeking an efficient reconstruction method that achieves fast and high-quality reconstruction remains a critical issue.

In order to address the above challenges, naturally, the main motivation of the paper is to improve the reconstruction speed and the reconstruction quality, and the main contributions are outlined as follows.

A dimensionality reduction-based multiple measurement vectors dynamic reconstruction model that utilizes the ECT measurement information and the dynamic evolution information of a dynamic object is proposed.

A robust sparse orthogonal projective non-negative matrix factorization (RSOPNMF) method is proposed for extracting the basis vectors from a set of snapshots, the original unknown variables are projected onto the subspaces spanned by the set of the basis vectors to obtain a low-dimensional model where the images are indirectly reconstructed by solving a low-dimensional RSOPNMF coefficient vector, and thus the dimensionality of the unknown variables is reduced, the underdetermined property is alleviated, the computational cost is decreased and the reconstruction speed is improved.

Based on the multiple measurement vectors and the dimensionality reduction model, a new objective functional that integrates the dynamic evolution information of the objects of interest, the ECT measurement information, the temporal constraint and the spatial constraint is proposed, in which the images are reconstructed by a batching pattern.

An iteration scheme that integrates the beneficial advantages of the SBI method and the forward–backward splitting (FBS) algorithm is developed for solving the proposed objective functional. Numerical simulations are implemented to validate the feasibility of the proposed algorithm. More importantly, this paper proposes a general framework for ECT image reconstruction problem, which may be applicable for other related tomography problems.

The rest of this paper is structured as follows. In the next section, a RSOPNMF method is proposed, then an iterative scheme is proposed for solving the RSOPNMF model, and finally a RSOPNMF-based dimensionality reduction model is proposed. A multiple measurement vector-based dimensionality reduction dynamic reconstruction model is presented, and an objective functional that fuses the dynamic evolution information of a dynamic object, the ECT measurement information, the temporal constraint and the spatial constraint is proposed. An iteration scheme, which integrates the merits of the SBI method and the FBS algorithm, is developed for solving the proposed objective functional. Numerical simulation results and discussions are illustrated and finally, conclusions are presented.

RSOPNMF-based dimensionality reduction method

In common reconstruction algorithms, including static and dynamic reconstruction algorithms, the original unknown variables are directly reconstructed, which is therefore not appropriate for imaging a dynamic object because of the time-consuming computation process. In this section, the non-negative matrix factorization (NMF) approach and the projective non-negative matrix factorization (PNMF) method are concisely introduced, a new RSOPNMF method is proposed for extracting the basis vectors from a set of snapshots, and the SBI method is developed for solving the RSOPNMF model. Finally, a procedure of extracting the basis vectors from the set of snapshots is outlined and a RSOPNMF-based dimensionality reduction method is presented.

Non-negative matrix factorization

The principle component analysis (PCA) method is an effective data processing approach. One of the important characteristics in the PCA method is that data representation is not purely additive, i.e. each principle component consists of both negative and positive entries. The linear combination in the PCA subspace mixes with both positive and negative weights, which are likely to cancel each other partially in the data representation. Finally, the partial cancellations lead to the loss of data in the PCA method (Han, 2010). Additionally, the negative components in the PCA method are hard to explain for specific applications, such as an ECT image.

The NMF method attempts to recover the hidden non-negative structures or patterns from the redundant data. Differing from common PCA method, the NMF method is a parts-based representation of an ensemble of signals or images. Far beyond a mathematical exploration, the philosophy underlying the NMF method, which tries to formulate a feasible model for learning object parts, is closely relevant to the perception mechanism. While the parts-based representation seems intuitive, it is indeed on the basis of physiological and psychological evidence: perception of the whole is based on the perception of its parts (Choi, 2008; Cichocki et al., 2009; Ding et al., 2006; Lee and Seung, 1999, 2001; Wang and Zhang, 2013). In fact, there are two complementary connotations in non-negativity (or non-negative component) and purely additive combination. On the one hand, the negative values of both observations and the latent components are physically meaningless in numerous real applications, such as an ECT image. On the other hand, objects of interest are most naturally characterized by the inventory of its parts, and the exclusively additive combination means that they can be reassembled by adding required parts together similar to identikits (Choi, 2008; Cichocki et al., 2009; Ding et al., 2006; Lee and Seung, 1999, 2001; Wang and Zhang, 2013).

Given an input data matrix $D$ , where its each column represents a candidate solution, the NMF method seeks the decomposition of a non-negative data matrix $D \geq 0$ , which is formulated by (Cichocki et al., 2009; Lee and Seung, 1999):

D \approx Φ H

(3)

where $Φ$ contains basis vectors in its columns; $H$ is considered a weighted matrix, and its each column describes how these basis vectors are combined fractionally to form the corresponding sample in $D$ . All entries in matrices $D$ , $Φ$ and $H$ are required to be non-negative.

In order to obtain the non-negative decompositions of a non-negative matrix, the NMF problem is cast as an optimization problem (Cichocki et al., 2009):

{\begin{matrix} min_{Φ, H} {‖ D - Φ H ‖_{F}^{2}} \\ s . t . Φ \geq 0, H \geq 0 \end{matrix}

(4)

RSOPNMF

The PNMF method is derived from the NMF technique, which is introduced to improve the locality of the part-based representations, and it can be formulated as (Yang et al., 2007; Yang and Oja, 2010):

min_{Φ \geq 0} {‖ D - Φ Φ^{T} D ‖_{F}^{2}}

(5)

where $Φ$ is a orthogonal matrix.

In real ECT applications, data matrix $D$ consists of a set of snapshots generated by measurement methods or numerical solutions or a combination of the two approaches, which is often corrupted by noises owing to the imperfect measurement or the deviations of the numerical computations. In common NMF models, previous studies often assume that the noises in the data matrix have bounded energy. Without such an assumption, unfortunately, it is hard to ensure a reliable estimation result. Obviously, satisfying such requirements is not a trivial task in real applications. Differing from common NMF models, following the ideas presented in Bao et al. (2012), Candès et al. (2011), Nguyen and Tran (2013), Wright et al. (2009) and Zhang et al. (2011a), the RPNMF method is introduced to deal with the outliers in data matrix, which can be expressed as:

D \approx Φ Φ^{T} D + R

(6)

where $R$ represents the noise, and will be solved iteratively.

In practice, imposing additional orthogonality and sparsity constraints is useful. In particular, an orthogonal matrix will form a basis of a subspace, and thus it will facilitate the geometric interpretation and the data reconstruction. Finally, a RSOPNMF method with the orthogonality and sparsity constraints on $Φ$ is proposed, and in a mathematical notation it can be formulated as the following optimization problem:

{\begin{matrix} min_{Φ \geq 0} {\frac{1}{2} ‖ D - Φ Φ^{T} D - R ‖_{F}^{2} + λ_{1} {‖ Φ ‖}_{1} + λ_{2} {‖ R ‖}_{1}} \\ s . t . Φ^{T} Φ = I \end{matrix}

(7)

where $‖ \cdot ‖_{1}$ represents the 1-norm for a matrix; $λ_{1} > 0$ and $λ_{2} > 0$ are called as the regularization parameters. It is worth mentioning that Yang et al. (2007) observed that although the optimization of the PNMF does not involve an explicit orthogonalization step, the resulting matrix $Φ$ has high sparsity, locality and orthogonality. Following the ideas presented in Zhang (2009) and Cichocki et al. (2009), however, in this paper the additive constraints are explicitly imposed to ensure the orthogonality and the sparsity of matrix $Φ$ , which is highly appropriate for real applications.

Proposed iterative scheme

Equation (7) involves two unknown matrix variables, $Φ$ and $R$ . For the sake of easy computation, following the computational strategies presented in Blondel et al. (2013), Chan and Wong (2000), Lv et al. (2013), Qin et al. (2013), Tseng (2001), Xiao et al. (2011) and Xu and Yin (2013), it is numerically appealing to decouple Equation (7) as follows:

min_{R} {\frac{1}{2} ‖ D - Φ Φ^{T} D - R ‖_{F}^{2} + λ_{2} {‖ R ‖}_{1}}

(8)

{\begin{matrix} min_{Φ \geq 0} {\frac{1}{2} ‖ D - Φ Φ^{T} D - R ‖_{F}^{2} + λ_{1} {‖ Φ ‖}_{1}} \\ s . t . Φ^{T} Φ = I \end{matrix}

(9)

Finally, it can be observed from Equations (8) and (9) that the algorithmic structure can be expressed as an iterative series ${R^{τ}, Φ^{τ}, R^{τ + 1}, Φ^{τ + 1}, R^{τ + 2}, Φ^{τ + 2}, \dots}$ . Equation (8) can be easily solved by the shrinkage algorithm, and Equation (9) can be rewritten as a constrained optimization problem by introducing an additional variable $E$ and an equality constraint $E = Φ$ :

{\begin{matrix} min_{Φ \geq 0, E \geq 0} {\frac{1}{2} ‖ D - Φ Φ^{T} D - R ‖_{F}^{2} + λ_{1} {‖ E ‖}_{1}} \\ s . t . Φ^{T} Φ = I, E = Φ \end{matrix}

(10)

According to the SBI method (Cai et al., 2009a, 2009b; Goldstein and Osher, 2009; Goldstein et al., 2010; Osher et al., 2010; Yin et al., 2008; Zhang et al., 2011b), Equation (10) is iteratively solved by:

(Φ^{τ + 1}, E^{τ + 1}) = {\begin{matrix} min_{Φ \geq 0, E \geq 0} {\frac{1}{2} ‖ D - Φ Φ^{T} D - R ‖_{F}^{2} + λ_{1} {‖ E ‖}_{1} + \frac{μ_{1}}{2} ‖ E - Φ - θ^{τ} ‖_{F}^{2}} \\ s . t . Φ^{T} Φ = I \end{matrix}

(11)

θ^{τ + 1} = θ^{τ} + (Φ^{τ + 1} - E^{τ + 1})

(12)

For the sake of easy computation, similarly, Equation (11) can be decoupled as the following two sub-problems (Blondel et al., 2013; Chan and Wong, 2000; Lv et al., 2013; Qin et al., 2013; Tseng, 2001; Xiao et al., 2011; Xu and Yin, 2013):

E^{τ + 1} = min_{E \geq 0} {λ_{1} {‖ E ‖}_{1} + \frac{μ_{1}}{2} ‖ E - Φ^{τ} - θ^{τ} ‖_{F}^{2}}

(13)

{\begin{matrix} Φ^{τ + 1} = min_{Φ \geq 0} {J = \frac{1}{2} ‖ D - Φ Φ^{T} D - R ‖_{F}^{2} + \frac{μ_{1}}{2} ‖ E^{τ + 1} - Φ - θ^{τ} ‖_{F}^{2}} \\ s . t . Φ^{T} Φ = I \end{matrix}

(14)

Commonly, Equation (14) can be solved using the standard Lagrangian multipliers method by introducing the orthogonality constraint to the objective functional. For the sake of easy computation, however, applying the natural gradient on the Stiefel manifold to impose the orthogonality constraint on matrix $Φ$ yields the following expression (Choi, 2008; Cichocki et al., 2009; Ding et al., 2006; Edelman et al., 1999; Nishimori and Akaho, 2005; Yoo and Choi, 2010):

{\bar{\nabla}}_{Φ} J = \nabla_{Φ} J - Φ {[\nabla_{Φ} J]}^{T} Φ

(15)

where $\nabla (\cdot)$ is a gradient operator.

Based on the multiplicative update rule (Choi, 2008; Cichocki et al., 2009; Ding et al., 2006; Edelman et al., 1999; Lee and Seung, 2001; Nishimori and Akaho, 2005; Yoo and Choi, 2010) and Equation (14), finally, an iteration scheme of solving matrix variable $Φ$ can be formulated as:

Φ_{ij}^{τ + 1} = Φ_{ij}^{τ} \frac{{(D D^{T} Φ + {(Φ Φ^{T} R D^{T} Φ)}_{+} + {(Φ Φ^{T} D R^{T} Φ)}_{+} + μ_{1} Φ {({(ϑ - E)}^{T})}_{+} Φ)}_{ij}^{τ}}{{({(R D^{T} Φ)}_{+} + {(D R^{T} Φ)}_{+} + μ_{1} {(ϑ - E)}_{+} + Φ Φ^{T} D D^{T} Φ)}_{ij}^{τ}}

(16)

where ${(x_{i, j})}_{+} = max (ε, x_{i, j})$ and $ε = 10^{- 8}$ .

According to the above discussions, an iterative scheme is proposed for solving Equation (7), which can be summarized as follows:

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

Step 2. A set of snapshots is generated by measurement methods or numerical solutions or a combination of the two approaches.

Step 3. Update $R$ by solving Equation (8) using the shrinkage algorithm.

Step 4. Update $E$ by solving Equation (13) using the shrinkage algorithm.

Step 5. Update $Φ$ using Equation (16).

Step 6. Update $θ$ using Equation (12).

Step 7. Loop to Step 3 until a predetermined stopping criterion is satisfied.

RSOPNMF-based dimensionality reduction method

Following the above discussions, if a basis matrix $Φ$ is provided, the original data matrix $X$ is approximated by:

X \approx Φ A

(17)

where $A$ is a low-dimensional matrix vector. Unlike other bases in the finite element method and the spectral method, the RSOPNMF basis is directly extracted from the physical properties of the systems of interest, and thus the RSOPNMF basis contains main characteristics of the expected solutions.

According to the discussions presented in previous sections, finally, a RSOPNMF-based dimensionality reduction approach can be outlined as follows:

Step 1. A set of snapshots is generated.

Step 2. Extract the basis vectors from a set of snapshots using the proposed iterative scheme presented above.

Step 3. Use the basis vectors to represent an unknown matrix, and thus a low-dimensional model is obtained.

Reconstruction model

The common reconstruction model, Equation (2), considers the inaccuracy of the capacitance measurement data. Studies indicate that the model approximation deviations derived from the approximation and discretization of a real problem and physically implementing an ECT sensor will undoubtedly introduce errors, and thus developing a model that simultaneously considers the inaccurate properties on the capacitance measurement data, the model approximation and the sensitivity matrix may be essential for improving the reconstruction quality. Differing from vector-based reconstruction models, in this paper a multiple measurement vector-based reconstruction model considering the above inaccurate properties can be formulated as:

(S + E) X + U = Y + N

(18)

where $U$ is an $m \times t$ dimensional matrix standing for the linearization approximation errors; $E$ defines an $m \times n$ dimensional matrix representing the inaccuracy of the sensitivity matrix derived from physically implementing an ECT sensor; $Y$ is an $m \times t$ dimensional matrix indicating the capacitance measurement data; $X = [x_{1}, x_{2}, \dots, x_{t}]$ represents an $n \times t$ dimensional matrix, and its each column stands for a specific permittivity distribution at a specific measurement instant; $N$ is an $m \times t$ dimensional matrix standing for the capacitance measurement noises; $t > 1$ defines the number of measurements or the measurement time window. It is worth emphasizing that in Equation (2) the unknown variable is a vector; on the contrary, the unknown variables in Equation (18) are three matrices, which is a distinct feature of differing from other vector-based reconstruction models.

Commonly, the dynamic evolution process of a dynamic object in real applications can be generally depicted by the following formula:

f (x_{k}, v_{k}) = 0

(19)

where $f$ describes the dynamic evolution process of the objects of interest, and it can be formulated as different expressions according to the different dynamic objects. In particular, $f$ can be specified as a set of the partial differential equations in the gas-solid two-phase flow measurements (Li et al., 2005); $v_{k}$ defines the uncertainties of the dynamic evolution process and $k$ stands for the discrete time index. In real applications, the dynamic evolution information can be in advance solved by means of numerical methods such as the finite element method and the finite difference method.

ECT imaging objects are often in a fast dynamic process. To visualize the detailed characteristics of the dynamic process, achieving fast reconstruction is essential, and thus introducing the linearization approximation to Equation (19) yields:

x_{k + 1} = F_{k} x_{k} + v_{k}

(20)

where $F_{k}$ represents a state-transition operator at time instant $k$ ; If set $F_{k} = I$ , where $I$ defines an identity matrix, Equation (20) is a purely random-walk evolution model. The model is often used in practices when a better temporal dynamic model is not known, and although it is simple, a dynamic reconstruction with a random-walk model may be superior to a static reconstruction model (Butala et al., 2010). Ideally, the non-linear dynamic model can be incorporated into the state transition operator $F_{k}$ . Like Butala et al. (2010) and Soleimani et al. (2007), for the sake of easy computation, in this paper the random-walk evolution model is used.

Finally, a dynamic reconstruction model that integrates the ECT measurement information and the dynamic evolution information of a dynamic object is formulated as:

{\begin{matrix} f (x_{k}, v_{k}) = 0 \\ (S + E) X + U = Y + N \end{matrix}

(21)

In common static and dynamic reconstruction methods, the number of the unknown variables is more than that of the independent equations, which restricts the improvement of the reconstruction speed. Differing from common reconstruction methods, in this paper the original unknown variables are projected onto the spaces spanned by the basis vectors extracted by the RSOPNMF method to reduce the dimensionality of the unknown variable, and a dimensionality reduction dynamic image reconstruction model can be formulated as:

{\begin{matrix} f (Φ A_{k}, v_{k}) = 0 \\ (S + E) Φ A + U = Y + N \end{matrix}

(22)

Objective functional

In order to ensure a stable numerical solution, the Tikhonov regularization method is employed to solve Equation (22), and thus the equation can be cast as the following optimization problem:

min {\begin{matrix} \frac{1}{2} ‖ (S + E) Φ A + U - Y ‖_{F}^{2} + \frac{α}{2} ‖ Φ A - Φ A_{c} ‖_{F}^{2} + \\ α_{1} Ω_{1} (A) + α_{2} Ω_{2} (U) + \frac{α_{3}}{2} ‖ Φ AW ‖_{F}^{2} + α_{4} Ω_{3} (E) \end{matrix}}

(23)

where $α > 0$ , $α_{1} > 0$ , $α_{2} > 0$ , $α_{3} > 0$ and $α_{4} > 0$ are the regularization parameters; term $‖ (S + E) Φ A + U - Y ‖_{F}^{2}$ is used to measure the data fidelity; term $Φ A_{c}$ can be solved by Equations (19) or (20); term $‖ Φ A - Φ A_{c} ‖_{F}^{2}$ is employed to fuse the dynamic evolution information of a dynamic object; $Ω_{1}$ , $Ω_{2}$ and $Ω_{3}$ represent the spatial constraints; $‖ Φ AW ‖_{F}^{2}$ defines a temporal constraint, which is introduced to impose the temporal correlations of a dynamic object; the weighted matrix $W$ is defined as:

W = {[\begin{matrix} - 1 & 0 & \dots & 0 \\ 1 & - 1 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ ⋮ & ⋮ & \dots & ⋮ \\ 0 & 0 & \dots & - 1 \\ 0 & 0 & \dots & 1 \end{matrix}]}_{t \times t - 1}

(24)

To improve the robustness of the least squares estimation, in this paper the absolute value function is used to replace the squared residuals (Cichocki and Amari, 2002; Huber 1981), and thus Equation (23) is reformulated by:

min {\begin{matrix} ‖ (S + E) Φ A + U - Y ‖_{1} + \frac{α}{2} ‖ Φ A - Φ A_{c} ‖_{F}^{2} + \\ α_{1} Ω_{1} (A) + α_{2} Ω_{2} (U) + \frac{α_{3}}{2} ‖ Φ AW ‖_{F}^{2} + α_{4} Ω_{3} (E) \end{matrix}}

(25)

where $‖ \cdot ‖_{1}$ represents the 1-norm for a matrix.

For the sake of easy computation, $Ω_{1}$ , $Ω_{2}$ and $Ω_{3}$ are defined as:

Ω_{1} = ‖ A ‖_{1}

(26)

Ω_{2} = ‖ U ‖_{1}

(27)

Ω_{3} = ‖ E ‖_{1}

(28)

Finally, an objective functional for ECT reconstruction is specified as:

min {\begin{matrix} ‖ (S + E) Φ A + U - Y ‖_{1} + \frac{α}{2} ‖ Φ A - Φ A_{c} ‖_{F}^{2} \\ + α_{1} ‖ A ‖_{1} + α_{2} ‖ U ‖_{1} + \frac{α_{3}}{2} ‖ Φ AW ‖_{F}^{2} + α_{4} ‖ E ‖_{1} \end{matrix}}

(29)

Equation (29) belongs to the framework of the Tikhonov regularization method, which has the following desirable properties:

Differing from common reconstruction algorithms, in which the unknown variables are directly reconstructed, the dimensionality of the unknown variables in Equation (29) is reduced by the RSOPNMF method, which will facilitate the improvement of the reconstruction speed.

The prior information from previous measurement or numerical simulations can be integrated to Equation (29) in solving the basis vectors, and thus the quantity of information is indirectly increased.

The dynamic evolution information of a dynamic object and the ECT measurement information are exploited in Equation (29), and thus the quantity of information in ECT image reconstruction is increased, and it will facilitate the improvement of the reconstruction quality.

Owing to the fact that it is hard to describe exactly the dynamic evolution information of the objects and to achieve exact capacitance measurements, the inaccurate properties of the dynamic evolution information and the measurement information are emphasized in Equation (29). In the measurement equation, particularly, the measurement noises, the model approximation deviations and the inaccuracy of the sensitivity matrix are simultaneously considered, which distinctly differs from common dynamic and static reconstruction algorithms.

The spatial constraint and the temporal constraint are simultaneously emphasized in Equation (29). In the case of the multiple measurement vectors, especially, the temporal constraint is introduced to impose the temporal correlations of a dynamic object.

The unknown variables in Equation (29) are three matrices, and it differs from other vector-based reconstruction algorithms, such as the ART method, the GVSPM algorithm, etc. In particular, the images are reconstructed by a batching pattern in Equation (29), which will facilitate the exploitation of the temporal correlations of a dynamic object and the improvement of the reconstruction speed.

Owing to the ill-posed nature of the inverse problem, ECT inversion solution is unstable, that is, a small fluctuation in the input data may bring a fatal effect on final results, which may make the inversion solution meaningless. As a result, methods that ensure a stable numerical solution while improving the reconstruction quality should be used. In Equation (29), the Tikhonov regularization technique is introduced to stabilize an inversion solution.

In real applications, measurement noises are ubiquitous and complicated. In Equation (29), the robust estimation is employed to improve the robustness of estimation.

Solving of the objective functional

Seeking an efficient algorithm to solve Equation (29) is crucial. In this section, the FBS algorithm is introduced, and an iteration scheme that integrates the merits of the SBI method and the FBS algorithm is proposed for solving Equation (29).

Forward–backward splitting method

Owing to the low computational complexity, the FBS algorithm is used to solve the proposed objective functional. In a mathematical notation, the FBS method aims at solving the following optimization problem (Combettes and Wajs, 2005; Duchi and Singer, 2009; Montefusco et al., 2011; Zhang et al., 2010):

min_{u} {μ J (u) + H (u)}

(30)

where $J$ and $H$ are the known functions.

After the corresponding deductions, the final FBS algorithm is formulated as:

u^{τ + 1} = Pro x_{δ μ J} (u^{τ} - δ \partial H (u^{τ}))

(31)

where the proximal operator $Pro x_{δ μ J} (v)$ is defined as:

\underset{u}{\arg min} {μ J (u) + \frac{1}{2 δ} {‖ u - v ‖}^{2}}

(32)

In the case of $H (u) = \frac{1}{2} ‖ ψ u - f ‖^{2}$ , the solution of Equation (30) is computed by the following two-step algorithm:

v^{τ + 1} = u^{τ} - δ ψ^{T} (ψ u^{τ} - f)

(33)

u^{τ + 1} = \underset{u}{\arg min} {μ J (u) + \frac{1}{2 δ} {‖ u - v^{τ + 1} ‖}^{2}}

(34)

Proposed iteration scheme

According to the optimization theory, Equation (29) is reformulated as a constrained optimization problem by introducing an equality constraint:

{\begin{matrix} min {{‖ d ‖}_{1} + \frac{α}{2} ‖ Φ A - Φ A_{c} ‖_{F}^{2} + α_{1} {‖ A ‖}_{1} + α_{2} {‖ U ‖}_{1} + \frac{α_{3}}{2} ‖ W Φ A ‖_{F}^{2} + α_{4} {‖ E ‖}_{1}} \\ s . t . d = (S + E) Φ A + U - Y \end{matrix}

(35)

Applying the SBI algorithm to Equation (35) yields:

(d^{τ + 1}, A^{τ + 1}, U^{τ + 1}, E^{τ + 1}) = min {\begin{matrix} ‖ d ‖_{1} + \frac{α}{2} ‖ Φ A - Φ A_{c} ‖_{F}^{2} + α_{1} ‖ A ‖_{1} + α_{2} ‖ U ‖_{1} + \frac{α_{3}}{2} ‖ Φ AW ‖_{F}^{2} \\ + \frac{μ}{2} ‖ d - ((S + E) Φ A + U - Y) - b^{τ} ‖_{F}^{2} + α_{4} ‖ E ‖_{1} \end{matrix}}

(36)

b^{τ + 1} = b^{τ} - (d^{τ + 1} - ((S + E^{τ + 1}) Φ A^{τ + 1} + U^{τ + 1} - Y))

(37)

There are four unknown variables in Equation (36), and efficiently solving the equation is crucial. For the sake of easy computation, following the computational strategies presented in Blondel et al. (2013), Chan and Wong (2000), Lv et al. (2013), Qin et al. (2013), Tseng (2001), Xiao et al. (2011) and Xu and Yin (2013), it is numerically appealing to decouple Equation (36) as:

U^{τ + 1} = min {α_{2} {‖ U ‖}_{1} + \frac{μ}{2} ‖ d^{τ} - ((S + E^{τ}) Φ A^{τ} + U - Y) - b^{τ} ‖_{F}^{2}}

(38)

d^{τ + 1} = min {{‖ d ‖}_{1} + \frac{μ}{2} ‖ d - ((S + E^{τ}) Φ A^{τ} + U^{τ + 1} - Y) - b^{τ} ‖_{F}^{2}}

(39)

A^{τ + 1} = min {\begin{matrix} \frac{α}{2} ‖ Φ A - Φ A_{c} ‖_{F}^{2} + α_{1} ‖ A ‖_{1} + \frac{α_{3}}{2} ‖ Φ AW ‖_{F}^{2} \\ + \frac{μ}{2} ‖ d^{τ + 1} - ((S + E^{τ}) Φ A + U^{τ + 1} - Y) - b^{τ} ‖_{F}^{2} \end{matrix}}

(40)

E^{τ + 1} = min {α_{4} {‖ E ‖}_{1} + \frac{μ}{2} ‖ d^{τ + 1} - ((S + E) Φ A^{τ + 1} + U^{τ + 1} - Y) - b^{τ} ‖_{F}^{2}}

(41)

Following the discussions presented in previous sections, an iteration scheme is proposed for solving Equation (29), which is outlined as follows:

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

Step 2. Obtain $Φ$ by solving Equation (7).

Step 3. Update variable $U$ by solving Equation (38) using the shrinkage algorithm.

Step 4. Update variable $d$ by solving Equation (39) using the shrinkage algorithm.

Step 5. Update variable $A$ by solving Equation (40) using the FBS algorithm.

Step 6. Update variable $E$ by solving Equation (41) using the FBS algorithm.

Step 7. Update variable $b$ according to Equation (37).

Step 8. Loop to Step 3 until a predetermined iteration stopping criterion is satisfied.

Step 9. Solve Equation (17) to obtain the estimation of the original unknown variable.

Numerical simulations and discussions

For simplicity, the proposed method is called as the RSOPNMF-based dynamic reconstruction (RSOPNMFDR) algorithm. In this section, the dynamic reconstruction cases with multiple measurement vectors are implemented to numerically evaluate the feasibility and effectiveness of the RSOPNMFDR algorithm, and the reconstruction quality and speed are compared with the projected Landweber iteration (PLI) method. The robustness of the RSOPNMFDR algorithm is evaluated using the noise-contaminated capacitance data. Additionally, the effect of the number of the basis vectors on the reconstruction quality is numerically exploited.

A 12-electrode square ECT sensor is selected for simulations and an image is visualized using 32×32 pixels. All algorithms are implemented using the MATLAB 7.0 software on a PC with a Pentium IV 2.4-GHz CPU and 4 Gbytes memory. In this paper, the image error is used to evaluate the quality of an inversion solution, which can be specified as:

η = \frac{‖ x_{Original} - x_{Computed} ‖}{‖ x_{Original} ‖} \times 100 %

(42)

where $η$ stands for the image error; $x_{Original}$ and $x_{Computed}$ represent the original and computed permittivity distributions, respectively.

Case 1

In this section, numerical simulations are implemented to evaluate the feasibility and effectiveness of the RSOPNMF method, and the initial value is solved by the singular values decomposition method. Figure 1 illustrates the basis vectors extracted by the RSOPNMF method from a set of snapshots and Figure 2 shows the orthogonality of the basis matrix.

Four basis vectors extracted by the RSOPNMF method are shown in Figure 1. It can be found from Figure 1 that differing from the common PCA method, the RSOPNMF method can ensure the sparseness and non-negativity of a basis vector, which indicates that the RSOPNMF algorithm is successful in exacting the basis vectors. It is worth emphasizing that the basis vectors can be in advance solved, which will not influence the reconstruction speed.

In practice, the orthogonality of the basis matrix is highly useful. In particular, an orthogonal matrix will form a basis of a subspace, and thus facilitates the geometric interpretation and the data reconstruction. Figure 2 illustrates the orthogonality of the basis matrix, that is, the image of the matrix $Φ^{T} Φ$ . In theory, the matrix $Φ$ is orthogonal if the $Φ^{T} Φ = I$ . It is found that the basis matrix solved by the proposed algorithm has a good orthogonality. It is worth emphasizing that for the sakes of easy computation and few algorithmic parameters, in this paper the orthogonalilty constraint is imposed by introducing the natural gradient on the Stiefel manifold and numerical results shows the feasibility of the method. In theory, the orthogonalilty constraint can be introduced by the standard Lagrangian multipliers method. However, the method proposed in this paper requires few algorithmic parameters, which is appropriate for real applications.

Case 2

ECT imaging objects are often in a dynamic process, a successful algorithm should be able to tackle such reconstruction tasks. In this section, a dynamic reconstruction case is simulated to evaluate the feasibility of the RSOPNMFDR algorithm. Subfigures (a), (b), (c) and (d) in Figure 3 stand for the different permittivity distributions at time instants $k$ , $k + 1$ , $k + 2$ and $k + 3$ , respectively. In Figure 3, the diameters of the cylinders are 20 mm, the permittivity of the black part is defined as 2.6, and the permittivity of the rest of reconstruction region is 1.0. In essence, the RSOPNMFDR algorithm belongs to the framework of the Tikhonov regularization method, in which the selection of the regularization parameter plays a crucial role in real applications. Currently, various methods, including the discrepancy principle method (Morozov, 1984), the generalized cross-validation method (Golub et al., 1979), the L-curve method (Hansen, 1992) and the Engl criteria method (Engl et al., 2000), are available for determining the regularization parameter, and more discussions on the selection of the regularization parameter can be traced back to Huang and Zhao (2005) and Wang (2007). Owing to the complexity and particularity of the problem, determining an optimal regularization parameter remains an open problem presently, and in the field of ECT image reconstruction the regularization parameter is often chosen empirically (Yang and Peng, 2003). For the sake of easy computation, in this paper the regularization parameters are empirically determined, which are specified as $α = 0.0004$ , $α_{1} = 0.004$ , $α_{2} = 10$ , $α_{3} = 0.0001$ , $α_{4} = 1$ , $μ = 1$ , and the number of iterations is 10. Initial solutions for the PLI method and the RSOPNMFDR algorithm are computed by the STR method. The images reconstructed by the PLI method and the RSOPNMFDR algorithm are illustrated in Figures 4 and 5, respectively. Table 1 lists the algorithmic parameters for the PLI algorithm. The image errors and the CPU times are listed in Tables 2 and 3, respectively.

Table 1.

Algorithmic parameters for the projected Landweber iteration (PLI) algorithm.

Algorithmic parameters	Figure 3(a)	Figure 3(b)	Figure 3(c)	Figure 3(d)
Relaxation factor	1	1	1	1
Number of iteration	596	500	530	411

Table 2.

Image errors (%).

Algorithms	Figure 3(a)	Figure 3(b)	Figure 3(c)	Figure 3(d)
PLI	12.62	15.09	18.84	14.61
RSOPNMFDR	0.71	1.49	3.33	0.29

PLI, projected Landweber iteration; RSOPNMFDR, robust sparse orthogonal projective non-negative matrix factorization-based dynamic reconstruction.

Table 3.

CPU times (seconds).

Figures	PLI	RSOPNMFDR
Figure 3(a)	0.88	0.21
Figure 3(b)	0.69
Figure 3(c)	0.71
Figure 3(d)	0.56

PLI, projected Landweber iteration; RSOPNMFDR, robust sparse orthogonal projective non-negative matrix factorization-based dynamic reconstruction.

Figure 4 shows the images reconstructed by the PLI algorithm. Owing to the fact that only gradient information of the objective functional is used in implementing image reconstruction, the numerical implementation of the PLI algorithm is easy and the computational complexity is low. In the past years, the PLI algorithm had found numerous applications in different fields of the inverse problems, including ECT image reconstruction. Like other static reconstruction algorithms, unfortunately, the PLI algorithm merely utilizes the ECT measurement information, and fails to consider the dynamic evolution information and the temporal correlations of a dynamic object, which is not optimal for the dynamic reconstruction tasks. In fact, it can be found from Figure 4 that in the case of the dynamic reconstruction, the quality of the images reconstructed by the PLI method is far from perfect, the distortions are serious, and the detailed edge information of the reconstructed images is blurred. It can be observed form Table 2 that the image errors of the PLI algorithm for dynamic reconstruction objects, Figures 3(a)–(d), are 12.62%, 15.09%, 18.84% and 14.61%, which are higher than that of the RSOPNMFDR algorithm.

Figure 5 is the images reconstructed by the RSOPNMFDR algorithm. Numerical simulation results indicate that the RSOPNMFDR algorithm can ensure the numerical stability of an inversion solution owing to the fact that the Tikhonov regularization technique is employed, and the computational cost is low due to the dimensionality reduction of the original unknown variables. In particular, the RSOPNMFDR algorithm simultaneously utilizes the ECT measurement information and the dynamic evolution information of a dynamic object, and thus the quantity of information is increased, which facilitates the improvement of the reconstruction quality. It can be seen from Figure 5 that the quality of the images reconstructed by the RSOPNMFDR algorithm is improved compared with the PLI method. Furthermore, it can be seen from Table 2 that the image errors of the RSOPNMFDR algorithm for dynamic reconstruction objects, Figures 3(a)–(d), are 0.71%, 1.49%, 3.33% and 0.29%, which are lower than that of the PLI algorithm.

The image errors for the PLI method and the RSOPNMFDR algorithm are listed in Table 2. As can be expected, the reconstruction quality of the RSOPNMFDR algorithm is superior to the PLI method. Furthermore, it can be observed from Table 2 that for the cases simulated in this section, the RSOPNMFDR algorithm gives the smallest image errors, which indicates that the RSOPNMFDR algorithm is a promising candidate for solving ECT image reconstruction problems.

ECT measurement objects are often in a time-vary process, and achieving fast reconstruction is crucial for real applications. Based on the multiple measurement vectors model and the dimensionality reduction of the unknown variables, in the RSOPNMFDR algorithm the images are indirectly reconstructed by solving the corresponding low-dimensional coefficient vector in a batching pattern, which differs from existing reconstruction algorithms. Owing to the dimensionality reduction of the original unknown variables, especially, the reconstruction speed is improved and the CPU times are reduced. In fact, it can be found from Table 3 that the CPU times of the RSOPNMFDR algorithm for reconstructing Figures 3(a)–(d) are about 0.21 s, which is significantly fewer than that of the PLI method, 2.84 s. This advantage is highly appropriate for real applications.

Case 3

The ECT image reconstruction process is an ill-posed problem, and small perturbations in input data may lead to a large fluctuation of an inversion solution, which may make the final solution meaningless. In practice, measurement noises are ubiquitous and complicated. A successful algorithm should be able to deal with such problem. In this section, the noise-contaminated capacitance data is used to evaluate the robustness of the RSOPNMFDR algorithm, and the definition of the noise level can be specified as:

ϖ = \frac{‖ Y_{Contaminated} - Y_{Original} ‖}{‖ Y_{Original} ‖} \times 100 %

(43)

where $ϖ$ is the noise level; $Y_{Original}$ stands for the true capacitance data; $Y_{Contaminated}$ is the noise-contaminated capacitance data; $Y_{Contaminated} = Y_{Original} + r$ , where $r = σ_{1} \cdot randn$ ; $σ_{1}$ represents the standard deviation and $randn$ stands for a Gaussian distribution random number with the mean of 0 and the standard deviation of 1, which can be achieved by the function ‘randn’ in the MATLAB software.

Algorithmic parameters for the RSOPNMFDR algorithm are the same as above. Figures 6 –8 illustrate the images reconstructed by the RSOPNMFDR algorithm under the noise levels of 5%, 16% and 30%, respectively. The image errors are listed in Table 4.

Table 4.

Image errors (%).

Noise levels	Figure 3(a)	Figure 3(b)	Figure 3(c)	Figure 3(d)
5%	0.78	1.97	3.60	0.78
16%	2.91	3.39	5.32	1.83
30%	3.17	5.38	8.51	4.72

Figures 6 –8 show the images reconstructed by the RSOPNMFDR algorithm under the noise levels of 5%, 16% and 30%, respectively. Owing to the fact that the M-estimation technique and the Tikhonov regularization method are introduced to the proposed objective functional, as can be expected, the RSOPNMFDR algorithm shows satisfactory robustness, and the reconstructed images under different noise levels are in a good agreement with the original dynamic objects. In particular, it can be observed from Table 4 that when the noise level is 30%, the image errors for the original reconstruction objects, Figures 3(a)–(d), are 3.17%, 5.38%, 8.51% and 4.72%, respectively. Additionally, owing to the fact that the dimensionality of the unknown variables in the RSOPNMFDR algorithm, the ill-posed property may be therefore alleviated, finally the sensitivity of the inversion solution to measurement noises is decreased.

Case 4

It is worth emphasizing that in real applications the sensitivity matrix is often inaccurate derived from physically implementing an imperfect ECT sensor and the approximation solution of the sensitivity matrix, and thus a successful algorithm should be able to deal with such inaccuracy. In this section, the influence of the inaccuracy of the sensitivity matrix on the reconstruction quality is evaluated, and the inaccuracy of the sensitivity matrix is defined as:

S_{contaminated} = S_{original} + E

(44)

where $E = σ_{2} \cdot randn$ , and $σ_{2}$ is the standard deviation; $randn$ represents a normal distribution random number with the mean of 0 and the standard deviation of 1, which can be achieved by the function ‘randn’ in the MATLAB software; $S_{original}$ and $S_{contaminated}$ stand for the original and noise-contaminated sensitivity matrices, respectively.

Algorithmic parameters for the RSOPNMFDR algorithm are the same as before. Figures 9 and 10 illustrate the images reconstructed by the RSOPNMFDR algorithm when the standard deviations of the sensitivity matrix are 0.002 and 0.004, respectively. The image errors are presented in Table 5.

Table 5.

Image errors (%).

Standard deviations	Figure 3(a)	Figure 3(b)	Figure 3(c)	Figure 3(d)
$σ_{2} = 0.002$	1.43	3.29	4.53	0.43
$σ_{2} = 0.004$	2.85	3.49	8.23	2.25

When the standard deviations of the sensitivity matrix are 0.002 and 0.004, the images reconstructed by the RSOPNMFDR algorithm are shown in Figures 9 and 10, respectively. It can be observed from Figures 9 and 10 that owing to the fact that the inaccurate property of the sensitivity matrix is considered in the proposed objective functional, the quality of the images reconstructed by the RSOPNMFDR algorithm under the conditions of the inaccuracy of the sensitivity matrix is satisfactory. Particularly, it can be seen from Table 5 that when the standard deviation of the sensitivity matrix is 0.004, the image errors for the original dynamic reconstruction objects, Figures 3(a)–(d), are 2.85%, 3.49%, 8.23% and 2.25%, which indicates that the RSOPNMFDR algorithm is robust to the inaccuracy of the sensitivity matrix.

Case 5

In order to further evaluate the robustness of the RSOPNMFDR algorithm, in this section the influences of the inaccuracies of the sensitivity matrix and the measurement data on the reconstruction quality are studied. In this case, algorithmic parameters for the RSOPNMFDR algorithm are the same as before. Figures 11 and 12 illustrate the images reconstructed by the RSOPNMFDR algorithm when the inaccurate properties on the sensitivity matrix and the measurement data, and the image errors are presented in Table 6.

Table 6.

Image errors (%).

Noise levels	Standard deviations	Figure 3(a)	Figure 3(b)	Figure 3(c)	Figure 3(d)
5 %	$σ_{2} = 0.001$	1.67	2.42	4.47	1.78
10 %	$σ_{2} = 0.003$	3.27	3.39	7.85	3.86

Figures 11 and 12 present the images reconstructed by the RSOPNMFDR algorithm when the inaccurate properties on the sensitivity matrix and the measurement data are simultaneously considered. It can be observed from Figures 11 and 12 that owing to the fact that the inaccurate properties on the sensitivity matrix, the measurement data and the reconstruction model are considered in the proposed objective functional and the ill-posed nature is illustrated via the dimensionality-reduction of the unknown variables, the quality of the images reconstructed by the RSOPNMFDR algorithm under the conditions of the inaccuracies of the sensitivity matrix and the measurement data is satisfactory. In particular, it can be seen from Table 6 that when the standard deviation of the sensitivity matrix is 0.003 and the noise level of the capacitance data is 10%, the image errors for the original dynamic reconstruction objects, Figures 3(a)–(d), are 3.27%, 3.39%, 7.85% and 3.86%, which confirms that the RSOPNMFDR algorithm is robust to the inaccuracies of the sensitivity matrix and the measurement data.

Case 6

In this section the noise-contaminated capacitance data is used to further evaluate the robustness of the RSOPNMFDR algorithm. Subfigures (a), (b) and (c) in Figure 13 stand for the different permittivity distributions at time instants $k$ , $k + 1$ and $k + 2$ , respectively. In Figure 13, the permittivity of the black part is defined as 2.6, and the permittivity of the rest of reconstruction region is 1.0. In the RSOPNMFDR algorithm, $α = 0.0004$ , $α_{1} = 1$ , $α_{2} = 10$ , $α_{3} = 0.1$ , $α_{4} = 1$ , $μ = 1$ , and the number of iterations is 10. The images reconstructed by the RSOPNMFDR algorithm under the noise levels of 3%, 8% and 15% are illustrated in Figures 14 –16, respectively. The image errors are listed in Table 7.

Table 7.

Image errors (%).

Noise levels	Figure 13(a)	Figure 13(b)	Figure 13(c)
3%	6.43	3.31	7.06
8%	7.58	3.34	7.61
15%	9.36	4.22	9.64

Figures 14 –16 show the images reconstructed by the RSOPNMFDR algorithm under the noise levels of 3%, 8% and 15%, respectively. Owing to the fact that the M-estimation technique and the Tikhonov regularization method are introduced to the proposed objective functional, as can be expected, the RSOPNMFDR algorithm shows satisfactory robustness, and the reconstructed images under different noise levels are in a good agreement with the original dynamic objects. It can be observed from Table 7 that the RSOPNMFDR algorithm is robust to the inaccuracies of the measurement data. However, it can be found from Table 7 that the image errors increase with the rising noise levels, which indicates that the measurement noises needs to be carefully tackled in real applications, and this issue should be further investigated in the future.

Case 7

In this section, the effect of the number of the basis vector on the reconstruction quality is numerically evaluated. Figure 17 shows the image errors of the RSOPNMFDR algorithm for original dynamic reconstruction objects under different numbers of the basis vectors.

Following the discussions presented in previous sections, it can be seen that differing from existing reconstruction methods, in the RSOPNMFDR algorithm the images are indirectly reconstructed by solving a low dimensional coefficient vector. It can be found from Figure 17 that for the dynamic reconstruction cases presented in this paper, the number of the basis vectors will affect the reconstruction quality, which indicates that the snapshots that are used to solve the basis vectors should be cautiously determined. In fact, previous studies indicate that the basis vectors from a poor set of snapshots may impede the improvement of the reconstruction quality.

The simulation results presented in previous sections indicate that the RSOPNMFDR algorithm is a promising candidate for imaging the dynamic reconstruction objects. However, it is worth emphasizing that this paper presents a general framework for solving inverse problems in dynamic image reconstruction tasks, which may be useful for the image reconstruction problems in other related fields.

Conclusions

ECT image reconstruction process is an ill-posed problem, in which the accuracy and speed of image reconstruction algorithms play a crucial role in real applications. In this paper, a multiple measurement vector-based dimensionality reduction dynamic reconstruction model that utilizes the ECT measurement information and the dynamic evolution information of a dynamic object is proposed. A RSOPNMF method is proposed for extracting the basis vectors from a set of snapshots, and the SBI method is used to solve the RSOPNMF model. The original unknown variables are projected onto the subspaces spanned by a set of the basis vectors extracted by the RSOPNMF method from a set of snapshots to obtain a low-dimensional model where the images are indirectly reconstructed by solving the corresponding coefficients, and thus the dimensionality of the unknown variables is reduced and the reconstruction speed is improved. Based on the multiple measurement vectors and the dimensionality reduction model, an objective functional that incorporates the ECT measurement information, the dynamic evolution information of a dynamic object, the spatial constraint and the temporal constraint is proposed, in which the unknown variables are solved in a batching pattern. An iterative scheme that integrates the beneficial advantages of the SBI method and the FBS algorithm is developed for solving the proposed objective functional. Numerical simulation results indicate that the proposed algorithm is advantageous on the improvement of the reconstruction quality and the reconstruction speed. As a result, a promising algorithm is provided for solving ECT inverse problems.

Applications indicate that different algorithm may obtain different numerical results when different reconstruction tasks are implemented, and the choice of a specific reconstruction algorithm depends mainly on the measurement requirements and the understanding of a reconstruction object. Our work provides a candidate method for solving ECT image reconstruction problems, which needs to be further validated by more cases and improved in the future.

Footnotes

Funding

The authors wish to thank the National Natural Science Foundation of China (grant number 51206048) and the Fundamental Research Funds for the Central Universities (grant number 13MS11) for supporting this research.

References

Banasiak

Soleimani

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

Bao

Liu

Yan

(2012) Inductive robust principal component analysis. IEEE Transactions on Image Processing 21(8): 3794–3800.

Blondel

Seki

Uehara

(2013) Block coordinate descent algorithms for large-scale sparse multiclass classification. Machine Learning 93(1): 31–52.

Butala

Hewett

Frazin

Kamalabadi

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

Cai

Osher

Shen

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

Cai

Osher

Shen

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

Candès

Wright

(2011) Robust principle component analysis?Journal of the ACM 58(3): 1–37.

Cao

Fan

Wang

(2011) Electrical capacitance tomography for sensors of square cross sections using Calderon’s method. IEEE Transactions on Instrumentation and Measurement 60(3): 900–907.

Chan

Wong

(2000) Convergence of the alternating minimization algorithm for blind deconvolution. Linear Algebra and its Applications 316(1–3): 259–285.

10.

Choi

(2008) Algorithms for orthogonal nonnegative matrix factorization. In: IEEE International Joint Conference on Neural Networks, Hong Kong, China, 1–8 June, 2008, pp. 1828–1932.

11.

Cichocki

Amari

(2002) Adaptive Blind Signal and Image Processing. New York: John Wiley & Sons.

12.

Cichocki

Zdunek

Phan

Amari

(2009) Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-way Data Analysis and Blind Source Separation. New York: John Wiley & Sons.

13.

Combettes

Wajs

(2005) Signal recovery by proximal forward–backward splitting. SIAM Journal on Multiscale Modeling and Simulation 4(4): 1168–1200.

14.

Ding

Peng

Park

(2006) Orthogonal nonnegative matrix tri-factorizations for clustering. In: Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Philadelphia, PA, 20-23 August 2006, pp. 126–135.

15.

Duchi

Singer

(2009) Efficient online and batch learning using forward backward splitting. Journal of Machine Learning Research 10(12): 2899–2934.

16.

Edelman

Arias

Smith

(1999) The geometry of algorithm with orthogonality constraints. SIAM Journal on Matrix Analysis and Applications 20(2): 303–353.

17.

Engl

Hanke

Neubauer

(2000) Regularization of Inverse Problems. Dordrecht: Kluwer Academic Publishers.

18.

Fang

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

19.

Goldstein

Osher

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

20.

Goldstein

Bresson

Osher

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

21.

Golub

Heath

Wahba

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

22.

Han

(2011) Computational Methods and Applications of Nonlinear Ill-posed Problems. Beijing: Science Press.

23.

Han

(2010) Nonnegative principle component analysis for cancer molecular pattern discovery. IEEE/ACM Transactions on Computational Biology and Bioinformatics 7(3): 537–549.

24.

Hansen

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

25.

Huang

Zhao

(2005) Inverse Problem and Applications in the Electromagnetic Field. Beijing: Science Press.

26.

Huber

(1981) Robust Statistics. New York: John Wiley & Sons.

27.

Jang

Lee

Kim

Choi

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

28.

Kirsch

(1996) An introduction to the Mathematical Inverse Problems. New York: Springer–Verlag.

29.

Kortschak

Wegleiter

Brandstatter

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

30.

Landweber

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

31.

Lee

Seung

(1999) Learning of the parts of objects by nonnegative matrix factorization. Nature 401: 788–791.

32.

Lee

Seung

(2001) Algorithms for non-negative matrix factorization, In: Advances in Neural Information Processing Systems, 13, pp. 556–562. Cambridge, MA: MIT Press.

33.

Lei

Liu

Wang

(2012) Dynamic inversion in electrical capacitance tomography using the ensemble Kalman filter. IET Science, Measurement and Technology 6(2): 63–77.

34.

Lei

Wang

Liu

(2013) Dynamic reconstruction algorithm based on split Bregman iteration for electrical capacitance tomography. Transactions of the Institute of Measurement and Control 35(4): 476–490.

35.

Ouyan

Gao

Yang

Song

(2005) Multiscale Simulation of the Complex System of the Particle Fluid. Beijing: Science Press.

36.

Yang

(2008) Image reconstruction by nonlinear Landweber iteration for complicated distributions. Measurement Science and Technology 19(9): 1–8.

37.

Liu

(2005) Regularization Methods of the Ill-posed Problems and its Applications. Beijing: Science Press.

38.

Liu

Yang

Wang

H.G.

Jiang

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

39.

Song

Wang

(2013) Image restoration with a high-order total variation minimization method. Applied Mathematical Modelling 37(16–17): 8210–8224.

40.

Montefusco

Lazzaro

Papi

Guerrini

(2011) A fast compressed sensing approach to 3D MR image reconstruction. IEEE Transactions on Medical Imaging 30(5): 1064–1075.

41.

Morozov

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

42.

Mou

Peng

Yao

Xiao

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

43.

Nguyen

Tran

(2013) Robust lasso with missing and grossly corrupted observations. IEEE Transactions on Information Theory 59(4): 2036–2058.

44.

Nishimori

Akaho

(2005) Learning algorithms utilizing quasi-geodesic flows on the Stiefel manifold. Neurocomputing 67(1): 106–135.

45.

Olmos

Botella

Castillo

Morales

Banqueri

García

(2012) A reconstruction method for electrical capacitance tomography based on image fusion techniques. Digital Signal Processing 22(6): 885–893.

46.

Ortiz-Aleman

Martin

Gamio

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

47.

Osher

Mao

Dong

Yin

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

48.

Qin

Scheinberg

Goldfarb

(2013) Efficient block-coordinate descent algorithms for the group lasso. Mathematical Programming Computation 5(2): 143–169.

49.

Soleimani

Lionheart

WRB

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

50.

Soleimani

Mitchell

Banasiak

Wajman

Adler

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

51.

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(11): 3287–3294.

52.

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(8): 2104–2112.

53.

Tarantola

(2005) Inverse Problem Theory and Methods for Model Parameter Estimation. Philadelphia, PA: SIAM.

54.

Tikhonov

Arsenin

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

55.

Tseng

(2001) Convergence of a block coordinate descent method for nondifferentiable minimization. Journal of Optimization Theory and Applications 109(3): 475–494.

56.

Vogel

(2002) Computational Methods for Inverse Problems. Philadelphia, PA: SIAM.

57.

Wang

Tang

Cao

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

58.

Wang

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

59.

Wang

Zhang

(2013) Nonnegative matrix factorization: a comprehensive review. IEEE Transactions on Knowledge and Data Engineering 25(6): 1336–1353.

60.

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(12): 2198–2210.

61.

Watzenig

Brandner

Steiner

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

62.

Wright

Peng

Ganesh

Rao

(2009) Robust principal component analysis: exact recovery of corrupted low-rank matrices via convex optimization. In: Proceedings of the 23rd Annual Conference on the Neural Information Processing Systems, Vancouver, British Columbia, Canada, 7–10 December 2009, pp. 2080–2088.

63.

Xiao

Zeng

(2011) Restoration of images corrupted by mixed Gaussian-impulse noise via L1–L0 minimization. Pattern Recognition 44(8): 1708–1720.

64.

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, Circuits, Devices and Systems 139(1): 89–98.

65.

Yin

(2013) A block coordinate descent method for regularized multi-convex optimization with applications to nonnegative tensor factorization and completion. SIAM Journal on Imaging Sciences 6(3): 1758–1789.

66.

Yang

Peng

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

67.

Yang

Spink

York

McCann

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

68.

Yang

Oja

(2010) Linear and nonlinear projective nonnegative matrix factorization. IEEE Transactions on Neural Networks 21(5): 734–749.

69.

Yang

Yuan

Laaksonen

(2007) Projective non-negative matrix factorization with applications to facial image processing. International Journal of Pattern Recognition and Artificial Intelligence 21(8): 1353–1362.

70.

Yin

Osher

Goldfarb

Darbon

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

71.

Yoo

Choi

(2010) Nonnegative matrix factorization with orthogonality constraints. Journal of Computing Science and Engineering 4(2): 97–109.

72.

Zhang

Chen

Zhang

(2011) Robust non-negative matrix factorization. Frontiers of Electrical and Electronic Engineering in China 6(2): 192–200.

73.

Zhang

Burger

Osher

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

74.

Zhang

Burger

Bresson

Osher

(2010) Bregmanized nonlocal regularization for deconvolution and sparse reconstruction. SIAM Journal on Imaging Sciences 3(3): 253–276.

75.

Zhang

(2009) Subspace-Based Face Recognition. Beijing: Tsinghua University.

76.

Zhang

Chen

Zhang

(2011a) Robust nonnegative matrix factorization. Frontiers of Electrical and Electronic Engineering in China 6(2): 192–200.

77.

Zhang

Burger

Osher

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

78.

Zhu

(2010) Algorithms and Applications in Particle Filter. Beijing: Science Press.