Statistical Inversion Based on Nonlinear Weighted Anisotropic Total Variational Model and Its Application in Electrical Impedance Tomography

Pengfei Qi

School of Mathematics and Statistics, Shandong Normal University, Jinan, China.

DOI: 10.4236/eng.2024.161001                

Abstract

Electrical impedance tomography (EIT) aims to reconstruct the conductivity distribution using the boundary measured voltage potential. Traditional regularization based method would suffer from error propagation due to the iteration process. The statistical inverse problem method uses statistical inference to estimate unknown parameters. In this article, we develop a nonlinear weighted anisotropic total variation (NWATV) prior density function based on the recently proposed NWATV regularization method. We calculate the corresponding posterior density function, i.e., the solution of the EIT inverse problem in the statistical sense, via a modified Markov chain Monte Carlo (MCMC) sampling. We do numerical experiment to validate the proposed approach.

Keywords

Statistical Inverse Problem, Electrical Impedance Tomography, NWATV Prior, Markov Chain Monte Carlo Sampling

1. Introduction

Electrical impedance tomography (EIT) [1] is an imaging modality that aims to reconstruct the conductivity distribution by injecting a current into the body through pairs of electrodes attached to the surface of the target and measuring the voltage data. EIT for medical imaging has the advantages of being noninvasive, real-time monitoring, portable, and low cost.

However, traditional regularization based method would suffer from error propagation due to the iteration process. The statistical inverse problem method uses statistical inference to estimate unknown parameters. Meanwhile, the posterior density function is sampled in solving to avoid the error propagation and accumulation caused by iteration.

The statistical inverse problem method for EIT is proposed in [2] . In [2] the resistivity sampling based on Markov chain Monte Carlo (MCMC) is given. It is shown in reference [3] that the Bayesian conditional mean (CM) estimation of the prior distribution of total variation (TV) cannot preserve the edge.

In reference [4] , the statistical inverse problem of the EIT is realized by selecting the TV regularization term as the prior density function, and the visualization image is obtained using the MCMC algorithm. In the literature [5] , the TV prior density function and MCMC sampling algorithm are also used to achieve 3D EIT image reconstruction. However, as in the non-statistical method, staircasing artifacts occur.

In Song [6] , the authors proposed having an edge-preserving effect nonlinear weighted anisotropic total variation (NWATV) regularization method to solve the EIT inverse problem. The reconstruction quality and edge-preserving effect are improved.

In this article, we develop a NWATV prior density function based on the NWATV regularization method. We estimate the corresponding posterior density function, i.e., the solution of the EIT inverse problem in the statistical sense, via a modified MCMC sampling method. Numerical results show that the proposed NWATV prior density function is feasible, and the edge-preserving effect is improved.

2. NWATV Priors and Posterior Distribution of EIT Problems

2.1. NWATV Priors Density Function

In Song [6] , the authors proposed having an edge-preserving effect NWATV regularization term ${\Vert \nabla \sigma /{\left|\nabla \sigma \right|}^2\Vert }_{L^1\left(\Omega \right)}$ , here $\sigma$ is conductivity. The conductivity is then reconstructed by solving the following minimization problem:

$\delta \sigma =\underset{\delta \sigma }{\arg \min }\left\{\frac{1}{2}{\Vert S\delta \sigma -\delta V\Vert }^2+\alpha {\Vert p\cdot D\left(\delta \sigma \right)\Vert }_{l_1}\right\}\mathrm{,}$ (1)

here, $\delta \sigma$ is the change of conductivity, $S$ is the sensitivity matrix, $p=\left(\zeta \left(\delta \sigma \right)\mathrm{<mo>\; ;\; </mo>}\zeta \left(\delta \sigma \right)\right)\in R^{2N}$ , $\zeta \left(r\right)\mathrm{<mo>\; :\; </mo>}=1/{\left|\nabla \delta \sigma \left(r\right)\right|}^2$ , $r=\left(x\mathrm{,}y\mathrm{,}z\right)$ , $D\left(\delta \sigma \right)=\left(D_x\delta \sigma \mathrm{<mo>\; ;\; </mo>}{D}_y\delta \sigma \right)\in R^{2N}$ , where $D_x\mathrm{,}{D}_y\in R^{N\times N}$ are respectively the first-order difference operators along the x and y directions.

Next, we will analyze the prior density function corresponding to NWATV regularization in the sense of statistical inference. Precisely, we consider the following form of prior density:

${\pi }_{\text{pr}}\left(\sigma \right)\propto {\pi }_{+}\left(\sigma \right){\tilde{\pi }}_{\text{pr}}\left(\sigma \right)\mathrm{,}$ (2)

where ${\tilde{\pi }}_{\text{pr}}\left(\sigma \right)$ is a regularization prior, ${\pi }_{+}\left(\sigma \right)$ is the positivity prior

${\pi }_{+}\left(\sigma \right)=\{\begin{array}{cc}1& \text{if}0\le \sigma \left(k\right)\le \infty \text{for}\text{all}k=\mathrm{1,2,}\cdots \mathrm{,}N\mathrm{,}\\ 0& \text{otherwise}\mathrm{.}\end{array}$

Here, N represents the number of components discretized by the finite element method for the conductivity distribution in the reconstructed region.

Considering the prior density function corresponding to NWATV regularization

$\text{NWATV}\left(\sigma \right)={\Vert p\cdot D\left(\sigma \right)\Vert }_{l_1}\mathrm{,}$ (3)

we obtain

${\tilde{\pi }}_{\text{pr}}\left(\sigma \right)\propto {\text{e}}^{-\alpha {\Vert p\cdot D\left(\sigma \right)\Vert }_{l_1}}\mathrm{.}$ (4)

Therefore, a nonlinear weighted anisotropic total variational prior density function is obtained

${\pi }_{\text{pr}}\left(\sigma \right)\propto {\pi }_{+}\left(\sigma \right){\text{e}}^{-\alpha {\Vert p\cdot D\left(\sigma \right)\Vert }_{l_1}}\mathrm{,}$ (5)

where $\alpha >0$ is a parameter related to the confidence of the regularized priors. When considering the change of conductivity $\delta \sigma$ , (5) can be rewritten as follows

${\pi }_{\text{pr}}\left(\delta \sigma \right)\propto {\text{e}}^{-\alpha {\Vert p\cdot D\left(\delta \sigma \right)\Vert }_{l_1}}\mathrm{.}$ (6)

2.2. The Solution of EIT Inverse Problem in Statistical Sense

Now, consider the EIT problem as a concrete form of the statistical inverse problem. The EIT statistical model can be expressed as follows

$\delta V=S\delta \sigma +W\mathrm{,}$ (7)

where, $\delta V$ represents the observed quantity, $\delta \sigma$ represents the change of conductivity and $W$ represents noise. We assume that $\delta \sigma$ and $W$ are independent.

Suppose that the noise vector $W$ is a Gaussian random vector with zero mean, and the covariance matrix $\Gamma$ is positive definite, i.e.,

$W~\mathcal{N}\left(\mathrm{0,}\Gamma \right)\mathrm{.}$ (8)

The likelihood density function is obtained as

$\pi \left(\delta v\mathrm{<mo>\; |\; </mo>}\delta \sigma \right)\propto \exp \left(-\frac{1}{2}{\left(S\delta \sigma -\delta V\right)}^{\text{T}}{\Gamma }^{-1}\left(S\delta \sigma -\delta V\right)\right)\mathrm{.}$ (9)

According to the Bayesian formula, the posterior density function of conductivity is

$\pi \left(\delta \sigma \mathrm{<mo>\; |\; </mo>}\delta v\right)\propto \exp \left(-\frac{1}{2}{\left(S\delta \sigma -\delta V\right)}^{\text{T}}{\Gamma }^{-1}\left(S\delta \sigma -\delta V\right)-\alpha {\Vert p\cdot D\left(\delta \sigma \right)\Vert }_{l_1}\right)\mathrm{.}$ (10)

In the experimental application, we assume that the noise covariance matrix is $\Gamma ={\beta }^{2}I$ . We obtain

$\pi \left(\delta \sigma \mathrm{<mo>\; |\; </mo>}\delta v\right)\propto \exp \left(-\frac{1}{2{\beta }^2}{\Vert S\delta \sigma -\delta V\Vert }_{l_2}^2-\alpha {\Vert p\cdot D\left(\delta \sigma \right)\Vert }_{l_1}\right)\mathrm{.}$ (11)

Next, we will estimate the posterior density function using the CM method. The formula for calculating the CM method is as follows

$\delta {\sigma }_{\text{CM}}={\displaystyle {\int }_{{ℝ}^n}}\delta \sigma \pi \left(\delta \sigma \mathrm{<mo>\; |\; </mo>}\delta v\right)\text{d}\delta \sigma \mathrm{.}$ (12)

Since an image contains a large number of pixels and hence the above integral dimension is huge. We use the MCMC algorithm to sample the posterior density function (11) to obtain the sample sequence $\delta {\sigma }^{\left(1\right)}\mathrm{,}\delta {\sigma }^{\left(2\right)}\mathrm{,}\cdots \mathrm{,}\delta {\sigma }^{\left(M\right)}$ .

The burn-in period of the probe posterior density function is assumed to have a sample size of $m_0$ . When the total sample quantity M is sufficiently large, the remaining quantity after removing the burn-in period samples is given by $M-m_0$ . Therefore, the above integral (12) can be approximated by the average of $\left(M-m_0\right)$ samples, i.e.,

$\delta {\sigma }_{\text{CM}}={\displaystyle {\int }_{{ℝ}^n}}\delta \sigma \pi \left(\delta \sigma \mathrm{<mo>\; |\; </mo>}\delta v\right)\text{d}\delta \sigma \approx \frac{1}{M-m_0}{\displaystyle \underset{m=m_0+1}{\overset{M}{\sum }}}\delta {\sigma }^{\left(m\right)}\mathrm{.}$ (13)

Specifically, we probed and sampled the posterior density function (11) using the Metropolis-Hastings [7] [8] method from MCMC. Thus, we obtained a distribution Q on ${ℝ}^n$ similar to the posterior density function (11).

Fix $1\le N_{\text{update}}\le n$ and $\kappa >0$ . Generate an alternative value $\delta \sigma \in {ℝ}^n$ from Q, satisfying $\delta \sigma =\delta {\sigma }^{\left(m\right)}+\epsilon$ , where

$\epsilon ={\left[\mathrm{0,}\cdots \mathrm{,0,}{{\epsilon }^{\prime }}_1\mathrm{,0,}\cdots \mathrm{,0,}{{\epsilon }^{\prime }}_2\mathrm{,0,}\cdots \mathrm{,0,}{{\epsilon }^{\prime }}_{N_{\text{update}}}\mathrm{,0,}\cdots \mathrm{,0}\right]}^{\text{T}},$

with ${{\epsilon }^{\prime }}_l~\mathcal{N}\left(\mathrm{0,}\kappa \right)$ , $l=\mathrm{1,2,}\cdots \mathrm{,}{N}_{\text{update}}$ . In addition, for the EIT inverse problem, it is generally possible to determine the conductivity range of the object. We also incorporate this prior information into our reconstruction algorithm. The sampling algorithm is as follows

Step 1: Set $m:=0$ , B, $N_{\text{update}}$ and initialize $\delta {\sigma }^{\left(0\right)}$ by e.g. $\delta {\sigma }^{\left(0\right)}\mathrm{<mo>\; :\; </mo>}={\left[\mathrm{0,}\cdots \mathrm{,0}\right]}^{\text{T}}$ .

Step 2: Set $\delta \sigma \mathrm{<mo>\; :\; </mo>}=\delta {\sigma }^{\left(m\right)}+\epsilon$ , where $\epsilon ={\left[\mathrm{0,}\cdots \mathrm{,0,}{{\epsilon }^{\prime }}_1\mathrm{,0,}\cdots \mathrm{,0,}{{\epsilon }^{\prime }}_2\mathrm{,0,}\cdots \mathrm{,0,}{{\epsilon }^{\prime }}_{N_{\text{update}}}\mathrm{,0,}\cdots \mathrm{,0}\right]}^{\text{T}}$ , with ${{\epsilon }^{\prime }}_l~\mathcal{N}\left(\mathrm{0,}\kappa \right)$ independent random numbers.

Step 3: If $\delta \sigma \left(i\right)<B$ , then perform n steps to retransfer again, such that $\delta \sigma \left(i\right)=\delta {\sigma }^{\left(m\right)}\left(i\right)+w\ge B$ , $i=1,\cdots ,N$ , $w~\mathcal{N}\left(\mathrm{0,}{\gamma }^2/2^n\right)$ , $n=\mathrm{1,2,3,}\cdots$ , where B is the background pixel value of the image.

Step 4: If $\pi \left(\delta \sigma \mathrm{<mo>\; |\; </mo>}\delta v\right)\ge \pi \left(\delta {\sigma }^{\left(m\right)}\mathrm{<mo>\; |\; </mo>}\delta v\right)$ then set $\delta {\sigma }^{\left(m+1\right)}=\delta \sigma$ and go to Step 5.

Step 5: Draw a random number s from the uniform distribution on $\left[\mathrm{0,1}\right]$ . If $s\le \frac{\pi \left(\delta \sigma \mathrm{<mo>\; |\; </mo>}\delta v\right)}{\pi \left(\delta {\sigma }^{\left(m\right)}\mathrm{<mo>\; |\; </mo>}\delta v\right)}$ then set $\delta {\sigma }^{\left(m+1\right)}\mathrm{<mo>\; :\; </mo>}=\delta \sigma$ ; else set $\delta {\sigma }^{\left(m+1\right)}\mathrm{<mo>\; :\; </mo>}=\delta {\sigma }^{\left(m\right)}$ .

Step 6: If $m=M$ then stop; else set $m\to m+1$ and go to Step 2.

3. Result

We used the EIDORS [9] software package to solve the forward problem. Experiments were conducted using adjacent injection current and voltage measurement modes with a maximum injection current of 1 mA. In the sampling process, we maintain the acceptance rate in the range 0.25 - 0.35.

We used the sample data to calculate and visualize the upper and lower bounds for 90% credibility interval of the imaging results to evaluate the reconstruction effect. The credibility interval is defined as follows

$\left[CL\left(i\right)\mathrm{,}CU\left(i\right)\right]=\left[\delta \hat{\sigma }\left(i\right)-a\gamma \left(i\right)\mathrm{,}\delta \hat{\sigma }\left(i\right)+a\gamma \left(i\right)\right]\text{for}i=\mathrm{1,}\cdots \mathrm{,}n\mathrm{.}$ (14)

Here CL represents the lower credibility bound, CU represents the upper credibility bound, a is the credibility coefficient (In the experiment, the value is 1.645.), $\delta \hat{\sigma }\left(i\right)$ and $\gamma \left(i\right)$ are the mean and variance of the sample data for the i-th component after fitting the probability distribution.

To verify the feasibility of the theory and obtain the visualization image, we used a modified MCMC algorithm to perform a 2D numerical experiment of EIT. The parameters of the numerical experiment were set to $B=0$ , $M=1\times {10}^6$ , $m_0=5\times {10}^5$ , $N_{\text{update}}=5$ , $\kappa =5\times {10}^{-3}$ , $\alpha =1\times {10}^{-12}$ , $\beta =1\times {10}^{-3}$ . The reconstruction results are shown in Figure 1. We also visualized the upper and lower bounds for 90% credibility interval to evaluate the effect of the CM estimate reconstruction, as shown in Figure 2. As a comparison, we also used the Tikhonov regularization [10] method to reconstruct the conductivity distribution, and the results are shown in Figure 1.

The experimental results show that the CM method is better than the Tikhonov regularization method regarding conductivity reconstruction quality and boundary preservation. In addition, we obtain more information about the solution, such as credibility intervals.

Figure 1. Reconstruction of the conductivity distribution of circular objects. The first line is the real conductivity distribution image, the second line is the Tikhonov reconstruction of the conductivity distribution, and the third line is the CM estimate of the conductivity distribution.

Figure 2. The upper and lower bounds for 90% credibility interval. The first line is the lower credibility bound, and the second is the upper credibility bound.

4. Conclusions

In this article, the traditional NWATV regularization term is extended to the NWATV prior density function, and the corresponding posterior density function is obtained. The modified MCMC sampling algorithm was used to sample the posterior density function, and the visualization image was obtained from the CM estimation. Meanwhile, the edge-preserving effect is improved.

Compared with the traditional Tikhonov iterative algorithm, the statistical sampling algorithm avoids error propagation and accumulation in the iterative process. Moreover, for the statistical inverse problem of EIT, except for obtaining a numerical solution, we obtain the posterior density function, which contains rich information about the solution.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Cheney, M., Isaacson, D. and Newell, J.C. (1999) Electrical Impedance Tomography. SIAM Review, 41, 85-101. https://doi.org/10.1137/S0036144598333613
[2]	Kaipio, J.P., Kolehmainen, V., Somersalo, E. and Vauhkonen, M. (2000) Statistical Inversion and Monte Carlo Sampling Methods in Electrical Impedance Tomography. Inverse Problems, 16, 1487-1522. https://doi.org/10.1088/0266-5611/16/5/321
[3]	Lassas, M. and Siltanen, S. (2004) Can One Use Total Variation Prior for Edge- Preserving Bayesian Inversion? Inverse Problems, 20, 1537-1563. https://doi.org/10.1088/0266-5611/20/5/013
[4]	Somersalo, E., Kaipio, J.P., Vauhkonen, M., Baroudi, D. and Jarvenpaa, S. (1997) Impedance Imaging and Markov Chain Monte Carlo Methods. Proc SPIE’s 42nd Annual Meeting, Computational, Experimental and Numerical Methods for Solving Illposed Inverse Imaging Problems: Medical and Nonmedical Applications, San Diego, USA, June 27-August 1, 1997, 175-185.
[5]	Kolehmainen, V., Somersalo, E., Vauhkonen, P.J., Vauhkonen, M. and Kaipio, J.P. (1998) A Bayesian Approach and Total Variation Priors in 3D Electrical Impedance Tomography. Proceedings of the 20th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Hong Kong, 1028-1031.
[6]	Song, Y., Wang, Y. and Liu, D. (2022) A Nonlinear Weighted Anisotropic Total Variation Regularization for Electrical Impedance Tomography. IEEE Transactions on Instrumentation and Measurement, 71, 1-13. https://doi.org/10.1109/TIM.2022.3220288
[7]	Gilks, W.R., Richardson, S. and Spiegelhalter, D.J. (1996) Markov Chain Monte Carlo in Practice. Chapman and Hall, London. https://doi.org/10.1201/b14835
[8]	Hastings, W.K. (1970) Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 57, 97-109. https://doi.org/10.1093/biomet/57.1.97
[9]	Adler, A. and Lionheart, W.R. (2006) Uses and Abuses of EIDORS: An Extensible Software Base for EIT. Physiological Measurement, 27, S25-S42. https://doi.org/10.1088/0967-3334/27/5/S03
[10]	Vauhkonen, M., Vadász, D., Karjalainen, P.A., Somersalo, E. and Kaipio, J.P. (1998) Tikhonov Regularization and Prior Information in Electrical Impedance Tomography. IEEE Transactions on Medical Imaging, 17, 285-293. https://doi.org/10.1109/42.700740