Dynamic imaging method using the low n-rank tensor for electrical capacitance tomography

Flow Measurement and Instrumentation 41 (2015) 104–114

Contents lists available at ScienceDirect

Flow Measurement and Instrumentation journal homepage: www.elsevier.com/locate/flowmeasinst

Dynamic imaging method using the low n-rank tensor for electrical capacitance tomography J. Lei a,n, W.Y. Liu a, S. Liu a, X.Y. Wang b a b

School of Energy, Power and Mechanical Engineering, North China Electric Power University, Changping District, Beijing 102206, China Institute of Engineering Thermophysics, Chinese Academy of Sciences, Haidian District, Beijing 100190, China

art ic l e i nf o

a b s t r a c t

Article history: Received 11 August 2014 Received in revised form 29 October 2014 Accepted 2 November 2014 Available online 11 November 2014

Imaging objects in electrical capacitance tomography (ECT) measurement are often in a dynamic evolution process, and exploiting the spatial–temporal properties of the dynamic reconstruction objects is crucial for the improvement of the reconstruction quality. Based on the multiple measurement vectors, in this paper a robust dynamic reconstruction model that incorporates the ECT measurement information and the dynamic evolution information of a dynamic object, in which a series of dynamic images is cast as a third-order tensor that the first two dimensions are space and the third is time, is proposed. Under the considerations of the two-dimensional spatial structure property of a difference image and the spatial–temporal property of a third-order image tensor, a new objective functional that fuses the ECT measurement information, the dynamic evolution information, the temporal constraint, the spatial constraint, the low rank constraint of a difference image and the low n-rank constraint of a third-order tensor is proposed, where the images are reconstructed by a batching pattern. The split Bregman iteration (SBI) algorithm is developed for solving the proposed objective functional. Numerical simulations are implemented to demonstrate the advantages of the proposed algorithm on improving the reconstruction quality and the robustness. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Electrical capacitance tomography Volume imaging Image tensor Low n-rank constraint Image reconstruction

1. Introduction

1.1. Dynamic behaviors of the measurement objects

As a promising noninvasive visualization measurement method, ECT technology attempts to reconstruct the permittivity distributions of the cross-section via an appropriate algorithm from the given capacitance measurement data, in which reconstructing high-quality images plays a crucial role in the understanding of the underlying physical or chemical mechanisms of the dynamic behaviors of a dynamic imaging object. Currently, various methods, including static reconstruction methods [1–20] and dynamic reconstruction algorithms [21–25], have been developed for improving the reconstruction quality. ECT imaging objects are often in a dynamic evolution process and the key difference of the static reconstruction and the dynamic reconstruction dependents mainly on the additional time dimension. Consequentially, the dynamic reconstruction has the following distinct features.

ECT measurement objects, such as the multiphase flows and the combustion flame, are in a dynamic evolution process, and thus considering the dynamic evolution information or corresponding prior information of a dynamic object will be crucial for the improvement of the reconstruction quality. Though static reconstruction methods are often employed to image a dynamic object, these methods naturally fail to exploit the dynamic evolution information of a dynamic object, which may not optimal for dynamic reconstruction tasks. 1.2. Temporal correlations In dynamic reconstruction tasks, the measurement images at adjacent frames present the temporal correlations. Obviously, exploiting the inter-frame correlations will be important for the improvement of the reconstruction quality. Unfortunately, the temporal correlations are omitted in static reconstruction algorithms. 1.3. Image volume property

n

Corresponding author. Tel.: þ 86 10 61772472; fax: þ 86 10 61772219. E-mail address: [email protected] (J. Lei).

http://dx.doi.org/10.1016/j.flowmeasinst.2014.11.001 0955-5986/& 2014 Elsevier Ltd. All rights reserved.

For a dynamic object, at a measurement cross-section, a series of two-dimensional dynamic images will form an image volume or a third-order tensor along the time direction at a specific

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

105

ECT image reconstruction task as an optimization problem is proposed. Section 6 presents an iteration scheme designed for solving the proposed objective functional. In Section 7, numerical simulations are implemented to evaluate the feasibility of the proposed algorithm and the detailed discussions on the numerical results are provided. Finally, Section 8 outlines the main conclusions.

2. Related work

Fig. 1. Image volume. (a) A series of dynamic images; (b) third-order tensor.

measurement time window due to the temporal consistency, and thus introducing effective data analysis methods to exploit such spatial–temporal properties will facilitate the improvement of the reconstruction quality. Owing to the above distinct features, naturally, in this paper the following key problem is studied: how to exploit the spatial– temporal properties of a dynamic object to improve the quality of the image reconstruction? The major contributions of this paper are outlined as follows. (1) The space–time volume [26–28], which stacks a set of dynamic images to form a three-dimensional data structure (see Fig. 1), is introduced to ECT image reconstruction. A series of dynamic images is represented by a three-dimensional function f ðx; y; tÞ, where ðx; yÞ denotes the coordinate in space and t defines the coordinate in time. Suppose that each image has M rows, N columns, and there are K frames, then, the dynamic images f ðx; y; tÞ for x ¼ 1; 2; …; M and y ¼ 1; 2; …; N form a third-order tensor of size M  N  K. It is found that common matrix or vector based data analysis methods may not be optimal in extracting the spatial– temporal information from a multidimensional perspective. On the other hand, the tensor-based data analysis method has shown that tensor models are able to take full advantage of the multilinear structures to provide better understanding of a dynamic object of interest. In particular, the tensor representation can better expose the spatial–temporal structure of a series of dynamic images. (2) A multiple measurement vectors based dynamic reconstruction model that integrates the dynamic evolution information of a dynamic object and ECT measurement information is proposed. Differing from common reconstruction methods, particularly, the inaccurate properties on the reconstruction model, the sensitivity matrix and the measurement data are simultaneously emphasized. (3) Differing from common vector-based reconstruction methods, in order to emphasize the temporal correlations of a dynamic object and the spatial structure property of a two-dimensional image, the low rank constraint of a difference image at adjacent frames is introduced in this paper. (4) By introducing the regularization functions along the spatial and temporal directions, both spatial and temporal constraints are imposed. In contrast to existing reconstruction methods, especially, the image volume is cast as a third-order tensor in this paper, and the low n-rank constraint of a tensor is introduced to reveal the spatial–temporal information of a dynamic object. (5) This paper presents a general framework for the improvement of the reconstruction quality in dynamic ECT image reconstruction tasks, which may be useful for image reconstruction problems in other related tomography fields. The rest of this paper is organized as follows. In Section 2, previous studies on ECT reconstruction algorithms are reviewed concisely. A multiple measurement vectors based robust dynamic reconstruction model is proposed in Section 3. Section 4 concisely introduces the tensor theory. In Section 5, a new objective functional that casts the

ECT image reconstruction process consists of two key procedures, the forward problem with a aim of computing the capacitance values from the known permittivity distributions and the inverse problem with a motivation of estimating the permittivity distribution from the given capacitance data. In a mathematical notation, common ECT image reconstruction model can be approximated as [7] Sx ¼ C þ r

ð1Þ

where S represents a matrix of dimension m  n; x is an n  1 dimensional unknown vector; C stands for a capacitance vector of dimension m  1 and r is an m  1 dimensional vector indicating the stochastic measurement noises. More details on ECT image reconstruction model can be traced back to [1,7,16]. Eq. (1) is a basic linearization model for ECT image reconstruction, and it can be modified according to other assumptions or prior information. The image reconstruction method plays a crucial role in successful applications of ECT visualization measurement technique. Presently, numerous algorithms, including static and dynamic reconstruction algorithms, had been developed for ECT image reconstruction. Common static reconstruction algorithms include the linear back-projection (LBP) method [1], the standard Tikhonov regularization (STR) method [2], the Landweber iteration algorithm and its variants [3–5], the offline iteration and online reconstruction (OIOR) algorithm [6], the truncated singular value decomposition (TSVD) method [7], the algebraic reconstruction technique (ART) [7], the simultaneous iterative reconstruction technique (SIRT) [8], the genetic algorithm (GA) [9], the generalized vector sampled pattern matching (GVSPM) method [10], the generalized Tikhonov regularization techniques [11–13], the least squares support vector machine based method [14], the simulated annealing (SA) method [15], the neural network algorithm [16], the level set method [17], the Calderon method [18], the nonlinear Landweber iteration algorithm [19] and the dbar method [20]. Popular dynamic reconstruction algorithms include the particle filter (PF) method [21], the Kalman filter (KF) algorithm [22], the four-dimensional imaging technique [23], the Ensemble Kalman filter (EnKF) method [24] and the generalized dynamic inversion algorithm [25]. Static reconstruction methods have been intensely studied in the past years, which may not be optimal for a dynamic object because the dynamic evolution information and temporal correlations of a dynamic object are omitted. With the increasing requirements of real applications, dynamic reconstruction algorithms had attracted increasing attention in recent years. Owing to the complexities and particularities of the dynamic reconstruction tasks, seeking an efficient algorithm remains an open problem. More importantly, existing algorithms, including static and dynamic reconstruction methods, fail to fully consider three-dimensional spatial–temporal properties of the image volume, and thus the improvement of the reconstruction quality may be restricted.

3. Multiple measurement vectors based dynamic reconstruction model For a dynamic object, the space-time volume at a specific time window can be written as X ¼ ½x1 ; x2 ; …; xt , in which the each

106

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

column of X represents the permittivity distribution at a specific instant. Owing to the limitations of the imperfect measurement conditions and the ubiquitous measurement noises, common reconstruction models often focus on the measurement noises. It is found that the model approximation deviations derived from the approximation and discretization of a real problem and physically implementing an imperfect ECT sensor in real applications and the approximation solution of the sensitivity matrix will undoubtedly introduce deviations. Especially, most of previous studies assume that measurement noises have bounded energy. Without such assumption, unfortunately, it is hard to ensure a reliable estimation result [29]. Owing to the restrictions of the imperfect measurement conditions, unfortunately, it is hard to satisfy such requirements in real applications. Differing from common reconstruction models, in this paper the outliers in the measurement data are emphasized with a motivation of satisfying the requirements of real measurement conditions. Naturally, in order to exploit the temporal correlations of a dynamic object of interest and the inaccurate properties on the sensitivity matrix, the reconstruction model and the measurement data, a new multiple measurement vectors based robust reconstruction model is proposed, which can be formulated as ðS þEÞX þ U ¼ Y þ R þ N

ð2Þ

where U is an m  t dimensional matrix defining the model approximate deviations; E defines an m  n dimensional matrix representing the inaccuracy of the sensitivity matrix; Y stands for an m  t dimensional matrix indicating the capacitance measurement data, and its each column defines the capacitance measurement vector at a specific measurement instant; N is an m  t dimensional matrix standing for the conventional stochastic measurement noises and R defines an m  t dimensional matrix representing the gross errors. Additionally, the similar studies on the inaccuracies of the input data can be found in [30,31]. In real applications, the dynamic evolution process of a dynamic object can be generally formulated as f ðxτ ; vτ Þ ¼ 0

ð3Þ

where f describes the dynamic evolution information of a dynamic object, and it can be specified as different equations according to different measurement objects; vτ depicts the uncertainties of the evolution equation and subscript τ represents the discrete time index. Achieving fast reconstruction is essential for the visualization of the dynamic behaviors of an object of interest, and thus Eq. (3) is approximated by: xτ þ 1 ¼ F τ xτ þ vτ

ð4Þ

where F τ is the state-transition operator at time instant τ; If set F τ ¼ I and I is an identity matrix, Eq. (4) is a purely random-walk evolution model. One of the main bottlenecks that impede the improvement of the reconstruction quality in ECT imaging is the lack of information. In the case of the dynamic reconstruction, common static reconstruction algorithms merely utilize ECT measurement information, and the temporal correlations and the dynamic evolution information of a dynamic object are omitted. Studies indicate that a lack of information cannot be remedied by any mathematical trickery [32,33]. For a dynamic object, in addition to increasing the number of the capacitance measurements, a natural approach of adding the amount of information is to fuse the dynamic evolution information of a dynamic object, which does not require the additional and costly devices. Owing to the limitations of the measurement conditions and the additional and costly devices, moreover, increasing the number of the capacitance measurement may be impractical. As a result, integrating the temporal

correlations and the dynamic evolution information of a dynamic object in ECT image reconstruction will be essential for the improvement of the reconstruction quality. In contrast to existing sequential dynamic reconstruction models, under the consideration of the multiple measurement vectors, a robust dynamic reconstruction model that incorporates the ECT measurement information and the dynamic evolution information of a dynamic object is finally formulated by ( f ðxτ ; vτ Þ ¼ 0 ð5Þ ðS þ EÞX þ U ¼ Y þ R þ N It is worth mentioning that differing from common sequential dynamic reconstruction methods such as the PF method and the KF algorithm, in which the images are sequentially reconstructed one by one, the images in Eq. (5) are reconstructed by means of a batching pattern at a specific time window.

4. Tensor background In this section, some essential nomenclatures and notations used in this paper are concisely introduced, and more details can be found in [34–40]. A tensor is a multi-way array or multidimensional matrix. The order of a tensor is the number of dimensions, also known as ways or modes. In a mathematical notation, a n-order tensor is defined as X A ℝn1 n2 ⋯nN , whose elements are denoted as xi1 ⋯ik ⋯iN , in which 1 r ik r nk and 1 r k r N. The matricization, also known as unfolding or flattening, is the process of reordering the elements of a tensor into a matrix. In this paper, X ðnÞ is defined as the mode-n unfolding of a tensor X . Specially, the mode-n unfolding maps the tensor elements ði1 ; i2 ; …; iN Þ to the matrix element ðin ; jÞ according to j ¼ 1þ

N

∑

ðik  1ÞV k

k ¼ 1;k a n

with

Vk ¼

k1

∏

m ¼ 1;m a n

nm

ð6Þ

That is, the unfold operation on a tensor X is defined as unfoldn ðX Þ ¼ X ðnÞ A ℝnn I n with I n ¼ ∏N k ¼ 1;k a n nk . The opposite operation fold is defined as foldðX ðnÞ Þ ¼ X . As an example, Fig. 2 shows the mode-1, mode-2 and mode-3 unfoldings of a thirdorder tensor. It is worth emphasizing that the latent information of the multi-way data structure can be better revealed by the unfolding operation of a tensor, which is highly applicable for ECT image reconstruction. The inner product of two same-sized tensors X , Y A ℝn1 n2 ⋯nN is the sum of the products of their entries, i.e. n1   X; Y ¼ ∑

n2

nN

∑ ⋯ ∑ xi1 i2 ⋯iN yi1 i2 ⋯iN

i1 ¼ 1 i2 ¼ 1

iN ¼ 1

ð7Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi X ; X . When the tensor X The Frobenius norm is ‖X ‖F ¼ reduces to the matrix X, and ‖X‖F is just the Frobenius norm of the matrix X. The n-mode product of a tensor X A ℝn1 n2 ⋯nN with a matrix ψ A ℝJnn is denoted by X n ψ and is of size n1  ⋯  nn  1 J  nn þ 1  ⋯  nN . Each mode-n fiber is multiplied by matrix ψ. This idea is compactly expressed in terms of the unfolded tensors [38]: Y ¼ X n ψ 3 Y ðnÞ ¼ ψX ðnÞ

ð8Þ

The n-rank of an n-order tensor X is the tuple of the ranks of the mode-n unfoldings, which can be formulated as [38] n  rankðX Þ ¼ ðrank X ð1Þ ; rank X ð2Þ ; …; rank X ðNÞ Þ

ð9Þ

The tensor rank is difficult to handle. However, the n-rank of a tensor on the other hand is easy to compute. Therefore, the tensor

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

difference operation of the images at adjoining frames emphasizes the temporal correlations and the low rank constraint imposes the spatial structure property of a two-dimensional image, and thus the spatial–temporal properties are simultaneously considered. (2) Exploiting the spatial–temporal correlations of a dynamic object is crucial for the improvement of the reconstruction quality. Differing from common reconstruction methods, in this paper a series of dynamic images is cast as a third-order tensor that the first two dimensions are space and the third is time. Studies indicate that the latent spatial–temporal information of the multi-way data structure in an image tensor is better revealed via the mode-n unfolding operation. It is found that owing to the temporal correlations and spatial redundancies of a dynamic object, the matrices of the unfoldings of an image tensor along each mode present the low rank property. In this paper, such prior information is introduced to the proposed objective functional. It is worth emphasizing that because a series of dynamic images is cast as a third-order tensor, the spatial–temporal properties of a dynamic object are better considered, and the underlying spatial–temporal information and details can be more adequately revealed as compared with common reconstruction methods. Following the above discussions, finally, a new objective functional for ECT image reconstruction is specified as n α min 12‖ðS þ EÞX þU  Y  R‖2F þ α1 ‖E‖1 þ α2 ‖U‖1 þ 23 ‖X  X c ‖2F E;X;R;U ! m n   þ α4 ∑ ∑ W ℓ ;ℓ X ℓ ;ℓ  þ α5 ‖R‖1

Fig. 2. Mode-1, mode-2 and mode-3 unfoldings of a third-order tensor.

n-rank minimization problem can be formulated as [38] min f ðn  rankðX ÞÞ s:t: ΑðX Þ ¼ b

107

ℓ2 ¼ 1

ð10Þ

where f ðn  rankðX ÞÞ ¼ f ðrank X ð1Þ ; rank X ð2Þ ; …; rank X ðNÞ Þ and A is a linear map. In real applications, for the sake of easy computation, objective function f can be replaced with the sum of the ranks of different unfoldings f ðn  rankðX ÞÞ ¼ ∑N i ¼ 1 rankðX ðiÞ Þ [38]. It is worth emphasized that function f can be generalized as f ðn  rankðX ÞÞ ¼ ∑N i ¼ 1 ωi rankðX ðiÞ Þ, where ωi can be regarded as a weighted value that indicates the preference towards different unfoldings [36,37]. A different parameter for each unfolding is assigned to make the nrank function more flexible, and thus the low rank assumptions on any mode can be discarded by easily setting the weighted value as zeros. The low n-rank tensor recovery problem is a difficult non-convex optimization problem. The nuclear norm is the greatest convex minorant of the rank function, and thus each rank term can be replaced by the norm for easy computation [38].

5. Design of the objective functional Directly solving Eq. (5) is challenging due to the ill-posed nature of the inverse problem, an alternative approach is to cast the problem as an optimization of an objective functional consisting of the data fidelity and the regularizer that incorporates the prior assumptions on the underlying images, by introducing the Tikhonov regularization technique to ensure a stable inversion solution. Theoretically, a Tikhonov regularization solution is a result of balancing the stability and accuracy of an inversion solution. In this paper, a new objective functional with the following distinct properties is proposed. (1) ECT measurement objects are often in a fast dynamic evolution process and the measurement images present the interframe correlations. When the multiple measurement vectors model is used, it is found that a two-dimensional difference image at adjacent frames presents the low rank property due to the spatial redundancy and the temporal correlations. More theoretical details and successful applications on the low rank representation can be found in [41–47]. In this paper, the low rank property of a twodimensional difference image at a specific time window is introduced to the objective functional. It is worth mentioning that the

ℓ1 ¼ 1

1

2

1

2

n

N

j¼2

i¼1

þ λ1 ∑ ‖X j  X j  1 ‖n þ βi ∑ ‖X ðiÞ ‖n

) ð11Þ

where operator ‖ U ‖F defines the Frobenius norm for a matrix; operator ‖U ‖1 represents the 1-norm for a matrix; operator ‖ U‖n defines the nuclear norm for a matrix, and it can be specified as ‖X j ‖n ¼ ∑ σ k , in which σ k is the singular values of the matrix X j k

that is formed by reshaping the jth column of matrix X; X ðiÞ stands for the unfoldings of a third-order image tensor X ; α1 40, α2 4 0, α3 40, α4 4 0, α5 4 0, λ1 40 and βi 4 0 are the regularization parameters; W is a weighted matrix; X c represents the dynamic evolution information, which can be in advance obtained by solving Eq. (3) or (4). In this paper, X c is computed using the random-walk evolution model [22]. It is worth mentioning that ∑N i ¼ 1 ‖X ðiÞ ‖n is employed to impose the low n-rank constraint of a third-order image tensor [37,38,40]. Eq. (11) belongs to the framework of the Tikhonov regularization method, which has the following appealing properties. (1) Differing from common reconstruction algorithms, in Eq. (11) a series of dynamic images is cast as a third-order tensor, in which the spatial–temporal properties are better presented. Especially, the low n-rank constraint of a tensor is introduced to Eq. (11) to exploit the spatial–temporal properties of a dynamic object, and it will facilitate the improvement of the reconstruction quality. (2) One of the main bottlenecks that restrict the improvement of the reconstruction quality is derived from the lack of information. Differing from existing static reconstruction models, Eq. (11) fuses the dynamic evolution information of a dynamic object and the ECT measurement information, and thus the quantity of information is increased and the reconstruction quality may be improved in theory. Particularly, in contrast to common sequential dynamic reconstruction methods such as the PF method, the KF algorithm and the EnKF method, in which the images are sequentially reconstructed one by one, the images in Eq. (11) are reconstructed by means of a batching pattern at a specific time window. It is worth mentioning that in sequential reconstruction methods, the optimization problem is implemented at a specific measurement instant; however, the optimization problem in this

108

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

proposed reconstruction method is implemented at a specific time window, and thus the inter-frame correlations can be better exploited. (3) The dynamic evolution information and the ECT measurement information are incorporated in Eq. (11). Furthermore, it should be pointed out that differing from common multimodal systems, incorporating the dynamic evolution information of a dynamic object in ECT image reconstruction does not require additive cost, and it is highly attractive for real industrial applications. (4) Studies indicate that it is hard for common static reconstruction algorithms to exploit the temporal correlations of a dynamic object. In Eq. (11), the multiple measurement vectors are used to model the dynamic image reconstruction problem, which will facilitate the improvement of the reconstruction quality and the consideration of the temporal correlations of a dynamic object. (5) Owing to the temporal correlations of a dynamic object, the difference images at adjacent frames will present the low rank property. In Eq. (11), the low rank constraint of a two-dimensional difference image is introduced to impose such prior information, and thus the spatial–temporal properties can be better exploited. (6) Exact descriptions of the dynamic evolution process of a dynamic object of interest and the acquirements of the exact capacitance measurements are impracticable in real applications. Eq. (11) simultaneously considers the inaccurate properties of the dynamic evolution information and the measurement information. In particular, measurement noises, model approximation deviations and sensitivity matrix errors are emphasized in Eq. (11). Differing from existing static and dynamic reconstruction models, more importantly, in Eq. (11) the gross errors in measurement data are also emphasized. (7) The Tikhonov regularization method is introduced to Eq. (11) to overcome the numerical instability derived from the illposed nature in ECT image reconstruction problems, and thus a stable inversion solution can be ensured.

6. Solving of the objective functional In previous sections, ECT image reconstruction task is cast as an optimization problem by employing the framework of the Tikhonov regularization method, and seeking a reliable method to solve the proposed objective functional is crucial for real applications. Owing to the distinct advantages, including easy implementation,

X k þ 1 ¼ min

8 m > > > > 12‖ðS þ E k þ 1 ÞX þ U k þ 1  Y  Rk ‖2F þ α23 ‖X  X c ‖2F þ α4 ∑ > < ℓ2 ¼ 1

X

kþ1

U k þ 1 ¼ min





1 ‖ðS þE k ÞX k þ U  Y Rk ‖2F þ α2 ‖U‖1 2



1 ‖ðS þ EÞX k þU k þ 1  Y  Rk ‖2F þ α1 ‖E‖1 2

ð12Þ  ð13Þ



1 ‖ðS þE k þ 1 ÞX k þ U k þ 1  Y  R‖2F þ α5 ‖R‖1 2

 ð14Þ

9 81 ‖ðS þ Ek þ 1 ÞX þ U k þ 1  Y  Rk þ 1 ‖2F þ α23 ‖X X c ‖2F > > > > = <2 ! m n n N   ¼ min   > þα W X ‖X X ‖ þβ ‖X ‖ þ λ ∑ ∑ ∑ ∑ ℓ1 ;ℓ2 ℓ1 ;ℓ2 4 1 j j1 n ðiÞ n > > > i ; : ℓ ¼1 ℓ ¼1 j¼2 i¼1 2

1

ð15Þ According to the optimization theory, Eq. (15) can be rewritten as the following constrained optimization problem: 8 9 81 kþ1 > ÞX þ U k þ 1  Y  Rk þ 1 ‖2F þ α23 ‖X  X c ‖2F > > > > > 2‖ðS þ E > = < > ! > > m n n N < min     > þ α ∑ ∑ W X ∑ ‖d ‖ þ β ∑ ‖d ‖ þ λ ℓ1 ;ℓ2 ℓ1 ;ℓ2 4 1 1;j n 2;i n > > > i ; : > ℓ2 ¼ 1 ℓ1 ¼ 1 j¼2 i¼1 > > > > > : s:t: d1;j ¼ X j  X j  1 ; d2;i ¼ X ðiÞ

ð16Þ Applying the SBI algorithm to (16) yields 9 8 kþ1 1 > > ÞX þ U k þ 1  Y  Rk þ 1 ‖2F þ α23 ‖X  X c ‖2F > > 2‖ðS þ E > > > > ! > > > > m n n N > >   > > < þ α4 ∑ ∑ W ℓ1 ;ℓ2 X ℓ1 ;ℓ2  þ λ1 ∑ ‖d1;j ‖n þ βi ∑ ‖d2;i ‖n = min ℓ2 ¼ 1 ℓ1 ¼ 1 j¼2 i¼1 > > > > > > > μ n μ2 N > > k 2 k 2 > 1 > > ∑ þ ∑ ‖d  ðX  X Þ b ‖ þ ‖d  X  b ‖ > > 1;j j j1 2;i ðiÞ F F 1;j 2;i > > 2 ; : 2 j¼2

i¼1

ð17Þ ¼ b1;j  ðd1;j  ðX j

kþ1

¼ b2;i  ðd2;i  X kðiÞþ 1 Þ

b2;i

kþ1

kþ1

kþ1

b1;j

k

kþ1

 X j  1 ÞÞ

ð18Þ

kþ1

k

ð19Þ

Similarly, following the computational strategies proposed in [56–61], Eq. (17) can be decoupled as ( ) n μ n k k kþ1 k d1;j ¼ min λ1 ∑ ‖d1;j ‖n þ 1 ∑ ‖d1;j  ðX j  X j  1 Þ  b1;j ‖2F 2 j¼2 j¼2 ð20Þ ( kþ1

d2;i

N

¼ min βi ∑ ‖d2;i ‖n þ i¼1

!9   > >   > ∑ W ℓ1 ;ℓ2 X ℓ1 ;ℓ2 > > = ℓ1 ¼ 1

μ2 N k ∑ ‖d  X ðiÞ  b2;i ‖2F 2 i ¼ 1 2;i

) ð21Þ

n

n > μ N > kþ1 k kþ1 k >  ðX j X j  1 Þ  b1;j ‖2F þ 2 ∑ ‖d2;i X ðiÞ  b2;i ‖2F þ μ1 ∑ ‖d > > : 2 j ¼ 2 1;j 2i¼1

high computational efficiency and low computational complexity, the SBI method has found wide applications in various fields, and more details on the method can be traced back to [48–55]. In this section, the SBI algorithm is introduced to solve Eq. (11). Essentially, Eq. (11) is an unconstrained optimization problem. For the sake of easy computation, following the computational strategies proposed in [56–61], it is numerically appealing to decouple Eq. (11) as:

E k þ 1 ¼ min

Rk þ 1 ¼ min

> > > > > ;

ð22Þ

Consequentially, following the discussions presented in previous sections, an iteration scheme is developed for solving the proposed objective functional, which can be outlined as follows. Step 1. Specify the algorithmic parameters and the initial values. Step 2. Update variable U k þ 1 by solving Eq. (12) using the shrinkage algorithm. Step 3. Update variable E k þ 1 by solving Eq. (13) using the IST algorithm [62]. Step 4. Update variable Rk þ 1 by solving Eq. (14) using the shrinkage algorithm. kþ1 Step 5. Update variable d1;j by solving Eq. (20) using the singular value shrinkage algorithm [63].

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114 kþ1

Step 6. Update variable d2;j by solving Eq. (21) using the singular value shrinkage algorithm. Step 7. Update variable X k þ 1 by solving Eq. (22). kþ1

Step 8. Update variable b1;j

according to Eq. (18).

kþ1

Step 9. Update variable b2;i according to Eq. (19). Step 10. Loop to Step 2 until a predetermined iteration stopping criterion is satisfied. In real applications, 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 iterative scheme: n o ð23Þ X k þ 1 ¼ Project X k þ 1 where, 8 > Θ1 ;  < Project Z i;j ¼ Z i;j ; > :Θ ; 2

if

Z i;j o Θ1

if

Θ1 rZ i;j r Θ2

if

Z i;j 4 Θ2

ð24Þ

109

7. Numerical simulations and discussions In previous sections, ECT image reconstruction task is cast as an optimization problem, and an iteration scheme is developed for the solving the proposed objective functional. For easy notation, the proposed method can be concisely called as the tensor based volume imaging (TVI) algorithm. In this section, numerical simulations are implemented to evaluate the feasibility and effectiveness of the TVI algorithm and the reconstruction quality is compared with the projected Landweber iteration (PLI) method and the KF method. The robustness of the TVI algorithm to the inaccuracies of the measurement data and the sensitivity matrix is numerically evaluated. A 12 electrodes square ECT sensor is used for simulations and an image is presented using 32  32 pixel. All algorithms are implemented using the MATLAB 7.0 software on a PC with a Pentium IV 2.4 GHz CPU, and 4 Gb memory. The image error is employed to evaluate the quality of an inversion solution, which is defined as γj ¼

‖X Original;j  X Computed;j ‖  100% ‖X Original;j ‖

ð25Þ

Fig. 3. Dynamic reconstruction objects.

Table 1 Algorithmic parameters for the PLI algorithm. Algorithmic parameters

Fig. 3(a)

Fig. 3(b)

Fig. 3(c)

Fig. 3(d)

Fig. 3(e)

Relaxation factor Number of iteration

1 596

1 535

1 748

1 401

1 448

Fig. 4. Reconstructed images by the PLI algorithm.

Fig. 5. Reconstructed images by the TVI algorithm.

110

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

where γ j represents the imager error at the jth time instant; X Original;j and X Computed;j stand for the original and reconstructed values at the jth time instant, respectively. In essence, Eq. (11) belongs to the framework of the Tikhonov regularization method, in which the selection of the regularization parameter plays a crucial role in real applications. Various methods, including the discrepancy principle method [64], the generalized cross-validation method [65], the generalized Arcangeli criteria method [66], the L-curve method [67] and the Engl criteria method [68], are available for determining the regularization parameter, and more discussions on the selection of the regularization parameter can be found in [69]. Owing to the complexity and particularity of the problem, determining an optimal regularization parameter remains a critical problem, especially when the objective functionals include multiple regularizers. In the field of ECT image reconstruction, the regularization parameter is often chosen empirically [7]. For the sake of easy computation, in this paper the regularization parameter empirically determined. However, it is worth mentioning that this issue should be further investigated in the future.

7.1. Case 1 In this section, a dynamic reconstruction case with five noisefree measurement vectors is implemented to evaluate the feasibility of the TVI algorithm and the reconstruction quality is compared with the PLI method. Fig. 3 presents the dynamic permittivity distributions at discrete time instants, in which 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. Table 2 Image errors (%). Algorithms

Fig. 3(a)

Fig. 3(b)

Fig. 3(c)

Fig. 3(d)

Fig. 3(e)

PLI TVI

12.62 3.13

15.77 4.56

18.39 5.62

18.03 7.14

17.50 6.68

For the sake of easy computation and fast reconstruction, in Table 1 the relaxation factors of the PLI method are set as one, and more discussions on the algorithm can be found in [69–71]. For fair comparisons, the algorithmic parameters in Table 1 correspond to the conditions of the minimum image errors. For the sake of easy computation, in the TVI algorithm the regularization parameters are empirically determined, which can be specified as α1 ¼ 1, α2 ¼ 1, α3 ¼ 0:001, α4 ¼ 0:02, α5 ¼ 1, λ1 ¼ 1, βi ¼ 1,

  p W ℓ ;ℓ ¼ 1= X ℓ ;ℓ  þ ε , p ¼ 1, ε ¼ 10  10 , and the number of 1

2

1

2

iterations is 25. The initial values for the PLI method and the TVI algorithm are calculated by the STR method. The images reconstructed by the PLI method and the TVI algorithm are shown in Figs. 4 and 5, respectively. The image errors are listed in Table 2. The images reconstructed by the PLI method are presented in Fig. 4. Numerical simulation results indicate that distinct advantages of the PLI method include easy numerical implementation and low computational complexity and cost owing to the fact that only gradient information of the objective functional is used, and the final inversion solution is not sensitive to the selection of the initial solutions. It is worth mentioning that the PLI method has found successful applications in the field of ECT image reconstruction, including static and dynamic reconstruction tasks. Like other static reconstruction algorithms, such as the TSVD method, the ART algorithm, the GVSPM technique, the SIRT method and the OIOR technique, the PLI algorithm fails to utilize the dynamic evolution information and the spatial–temporal correlations of a dynamic object, and thus the improvement of the reconstruction quality may be restricted. Furthermore, in the PLI algorithm the images are independently reconstructed one by one, which may not be an optimal selection for a dynamic object. In fact, it can be observed from Fig. 4 and Table 2 that for the dynamic reconstruction cases simulated in this section, the quality of the images reconstructed by the PLI method is far from perfect, the reconstructed images are not in a good agreement with original reconstruction objects, and the image errors for original dynamic reconstruction objects, Fig. 3(a)–(e), are 12.62%, 15.77%, 18.39%, 18.03% and 17.50%. The images reconstructed by the TVI algorithm are illustrated in Fig. 5. Numerical simulation results confirm that the TVI

Fig. 6. Reconstructed image by the TVI algorithm when the noise level is 5%.

Fig. 7. Reconstructed image by the TVI algorithm when the noise level is 16%.

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

algorithm can ensure the numerical stability of an inversion solution by introducing the Tikhonov regularization technique to the objective functional. Under the condition of the multiple measurement vectors, the temporal and spatial constraints are better considered and the ECT measurement information and the dynamic evolution information of a dynamic object are integrated in the TVI algorithm, which facilitates the improvement of the reconstruction quality. Additionally, it is worth emphasizing that in the TVI algorithm the time–space image volume is cast as a third-order tensor, the low rank properties of a two-dimensional difference image and the unfoldings of a tensor are simultaneously utilized, the spatial–temporal properties of a dynamic object are exploited, and thus the reconstruction quality can be improved. As can be expected, it can be seen from Fig. 5 and Table 2 that the quality of the images reconstructed by the TVI algorithm is improved as compared with the PLI method, and the reconstructed images are in a good agreement with the original reconstruction objects, which indicates that the TVI algorithm is successful in solving ECT inverse problems. Acquiring high-quality images is highly desired for real applications. It can be found from Table 2 that for the cases simulated in this section, the image errors of the TVI algorithm for original dynamic reconstruction objects, Fig. 3(a)–(e), are 3.13%, 4.56%, 5.62%, 7.14% and 6.68%, which is smaller than that of the PLI method, 12.62%, 15.77%, 18.39%, 18.03% and 17.50%. These results confirm that the TVI algorithm is a promising candidate for solving ECT image reconstruction problems. In the TVI algorithm, the size of the time window should be in advance determined. In order to ensure the reconstruction speed, the size of the time window should not be too large. In fact, from the viewpoint of the temporal correlations, too large time window may not be suitable for a dynamic object owing to the fact that the inter-frame correlations may be weak when the interval of the Table 3 Image errors (%). Noise levels (%)

Fig. 3(a)

Fig. 3(b)

Fig. 3(c)

Fig. 3(d)

Fig. 3(e)

5 16

4.03 5.37

4.94 14.79

8.08 14.58

7.36 17.69

9.57 13.77

111

measurement instants is relatively large. In the case of twodimensional image reconstruction tasks with the M  N pixels, empirically selecting the size of the time window as M=4–M=2 may be appropriate. After a time window is determined, the time window is slided along the time dimension to achieve image reconstruction. 7.2. Case 2 The final inversion solution is sensitive to the inaccurate properties on the input data owing to the ill-posed nature in ECT image reconstruction. In practice, measurement noises are ubiquitous and complicated. From the viewpoint of applications, a successful algorithm should be robust to the capacitance measurement noises. In this section, the noise-contaminated capacitance data is used to evaluate the robustness of the TVI algorithm, and the noise level is defined as: ηj ¼

‖Y Original;j  Y Contaminated;j ‖  100% ‖Y Original;j ‖

ð26Þ

where ηj represents the noise level at the jth instant; Y Original;j and Y Contaminated;j stand for the original and contaminated capacitance data at the jth time instant, respectively; Y Contaminated;j ¼ Y Original;j þ σω, σ represents the standard deviation and ω stands for a normally distributed random number with mean 0 and standard deviation 1. The algorithmic parameters for the TVI algorithm are the same as Section 7.1. The initial values for the TVI algorithm are solved by the STR method. Figs. 6 and 7 are the images reconstructed by the TVI algorithm under the noise levels of 5% and 16%, respectively. The image errors are illustrated in Table 3. Figs. 6 and 7 are the images reconstructed by the TVI algorithm when the noise levels of the capacitance data are 5% and 16%, respectively. Owing to the fact that the inaccurate properties on the sensitivity matrix, the reconstruction model and the measurement data are considered in the proposed objective functional, as can be expected, it can be found from Figs. 6 and 7 that the TVI algorithm shows satisfactory robustness, and the quality of the images reconstructed by the TVI algorithm under different noise levels is satisfactory. Particularly, it can be seen from Table 3 that when the noise level is 16%, the image errors are 5.37%, 14.97%, 14.58%, 17.69%

Fig. 8. Reconstructed images when the standard deviation of the sensitivity matrix and the noise level of the capacitance data are 0.001 and 6%.

Fig. 9. Reconstructed images when the standard deviation of the sensitivity matrix and the noise level of the capacitance data are 0.004 and 8%.

112

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

and 13.77%, which indicates that the TVI algorithm is successful in treating with the inaccuracy of the measurement data.

and thus it can be concluded that the TVI algorithm is robust to the inaccuracies of the sensitivity matrix and the measurement data. 7.4. Case 4

7.3. Case 3 In order to further evaluate the robustness of the TVI algorithm, the influences of the inaccuracies of the sensitivity matrix and the measurement data on the reconstruction quality are studied in this section, and the inaccuracy of the sensitivity matrix is defined as: S contaminated ¼ S original þE

ð27Þ

where E ¼ σ Urandn; S original and S contaminated stand for the original and noise-contaminated sensitivity matrices, respectively. In this case, algorithmic parameters for the TVI algorithm are the same as Section 7.1. The initial values for the TVI algorithm are computed by the STR method. Figs. 8 and 9 illustrate the images reconstructed by the TVI algorithm when the inaccurate properties on the sensitivity matrix and the measurement data, and the image errors are presented in Table 4. Figs. 8 and 9 illustrate the images reconstructed by the TVI algorithm when the inaccurate properties of the sensitivity matrix and the measurement data are considered. It can be observed from Figs. 8 and 9 that owing to the fact that the inaccurate properties on the sensitivity matrix, the measurement data and the reconstruction model are emphasized in the proposed objective functional, the quality of the images reconstructed by the TVI algorithm is satisfactory. Especially, it can be seen from Table 4 that when the standard deviation of the sensitivity matrix is 0.004 and the noise level of the capacitance data is 8%, the image errors for the original dynamic reconstruction objects are 10.63%, 11.87%, 17.53%, 17.28% and 17.65%,

In this section, another dynamic reconstruction case is implemented to evaluate the feasibility of the TVI algorithm, and the reconstruction quality is compared with the KF method. Fig. 10 presents the dynamic permittivity distributions at different time instants, in which 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. The algorithmic parameters of the TVI method are the same as Section 7.1, and the standard deviation of the capacitance measurement noises is 10  3. The initial values for the KF method and the TVI algorithm are calculated by the STR method, and the random-walk evolution model is employed. The images reconstructed by the KF method and the TVI algorithm are shown in Figs. 11 and 12, respectively. The image errors for the KF method and the TVI algorithm are listed in Table 5. The images reconstructed by the KF method and the TVI algorithm are presented in Figs. 11 and 12, respectively. Numerical simulation results indicate that the distinct advantages of the KF method include easy numerical implementation and low computational complexity and cost. It is worth mentioning that the dynamic evolution information and the ECT measurement information are incorporated in the KF method, which is appropriate for dynamic reconstruction problems. In the past years, the KF method has found numerous successful applications in various fields. For the cases simulated in this section, however, it can be observed from Fig. 11 and Table 5 that the reconstruction quality of the KF method is not satisfactory, and the reconstructed images are blurred owing to the smoothness effect and the ill-posed nature in ECT image

Table 4 Image errors (%). Noise levels (%)

Standard deviations

Fig. 3(a)

Fig. 3(b)

Fig. 3(c)

Fig. 3(d)

Fig. 3(e)

6 8

σ ¼ 0:001 σ ¼ 0:004

5.07 10.63

6.34 11.87

10.85 17.53

9.96 17.28

12.61 17.65

Fig. 10. Dynamic reconstruction objects.

Fig. 11. Reconstructed images by the KF method.

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

113

Fig. 12. Reconstructed images by the TVI algorithm.

Table 5 Image errors (%). Algorithms

Fig. 10(a)

Fig. 10(b)

Fig. 10(c)

Fig. 10(d)

Fig. 10(e)

KF TVI

25.19 11.82

25.51 12.35

18.83 14.16

24.25 12.13

19.13 7.63

reconstruction problems. The image errors for original dynamic reconstruction objects, Fig. 10(a)–(e), are 25.19%, 25.51%, 18.83%, 24.25% and 19.13%. Additionally, it can be observed from Fig. 12 and Table 5 that the quality of the images reconstructed by the TVI algorithm is improved as compared with the KF method, which indicates that the TVI algorithm is successful in solving ECT image reconstruction problems. Acquiring high-quality images is crucial for real applications of ECT visualization measurement method. The simulation results presented in previous sections indicate that the TVI algorithm is a promising candidate for the dynamic reconstruction tasks. Particularly, it is worth emphasizing that this paper provides a general framework for solving dynamic image reconstruction problems in the field of ECT, which may be useful for the image reconstruction problems in other related fields.

8. Conclusions Reconstructing high-quality images is crucial for real applications of ECT visualization measurement method. In this paper, a robust dynamic reconstruction model that incorporates the ECT measurement information and the dynamic evolution information of a dynamic object is proposed, in which a series of dynamic images is cast as a third-order tensor that the first two dimensions are space and the third is time. Under the considerations of the two-dimensional spatial structure property of a difference image matrix and the spatial–temporal property of an image tensor, a new objective functional that fuses the ECT measurement information, the dynamic evolution information, the temporal constraint, the spatial constraint, the low rank constraint of a difference image and the low n-rank constraint of a third-order tensor is proposed, where the images are reconstructed by a batching pattern. The SBI algorithm is used to solve the proposed objective functional. Numerical simulations are implemented to validate the feasibility of the proposed algorithm. For the cases compared in this paper, the reconstruction quality and the robustness are improved, which indicates that the proposed algorithm is successful in solving ECT image reconstruction problems. Additionally, it is worth mentioning that this paper provides a general framework for solving ECT inverse problems, which may be applied to the solving of inverse problems in other related fields. In real applications, the improvement of the reconstruction quality depends mainly on the measurement conditions, the design of the sensor, the calibration and filter of the measurement

data, the image reconstruction algorithms, the understanding of the measurement objects, the exploitations of the prior information, etc. Our work provides an alternative approach for solving ECT image reconstruction problems, which needs to be further validated by more cases in the future and be further studied on the issues such as the adaptive selection of the algorithmic parameters and the efficient implementation of the numerical algorithm.

Acknowledgements The authors wish to thank the National Natural Science Foundation of China (nos. 51206048 and 51276059) and the Fundamental Research Funds for the Central Universities (no. 13MS11) for supporting this research. References [1] Xie CG, Huang SM, Hoyle BS, Thorn R, Lenn C, Snowden D, et al. Electrical capacitance for flow imaging: system model for development of image reconstruction algorithms and design of primary sensors. IEE Proc: G Circuits Devices Syst 1992;139:89–98. [2] Tikhonov AN, Arsenin VY. Solution of Ill-posed problems. New York: V.H. Winston & Sons; 1977. [3] Landweber L. An iteration formula for fredholm integral equations of the first kind. Am J Math 1951;73:615–24. [4] Yang WQ, Spink DM, York TA, McCann H. An image reconstruction algorithm based on Landweber’s iteration method for electrical capacitance tomography. Meas Sci Technol 1999;10:1065–9. [5] Jang JD, Lee SH, Kim KY, Choi BY. Modified iterative Landweber method in electrical capacitance tomography. Meas Sci Technol 2006;17:1909–17. [6] Liu S, Fu L, Yang WQ, Wang HG, Jiang F. Prior-online iteration for image reconstruction with electrical capacitance tomography. IEE Proc: Sci Meas Technol 2004;151:195–200. [7] Yang WQ, Peng LH. Image reconstruction algorithms for electrical capacitance tomography. Meas Sci Technol 2003;14:R1–13. [8] Su BL, Zhang YH, Peng LH, Yao DY, Zhang BF. The use of simultaneous iterative reconstruction technique for electrical capacitance tomography. Chem Eng J 2000;77:37–41. [9] Mou CH, Peng LH, Yao DY, Xiao DY. Image reconstruction using a genetic algorithm for electrical capacitance tomography. Tsinghua Sci Technol 2005;10:587–92. [10] Takei M. GVSPM image reconstruction for capacitance CT images of particles in a vertical pipe and comparison with the conventional method. Meas Sci Technol 2006;17:2104–12. [11] Soleimani M, Lionheart WRB. Nonlinear image reconstruction for electrical capacitance tomography using experimental data. Meas Sci Technol 2005;16:1987–96. [12] Wang HX, Tang L, Cao Z. An image reconstruction algorithm based on total variation with adaptive mesh refinement for ECT. Flow Meas Instrum 2007;18:262–7. [13] Fang WF. A nonlinear image reconstruction algorithm for electrical capacitance tomography. Meas Sci Technol 2004;15:2124–32. [14] Wang HM, Hu HL, Wang LJ, Wang HX. Image reconstruction for an electrical capacitance tomography (ECT) system based on a least squares support vector machine and bacterial colony chemotaxis algorithm. Flow Meas Instrum 2012;27:59–66. [15] Ortiz-Aleman C, Martin R, Gamio JC. Reconstruction of permittivity images from capacitance tomography data by using very fast simulated annealing. Meas Sci Technol 2004;15:1382–90. [16] Warsito W, Fan LS. Neural network based multi-criterion optimization image reconstruction technique for imaging two-and three-phase flow systems using electrical capacitance tomography. Meas Sci Technol 2001;12:2198–210.

114

J. Lei et al. / Flow Measurement and Instrumentation 41 (2015) 104–114

[17] Banasiak R, Soleimani M. Shape based reconstruction of experimental data in 3D electrical capacitance tomography. NDT & E Int 2010;43:241–9. [18] Cao Z, Xu LJ, Fan WR, Wang HX. Electrical capacitance tomography for sensors of square cross sections using Calderon’s method. IEEE Trans Instrum Meas 2011;60:900–7. [19] Li Y, Yang WQ. Image reconstruction by nonlinear Landweber iteration for complicated distributions. Meas Sci Technol 2008;19:1–8. [20] Cao Z, Xu LJ, Wang HX. Electrical capacitance tomography with a non-circular sensor using the dbar method. Meas Sci Technol 2010;21:1–6. [21] Watzenig D, Brandner M, Steiner G. A particle filter approach for tomographic imaging based on different state-space representations. Meas Sci Technol 2007;18:30–40. [22] Soleimani M, Vauhkonen M, Yang WQ, Peyton A, Kim BS, Ma XD. Dynamic imaging in electrical capacitance tomography and electromagnetic induction tomography using a Kalman filter. Meas Sci Technol 2007;18:3287–94. [23] Soleimani M, Mitchell CN, Banasiak R, Wajman R, Adler A. Four-dimensional electrical capacitance tomography imaging using experimental data. Prog Electromagnet Res 2009;90:171–86. [24] Lei J, Liu S, Wang XY. Dynamic inversion in electrical capacitance tomography using the ensemble Kalman filter. IET Sci Meas Technol 2012;6:63–77. [25] Lei J, Liu S. Dynamic inversion approach for electrical capacitance tomography. IEEE Trans Instrum Meas 2013;62:3035–49. [26] Jähne B. Spatio–temporal image processing: theory and scientific applications. New York, NY: Springer-Verlag; 1993. [27] Wexler Y, Shechtman E, Irani M. Space–time completion of video. IEEE Trans Pattern Anal Mach Intell 2007;29:463–76. [28] Shechtman E, Caspi Y, Irani M. Space–time super-resolution. IEEE Trans Pattern Anal Mach Intell 2005;27:531–45. [29] Nguyen NH, Tran TD. Robust lasso with missing and grossly corrupted observations. IEEE Trans Inf Theory 2013;59:2036–58. [30] Nissinen A, Kolehmainen V, Kaipio JP. Reconstruction of domain boundary and conductivity in electrical impedance tomography using the approximation error approach. Int J Uncertainty Quant 2011;1:203–22. [31] Kaipio JP, Somersalo E. Statistical and computational inverse problems. Springer-Verlag; 2005. [32] Lanczos C. Linear differential operators. New York, NY: Van Nostrand; 1961. [33] Kirsch A. An introduction to the mathematical inverse problems. New York, NY: Springer-Verlag; 1996. [34] Kolda TG, Bader BW. Tensor decompositions and applications. SIAM Rev 2009;51:455–500. [35] Tan H, Feng J, Feng G, Wang W, Zhang YJ. Traffic volume data outlier recovery via tensor model. Math Prob Eng 2013;2013:1–9. [36] Liu J, Musialski P, Wonka P, Ye J, Tensor completion for estimating missing values in visual data. In: IEEE 12th international conference on computer vision. Kyoto, Japan; September 27–October 4 2009. p. 2114–2121. [37] Liu J, Msialski P, Wonka P, Ye J. Tensor completion for estimating missing values in visual data. IEEE Trans Pattern Anal Mach Intell 2013;35:208–20. [38] Gandy S, Recht B, Yamada I. Tensor completion and low-n-rank tensor recovery via convex optimization. Inverse Prob 2011;27:1–20. [39] Tan H, Cheng B, Feng J, Feng G, Wang W, Zhang YJ. Low-n-rank tensor recovery based on multi-linear augmented Langrange multiplier method. Neurocomputing 2013;119:144–52. [40] Yang L, Huang ZH, Shi X. A fixed point iterative method for low n-rank tensor pursuit. IEEE Trans Signal Process 2013;61:2952–62. [41] Majumdar A. Improved dynamic MRI reconstruction by exploiting sparsity and rank-deficiency. Magn Reson Imaging 2013;31:789–95. [42] Candes EJ, Li X, Ma Y, Wright J. Robust principle component analysis? J ACM 2011;58:1–37. [43] Wright J, Peng Y, Ma Y, Ganesh A, Rao S. Robust principal component analysis: exact recovery of corrupted low-rank matrices via convex optimization. In: Bengio Y, Schuurmans D, Lafferty J, Williams C, editors. Advances in neural information processing systems; 2009. p. 2080–8.

[44] Liu YY, Jiao LC, Shang FH. An efficient matrix factorization based low-rank representation for subspace clustering. Pattern Recognit 2013;46:284–92. [45] Zheng YG, Zhang XR, Yang SY, Jiao LC. Low-rank representation with local constraint for graph construction. Neurocomputing 2013;122:398–405. [46] Liu G, Lin Z, Yan S, Sun J, Yu Y, Ma Y. Robust recovery of subspace structures by low-rank representation. IEEE Trans Pattern Anal Mach Intell 2013;35:171–84. [47] Hu Y, Lingala SG, Jacob M. A fast majorize–minimize algorithm for the recovery of sparse and low-rank matrices. IEEE Trans Image Process 2012;21:742–53. [48] Goldstein T, Osher S. The split Bregman method for L1-regularized problems. SIAM J Imag Sci 2009;2:323–43. [49] Yin W, Osher S, Goldfarb D, Darbon J. Bregman iterative algorithms for L1 minimization with applications to compressed sensing. SIAM J Imag Sci 2008;1:143–68. [50] Cai JF, Osher S, Shen Z. Split Bregman methods and frame based image restoration. Multiscale Model Simul 2009;8:337–69. [51] Cai JF, Osher S, Shen Z. Linearized Bregman iterations for compressed sensing. Math Comput 2009;78:1515–36. [52] Zhang X, Burger M, Osher S. A unified primal–dual algorithm framework based on Bregman iteration. J Sci Comput 2011;46:20–46. [53] Goldstein T, Bresson X, Osher S. Geometric applications of the split Bregman method: segmentation and surface reconstruction. J Sci Comput 2010;45: 272–93. [54] Osher S, Mao Y, Dong B, Yin W. Fast linearized Bregman iteration for compressive sensing and sparse denoising. Commun Math Sci 2010;8:93–111. [55] Gao H, Cai JF, Shen ZW, Zhao HK. Robust principal component analysis-based four-dimensional computed tomography. Phys Med Biol 2011;56:3181–98. [56] Chan TF, Wong CK. Convergence of the alternating minimization algorithm for blind deconvolution. Linear Algebra Appl 2000;316:259–85. [57] Xiao Y, Zeng T, Yu J, Ng MK. Restoration of images corrupted by mixed Gaussian-impulse noise via L1–L0 minimization. Pattern Recognit 2011;44: 1708–20. [58] Lv XG, Song YZ, Wang SX, Le J. Image restoration with a high-order total variation minimization method. Appl Math Modell 2013;37:8210–24. [59] Blondel M, Seki K, Uehara K. Block coordinate descent algorithms for largescale sparse multiclass classification. Mach Learn 2013;93:31–52. [60] Qin Z, Scheinberg K, Goldfarb D. Efficient block-coordinate descent algorithms for the group lasso. Math Program Comput 2013;5:143–69. [61] Tseng P. Convergence of a block coordinate descent method for nondifferentiable minimization. J Optim Theory Appl 2001;109:475–94. [62] Beck A, Tebouule M. A fast iteration shrinkage-thresholding algorithm for linear inverse problems. SIAM J Imag Sci 2009;2:183–202. [63] Cai JF, Candes EJ, Shen Z. A singular value thresholding algorithm for matrix completion. SIAM J Optim 2008;20:1956–82. [64] Morozov VA. Methods for solving incorrectly posed problems. New York, NY: Springer-Verlag; 1984. [65] Golub G, Heath M, Wahba G. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 1979;21:215–23. [66] Engl HW. Discrepancy principle for Tikhonov regularization of ill-posed problems leading to optimal convergence rates. J Optim Theory Appl 1987;52:209–15. [67] Hansen PC. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev 1992;34:561–80. [68] Engl HW, Hanke M, Neubauer A. Regularization of inverse problems. Dordrecht: Kluwer Academic Publishing; 1996. [69] Wang YF. Computational methods for inverse problems and their applications. Beijing: Higher Education Press; 2007. [70] Xiao TY, Yu SG, Wang YF. Numerical methods of inverse problems. Beijing: Science Press; 2003. [71] Liu S, Fu L, Yang WQ. Optimization of an iterative image reconstruction algorithm for electrical capacitance tomography. Meas Sci Technol 1999;10: L37–L39.