This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/3.0/).

Electrical capacitance tomography (ECT) attempts to reconstruct the permittivity distribution of the cross-section of measurement objects from the capacitance measurement data, in which reconstruction algorithms play a crucial role in real applications. Based on the robust principal component analysis (RPCA) method, a dynamic reconstruction model that utilizes the multiple measurement vectors is presented in this paper, in which the evolution process of a dynamic object is considered as a sequence of images with different temporal sparse deviations from a common background. An objective functional that simultaneously considers the temporal constraint and the spatial constraint is proposed, where the images are reconstructed by a batching pattern. An iteration scheme that integrates the advantages of the alternating direction iteration optimization (ADIO) method and the forward-backward splitting (FBS) technique is developed for solving the proposed objective functional. Numerical simulations are implemented to validate the feasibility of the proposed algorithm.

Acquiring the spatial distribution information of materials is vital for improving the system efficiency and reducing pollution emission in chemical reactors or multiphase flow units. ECT is a noninvasive imaging technique, which is used to acquire spatial distribution information from inaccessible objects in order to monitor industrial processes. Owing to its distinct advantages such as the non-intrusive sensing, radiation-free nature, high temporal resolution, affordable measuring device and easy implementation, ECT is proven to be useful in industrial process monitoring, multiphase flow measurements, the visualization of combustion flames in porous media and the identification of two-phase flow patterns [1–10].

ECT technology attempts to reconstruct the permittivity distribution of the cross-section via an appropriate reconstruction algorithm from the capacitance measurement data, where reconstructing high-quality images plays a crucial role in real applications. Due to the ill-posed nature of the inverse problem, the ‘soft-field’ effect and the underdetermined problem in ECT image reconstruction, achieving high-accuracy reconstruction of a dynamic object is challenging. The key issue for improving the reconstruction quality has attracted intensive attention, and thus various algorithms, which can be approximately divided into two categories, static and dynamic reconstruction algorithms, had been developed for ECT image reconstruction. Common static reconstruction algorithms include the linear back-projection (LBP) method [11], the Tikhonov regularization method [12], the Landweber iteration algorithm [13–15], the offline iteration and online reconstruction (OIOR) technique [16], the truncated singular value decomposition (TSVD) method [17], the algebraic reconstruction technique (ART) [17], the simultaneous iterative reconstruction technique (SIRT) [17], the genetic algorithm (GA) [18], the generalized vector sampled pattern matching (GVSPM) method [19], the generalized Tikhonov regularization methods [20–23], the simulated annealing (SA) algorithm [24], the neural network algorithm [25], the level set method [26,27]. Detailed discussions about the numerical performance of other reconstruction algorithms can be found in [17,28].

The above-mentioned algorithms have played an important role in promoting the development of ECT technology and found numerous successful applications. It is worth mentioning that static reconstruction algorithms are often used to image a dynamic object [4,8]. However, these approaches exploit only the spatial relationship of the objects of interest, without using any temporal dynamics of the underlying process, which are not optimal for reconstructing a dynamic object unless the inversion solution is temporally uncorrelated. ECT measurement tasks often involve time-varying objects, and will be more applicable to image a dynamic object using a dynamic reconstruction algorithm that considers the temporal correlations of a dynamic object. In the field of ECT image reconstruction, dynamic reconstruction algorithms do not attract enough attention at present. Fortunately, several algorithms, such as the particle filter (PF) technique [29], the Kalman filter (KF) method [30] and the four-dimensional imaging algorithm [31], had been proposed for tackling the dynamic reconstruction tasks. Overall, the investigations of the dynamic reconstruction algorithms in the field of ECT are far from perfect, and finding an efficient dynamic reconstruction algorithm remains a critical issue.

Based on the RPCA method, a dynamic reconstruction model that utilizes the multiple measurement vectors is presented in this paper, where the evolution process of a dynamic object is regarded as a sequence of 2-D images with different temporal sparse deviations from a common background. An objective functional that simultaneously considers the temporal constraint and the spatial constraint is proposed, in which the images are reconstructed in a batching pattern. An iteration scheme that integrates the merits of the ADIO method and the FBS technique is developed for solving the established objective functional. Numerical simulations are implemented to validate the feasibility of the proposed algorithm.

The rest of this paper is organized as follows: based on the RPCA method, a reconstruction model that utilizes the multiple measurement vectors is proposed in Section 2. The original image reconstruction model is formulated into an optimization problem, and a new objective functional is established in Section 3. In Section 4, an iteration scheme that integrates the advantages of the ADIO method and the FBS algorithm is developed for solving the proposed objective functional. In Section 5, numerical simulations are implemented to evaluate the feasibility of the proposed algorithm, and a concise discussion on the numerical results is provided. Finally, conclusions are presented in Section 6.

2.Model Representation2.1.Static Reconstruction Model

The ECT image reconstruction process involves two key phases: the forward problem and the inverse problem. The forward problem solves the capacitance values from a given permittivity distribution. It is worth mentioning that the forward problem is a well-posed problem, and it can be easily solved by numerical methods such as the finite element method or the finite difference method. The relationship between capacitance and the permittivity distribution can be formulated by [17]:
(1)C=QV=−1V∬Γε(x,y)∇ϕ(x,y)dΓwhere Q is the electric charge; V represents the potential difference between two electrodes forming the capacitance; ε (x, y) and ϕ (x, y) indicate the permittivity and electrical potential distributions, respectively; Г stands for the electrode surface.

The inverse problem attempts to estimate the permittivity distribution from the given capacitance data. In real applications, the static linearization image reconstruction model can be simplified as [17]:
(2)SG=C+rwhere G is an n×1 dimensional vector standing for the normalized permittivity distributions; C represents an m×1 dimensional vector indicating the normalized capacitance values; r is an m×1 dimensional vector representing the capacitance measurement noises; S stands for a matrix of dimension m×n, and it is called as the sensitivity matrix in the field of ECT image reconstruction, which can be formulated by [32,33]:
(3)Si,j(x,y)=−∫p(x,y)Ei(x,y)Vi⋅Ej(x,y)Vjdxdywhere Si,j (x, y) defines the sensitivity between the ith electrode and the jth electrode at p(x, y); Ei(x, y)stands for the electric field distribution inside the sensing domain when the ith electrode is activated as an excitation electrode by applying a voltage Vi to it.

2.2.Multiple Measurement Vectors Dynamic Reconstruction Model

Equation (2) only considers the instantaneous measurement information, and uses single measurement data to implement image reconstruction without any considerations of the temporal dynamics of the underlying process, which is not optimal for reconstructing a dynamic object. It is well known that ECT reconstruction objects are often in a dynamic evolution process, and the measurement results at different time instants have a close correlation [4]. Therefore, considering such information may be important for imaging a dynamic object. In this paper, we propose a multiple measurement vectors dynamic reconstruction model, which can be formulated as:
(4)SX+U=Y+Nwhere U is an m×t dimensional matrix representing the model distortions derived from the facts [17,23,34], such as: (1) the simplification of a true physical process, (2) the linearization approximation distortions of the reconstruction model, and (3) physically implementing an ECT sensor; t >1 defines the number of measurements or the measurement time window; Y stands for an m×t dimensional matrix indicating the capacitance measurement data; X represents an n×t dimensional matrix, and each column of X stands for the permittivity distributions at the measurement time window t; N is an m×t dimensional matrix standing for the capacitance measurement noises.

In the static reconstruction model, the solution merely reflects an instantaneous measurement without any considerations of temporal dynamics of the underlying process, whereas in the case of the dynamic reconstruction model it reflects a sequence of temporally successive measurement, such that the temporal correlations of a dynamic object of interest should be imposed. In other word, the dependence of the capacitance measurement on the evolution of the permittivity distribution is explicitly considered in the dynamical reconstruction model. If the evolution of the permittivity distribution does not follow any dynamics, the dynamical reconstruction model reduces to the static reconstruction model. Obviously, the static reconstruction model is a special case of the dynamic reconstruction model.

The PCA method is an efficient data processing technique, which have enjoyed wide popularity in various fields. Unfortunately, the performances of the PCA technique suffer from the outliers in the data matrix, and thus different approaches had been developed for improving the PCA method. The RPCA method provides a new insight for modern data analysis approaches, which alleviates the deficiencies of the PCA method by applying the ℓ1-regularization and the nuclear norm on the matrix entries. Therefore, the RPCA method is robust to grossly corrupted observations of the underlying low-rank matrix. In a word, the RPCA method tries to recover principal component A (modeled by a low-rank matrix) from data matrix D with outliers E (modeled by a sparse matrix), which can be formulated into the following minimization problem under the certain conditions [35,36]:
(5)(A,E)=min{‖A‖∗+λ‖E‖1}subject toD=A+Ewhere, ‖·‖* defines the nuclear norm for a matrix, and it can be specified as
‖A‖∗=∑kσk, and σk are the singular values of matrix A; ‖· ‖1 represents the ℓ1-norm for a matrix, which can be defined as
‖E‖1=∑j(∑i|Ei,j|); λ > 0 is a regularization parameter. It is worth mentioning that term ‖E‖1 is introduced to satisfy the sparsity assumption of matrix E.

Studies indicate that the evolution process of a dynamic object can be regarded as a sequence of images with different temporal sparse deviations from a common background. Motivated by this observation, the following low-rank and sparse decomposition of X can be obtained [37]:
(6)X=X1+X2where X1 is the low-rank matrix component for modeling the ‘background components’ of X. It is worth mentioning that X1 is assumed to resemble each other rather than to be constant in time, which can be described as a low-rank matrix in mathematics under the RPCA framework [37]. On the other hand, X2 is the sparse matrix component for modeling the sparse deviation from the background X1. Submitting Equation (6) to Equation (4) yields the following expression:
(7)S(X1+X2)+U=Y+N

3.Design of the Objective Functional

Equation (7) is a matrix-based reconstruction model, and it includes three unknown matrix variables. Therefore, it is hard to directly solve the equation. ECT image reconstruction process is a typical ill-posed problem, methods that ensure the numerical stability while improving the quality of an inversion solution should be employed. The Tikhonov regularization technique is an efficient method for solving the ill-posed problems, which has enjoyed wide popularity in various fields. A Tikhonov regularization solution is essentially a result of balancing the accuracy and stability of an inversion solution. According to the Tikhonov regularization method and the optimization theory, Equation (7) can be reformulated as:
(8)min{12‖S(X1+X2)+U−Y‖F2+α1‖X1‖∗+α2‖X2‖1+α3Ω(U)+α42‖(X1+X2)W‖F2}where operator ‖·‖F defines the Frobenius norm for a matrix;
‖S(X1+X2)+U−Y‖F2=Tr((S(X1+X2)+U−Y)T(S(X1+X2)+U−Y)), and Tr(·) represents the trace of a matrix; α1 > 0, α2 > 0, α3 > 0 and α4 > 0 are the regularization parameters; Ω represents a stabilizing functional;
‖(X1+X2)W‖F2=Tr(((X1+X2)W)T((X1+X2)W)) stands for a temporal constraint, where W is defined as:
(9)W=[−10⋯01−1⋯001⋯0⋮⋮⋯⋮00⋯−100⋯1]t×t−1

Equation (8) can be called as a generalized Tikhonov functional, where
‖S(X1+X2)+U−Y‖F2 measures the accuracy of an inversion solution, and ‖X1‖*, ‖X2‖1, Ω(U) and
‖(X1+X2)W‖F2 achieve the numerical stability. Particularly, terms ‖X1‖* and ‖X2‖1 are used to model the background component and the sparse deviation from the background of a dynamic process, respectively. It is worth mentioning that
‖(X1+X2)W‖F2 can be considered as a temporal constraint, which is introduced to impose the temporal correlation of a dynamic object.

The design of function Ω is vital for Equation (8). For simplicity, function Ω is defined as:
(10)Ω(U)=‖U‖1

Following the discussions presented in previous sections, an objective functional for ECT image reconstruction can be specified by:
(11)min{12‖S(X1+X2)+U−Y‖F2+α1‖X1‖∗+α2‖X2‖1+α3‖U‖1+α42‖(X1+X2)W‖F2}

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

In traditional reconstruction models, only single measurement data is used to independently implement image reconstruction; in Equation (11), however, multiple measurement vectors are employed to image a dynamic object and the images are reconstructed by a batching pattern, which differs from other vector-based reconstruction algorithms such as the ART method and the GVSPM algorithm.

In Equation (11), the image series are divided into two matrix-based components by the RPCA method, and the evolution process of a dynamic object is considered as a sequence of images with different temporal sparse deviations from a common background. It is worth mentioning that in real applications different constraints can be imposed on the two matrix-based components according to the prior information of a dynamic object, which will facilitate the improvement of the reconstruction quality.

In Equation (11), the temporal constraint and the spatial constraint are simultaneously considered. Particularly, the temporal constraint is introduced to impose the temporal correlations of a dynamic object. Furthermore, the nuclear norm and the ℓ1-norm are used to model the background component and the sparse component of a dynamic evolution process, respectively.

The measurement noises and the model approximation distortions are simultaneously considered in Equation (11), which is highly suitable for real applications. In particular, the accuracy measure of an inversion solution is considered under a measurement time window, which differs from single measurement vector-based methods.

ECT image reconstruction process is a typical ill-posed problem, methods that ensure the numerical stability while improving the quality of an inversion solution should be employed. In Equation (11), the Tikhonov regularization technique is introduced to ensure the numerical stability of an inversion solution.

The unknown variables in Equation (11) are three matrices, and the computational approaches and strategies will be distinctly different from other vector-based reconstruction algorithms such as the ART method and the GVSPM algorithm.

4.Solving of the Objective Functional

Developing an efficient algorithm to solve Equation (11) is crucial for real applications of the proposed dynamic reconstruction model. In this section, the ADIO method and the FBS technique are concisely introduced, and an iteration scheme that integrates the advantages of the both methods is developed for solving Equation (11).

4.1.Alternating Direction Iteration Optimization Method

Equation (11) includes three unknown matrix variables U, X1 and X2, and directly solving the equation is challenging. In the ADIO method, different unknown variables can be alternately solved [38], and thus Equation (11) can be reformulated into the following three sub-problems:
(12)Uk+1=minU{12‖S(X1k+X2k)+U−Y‖F2+α3‖U‖1}where
‖S(X1k+X2k)+U−Y‖F2=Tr((S(X1k+X2k)+U−Y)T(S(X1k+X2k)+U−Y)).

According to the above discussions, Equation (11) is divided into three sub-problems, in which Equation (12) can be solved by the shrinkage algorithm [39]. Obviously, solving Equations (13) and (14) plays a crucial role in real applications.

4.2.Forward-Backward Splitting Technique

The FBS technique is originally proposed for solving the following optimization problem [40–44]:
(15)minu{μJ(u)+H(u)}where J and H are the known functions.

According to the corresponding deductions, the resulting FBS algorithm can be formulated as:
(16)uk+1=ProxδμJ(uk−δ∂H(uk))where the proximal operator ProxδμJ (ν) is defined as:
(17)argminu{μJ(u)+12δ‖u−v‖2}

In the case of
H(u)=12‖Au−f‖2, we can obtain ∂H(u) = AT(Au − f). Therefore, Equation (15) can be solved by the following two-step algorithm:
(18)vk+1=uk−δAT(Auk−f)(19)uk+1=argminu{μJ(u)+12δ‖u−vk+1‖2}

4.3.Proposed Iteration Scheme

Following the above discussions, an iteration scheme can be developed for solving Equation (11), which can be summarized as follows:

Additionally, it can be known in advance that the inversion solution belongs to the range [Θ1, Θ2 ], therefore a projected operator is introduced to the iteration scheme:
(20)Xk+1=Project{Xk+1}where:
(21)Project[Di,j]={Θ1,ifDi,j<Θ1f(x),ifΘ1≤Di,j≤Θ2Θ2ifDi,j>Θ2

5.Numerical Simulations and Discussions

According to the above discussions, the proposed reconstruction technique can be concisely called as the multiple measurement vectors dynamic reconstruction (MMVDR) algorithm. In this section, dynamic reconstruction cases are implemented to evaluate the feasibility of the MMVDR algorithm, and the reconstruction quality is compared with the projected Landweber iteration (PLI) method. The initial values are computed by the standard Tikhonov regularization method. All algorithms are implemented using the MATLAB 7.0 software on a PC with a Pentium IV 2.4 G Hz CPU, and 4 G bytes memory. The image error is used to evaluate the reconstruction quality, which is defined as [17]:
(22)η=‖GTrue−GEstimated‖‖GTrue‖×100%where, GTrue and GEstimated represent the true and estimated permittivity distributions, respectively. A 12 electrodes square ECT sensor, which is present in Figure 1, is selected for simulations and an image is presented using 32 × 32 pixels.

Figure 1 illustrates a square ECT sensor, which contains an array of 12 electrodes around the frame. The sensor is enclosed by an earthed shielding to protect the sensor against the effects of external charged objects. The size of the frame is 80 × 80 mm, and the width of the frame is 6 mm. The length and width of the electrodes are 18 mm and 1 mm, respectively. The length and width of the axial guards are 4 mm and 1 mm, respectively. The size of shielding is 108 × 108 mm.

5.1.Case 1

In the section, the image reconstructions with three measurement vectors (t = 3) are firstly implemented to evaluate the feasibility of the MMVDR algorithm and the reconstruction quality is compared with the PLI method. Subfigures a, b and c in Figure 2, where the black color stands for the high permittivity materials with a value of 2.6 and the white color represents the low permittivity materials with a value of 1.0, stand for the dynamic permittivity distributions of a dynamic object at time instants t, t + 1 and t +2. The diameters of the cylinders in Figure 2 are 20 mm.

Table 1 lists the algorithmic parameters for the PLI method. In the MMVDR algorithm, α1 = 1.5, α2 = 0.05, α3 = 0.07, α4 = 0.001, and the number of iterations is 130. Figures 3 and 4 illustrate the images reconstructed by the PLI method and the MMVDR algorithm, respectively. The image errors are shown in Table 2.

The images reconstructed by the PLI method are presented in Figure 3. Numerical simulation results indicate that distinct advantages of the PLI method involve easy numerical implementation and low computational complexity and cost owing to the fact that only gradient information of the objective functional is used. However, the PLI fails to consider the temporal correlations of a dynamic object, and the improvement of the reconstruction quality is therefore restricted. For the cases simulated in this section, the quality of the images reconstructed by the PLI method is far from perfect and the distortions are relatively serious.

The images reconstructed by the MMVDR algorithm are illustrated in Figure 4. Numerical simulation results indicate that the MMVDR algorithm can ensure the numerical stability of an inversion solution owing to the fact that the Tikhonov regularization technique is used to stabilize the numerical solution. Particularly, the computational complexity and cost of the MMVDR algorithm are low, and the numerical implementation is easy. As can be expected, it can be seen from Figure 4 that owing to the fact that the temporal constraint and spatial constraint are simultaneously considered, the quality of the images reconstructed by the MMVDR algorithm is improved as compared to the PLI method, which indicates that the MMVDR algorithm is competent in solving ECT inverse problems.

Meanwhile, it can be found that in the MMVDR algorithm, the images are reconstructed by a batching pattern, and thus the temporal correlations of a dynamic object can be better considered, which differs from other vector-based reconstruction algorithms. In addition, it can be found that the MMVDR method is an iterative algorithm, and it is hard to achieve the online reconstruction presently. In the future, more investigations on improving reconstruction speed should be implemented.

The image errors are listed in Table 2. It can be found that the quality of the images reconstructed by the MMVDR algorithm is higher than with the PLI method, which confirms that the MMVDR algorithm is a promising candidate for solving ECT image reconstruction problems.

5.2.Case 2

In this section, the image reconstructions with four measurement vectors (t = 4) are implemented to further evaluate the feasibility of the MMVDR algorithm, and the reconstruction quality is compared with the PLI method.

Figure 5 simulates the dynamic process of the objects of interest under different time instants. Subfigures a, b, c and d in Figure 5, where the black color stands for the high permittivity materials with a value of 2.6 and the white color represents the low permittivity materials with a value of 1.0, represent the dynamic permittivity distributions of a dynamic object at time instants t, t + 1, t + 2 and t + 3. The diameters of the cylinders in the subfigures a, b and d in Figure 5 are 20 mm, and the diameter of the cylinder in Figure 5(c) are 30 mm.

Algorithmic parameters for the PLI method are listed in Table 3. The algorithmic parameters for the MMVDR algorithm are the same as Section 5.1. The images reconstructed by the PLI method and the MMVDR algorithm are illustrated in Figures 6 and 7, respectively. Table 4 lists the image errors for the PLI algorithm and the MMVDR algorithm.

The images reconstructed by the PLI method are presented in Figure 6. It can be seen that for the dynamic reconstruction case simulated in this section, the quality of the images reconstructed by the PLI method is far from perfect and the distortions are relatively large.

Figure 7 illustrates the images reconstructed by the MMVDR algorithm. As can be expected, it can be seen from Figure 7 that the quality of the images reconstructed by the MMVDR algorithm is improved as compared to the PLI method. At the same time, it can be observed from Table 4 that for the case simulated in this section, the MMVDR algorithm gives the smallest image errors, which indicates that the MMVDR algorithm is successful in solving ECT inverse problems.

5.3.Case 3

In this section, the noise-contaminated capacitance data is used to evaluate the robustness of the MMVDR algorithm. In this case, two measurement vectors (t = 2) are used to implement the image reconstruction. The noise level is defined by [17]:
(23)γ=‖CContaminated−CTrue‖‖CTrue‖×100%

Where γ is the noise level; CTrue stands for the true capacitance data; CContaminated is the noise-contaminated capacitance data; CContraminated = CTrue + r, where r = σ1 · 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

Figure 8 shows the dynamic process of the objects of interest, in which subfigures a and b represent the dynamic permittivity distributions of a dynamic object at time instants t and t + 1. In Figures 8(a) and b, the diameters of the cylinders are 20 mm, the permittivity of the cylinder is 2.6, and the permittivity of the rest of the reconstruction region is 1.0.

In the MMVDR algorithm, the algorithmic parameters are the same as Section 5.1. Figures 9–12 show the images reconstructed by the MMVDR algorithm under the noise levels of 0, 9%, 24% and 33%. The image errors are listed in Table 7.

Figures 9–12 are the images reconstructed by the MMVDR algorithm under the noise levels of 0, 9%, 24% and 33%, respectively. As can be expected, it can be seen from Figures 9–12 that the MMVDR algorithm shows satisfactory robustness, and the quality of the images reconstructed under different noise levels is acceptable, which is highly desirable for real applications owing to the fact that the measurement noises are inevitable in practice. When the noise level is 33%, especially, the image errors for the dynamic reconstruction objects, Figure 8(a,b), are 13.81% and 10.20%, which indicates that the MMVDR algorithm is successful in treating with the measurement noises. Additionally, it can be seen from Figures 9–12 and Table 5 that the image errors increase with the increment of the noise levels, which indicates that the inaccuracy of the capacitance measurement data should be seriously tackled and this issue should be further studied in the future.

6.Conclusions

ECT is considered a promising visualization measurement technology, in which reconstructing high quality images is highly desirable for real applications. In this paper, based on the RPCA technique, a dynamic reconstruction model that utilizes the multiple measurement vectors is presented, in which the evolution process of a dynamic object is considered as a sequence of images with different temporal sparse deviations from a common background. An objective functional that simultaneously considers the temporal constraint and the spatial constraint is proposed, where the images are reconstructed in a batching pattern. An iteration scheme that integrates the advantages of the ADIO method and the FBS technique is developed for solving the proposed objective functional. Numerical simulation results indicate that the proposed algorithm can ensure a stable numerical solution. For the cases simulated in this paper, the quality of the images reconstructed by the proposed algorithm is improved, which indicates that the proposed algorithm is successful in solving ECT inverse problems. As a result, a promising algorithm is introduced for ECT image reconstruction.

Applications indicate that each algorithm may show different numerical performances to different reconstruction tasks. In practice, the selection of a 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 validated by more cases in the future. At the same time, more investigations on the improvement of the reconstruction speed should be undertaken.

The authors wish to thank the National Natural Science Foundation of China (Nos. 51206048, 61072005 and 51006106) and the Program for the Changjiang Scholars and Innovative Research Team in University (No.IRT0952) for supporting this research.