Comsol Electrical Impedance Tomography

Shengheng Liu

Abstract— Electrical impedance tomography (EIT) is developed to investigate the internal conductivity changes of an object through a series of boundary electrodes, and has become increasingly attractive in a broad spectrum of applications. However, the design of optimal tomogra- phy image reconstruction algorithms has not achieved the adequate level of progress and matureness. In this paper, we propose an efficient and high-resolution EIT image reconstruction method in the framework of sparse Bayesian learning. Significant performance improvement is achieved by imposing structure-aware priors on the learning process to incorporate the prior knowledge that practical conduc- tivity distribution maps exhibit clustered sparsity and intra- cluster continuity. The proposed method not only achieves high-resolution estimation and preserves the shape infor- mation even in low signal-to-noise ratio scenarios but also avoids the time-consuming parameter tuning process. The effectivenessof the proposedalgorithmis validatedthrough comparisons with state-of-the-art techniques using exten- sive numerical simulation and phantom experiment results.

Index Terms— Inverse problem, electrical impedance tomography (EIT), sparse Bayesian learning (SBL), image reconstruction, maximum a posteriori (MAP) estimation.

Continuous

real-time visualization of the dynamic electric conductivity distribution of the interior of a body. To perform EIT, we apply weak low-frequency alternating currents (the typical frequency range is 1–100 kHz and the magnitude is 1–5 mA) in multiple manners and measure the corresponding peripheral voltages Manuscript received January 18, 2018; revised March 7, 2018; accepted March 12, 2018. Date of publication March 29, 2018; date of current version August 30, 2018. The work of S. Liu, J. Jia, and Y. Yang was supported by the U.K. EPSRC under Grant EP/P006833/1.

on Imaging Systems and Techniques (IST 2017), Beijing, China .

(Corresponding Author: Jiabin Jia.)

S. Liu, J. Jia, and Y. Yang are with the Agile Tomography Group, School

Pa 19122 Usa.

Color versions of one or more of the figures in this paper are available

Digital Object Identifier 10.1109/Tmi.2018.2816739

through an array of electrodes attached on the surface of an object. Then, the cross-sectional image of the internal spatial conductivity distribution is recovered from the resulting boundary voltage measurements. Despite its relatively low spatial resolution (∼10% of the sensor diameter) that hinders its general applicability, EIT has become a well-accepted tomographic imaging technique because of its desirable prop- erties of portability, low cost, no known hazards, and high temporal resolution in comparison with other techniques avail- able for probing internal dynamics. As a result, EIT finds broad applications in a number of fields, such as industrial

Process Monitoring –, Geophysical Exploration –,

and biomedical diagnosis –. Mathematically, the reconstruction of conductivity maps in EIT amounts to solving a nonlinear ill-posed inverse problem from noisy data. Regularization techniques can be employed to mitigate the instability of the solutions. One of the most widely used families is the one-step Gauss–Newton (GN) reconstruction approach , which allows the use of sophis- ticated regularized models to describe the EIT inverse problem through a heuristically determined prior , . Landweber iteration is a variation of the steepest gradient descent method

And Is Also Widely Used In Eit , . The Algebraic

reconstruction technique (ART) is a simple and effective image reconstruction method for computerized tomography that can be applied to EIT , . Other important methods include regularization via the total variation (TV) functional , , which allows image reconstruction with edge preservation.

As EIT solutions generally manifest themselves as sparse vectors, a number of sparse regularization approaches have been proposed to stabilize the inversion. For example, in –, an ℓ1 penalty was incorporated into the reg- ularization to promote sparsity. It is shown in that, by clustering relevant signal components together and tak- ing the dependencies between them into account, a superior performance is achieved through the exploitation of structure- based approaches in various practical scenarios. On this basis, sparse signal recovery algorithms have been developed in the literature to exploit the cluster structure of the signals.

These methods generally require certain knowledge about the signals, such as the size and distribution of the partition, and the number of nonzero entries of the signal to be recovered.

LIU et al.: IMAGE RECONSTRUCTION IN EIT BASED ON SA-SBL

2091

In many applications, however, such knowledge is unavail- able. In , an adaptive group sparsity (AGS) constraint is proposed to obtain enhanced image quality, where the required a priori knowledge is first obtained through a coarse image reconstruction, and images with a finer quality are then obtained.

A natural mechanism to incorporate the structure knowl- edge in solving the inverse problem is to employ Bayesian approaches , , which are aimed to characterize the pos- terior distribution, e.g., computing posterior moments or other posterior expectations. Sampling via Markov chain Monte Carlo (MCMC) techniques is among the most commonly used paradigms for Bayesian inversion. Although various efforts have been made to accelerate the MCMC approach by exploiting, for example, the gradient and Hessian infor- mation of the posterior density , or multilevel sampling strategies , generating sufficient MCMC samples is still computationally costly . Alternatively, we can use the fil- tering and variational methods, which provide approximations of the exact posterior by finding the best representative from a class of simpler distributions. Readers are referred to and the references therein for an overview.

The past few years have witnessed extensive research and a significant progress on Bayesian inversion in terms, e.g., of the design of correlation priors and the use of hierarchical models.

To name a few, in , a Bayesian hierarchical model with conditionally Gaussian priors was used to resolve piecewise smooth structures in noisy and blurred images. In , a uni- fied and efficient Bayesian method based on Krylov subspace iterative techniques and preconditioners was further derived to compute the approximated maximum a posteriori (MAP) estimate of the image and the prior variance. The algorithm proposed in was later successfully applied to cerebral source localization , where the variances are assumed to be guided by a hyperprior with the form of generalized Gamma family, so as to favor small values while permitting rare large outliers which correspond to high source amplitudes.

Besov Space Priors Defined By A Wavelet Expansion

with random coefficients have attracted broad attention in the medical imaging community due to their merits such as edge- preserving, sparsity-promoting, and discretization invariant, with imaging results similar to TV regularization –.

In , the smoothness of unknown targets is modelled by Whittle-Matérn fields, which are essentially stationary Gaussian random fields specified by a correlation function to control the anisotropic properties of the prior distribution.

A finite-dimensional approximation of the Whittle-Matérn pri- ors is derived from sparse inverse covariance matrices by using a stochastic partial differential equation. The problem consid- ered in is to reconstruct the conductivity consisting of well-defined inclusion-type targets in an approximately homo- geneous background. The proposed iterative algorithm is based on the use of a nonlinear edge-preferring prior density and the minimization of the corresponding Tikhonov functional by efficiently solving an approximate sequence of linearized problems with the help of prior-conditioning and least squares with QR factorization (LSQR). In , the inverse problem of magnetoencephalography (MEG), which is very similar to EIT, is solved by postulating a hierarchical conditionally Gaussian prior model, where an anatomical prior is introduced to reflect the direction preference of each dipole based on the a priori magnetic resonance imaging (MRI) information.

The hyperparameter vector consists of prior variances of the dipole moments, which are assumed to follow a non- conjugate gamma distribution with variable scaling and shape parameters. By combining the iterative alternating sequential algorithm with Krylov subspace iterative solver, satisfactory sparsity control and convergence rate are both achieved.

The recently proposed sparse Bayesian learning (SBL)

Framework –, Which Is Closely Related To The Term

automatic relevant determination (ARD) in the context of neural network, can adaptively and flexibly explore and exploit signal structures such as clustering or continuity without any a priori information. In addition, SBL is more advantageous than other families of aforementioned algorithms in the sense that it is more robust in noisy environments, and offers better performance when the columns of the dictionary matrix are highly correlated and/or the image to be reconstructed is not highly sparse . In this paper, therefore, we introduce the concept of structure-aware sparse Bayesian learning (SA-SBL) to perform EIT imaging, and both the clustered sparsity and intra-cluster correlation are utilized in the image reconstruction to achieve improved reconstruction accuracy. Compared with the state-of-the-art methods, the proposed approach is advan- tageous because its structure-aware modeling capability pro- motes clustered sparsity and eliminates irrelevant components.

In addition, the proposed method yields a better approximation to the ℓ0-norm sparsity measure and, as a result, achieves enhanced EIT imaging with a higher spatial resolution and improved robustness against Gaussian noise.

The paper is organized as follows. The forward and inverse models are briefly described in Section II. In Section III, we first propose the SA-SBL algorithm for EIT image recon- struction. In Section IV, numerical and real data experiments are designed to simulate the challenging medical application scenarios, and the performance of the proposed algorithms is compared with other state-of-the-art approaches. The paper is concluded in Section V.

Notations: Lower-case (/upper-case) bold characters are used to denote vectors (/matrices).

¯(·) And |·| Respectively

return the average of a given vector and the modulus of a given complex number. ∇is the Nabla symbol. diag{A} returns a column vector consisting of the main diagonal entries, whereas and diag{A, 1} returns one corresponding to the first-diagonal entries above the main diagonal. IN denotes an N × N identity matrix. tr(·) and (·)T respectively represent the trace and transpose operation of a matrix. ∥·∥p represents the ℓp-norm of a vector. E(·) returns the expected value of a discrete random variable. p(·) denotes the probability density function. N(·) denotes Gaussian distribution. R is the set of real numbers.

A. Forward Model

In EIT, the computation of the voltage measurements from the known currents and conductivity distribution is referred to

2092

as the forward problem. The EIT forward model is mathemat- ically established from a low-frequency approximation of a subset of the Maxwell’s electromagnetic equations with some mixed Dirichlet/Neumann boundary conditions . Consider a bounded domain  with piecewise smooth boundary ∂.

Let σ denote the real-valued conductivity distribution in . Assume that the boundary electrodes are perfectly conductive. Then, applying a current flux through these electrodes results in the scalar interior electric potential u characterized as the

(1)

where we assume that the static excitation condition is satisfied and no current sources/sinks exist within . The complete electrode model (CEM) is the most accurate model for EIT since it accounts for the effects of the electrodes and the contact impedances between the electrolyte and the electrodes . Assume that a total of L electrodes are attached to the boundary ∂, whose locations are denoted as ∂l ⊂∂, l = 1, 2, . . . , L. Denote their contact impedances as zl. Then, the boundary conditions of the CEM is expressed

(4)

where⃗e is the outer normal vector on ∂, and Il denotes the injected current at the l-th electrode.

B. Inverse Model

The inverse problem in EIT is to determine the conductivity in the interior of a object based on the simultaneous measure- ments of alternating currents and voltages at the boundary.

In this paper, linearized difference imaging is considered. By solving the Laplacian elliptic partial differential equation given in (1), we can readily relate the map of internal conduc- tivity perturbation δσ ∈RN×1 within the region of interest to its corresponding boundary voltage variation δv ∈RM×1

(5)

where n is the measurement noise vector. In this paper, we assume that the noise vector is additive and follows the zero-mean Gaussian distribution with n ∼N (0, γ0I). Note that, in practice, EIT systems are susceptible to other types of noise and modeling errors. Sources of such noise and errors include multiplicative noise and bias due to amplifier gain distortions and offsets, inter-channel measurement noise due to cross-talk between measurement channels, linearization errors when reference conductive distribution does not match the reference voltages used to compute the Jacobian, and errors in electrode position, domain shape/size and FEM discretization level. However, given the space limitation, it is difficult to examine different types of noise and errors in this work.

As a result, we follow the common practice in EIT image reconstruction and assume the additive noise model in (5) with Gaussian distributed entries (e.g., , –). Still, the above-mentioned annoying noise and errors may co-exist in the real-data experiments in Section IV-B. As shown later, the proposed method performs better in such cases.

The sensitivity map J ∈RM×N, also frequently referred to as the Jacobian matrix, is determined by mesh, electrode posi- tions, and current injection/measurement protocol. By integrat- ing over the k-th (k = 1, 2, . . . , N) simplex k from the scalar product of the gradients of the potential fields u that are induced by driving current pattern Id and measurement pattern Im with the finite element method (FEM), J can be

(6)

where d = m ∈{1, 2, . . ., M} represents the sample index of selected electrode pairs in the m-th measurement under the d-th current pattern. For notational convenience, in the follow- ing discussion, we simplify the notations δσ and δv as σ and v, respectively. In this case, (5) is simplified as v = Jσ + n.

(7)

This is a typical example of inverse problems, where the task is to reconstruct a signal from observations that are subject to the forward model. Like many other inverse problems, EIT image reconstruction is severely ill-posed and thus is difficult to obtain stable and reliable solutions. Therefore, it is often formulated as the following regularized optimization

(8)

where d (v, Jσ) is the data fidelity term that enforces σ to satisfy the observations v, and R(σ) is a stable operator, or a regularization term which restricts the solution to comply with a predefined model over σ. The minimization in (8) represents a tradeoff between fitting the data exactly and stabilizing the solution of σ, which is controlled by the regularization parameter δ. The major problem in this framework is that it is difficult to determine a proper R(σ) such that it represents a most appropriate model for EIT image reconstruction and, at the same time, its optimization procedure can be easily implemented.

In this work, equation (8) is interpreted in a Bayes perspec- tive such that finding the solution amounts to performing the MAP estimation. In this case, d (v, Jσ) takes the form of a log- likelihood whereas R(σ) takes the form of a parametric log- prior log p(σ; ) over variable σ. As such, the minimization

Min

σ∈RN×1 log p(v |σ ) + δ log p(σ; ).

(9)

This allows for more adequately modeling the a priori knowl- edge about the objective. Compared to the conventional regu- larization approach based on functional analysis, the Bayesian probabilistic approach provides a more attractive and rational framework for the parameter recovery, as it allows for the LIU et al.: IMAGE RECONSTRUCTION IN EIT BASED ON SA-SBL

2093

quantification of the uncertainty in the recovery, while natu- rally accommodating different types of data and rich models of a priori information.

Iii. Image Reconstruction Based On Sa-Sbl

In medical imaging applications, generally, the underlying conductivity distribution change σ can be viewed as a con- catenation of several clusters due to its inherent clustering structure. In practice, however, the actual clustering parti- tion pattern is more likely to be unknown. In this context, we consider a general case in which all the clusters overlap each other with an equal size h and the nonzero entries are arbitrarily distributed. It was proved in that, since real clustering partition can be learned during the SBL process by revoking and merging the preset clusters, the reconstruction performance is guaranteed to approach the result obtained from the known clustering partition counterparts. To facilitate the utilization of SA-SBL framework, we factorize σ as

Of

clusters.

∈Rn×H. The

underlying linear model in (7) can then be rewritten as

∈

RM×h. We assume that the prior of the weights x follows a zero-

(12)

where the stretched covariance matrix is expressed as

0

...

⎤

⎥⎦∈Rgh×gh.

(13)

The structured sparsity in x is determined by γi, i = 1, . . . , g, in 0. If γi = 0, the values of the entries in the associated i-th cluster of x becomes zero. By considering the fact that n ∼N (0, γ0I), the a posteriori belief for the weights x is subject to the following Gaussian distribution

(16)

with v ≜γ0I + 0T. The MAP probability estimate of x is obtained from the posterior mean μx, prior to which, the hyperparameters must be estimated first. The expectation-maximization (EM) method is employed to maximize p (v; ), which is equivalent to minimizing −log p (v; ). Thus we obtain the following

Dx

= log |v| + vTv−1v.

(17)

The correlation structure matrix Bi can be updated using the gradient of the cost function with respect to Bi. To avoid the overfitting problem, we introduce Bi as an intermediate

:

ih, (i −1)h + 1 : ih). The estimation of Bi is then obtained by constraining it to the following Toeplitz form,

I

... ... ...

...

Diag(Bi)

.

(21)

Note that the averaging of ri in the above computation effectively prevents against overfitting. An computationally efficient learning rule for the hyperpa- rameter γi can be derived by minimizing the cost function via the majorization-minimization approach , , which can

, (22)

where β ∈[0, 1] models the pattern coupling between the hyperparameter γi of cluster Bi and the hyperparameters {γi+, γi−} of its neighboring clusters, and subscripts i+ and i−respectively indicate the neighboring clusters of the i-th cluster with larger and smaller indices. The mapping from vector σ to the actual geometry and the neighboring clusters are illustrated in Fig. 1. By introducing parameter β into the learning process, the relevance of the neighboring clus- ters, which are not adjacent to the underlying cluster within vector σ, is also exploited. Note that, while the update formu- lae provide good numerical behavior in all the experiments, there is no theoretical analysis on the stability of updating rules of Bi and γi. Nevertheless, as EIT has the merit of high time resolution, we can discard the abnormal frames.

2094

Fig. 1. Illustration of the structure in the conductivity distribution σ. By treating x as a latent parameter, and setting the derivative of the cost function over γ0 to 0, the learning rule for γ0 is

. (23)

A pseudo-code implementation of the proposed algorithm for EIT image reconstruction is provided in Algorithm 1. To update the parameters through the EM approach, the com- putational complexity of the proposed method is O(M2gh) for each iteration. It is worth emphasizing that, the time- consuming parameter tuning process is averted in the proposed algorithm. The concept of pattern coupling was studied in , where the recovery performance is shown to be insensitive to the choice of β as long as β > 0. It was also proved that taking different values of h will lead to negligible difference .

For simplicity, we set β = 0.25 and h = 4 throughout this paper. The other two input parameters min and ϑmax are related to the algorithm precision and are selected according to the required performance and affordable computational complexity. As such, none of them requires a time-consuming parameter tuning process.

Remarks: The proposed hierarchical Bayesian model allows the estimation of prior parameters in an unsupervised manner, i.e., the proposed algorithm does not require any information regarding either the sparsity or the clustering prior. Rather, improved performance is achieved by adaptively and flexibly exploring and exploiting such signal structures.

Iv. Results And Discussions

In this section, we evaluate the performance of the proposed SA-SBL based method through extensive numerical simula- tions and phantom experiments. To demonstrate the superiority of the proposed method, comparisons with existing state-of- the-art approaches, such as ART , TV regularization , ℓ1 regularization , Nissinen’s Bayesian method , and AGS constraint method , are provided. Note that since modeling errors are not considered in this work, the term “Nissinen’s Bayesian method” here refers to the conven- tional MAP estimate without the approximation error method Algorithm 1 EIT Image Reconstruction Based on SA-SBL

: V, J, H, Β, Min, Θmax

Initialize : Set = 1, ϑ = 0, μx = 0gh×1, x = 0gh×gh,

|Vi −¯V|2,

Bi = Toeplitz([0.90, . . . , 0.9h−1]).

8

ϑ = ϑ + 1.

: Σ = Μx

(see [49, Sec. 2.1] or [50, Sec. 2.2.3]). Note that in terms of the spatial resolution of the reconstructed images, exist- ing Bayesian inversion based methods do not show signifi- cant improvement when compared to traditional determinis- tic approaches. For instance, the reconstructed results using Bayesian inversion with Whittle-Matérn priors are comparable to the solution of the gradient-based smoothness methods .

It was also pointed out in that, the TV and Perona- Malik regularization schemes can both be seen as MAP estimations of Bayesian hierarchical models proposed therein.

As discussed in and will be presented later, Nissinen’s Bayesian method demonstrates similar performance with the Tikhonov regularization.

A. Synthetic Data Experiments

We present several synthetic data experiments in this sub- section to demonstrate the effectiveness of the proposed algo- rithm in comparison with other state-of-the-art approaches.

In the COMSOL Multiphysics® environment, we design a 16-electrode EIT sensor as shown in Fig. 2(a). The diameter of the sensor is set as 95 mm. As the neighboring bipolar pattern is adopted, the degrees of freedom of the measurement is M = 16 ×(16 −3)/2 = 104 (dimension of the voltage data vector).

Based on the sensor configuration, the sensitivity map J can be computed. The scaled summation of each row of J is plotted in Fig. 2(b). The sensitivity map is essentially the superposition of all the sensitivity distributions correspond- ing to each measurement, which is nonuniform as shown in Fig. 2(b). The matrix J has the highest sensitivity in the pix- els near the measurement surface, and becomes less sensitive in the center area. In the linearized algorithm, the sensitivity map directly affects the quality of the reconstructed image, as the objects near the excitations will be greatly augmented.

This will be demonstrated in the following experiments. The mutual coherence , (J), of the matrix J is defined as the LIU et al.: IMAGE RECONSTRUCTION IN EIT BASED ON SA-SBL

2095

Fig. 2. Synthetic data experimental settings. (a) Schematic illustration of the designed EIT sensor system in the COMSOL environment.

(b) Normalized sensitivity map of the EIT sensor configuration. maximum absolute value of the inner product between any

∥Ji∥2∥Jk∥2

.

(24)

The mutual coherence (J) is an important measure in the field of sparse representations that characterizes the spread of the columns of matrix J. (J) takes values within the range of 0 ≤ ≤1, with a lower value indicating a large spread and a low coherence. For the data being simulated, the mutual coherence of the sensitivity matrix is (J) = 0.9995, which indicates a high coherence between the columns of matrix J.

Under this condition, the SBL algorithm in has been proved effective. The SA-SBL algorithm proposed in the paper is devised through modifications of the SBL algorithm in and provides reliable image reconstruction even when the sensitivity matrix J is highly coherent.

In this paper, we assume that the voltage data is corrupted

2

.

(25)

In this subsection, since we have the ground truth of the conductivity distribution for the synthetic data, the correlation coefficient and the relative reconstruction error can be defined respectively as follows to quantitatively evaluate the recon-

(27)

where ς represents the true conductivity variation, and the number of simplexes in the inverse FEM mesh is N = 3228 in our phantom test. Note that to make the simulation more realistic and to avoid committing an inverse crime, we utilize a forward mesh consisting of 6570 domain elements and 304 boundary elements. The adopted forward and inverse FEM meshes are shown in Fig. 3. It can be observed from Fig. 3 Fig. 3.

Adopted forward and inverse FEM meshes. (a) Forward mesh consisting of 6570 domain elements and 304 boundary elements. (b) Inverse mesh consisting of 3228 elements.

Fig. 4. Ground truth of the generated multiple small objects. that, the inverse mesh consists of square simplices, while the forward mesh consists of triangular simplices. We construct the EIT forward model using the dense mesh, where the electric field potential of each triangular simplex is derived from the values of its corresponding vertices. On the other hand, the conductivity of each square simplex is determined by the conductivity at the center of the square simplex. Thus, the inverse mesh used when solving the inverse problem is totally different from the forward mesh used when generating the measurement data.

First, an illustrative example is presented to evaluate the performance of the proposed method when resolving multiple small objects. As shown in Fig. 4, we design a phantom consisting of three small titanium beta-21S squares with a conductivity value of 7.407 × 105 S/m and a side length of 10 mm, and the conductivity of the background saline is 0.05 S/m. Other simulation settings include an excitation frequency of 10 kHz and a peak-to-peak amplitude of 1.5 mA.

The iteration termination conditions for the iterative methods are set to min = 1 × 10−5 and ϑmax = 200. The target objective parameter in the ℓ1 regularization is set as 0.1. The iteration step of TV regularization is set as 0.01. For AGS constraint method, the maximum diameter of the group is set



. All the above simulation settings remains unchanged unless otherwise stated.

The reconstruction results are given in Table I, where the plots in each row show the results for a different method, whereas each column depicts the results corresponding to different input SNR values. Concretely, in the first three rows of Table I, the boundary between the two closer squares are not clearly resolved. However, compared with the other two

Related Journal Articles & DOI Links

Selected peer-reviewed publications relevant to 12 Lead ECG Acquisition. Click the DOI to access the full paper (may require institutional access).

Why Choose Us?

Bangalore guidance for robotics, Spectre and autonomous systems projects.