## Abstract

Electrical capacitance tomography (ECT) is a non-destructive detection technique for imaging the permittivity distributions inside an observed domain from the capacitances measurements on its boundary. Owing to its advantages of non-contact, non-radiation, high speed and low cost, ECT is promising in the measurements of many industrial or biological processes. However, in the practical industrial or biological systems, a deposit is normally seen in the inner wall of its pipe or vessel. As the actual region of interest (ROI) of ECT is surrounded by the deposit layer, the capacitance measurements become weakly sensitive to the permittivity perturbation occurring at the ROI. When there is a major permittivity difference between the deposit and the ROI, this kind of shielding effect is significant, and the permittivity reconstruction becomes challenging. To deal with the issue, an interface and permittivity simultaneous reconstruction approach is proposed. Both the permittivity at the ROI and the geometry of the deposit layer are recovered using the block coordinate descent method. The boundary and finite-elements coupling method is employed to improve the computational efficiency. The performance of the proposed method is evaluated with the simulation tests.

This article is part of the themed issue ‘Supersensing through industrial process tomography’.

## 1. Introduction

Electrical capacitance tomography (ECT) is a relatively new non-destructive detection technique. In ECT, a set of electrodes is attached at the periphery of an observed domain. The mutual capacitances between the electrodes are measured and used for reconstructing the permittivity distributions inside the observed domain. Owing to its advantages of non-contact, non-radiation, high-temporal resolution and low cost, ECT has been widely researched, and is promising in the measurements of many industrial and biological processes, e.g. oil transportation [1], pneumatic conveyors [2] and fluidized beds [3]. However, the permittivity reconstruction problem in ECT is nonlinear and highly ill-posed. Its spatial resolution is relative low.

To deal with the problem, many reconstruction algorithms have been proposed. The most commonly used algorithms are pixel based [4]. The permittivity inside the region of interest (ROI) of ECT is discretized into finite-volume elements, then reconstructed from the measured capacitances using a Tikhonov regularized least-squares formulation. This kind of algorithm is simple to implement, and is resistant to the measurement noise and model errors. However, the Tikhonov regularization method assumes the target permittivity is smoothing. When the target permittivity is piecewise constant, the smoothing assumption will conflict with the permittivity jumps occurring at the boundaries of the piecewise regions and blur the reconstructed images. Although the edge-enhanced images can be obtained from the iterative linearization method [5], the directly nonlinear method [6,7], the limited region method [8] or the boundary penalty operators [9,10], additional post-processing methods are required before applying the pixel-based images for quantitative analysis.

To directly obtain the quantitative results, some shape-based reconstruction algorithms are proposed. In these methods, the boundary of the target inclusion is parametrized by the analytical functions or level sets. The optimal boundary shape is estimated by minimizing the difference between the measured and predicted capacitances. Comparing with the pixel/voxel-based methods, the shape-based approaches mainly have two advantages. First, the prior information of the discontinuity of the permittivity distribution, such as the permittivity values of the target inclusions, can be explicitly incorporated into the shape reconstruction process. Secondly, as the inclusion boundary is reconstructed, the reconstruction results can be directly used for the quantitative analysis. However, as a kind of iterative linearization method, the shape-based methods are sensitive to measurement noise and model errors. Examples of the shape-based reconstruction methods are boundary element method (BEM)-based methods [11,12], the level-set methods [13,14], the integral equation method [15] and the boundary distributed source method [16]. Furthermore, shape-based approaches are also useful in the reconstruction of the outer boundary of the observed domain [17–20].

In the paper, we focus on the application of ECT on the reconstruction of permittivity at the ROI surrounded by the deposit layer. Deposit on the inner wall of its pipe or vessel is normally seen in industrial or biological systems, such as the oil transportation pipeline [21], heat exchanger network [22] and the gas transmission system [23]. Existence of the deposit layer will not only decrease the production efficiency, but will also pose some negative influences on the system safety. The associated clearance maintenance is always extremely expensive, and should be scheduled with the knowledge of the deposit thickness. Furthermore, as the actual ROI of ECT is surrounded by the deposit, the ECT measurements becomes weakly sensitive to permittivity changes occurring at the ROI. This kind of shielding effect will be significant, when there is a major permittivity contrast between the deposit layer and the ROI. As a result, the reconstruction of the permittivity at the ROI becomes challenging. Neither the pixel-based nor the shape-based method alone can provide reliable results.

To deal with the problem, an interface and permittivity simultaneous reconstruction approach is proposed. The deposit layer is assumed to be occupied by a homogeneous permittivity, while the actual ROI of ECT is occupied by an inhomogeneous permittivity. The permittivity distribution at the ROI is parametrized by the pixel-based approach, while the geometry of the deposit layer is parametrized by the shape-based approach. The block coordinate method is used to estimate the unknown permittivity and shape parameters from the boundary capacitance measurements. The boundary and finite-elements coupling method is used to simplify the numerical implementation of the simultaneous reconstruction approach. The simulation results show that the proposed method can mitigate the shielding effect of the deposit layer on the ROI, and significantly improve the permittivity reconstruction results. The remainder of the paper is organized as follow: in §2a, the measurement principles and governing equations of ECT are briefly described, in §2b, the target interface and permittivity distribution are represented by a set of parameters and in §2c, the unknown parameters are estimated by using the block coordinate descent method; in §3, the boundary and finite-elements coupling method in ECT is described; in §4, several numerical tests are conducted to evaluate the performance of the proposed method.

## 2. Methodology

### (a) Governing equations

Figure 1 gives a schematic of the capacitance tomography problem considered in the paper. The entire observed domain *Ω* consists of the electrode mounting layer *Ω*_{p}, the deposit layer *Ω*_{d} and the inhomogeneous permittivity section *Ω*_{m} (the ROI of the ECT). A set of electrodes *E* is mounted around the periphery of the observed domain. Each time, one of the electrodes is stimulated by a voltage signal (used as the drive electrode), while the others are at the virtual earth potential (used as the detecting electrodes). Mutual capacitances between the drive and detecting electrodes are sequentially measured by a measuring circuit.

Following Maxwell's theory, the electrical potential *u* at the entire observed region *Ω* is governed by the elliptic partial differential equation
2.1where *x*∈*R*^{k} is the *k*-dimensional spatial variable, *ε*∈*C*^{k} is the permittivity distribution. In the paper, we assume the electrode mounting layer and the deposit layer are occupied by the homogeneous permittivity. As the *ε* is constant in these regions, the associated governing equations are simplified to the Laplace equations
2.2

The boundary conditions at the exterior boundary ∂*Ω* are determined by the data collection strategy. During *i*th voltage stimulation, the electric potentials at the electrodes follow the Dirichlet condition
2.3where *δ*(*i*,*j*) is Kronecker delta equalling zero except when *i*=*j*, *V* _{c} is the drive voltage amplitude, and *E*_{j} denotes the *j*th electrode. The electric charge densities at the inter-electrode gaps follow the Neumann condition
2.4where *q*=∂*u*/∂*ν*, *ν* is the outward unit normal to the counterclockwise-oriented boundary, and is the combination of the inter-electrode gaps.

The boundary conditions at the interior interfaces are drawn from the integral formulation of Maxwell's equations
2.5and
2.6where *Γ*_{p}=*Ω*_{p}∩*Ω*_{d} is the interface between the electrode mounting layer and the deposit layer, and *Γ*_{d}=*Ω*_{d}∩*Ω*_{m} is the interface between the deposit layer and the inhomogeneous permittivity section. The subscripts *p*, *d* and *f* stand for the electric quantities calculated at the side of the electrode mounting layer, the deposit layer and the inhomogeneous permittivity section, respectively.

The boundary-value problem described above can be solved by the numerical computation method, such as the finite-element method (FEM). Then, the capacitance measurements are calculated by
2.7where *c*_{i,j} is the mutual capacitance between the drive electrode *E*_{i} and detecting electrode *E*_{j}. To emphasize the contribution of the interested permittivity and interface to the ECT measurements, we denote the capacitance vector with *C*(*Γ*_{d},*ε*_{d},*ε*_{m}). Correspondingly, the forward problem of ECT is to calculate the capacitances from given interface and permittivity configurations, while the inverse problem is to recover the interface and permittivity configurations from given capacitance measurements.

### (b) Permittivity and interface representation

To solve the ECT forward and inverse problems, *ε*_{m} and *Γ*_{d} must be represented by a set of parameters. One of the most popular methods is by using FEM mesh. The objective region *Ω*_{m} is divided into a finite number of small elements ∪△_{k}. The permittivity at *Ω*_{m} is approximated by
2.8where *M* is the total number of discrete elements, *ε*^{k} is a positive real number representing the permittivity value at *k*th element, and *δ*^{k} is a Kronecker delta equaling to zero except at △_{k}. To represent a two-dimensional closed interface, one can use the geometric nodes method [11] or the level sets method [14]. Here, we assume the deposit interface is smooth and star-shaped. One of the powerful methods to parametrize this kind of interfaces is the Fourier shape decomposition method [24], in which the points located at the target interface are defined by
2.9where *r* and *θ* are the radius and angle in polar system, *f*(*k*,*θ*) is the *k*th Fourier basis function, *γ*^{k} is the weighted coefficients and *N* is the total number of shape components. The method can not only significantly reduce the number of unknowns involved in the inverse problems, but also implicitly incorporate some smoothing constraints into the interface. If *N* is large enough, the approximated interface *Γ*_{d}(*γ*) will be close to the target. To simplify the notation, we denote *ε*_{m}=[*ε*^{1},*ε*^{2},…,*ε*^{M}] and *γ*_{d}=[*γ*^{1},*γ*^{2},…,*γ*^{N}], the vectors of the permittivity and shape coefficients, respectively.

### (c) Parameters estimation

The inverse problem considered here is to reconstruct the unknown permittivity and interface parameters from the boundary-capacitance measurements. The most popular solution to the problem is the output least-squares approach, in which the unknown parameters are estimated by minimizing the following cost function
2.10where *C* and *C** are the simulated and measured capacitances, respectively. As the inverse problem of ECT is ill-posed, directly minimizing (2.10) will lead to a unreliable solution. The generalized Tikhonov regularization is incorporated into the objective function
2.11where *L* is a positive definite matrix , and , and are the prior interface or permittivity guesses. Different *L* leads to different penalty manners. If *L* is a identity matrix, it will lead to a standard Tikhonov regularization procedure. If *L* is a diagonal weighted from *J*^{T}*J*, where *J* is the Jacobian matrix of *C*, it will lead to Newton's One-Step Error Reconstructor (NOSER) [25] or Levenberg–Marquardt (LM) method [26].

As it is known, the ECT measurements are more sensitive to the geometry and permittivity changes at the deposit layer than to the permittivity changes at the inhomogeneous permittivity section. Instead of updating all unknowns at the same time, an alternating minimization (Block Coordinate Descent) schedule is used. The unknown parameters are naturally partitioned into three blocks *γ*_{d}, *ε*_{d} and *ε*_{m}. Each block of variables is alternatively updated by solving the regularized least-square problem with respect to one block at a time while all the others are fixed. When the cost function is convex and differentiable, this kind of alternating minimization approach will converge to a reliable solution [27]. To simplify notification, when the permittivity in *Ω*_{m} is homogeneous, we will use (a positive real number) to denote the homogeneous property. The flowchart of the alternating minimization schedule is shown in figure 2.

Assume the permittivity at

*Ω*_{m}is homogeneous, and choose the initial guesses of*γ*_{d},*ε*_{d}and based on prior knowledge.Fix

*ε*_{d}and , recover the unknown interface parameters*γ*_{d}by minimizing the cost function*Ξ*using LM method.Update

*γ*_{d}with the results from the last step, fix*γ*_{d}, and estimate the permittivity values*ε*_{d}and by minimizing*Ξ*using LM method.Update

*ε*_{d}and ; if the capacitance residual or its variance rate is smaller than the tolerance, go to the next step; else, back to the second step.Fix

*γ*_{d}and*ε*_{d}, choose as the background permittivity distribution, and reconstruct*ε*_{m}using NOSER.

It is worth noting that the initial interface guess influences the reconstruction results slightly, while the initial permittivity guesses and influence the results significantly. A large tolerance *τ* in the stop condition will lead to more iterations. The trial-and-error method [28] is a reliable method to choice these parameters. NOSER and LM method are not the only choice in the alternating minimization schedule. They are used here because they are good enough to provide accurate results. When the deposit thickness is small, the constrained least-square approach is efficient.

To solve the nonlinear least-square problems by gradient type method, the Jacobian of *C* with respect to *γ*_{d}, *ε*_{d}, and *ε*_{m} must be calculated. Here, the perturbation method [12] is applied. The systematic details of the method can be found in [29,30]. The derivations of the boundary capacitances with respect to the shape coefficients are calculated by the boundary integral
2.12where *h*(*x*) is the displacement field of *Γ*_{d} due to the interface perturbation, and *u*_{i} and *u*_{j} are the potentials at *i*th and *j*th voltage drives, respectively. Based on (2.9), one further has
2.13where ∇_{τ}*u* denotes the tangential components of the electrical field. It is worth noting that the electrical variables involved in (2.13) are defined at the side of *Ω*_{m}. The Jacobian of the boundary capacitances with respect to discretized permittivity pixels are calculated by
2.14where △_{k} denotes *k*th discrete permittivity element. The Jacobian of the capacitances with respect to the constant permittivity values *ε*_{d} or can be calculated by substituting △_{k} in (2.14) with *Ω*_{d} or *Ω*_{m}. Based on the Gauss divergence theorem, one can obtain a more convenient formula
2.15Comparing with (2.14), the gradient operator is not involved in (2.15). Thus, the computation resource is saved, and it is also easy to implement in the boundary and finite-elements coupling method.

## 3. Numerical implementation

### (a) Boundary elements for homogeneous permittivity sections

In the BEM, the given boundary conditions are used to fit the values at the boundary of the space defined by a partial differential equation [31]. For the Laplace boundary-value problem, the boundary integral equation is derived from Green–Gauss theorem
3.1where ∂*Ω* is the boundary of *Ω*, *α* is the geometric coefficient and is calculated by the internal angle divided by 2*π*, and *x*′ and *x* are the source and filed points, respectively. The *u** and *q** in (3.1) are the fundamental solution for Laplace's equation and its normal derivation, respectively.

To formalize the forward problem of ECT by the BEM, the boundary conditions (2.3) and (2.4) are incorporated into (3.1). The boundary integral equation for the electric field at the electrode mounting layer *Ω*_{p} is governed by
3.2where *κ*_{d,p}=*ε*_{d}/*ε*_{p} is the permittivity contrast between the electrode mounting layer and the deposit layer. The minus before the third boundary integral terms in (3.2) is because all the interfaces are orientated counterclockwise. Similarly, the boundary integral equation at the deposit layer *Ω*_{d} is
3.3where *κ*_{m,d}=*ε*_{m}/*ε*_{d} is the permittivity contrast between the *Ω*_{d} and *Ω*_{m}.

To calculate the boundary integral equations derived from the BEM, the associated boundaries are divided into a finite set of small segments *Γ*=∪_{k}*l*_{k}. These segments are collected together to form a boundary-element mesh. And the potentials and fluxes on the boundaries are approximated by
3.4where *N* is the number of collocation points, *ϕ*^{k} is the nodal basis shape functions, *u*^{k} and *q*^{k} are the potential and flux at *k*th collocation point, respectively. Substituting (3.4) into (3.3), a linear algebraic system can be obtained.

### (b) Finite elements for inhomogeneous permittivity section

The FEM solution of ECT problem starts from the discretization of the objective domain by using the finite-elements mesh [28]. Let us suppose that there are *M* vertexes in the finite-element mesh. The electric potentials and fluxes at the inhomogeneous permittivity section is approximated by
3.5where *ψ*^{k} is a series of shape basis functions depending on the geometric configuration and interpretation modality of the FEM mesh.

The approximated electric quantities in (3.5) does not directly satisfy the differential equation (2.1). Noting that there is no internal electrical source inside *Ω*_{m}, the weak form of (2.1) is defined by
3.6where *v*(*x*) is the testing function. Let the *k*th testing function *v*^{k} equals to *ψ*^{k}, and substitute (3.5) into (3.6). The discrete FEM equation is obtained
3.7It is worth noting that only *q*_{m} at *Γ*_{d} is involved in (3.7). It will simplify the coupling process of the FEM and BEM.

### (c) Coupled system equation

To derive the boundary- and finite-elements coupling formulation, the continuous condition of the potentials and fluxes at *Γ*_{d} should be used. Let us denote the relationship between the electric variables in (3.4) and (3.5) with
3.8where *K* is a linear mapping from the electric variables in the BEM to those in FEM. To obtain a stable solution, *K* should be positive definite. A straightforward way is to make sure the discrete elements in the FEM and BEM at *Γ*_{d} have the same geometry configurations and shape functions. This will lead to an identity *K* and simplify the coding process. In the paper, we use the first-order triangle elements in FEM, and use first-order line elements in the BEM.

As the FEM and BEM systems are already formulated, the boundary and finite-elements coupling system is derived by combining (3.4) and (3.7) together
3.9where *F*, *G* and *H* are the coefficient matrix derived from the BEM, *B* and *A* is the coefficients matrix derived form FEM. By using the coupled algebraic system, the capacitance vector *C* can be computed with the given *γ*_{d}, *ε*_{d} and *ε*_{m}. The shape and permittivity Jacobians are calculated by some vector–matrix multiply operators. As the electrical potential and its derivatives are already involved at the forward solution, the Jacobian calculation is with high efficiency.

## 4. Simulation results and discussions

The proposed method is feasible in two and three dimensions. Considering that ECT sensor length is relatively smaller than pipe length, measurement setting in the numerical experiments is in two dimensions. Twelve electrodes are mounted at the periphery of an insulating pipe. The pipe wall thickness is 5% of the pipe diameter. The electrode mounting layer *Ω*_{p} is with a permittivity of 3.0. The deposit layer *Ω*_{d} is with a permittivity of 4.0. The background permittivity in the inhomogeneous permittivity *Ω*_{m} section is 2.0. The maximum permittivity contrast at *Ω*_{m} is 2:1. Noise components in the capacitance measurements are picked from zeros mean Gaussian random variables. Signal-to-noise ratio is defined by , where *σ*_{c} and *σ*_{n} are the variances of capacitances and noise, receptively. As the forward models are updated in each iteration, inverse crimes are naturally avoided.

### (a) Interface recovery with known permittivity

In the first group of tests, we assume the permittivity at *Ω*_{m} is homogeneous, and there is a exact knowledge of the permittivity values *ε*_{d} and . The inverse algorithm is to recover the deposit interface *Γ*_{d} from the capacitance measurements *C*. Figure 3 shows the interface reconstruction results with different interface shapes. The artificial measurements are contaminated with 20 dB, 30 dB and 40 dB noise, respectively. The target interfaces are controlled by 7 degree Fourier functions. As can be seen, after starting from the initial guesses, the estimated interfaces tend to the targets in a fast and smooth way. All the estimations converge to the reliable results within seven iterations.

The algorithm was run with several different targets. The results are quite similar to those illustrated in figure 3. The residuals (Res) between the measured and estimated capacitances, the Hausdorff distance (Hd) between the actual and recovered interfaces, and the symmetric difference area (Aerr) between the target and estimated deposit sections are used to quantitatively analyse the reconstruction performance. Fifty independent tests are conducted. The mean values (denoted by *_{e}) and the relative standard deviations (denoted by *_{v}) of these indicators are calculated. The results are shown in table 1. With increasing of the measurement noise, the reconstruction performances deteriorate a little. However, even with 20 dB noise, the results are still reliable. The Hausdorff distances are under 0.08 (4% of the diameter of the observed pipe). The area errors are below 0.15 (5% of the cross-sectional area of the observed pipe). Standard deviations of the indicators are quite similar in different noise cases, which means that the interface reconstructions are resistant to noise.

### (b) Interface recovery with incorrect prior permittivity values

In the second group of tests, we further assume the prior knowledge of *ε*_{d} and is incorrect. The inverse algorithm is to recover the *Γ*_{d}, *ε*_{d} and , simultaneously. Our solution starts with one stage of interface recovery, and another stage of permittivity value recovery. Then, the interface and permittivity estimations will be conducted alternately. Several tests are conducted. The artificial measurements are contaminated by 30 dB noise. The relative errors between the initial and actual permittivity values are within 10%.

Table 2 shows six groups of reconstruction results. The test number are shown in the first column. The initial and estimated permittivity values are shown in the next four columns. The quantitative indicators and the total number of iterations are shown in the last four columns. As can be seen, the permittivity value estimations for the deposit layer are accurate. The relative errors reduce from more than 6% to less than 3%. The reconstructed deposit layer interfaces are also reliable. The Hausdorff distances are under 0.04 (2% of the observed domain diameter), and the area errors are below 0.15 (5% of the cross-sectional area of the observed domain). The permittivity value estimations for the inhomogeneous permittivity section are not good. It is because the shielding effect caused by the deposit layer makes the ECT data less sensitive to the permittivity changes occurring at *Ω*_{m}. Figure 4 visualizes the interface reconstructions in table 2. As can be seen, the incorrect initial permittivity values will worsen the interface reconstruction results. One stage of shape recovery without permittivity correcting is always with low accuracy. The two-stage method with the permittivity correcting can improve the reconstruction accuracy.

### (c) Interface and permittivity simultaneous recovery

In the third group of tests, the permittivity at *Ω*_{m} is inhomogeneous. We assume there is a exact knowledge of *ε*_{d} and . The aim of our method is to estimate *Γ*_{d} and *ε*_{m} at the same time. The algorithm starts with one stage of interface recovery with the LM method, then conducts another stage of permittivity distribution estimation with the NOSER method. The artificial measurements are contaminated by 30 dB noise.

Figure 5 visualizes the results from three groups of tests. The first row shows the target interfaces and permittivity distributions. The second row shows the interface and permittivity reconstruction results from the proposed simultaneous reconstruction method (SRM). The third row shows the results from the one-stage permittivity reconstruction method (PRM), in which the permittivity images consist of *Ω*_{d} and *Ω*_{m}. As can be seen, one-step permittivity reconstruction can roughly image the deposit layer, while the reconstruction of the permittivity distribution at *Ω*_{m} is difficult. On the other hand, the proposed interface and permittivity simultaneous reconstruction method have a ability of imaging the geometry of the deposit layer and the permittivity distribution at *Ω*_{m} both in a good accuracy. Comparing the results shown in figures 5 and 3, the inhomogeneous permittivity at *Ω*_{m} only reduces the accuracy of the deposit interface reconstruction slightly.

### (d) Interface and permittivity recovery with incorrect prior knowledge

In the last group of tests, the prior knowledge of *ε*_{d} and is incorrect. The inverse problem is to estimate *Γ*_{d}, *ε*_{d}, and *ε*_{m} simultaneously. Our algorithm starts with alternate estimations of *Γ*_{d}, *ε*_{d} and by the LM method (denoted with electronic supplementary material, S1), then conducts one-step permittivity recovery for *ε*_{m} by using NOSER (see electronic supplementary material, S2). In the tests, the artificial capacitance measurements are contaminated by 30 dB noise, and the relative errors between the initial and target permittivity values are within 10%.

Results from six groups of tests are shown in table 3 and visualized in figure 6. As can be seen, the estimations of *ε*_{d} are reliable, while the estimations of are not very good. The inaccurate permittivity values only sightly influence the interface and permittivity distribution reconstructions in the electronic supplementary material, S2. The Hausdorff distances between the target and reconstructed deposit interfaces are smaller than 0.08 (4% of the diameter of the observed domain), and area errors between the target and estimated deposit layers are less than 0.19 (6% of the cross-sectional area of the observed domain). After a comparison between figures 5 and 6, one can find that the alternative minimization strategy can correct the incorrect initial permittivity guesses, and mitigate their influences on the following interface and permittivity reconstructions.

## 5. Conclusion

In the paper, an interface and permittivity simultaneous reconstruction approach is proposed. The method is based on the boundary and finite-elements coupling method and block coordinate descent algorithm. When the actual ROI of ECT is surrounded by the deposit layer, the approach can mitigate the shielding effect of the high dielectric deposit on the ROI. Furthermore, the numerical tests show that the proposed method can reconstruct the deposit layer profiles with high accuracy, and automatically correct the incorrect initial permittivity guesses. The proposed approach is also useful in the applications of electrical tomography on imaging of biological tissues, such as human brain. As a kind of nonlinear iterative algorithms, the proposed method is sensitive to the model errors. Further research will focus on the elimination of the model error and the application of the method with real data.

## Competing interests

We declare we have no competing interests

## Funding

This study was funded by the National Natural Science Foundation of China, nos. 61401304, 61571321.

## Footnotes

One contribution of 10 to a theme issue ‘Super-sensing through industrial process tomography’.

- Accepted March 7, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.