Detecting stochastic inclusions in electrical impedance tomography

TL;DR

This paper investigates detecting stochastic inclusions in electrical impedance tomography using the Factorization and Monotonicity Methods on mean Neumann-to-Dirichlet maps, demonstrating effectiveness with high-contrast anomalies through theoretical and numerical analysis.

Contribution

It extends existing methods to stochastic conductivities, showing that anomalies can be identified via mean maps under high contrast conditions.

Findings

Successful detection of stochastic inclusions using mean Neumann-to-Dirichlet maps.

Theoretical validation of methods for high-contrast anomalies.

Numerical examples in 2D confirm the approach's effectiveness.

Abstract

This work considers the inclusion detection problem of electrical impedance tomography with stochastic conductivities. It is shown that a conductivity anomaly with a random conductivity can be identified by applying the Factorization Method or the Monotonicity Method to the mean value of the corresponding Neumann-to-Dirichlet map provided that the anomaly has high enough contrast in the sense of expectation. The theoretical results are complemented by numerical examples in two spatial dimensions.

Figures

Equations

$$H_{\diamond }^1(D)$$

$$L_{\diamond }^2(\partial D)$$

$$\sigma :D\to \text{mathbbm}R$$

$$L_{+}^{\infty }(D)\ni \sigma \mapsto \Lambda (\sigma )\in \mathcal{L}(L_{\diamond }^2(\partial D)),$$

$$\Lambda (\sigma ):g\mapsto u_{\sigma }^g{|}_{\partial D}$$

$$\nabla \cdot \sigma \nabla u_{\sigma }^g=0\mathrm{in}D,\sigma {\partial }_{\nu }{u}_{\sigma }^g=g\mathrm{on}\partial D.$$

$${\int }_{\partial D}g\Lambda (\sigma )h\mathrm{d}s={\int }_D\sigma \nabla u_{\sigma }^g\cdot \nabla u_{\sigma }^h\mathrm{d}x.$$

$${\int }_D({\sigma }_1-{\sigma }_0)|\nabla u_0{|}^2\mathrm{d}x$$

$$\ge {\int }_D{\sigma }_0^2\left(\frac{1}{{\sigma }_0}-\frac{1}{{\sigma }_1}\right)|\nabla u_0{|}^2\mathrm{d}x,$$

$$A\ge B$$

$${\int }_{\partial D}g\left(A-B\right)g\mathrm{d}s\ge 0\text{mbox}forallg\in L_{\diamond }^2(\partial D).$$

$${\sigma }_0\le {\sigma }_1\text{ implies }\Lambda ({\sigma }_1)\le \Lambda ({\sigma }_0),$$

$$\sigma :\Omega \to L_{+}^{\infty }(D),$$

$$\Lambda (\sigma )(\omega ):=\Lambda (\sigma (\omega ))$$

$$f:\Omega \to X$$

$$\Vert f\Vert {}_{L^p(\Omega ;X)}:={\left({\int }_{\Omega }\Vert f(\omega ){\Vert }_X^p\mathrm{d}P\right)}^{1/p}<\infty .$$

$$\text{mathbbm}E(\sigma )$$

$$\text{mathbbm}E({\sigma }^{-1})$$

$$\text{mathbbm}E(\Lambda (\sigma ))$$

$$\frac{\Vert \Lambda (\sigma (\omega ))\Vert {}_{\mathcal{L}(L_{\diamond }^2(\partial D))}}{\Vert \Lambda (1)\Vert {}_{\mathcal{L}(L_{\diamond }^2(\partial D))}}$$

$$=1+\frac{{\sup }_{\Vert g\Vert {}_{L_{\diamond }^2(\partial D)}=1}\left|{\int }_{\partial D}g\left(\Lambda (\sigma )-\Lambda (1)\right)g\mathrm{d}s\right|}{{\sup }_{\Vert g\Vert {}_{L_{\diamond }^2(\partial D)}=1}\left|{\int }_{\partial D}g\Lambda (1)g\mathrm{d}s\right|}$$

$$\le 1+\max \left\{{\Vert 1-\sigma (\omega )\Vert }_{L^{\infty }(D)},{\Vert 1-\sigma (\omega {)}^{-1}\Vert }_{L^{\infty }(D)}\right\}$$

$$\le 2+\max \left\{{\Vert \sigma (\omega )\Vert }_{L^{\infty }(D)},{\Vert \sigma (\omega {)}^{-1}\Vert }_{L^{\infty }(D)}\right\}.$$

$$\sigma (\omega )=1+\sum _{q=1}^Q{\theta }_q(\omega ){\chi }_q,\omega \in \Omega ,$$

$${\theta }_q:\Omega \to [a_q,b_q]\subset (-1,\infty ),q=1,\dots ,Q,$$

$$\text{mathbbm}E(\sigma )$$

$$=1+\sum _{q=1}^Q\frac{1}{b_q-a_q}{\int }_{a_q}^{b_q}x\mathrm{d}x{\chi }_q=1+\sum _{q=1}^Q\frac{a_q+b_q}{2}{\chi }_q,$$

$$\text{mathbbm}E({\sigma }^{-1})$$

$$\sigma (\omega )=1+{\chi }_{\theta (\omega )},\omega \in \Omega ,$$

$$\Vert \sigma ({\omega }_1)-\sigma ({\omega }_2)\Vert {}_{L^{\infty }(D)}=\Vert {\chi }_{\theta ({\omega }_1)}-{\chi }_{\theta ({\omega }_2)}\Vert {}_{L^{\infty }(D)}=1,$$

$${\int }_D(\text{mathbbm}E(\sigma )-{\sigma }_0))|\nabla u_0{|}^2\mathrm{d}x$$

$$\ge {\int }_D{\sigma }_0^2\left({\sigma }_0^{-1}-\text{mathbbm}E({\sigma }^{-1})\right)|\nabla u_0{|}^2\mathrm{d}x,$$

$$\displaystyle\Lambda(\mathbbm{E}(\sigma))\leq\mathbbm{E}(\Lambda(\sigma)))\leq\Lambda\big{(}\mathbbm{E}(\sigma^{-1})^{-1}\big{)}.$$

$$1+\alpha$$

$$\sigma(\omega,x):=\left\{\begin{array}[]{l l }1&\quad\text{ for }x\in D\setminus A,\\ \sigma_{A}(\omega,x)&\quad\text{ for }x\in A,\end{array}\right.$$

Full text

1 Introduction

Electrical impedance tomography (EIT) is an imaging modality for finding information about the internal conductivity of an examined physical body from boundary measurements of current and potential. The reconstruction task of EIT is nonlinear and highly illposed. EIT has potential applications in, e.g., medical imaging, geophysics, nondestructive testing of materials and monitoring of industrial processes. For general information on EIT, we refer to [1, 4, 5, 6, 11, 28, 41, 42, 51] and the references therein.

This work considers detection of stochastic inclusions by EIT. The conductivity inside the examined object is assumed to consist of a known deterministic background with embedded inclusions within which the conductivity is a random field. The inclusion shapes are deterministic. The available measurement is assumed to be the expected value of the corresponding Neumann-to-Dirichlet boundary operator, i.e., the expectation of the idealized current-to-voltage map. We demonstrate that the Factorization Method and the Monotonicity Method can be utilized to reconstruct the supports of the stochastic inhomogeneities essentially in the same way as in the deterministic case if the inhomogeneities exhibit high enough conductivity contrast in the sense of expectation. Detecting inclusions with random conductivities from expected boundary measurements can be regarded an interesting problem setting from a mere mathematical standpoint, but it also has potential applications: consider, e.g., temporally averaged boundary measurements corresponding to an inhomogeneity whose conductivity varies in time (e.g., a cancerous lesion subject to blood flow).

To the best of our knowledge, this is the first work that considers reconstructing properties of a random conductivity from the expectation of the boundary measurements. With that being said, it should be mentioned that for inverse scattering in the half-plane it has been shown that a single realization of the backscatter data at all frequencies uniquely defines the principal symbol of the covariance operator for a random scattering potential or a Robin parameter in the boundary condition [27, 38]. Moreover, time reversal related inverse scattering techniques have been used extensively for the detection of inclusions embedded in random media; see, e.g., [2, 3, 7] and the references therein.

The Factorization Method was originally introduced in the framework of inverse obstacle scattering by Kirsch [35], and it was subsequently carried over to EIT by Brühl and Hanke [8, 9]. The Factorization Method has attracted a considerable amount of attention in the scientific community over the past twenty years; in the context of EIT we refer the reader to [10, 13, 14, 15, 16, 17, 18, 20, 22, 24, 29, 31, 36, 40, 44, 45, 46] and the overview articles [19, 21, 37] for more information on the developments since the seminal article of Kirsch. The Monotonicity Method [25] is a more recent inclusion detection technique that can be considered natural for EIT as it takes explicitly advantage of the coercivity of the underlying conductivity equation. In particular, unlike the Factorization Method, a variant of the Monotonicity Method can handle indefinite inclusions, i.e., it stays functional even if some parts of the inhomogeneities are more and some other parts less conductive than the background. Monotonicity-based methods were first formulated and numerically tested by Tamburrino and Rubinacci [49, 50]. Their rigorous justification is given in [25] based on the theory of localized potentials [15].

For simplicity, we employ the idealized continuum model for EIT, that is, we unrealistically assume to be able to drive any square-integrable current density into the object of interest and to measure the resulting potential everywhere on the boundary. Consult [30] for analysis of the relationship between the continuum model and the complete electrode model, which is the most accurate model for real-world EIT measurements [12, 48]. Observe also that the Factorization Method has, in fact, been demonstrated to function also for realistic electrode measurements [10, 24, 40], and that the monotonicity method allows to characterize the achievable resolution in realistic measurement settings [26].

This text is organized as follows. In Section 2, we introduce the stochastic Neumann-to-Dirichlet boundary map corresponding to a random conductivity field and examine how the fundamental monotonicity results for the related quadratic forms are carried over from the deterministic case. Subsequently, Section 3 formulates and proves the theorems that demonstrate the functionality of the Factorization and Monotonicity Methods for detecting stochastic inclusions. The two-dimensional numerical examples for the Factorization Method are presented in Section 4; the stochastic measurement data are simulated by applying the stochastic finite element method (cf., e.g., [47]) to the conductivity equation with a random coefficient.

2 The setting and monotonicity results

In this section, we first review certain monotonicity properties of the deterministic Neumann-to-Dirichlet map associated with the conductivity equation. These inequalities are subsequently interpreted from the viewpoint of stochastic conductivities.

2.1 The deterministic Neumann-to-Dirichlet operator

Let $D\subset\mathbbm{R}^{n}$, $n\geq 2$ be a bounded domain with a smooth boundary $\partial D$ and outer unit normal vector $\nu$. The subset of $L^{\infty}(D)$-functions with positive essential infima is denoted by $L_{+}^{\infty}(D)$. Moreover, $H_{\diamond}^{1}(D)$ and $L_{\diamond}^{2}(\partial D)$ are defined to be the subspaces of $H^{1}(D)$- and $L^{2}(\partial D)$-functions, respectively, with vanishing integral means over $\partial D$, i.e.,

[TABLE]

In the context of EIT, $D$ describes the imaging subject. A deterministic conductivity inside the subject is modeled by a function

[TABLE]

such that $\sigma\in L_{+}^{\infty}(D)$. The deterministic nonlinear forward operator of EIT is defined by

[TABLE]

where

[TABLE]

and $u_{\sigma}^{g}\in H_{\diamond}^{1}(D)$ solves

[TABLE]

The bilinear form associated with $\Lambda(\sigma)$ is

[TABLE]

It is well-known that $\Lambda(\,\cdot\,)$ is continuous and for each $\sigma\in L^{\infty}_{+}(D)$, the linear operator $\Lambda(\sigma):L_{\diamond}^{2}(\partial D)\to L_{\diamond}^{2}(\partial D)$ is self-adjoint and compact.

A useful tool in studying $\Lambda$ is the following monotonicity result from Ikehata, Kang, Seo, and Sheen [33, 34]. It has frequently been used in deterministic inclusion detection, cf., e.g., [21, 22, 23, 25, 32, 36].

Lemma 1**.**

Let $\sigma_{0},\sigma_{1}\in L^{\infty}_{+}(D)$. Then, for all $g\in L_{\diamond}^{2}(\partial D)$,

[TABLE]

with $u_{0}:=u^{g}_{\sigma_{0}}\in H^{1}_{\diamond}(D)$ solving (1).

Proof.

A short proof can be found, e.g., in [25, Lemma 3.1]. ∎

For two bounded self-adjoint operators $A,B:\ L_{\diamond}^{2}(\partial D)\to L_{\diamond}^{2}(\partial D)$ we write

[TABLE]

if the inequality holds in the sense of quadratic forms, i.e., if

[TABLE]

With this convention, Lemma 1 yields that for all $\sigma_{0},\sigma_{1}\in L^{\infty}_{+}(D)$,

[TABLE]

which is why we call Lemma 1 a monotonicity estimate.

2.2 The stochastic Neumann-to-Dirichlet operator

We model stochastic objects by treating the conductivity as a random field

[TABLE]

where $(\Omega,\Sigma,P)$ is a probability space. By pointwise concatenation

[TABLE]

the Neumann-to-Dirichlet operator becomes random as well.

To be able to rigorously introduce moments (such as expectation or variance), we recall the definition of the Lebesgue–Bochner space $L^{p}(\Omega;X)$:

Definition 2**.**

Given a Banach space $X$, a function

[TABLE]

is called strongly measurable if it is measurable (with respect to $\Sigma$ and the Borel algebra on $X$) and essentially separably valued. The space $L^{p}(\Omega;X)$, $1\leq p<\infty$ is the set of (equivalence classes with respect to almost everywhere equality of) strongly measurable functions satisfying

[TABLE]

For $Y\subseteq X$, $L^{p}(\Omega;Y)$ is defined as the subset of those $f\in L^{p}(\Omega;X)$ for which $f(\omega)\in Y$ almost surely.

The following lemma shows that the expected value of the Neumann-to-Dirichlet operator $\Lambda(\sigma)$ is well defined if both the conductivity $\sigma$ and its reciprocal, the resistivity $\sigma^{-1}$, have well defined expectations.

Lemma 3**.**

Let $\sigma,\sigma^{-1}\in L^{1}(\Omega;L^{\infty}_{+}(D))\subseteq L^{1}(\Omega;L^{\infty}(D))$.

Then $\Lambda(\sigma)\in L^{1}(\Omega;\mathcal{L}(L_{\diamond}^{2}(\partial D))),$ and the following expectations are well defined:

[TABLE]

Furthermore, $\mathbbm{E}(\Lambda(\sigma)):L_{\diamond}^{2}(\partial D)\to L_{\diamond}^{2}(\partial D)$ is self-adjoint and compact.

Proof.

The fact that $\mathbbm{E}(\sigma)$ and $\mathbbm{E}(\sigma^{-1})$ are well defined and belong to $L^{\infty}_{+}(D)$ follows directly from the assumption $\sigma,\sigma^{-1}\in L^{1}(\Omega;L^{\infty}_{+}(D))$. Furthermore, $\Lambda(\sigma):\ \Omega\to\mathcal{L}(L_{\diamond}^{2}(\partial D))$ is measurable since it is a concatenation of a measurable and a continuous function. Since continuous images of separable sets are separable, the same argument shows that $\Lambda(\sigma)$ is also strongly measurable.

To show that $\Lambda(\sigma)$ has a finite expectation, we fix the constant deterministic conductivity $\sigma_{0}\equiv 1\in L^{\infty}_{+}(\Omega)$, and employ the monotonicity lemma 1 and (2). It follows that for almost all $\omega\in\Omega$,

[TABLE]

Since $\sigma,\sigma^{-1}\in L^{1}(\Omega;L^{\infty}_{+}(D))$, this implies that $\Lambda(\sigma)\in L^{1}(\Omega;\mathcal{L}(L_{\diamond}^{2}(\partial D)))$.

The self-adjointness and compactness of $\mathbbm{E}(\Lambda(\sigma))$ follow as in the deterministic case. ∎

Our stochastic framework covers an anomaly of deterministic size with a random conductivity value but not an anomaly with random size. We will demonstrate this in the following example.

Example 4**.**

We consider a domain $D$ with homogeneous unit background conductivity and an embedded anomaly.

1. (a)

First assume that the anomaly consists of finitely many, mutually disjoint pixels, where the conductivity has a constant, uniformly distributed, stochastic value. To be more precise,

[TABLE]

*where $\chi_{q}$ is the characteristic function of the *$q$th pixel and

[TABLE]

are mutually independent, uniformly distributed random variables. Hence, each $\theta_{q}$ has a probability density function given by $(b_{q}-a_{q})^{-1}\chi_{[a_{q},b_{q}]}$. We then have, using the linearity of the expectation operator and the disjointness of the pixels

[TABLE] 2. (b)

Now assume that the domain $D=B_{1}(0)\subseteq\mathbbm{R}^{2}$ is the unit disk that contains a circular anomaly with stochastic radius, say,

[TABLE]

where $\theta:\Omega\to[0,1/2]$ is uniformly distributed and $\chi_{\theta}$ is the characteristic function of a disk with radius $\theta$ centered at the origin. Then for $\omega_{1},\omega_{2}\in\Omega$ such that $\theta(\omega_{1})\not=\theta(\omega_{2})$, it holds that

[TABLE]

which shows that $\sigma:\ \Omega\to L^{\infty}_{+}(D)$ cannot be essentially separably valued.

By replacing $\sigma_{1}$ in Lemma 1 by a random conductivity $\sigma$ and taking expectations, we immediately obtain the following stochastic monotonicity estimate.

Corollary 5**.**

Consider a deterministic conductivity $\sigma_{0}$ and a stochastic conductivity $\sigma$, i.e., $\sigma_{0}\in L^{\infty}_{+}(D)$ and $\sigma,\sigma^{-1}\in L^{1}(\Omega;L^{\infty}_{+}(D))$. Then, for all $g\in L_{\diamond}^{2}(\partial D)$,

[TABLE]

where $u_{0}:=u^{g}_{\sigma_{0}}\in H^{1}_{\diamond}(D)$ solves (1). In particular, by choosing in turns $\sigma_{0}=\mathbbm{E}(\sigma)$ and $\sigma_{0}=\mathbbm{E}(\sigma^{-1})^{-1}$, it follows that

[TABLE]

3 Detecting stochastic inclusions

We consider a domain $D$ that has known deterministic background properties but contains an anomaly $A$ with random conductivity. For simplicity, we derive our results only for the case of a homogeneous unit background conductivity and assume that $A$ is an open set for which $\overline{A}\subset D$ has a connected complement. Following the line of reasoning in [21, 25], it should be straightforward to generalize the results for (known and deterministic) piecewise analytic background conductivities and a more general subset $A$. In the rest of this section, $u_{0}:=u^{g}_{\sigma_{0}}\in H^{1}_{\diamond}(D)$ solves (1) for the considered background conductivity $\sigma_{0}\equiv 1$ and a current density $g\in L^{2}_{\diamond}(\partial D)$ that should be clear from the context.

In the deterministic case, rigorous noniterative inclusion detection methods such as the Factorization Method or the Monotonicity Method are known to correctly recover the shape of an unknown anomaly whenever there exists $\alpha>0$ such that the conductivity in the anomaly is either

[TABLE]

As noted above, here we consider a domain $D$ containing a stochastic anomaly $A$, that is,

[TABLE]

where $\sigma_{A}:\ \Omega\to L^{\infty}_{+}(A)$ is a random field with $\sigma_{A},\sigma_{A}^{-1}\in L^{1}(\Omega;L^{\infty}_{+}(A))$. It is tempting to try to detect $A$ from the stochastic boundary measurements by applying the deterministic Factorization or Monotonicity Method to the expected value $\mathbbm{E}(\Lambda(\sigma))$ of the stochastic Neumann-to-Dirichlet boundary operator. We will show that these methods indeed correctly recover the anomaly $A$ if there exists $\alpha>0$ such that, almost everywhere in $A$, either

[TABLE]

Note that, by Jensen’s inequality, $\mathbbm{E}(\sigma_{A})\geq\mathbbm{E}(\sigma_{A}^{-1})^{-1}$, meaning that (a) also implies that $\mathbbm{E}(\sigma_{A})\geq 1+\alpha$. Analogously, (b) implies that $\mathbbm{E}(\sigma_{A}^{-1})^{-1}\leq 1-\alpha$.

We start with the Factorization Method.

Theorem 6** (Factorization Method).**

Assume that there exists $\alpha>0$ such that one of the two conditions in (5) is valid. Then for any nonvanishing dipole moment $d\in\mathbbm{R}^{n}$ and every point $z\in D\setminus\partial A$, it holds that

[TABLE]

where $\Phi_{z,d}\in C^{\infty}(\overline{D}\setminus\{z\})\cap H^{-n/2+1-\varepsilon}(D)$, $\varepsilon>0$, is the unique solution of

[TABLE]

i.e., the potential created by a dipole at $z$ corresponding to the moment $d$, the unit conductivity and vanishing current density on $\partial D$.

Proof.

In our current setting, the monotonicity result of Corollary 5 reads

[TABLE]

for all $g\in L_{\diamond}^{2}(\partial D)$. In case (a), we have $\mathbbm{E}(\sigma_{A}^{-1})\leq(1+\alpha)^{-1}$, and thus for all $g\in L_{\diamond}^{2}(\partial D)$,

[TABLE]

In case (b), we use $\mathbbm{E}(\sigma_{A})\leq 1-\alpha$ to obtain

[TABLE]

for all $g\in L_{\diamond}^{2}(\partial D)$. The assertion now follows from the previous two estimates by the same arguments as in the proof of the deterministic Factorization Method in [21, Theorem 5]. ∎

Then it is the turn of the Monotonicity Method.

Theorem 7** (Monotonicity Method).**

1. (a)

For every constant $\alpha>0$ such that $\mathbbm{E}(\sigma_{A}^{-1})^{-1}\geq 1+\alpha$ almost everywhere on $A$, and for every open ball $B\subset D$,

[TABLE] 2. (b)

For every $0<\alpha<1$ such that $\mathbbm{E}(\sigma_{A})\leq 1-\alpha$ almost everywhere on $A$, and for every open ball $B\subset D$,

[TABLE]

Proof.

We prove the two cases separately.

1. (a)

If $B\subseteq A$, then $1+\alpha\chi_{B}\leq\mathbbm{E}(\sigma^{-1})^{-1}$, and we obtain from (3) and (4) that

[TABLE]

Now, let $B\not\subseteq A$. We will show that $\Lambda(1+\alpha\chi_{B})\not\geq\mathbbm{E}(\Lambda(\sigma))$. Since (3) proves that shrinking $B$ enlarges $\Lambda(1+\alpha\chi_{B})$, we may assume without loss of generality that $B\subseteq D\setminus\overline{A}$. Using the monotonicity result in Corollary 5, we obtain that

[TABLE]

for all $g\in L_{\diamond}^{2}(\partial D)$. Using a localized potential with large energy on $B$ and small energy on $A$ (see [15, 25]), it follows that there exists $g\in L_{\diamond}^{2}(\partial D)$ satisfying

[TABLE]

so that $\Lambda(1+\alpha\chi_{B})-\mathbbm{E}(\Lambda(\sigma))\not\geq 0$. 2. (b)

If $B\subseteq A$, then $1-\alpha\chi_{B}\geq\mathbbm{E}(\sigma)$, and we obtain from (3) and (4) that

[TABLE]

Now, let $B\not\subseteq A$. To show that $\Lambda(1-\alpha\chi_{B})\not\leq\mathbbm{E}(\Lambda(\sigma))$, we can again assume without loss of generality that $B\subseteq D\setminus\overline{A}$. Employing Corollary 5 with $\sigma_{0}=1+\alpha\chi_{B}$, we obtain this time around that

[TABLE]

for all $g\in L_{\diamond}^{2}(\partial D)$. As in case (a), we can use a localized potential with large energy on $B$ and small energy on $A$ to conclude that there exists $g\in L_{\diamond}^{2}(\partial D)$ satisfying

[TABLE]

so that $\Lambda(1-\alpha\chi_{B})-\mathbbm{E}(\Lambda(\sigma))\not\leq 0$.

This completes the proof. ∎

Theorem 7 provides an explicit characterization of the anomaly: If a suitable threshold value $\alpha>0$ is available, one can test whether any given open ball is enclosed by $A$ (assuming infinite measurement/numerical precision). The following theorem is essential in the numerical implementation of the Monotonicity Method. It demonstrates that, as in the deterministic case, the monotonicity test can be performed using the linearized Neumann-to-Dirichlet map, which considerably reduces the computational cost as one does not need to simulate a new nonlinear Neumann-to-Dirichlet operator for testing each new ball and threshold value. To get set, we define $\Lambda^{\prime}(1)\chi_{B}\in\mathcal{L}(L_{\diamond}^{2}(\partial D))$ for a ball $B\subseteq D$ via its quadratic form

[TABLE]

where once again $u_{0}\in H^{1}_{\diamond}(D)$ is the reference potential corresponding to the unit background conductivity and the current density $g$.

Theorem 8** (Linearized Monotonicity Method).**

1. (a)

If there exists $\alpha>0$ such that $\mathbbm{E}(\sigma_{A}^{-1})^{-1}\geq 1+\alpha$ almost everywhere on $A$, then for every open ball $B\subset D$ and for every $0<\beta\leq\frac{\alpha}{1+\alpha}$,

[TABLE] 2. (b)

If there exists $\alpha>0$ such that $\mathbbm{E}(\sigma_{A})\leq 1-\alpha$ almost everywhere on $A$, then for every open ball $B\subset D$ and for every $0<\beta<\alpha$,

[TABLE]

Proof.

1. (a)

For every $\alpha>0$ and $0<\beta\leq\frac{\alpha}{1+\alpha}$, it follows from [25, Theorem 4.3] that

[TABLE]

Hence, theorem 7 shows that

[TABLE]

The converse results from a similar combination of monotonicity and localized potentials arguments as above. 2. (b)

The claim follows as in part (a) by using [25, Remark 4.5].

∎

Remark 9**.**

In the next section we will numerically test the methods on the setting described in example 4(a) where the considered anomalies consist of a finite number of subregions, or pixels whose (constant) conductivity levels are mutually independent uniformly distributed random variables. However, it should be noted that our theoretic results hold for general conductivity distributions, as long as $\sigma,\sigma^{-1}\in L^{1}(\Omega;L^{\infty}_{+}(D))\subseteq L^{1}(\Omega;L^{\infty}(D))$ and (5) is fulfilled. In particular, this includes log-normal distributed conductivities.

4 Numerical examples

This section presents the numerical examples on the detection of stochastic inclusions from the mean of the Neumann-to-Dirichlet operator. We concentrate on the Factorization Method since it has thus far been more widely studied within the inverse problems community, making the comparison of our results to those achievable in the deterministic framework more straightforward. However, there is no reason to expect that the Monotonicity Method would perform any worse.

4.1 Simulation of data

Let us first discuss how the expectation $\mathbbm{E}(\Lambda(\sigma))$ of the Neumann-to-Dirichlet operator can be approximated for the numerical experiments that follow. The idea is to compute the expected boundary voltages $\mathbbm{E}(u_{\sigma}^{g}|_{\partial D})=\mathbbm{E}(u_{\sigma}^{g})|_{\partial D}$ separately for a (finite) number of boundary currents $g$.

Once the distribution of the conductivity is given, the expectations $\mathbbm{E}(u_{\sigma}^{g})$ can be obtained by Monte Carlo simulation, i.e., by sampling from the distribution and solving a deterministic problem (1) for each realization. Another approach, which we employ here, is known as the stochastic finite element method (SFEM). We briefly sketch the main idea of SFEM under the assumption that $\sigma\in L^{2}(\Omega;L^{\infty}_{+}(D))$, which will be the case in all our examples considered below. The forward problem describing EIT with a stochastic conductivity is to find $u_{\sigma(\omega)}^{g}\in H_{\diamond}^{1}(D)$ that solves

[TABLE]

$P$-almost surely for a given current density $g\in L_{\diamond}^{2}(\partial D)$, which is equivalent to the standard variational formulation

[TABLE]

$P$-almost surely. This is in turn equivalent to the extended variational formulation of finding $u_{\sigma}^{g}\in L^{2}(\Omega;H_{\diamond}^{1}(D))$ such that

[TABLE]

for all $\phi\in L^{2}(\Omega;H_{\diamond}^{1}(D))$. The key idea of the stochastic finite element method is to approximate statistical properties, such as moments, of the solution to (6) by applying the Galerkin method to the extended variational formulation (7). Hence, we choose finite dimensional subspaces $V_{\diamond}^{h}\subset H_{\diamond}^{1}(D)$ and $S_{m}\subset L^{2}(\Omega)$ and determine the function $\tilde{u}_{\sigma}^{g}\in V_{\diamond}^{h}\otimes S_{m}\subset L^{2}(\Omega;H_{\diamond}^{1}(D))$ that solves (7) for all $\phi\in V_{\diamond}^{h}\otimes S_{m}\subset L^{2}(\Omega;H_{\diamond}^{1}(D))$.

As $V^{h}$ we use a standard piecewise linear finite element space. In all our examples below, the probability space $(\Omega,\Sigma,P)$ can be parametrized by a finite number of independent random variables $\theta_{q}:\Omega\rightarrow\Theta_{q}\subseteq\mathbbm{R}$, $q=1,\ldots,Q$. Accordingly, we choose $S_{m}$ to be the space of all $Q$-variate polynomials with total degree less than or equal to $m$; the corresponding basis functions are chosen to be orthonormal, i.e., suitably scaled multivariate Legendre polynomials in our setting. Take note that the dimension of $S_{m}$ is

[TABLE]

For more general conductivities one can consider, e.g., truncated Karhunen–Loève expansions of the random conductivity field in question.

Both the Factorization and the Monotonicity Method require only the difference operator $\mathbbm{E}(\Lambda(\sigma))-\Lambda(1)$ as the input. Instead of approximating $u_{\sigma}^{g}$ as described above in order to calculate the expectation of $\Lambda(\sigma)$ and then subtracting the deterministic Neumann-to-Dirichlet map corresponding to the unit conductivity, a numerically more stable method is to directly compute the expectation $\mathbbm{E}(\Lambda(\sigma)-\Lambda(1))$. This can be done by considering the boundary value problem satisfied by the difference potential (cf., e.g., [16, Section 3])

[TABLE]

where $u_{0}^{g}:=u_{\sigma_{0}}^{g}$ is the solution to (1) with background conductivity $\sigma_{0}\equiv 1$.

Our numerical experiments are conducted in the unit disk $D\subset\mathbbm{R}^{2}$, for which the most natural $L^{2}(\partial D)$-normalized boundary currents are the Fourier basis functions

[TABLE]

for some $T\in\mathbbm{N}$ and with $\gamma$ denoting the angular coordinate. For these current densities, the gradient of the background solution $u_{0}^{g}$ can be presented in a closed form, which can be used explicitly on the right hand side of the variational formulation of (8). Using a stochastic finite element method as explained above, we solve (8) for all current patterns in (9), and then expand the expectation value of the correspondingly approximated difference boundary data $\mathbbm{E}(\tilde{w}_{\sigma}^{g_{t}})$ in the basis (9) by using FFT. This results in an approximate representation $L\in\mathbbm{R}^{2T\times 2T}$ of $\mathbbm{E}(\Lambda(\sigma))-\Lambda(1):L^{2}_{\diamond}(\partial D)\to L^{2}_{\diamond}(\partial D)$ with respect to the Fourier basis (9).

The considered inclusion geometries are illustrated in Figure 1. In each case, the anomalies consist of a finite number of subregions, or pixels whose (constant) conductivity levels are mutually independent uniformly distributed random variables. Thus, the random conductivity can be written as

[TABLE]

where $\chi_{q}$ is the indicator function of the $q$th pixel. The number of pixels $Q$ varies between one and nine depending on the configuration.

In all numerical tests we use third order stochastic basis polynomials, i.e., we choose $m=3$, and the deterministic FEM meshes are such that the dimension of $V_{h}$ is approximately $33000$. The resulting SFEM system matrices are very sparse [43]. The number of employed Fourier basis functions is $2T=100$ in all tests.

In addition to considering the unperturbed simulated measurement matrix $L\in\mathbbm{R}^{2T\times 2T}$ that contains only numerical errors, we also test the Factorization Method after adding $0.1\%$ noise to $L$. More precisely, we generate a random matrix $E\in\mathbbm{R}^{2T\times 2T}$ with uniformly distributed elements in the interval $[-1,1]$ and replace $L$ by its noisy version

[TABLE]

where $\|\cdot\|_{2}$ denotes the spectral norm. This is the same noise model as in [16] for the deterministic Factorization Method.

4.2 Factorization Method

Recall that Theorem 6 provides a binary test for deciding whether a point is inside the anomaly $A$ or not: Assuming that (5) holds for some $\alpha>0$, then for any $z\in D\setminus\partial A$,

[TABLE]

To numerically implement this test, we follow [16]. The boundary potential $\Phi_{z,d}|_{\partial D}$ belongs to $\mathcal{R}(|\mathbbm{E}(\Lambda(\sigma))-\Lambda(1)|^{1/2})$ if and only if it satisfies the so-called Picard criterion, which compares the decay of the scalar projections of $\Phi_{z,d}|_{\partial D}$ on the eigenfunctions of the compact and self-adjoint operator $|\mathbbm{E}(\Lambda(\sigma))-\Lambda(1)|^{1/2}$ to the decay of the corresponding eigenvalues. The leading idea is to write an approximate version of the Picard criterion by using the singular values and vectors of $L\in\mathbbm{R}^{2T\times 2T}$ in place of the eigensystem of $|\mathbbm{E}(\Lambda(\sigma))-\Lambda(1)|$. We refer to [39] for a proof that this approach can be formulated as a regularization strategy.

Let us introduce a singular value decomposition for the discretized boundary operator $L\in\mathbbm{R}^{2T\times 2T}$, that is,

[TABLE]

with nonnegative singular values $\{\lambda_{t}\}\subset\mathbbm{R}$, which are sorted in decreasing order, and orthonormal bases $\{u_{t}\},\{v_{t}\}\subset\mathbbm{R}^{2T}$. Denoting by $F_{z,d}\in\mathbbm{R}^{2T}$ the Fourier coefficients of $\Phi_{z,d}$, we introduce an indicator function

[TABLE]

where $1\leq\tau\leq 2T$ should be chosen so that $\lambda_{\tau+1}$ is the first singular value below the expected measurement/numerical error. Note that the formation of $F_{z,d}$ is trivial because the functional form of $\Phi_{z,d}$ is known explicitly when $D$ is the unit disk (cf., e.g., [9]). It is to be expected that ${\rm Ind}:D\to\mathbbm{R}_{+}$ takes larger values inside than outside the anomaly $A$ (cf., e.g., [16]), meaning that plotting ${\rm Ind}$ should provide information on the whereabouts of $A$. In our numerical tests, we plot the indicator function on an equidistant grid that is chosen independently of the finite element meshes used for solving the forward problems. The employed dipole moments are $d_{r}=(\cos(2\pi r/3),\sin(2\pi r/3))\in\mathbbm{R}^{2}$ for $r\in\{1,2,3\}$ and the indicator function used in reconstructions is actually the sum of the indicator functions ${\rm Ind}(y;d_{r},\tau)$ corresponding to those three dipole moments. For all noiseless experiments, we use the maximal spectral cut-off index $\tau=2T=100$. (We remind the reader that the Factorization Method is not designed to produce information about the conductivity of the inhomogeneity. In other words, only the ‘support’ of the indicator function is relevant.)

Test 1. Our first and simplest test case considers one disk-shaped inclusion, as depicted in the top left image of Figure 1. We choose $\Theta_{1}=[8,10]\subset\mathbbm{R}$ so that all realizations of the stochastic inclusion significantly differ from the unit background conductivity. In this case, $\mathbbm{E}(\sigma_{A}^{-1})^{-1}\approx 9.97$, which obviously satisfies part (a) of (5). Figure 2 shows the indicator function for both noiseless and noisy ($\tau=6$) Neumann-to-Dirichlet data; according to a visual inspection, the support of the indicator function matches well with the inclusion in both cases. This is not a big surprise: All possible realizations of the random conductivity satisfy the requirements of the deterministic Factorization Method, and thus one would expect reconstructions that are comparable to those in the deterministic setting (cf., e.g., [16]).

Test 2. The second experiment considers the nonconvex inclusion in the top right image of Figure 1. The expected conductivity values of the three pixels constituting the anomaly are now strictly less than the unit background. More precisely, the pixelwise values vary according to $\theta_{1}\in[-0.99,-0.41]$, $\theta_{2}\in[-0.99,-0.01]$ and $\theta_{3}\in[-0.99,0.39]$; see (10). Note that the condition (b) of (5) is satisfied, even though one of the three pixels can have either positive or negative contrast. The reconstructions produced by the Factorization Method are presented in Figure 3. Even without noise, the nonconvex edge of the inclusion is clearly smoothened, but otherwise the noiseless reconstruction is accurate. In the noisy case, the indicator function ($\tau=7$) contains practically no information on the nonconvexity of the inclusion, but the size and the location of the anomaly are still reproduced accurately.

Test 3. In the third example, we have two separate inclusions and altogether nine random variables as shown in the bottom left image of Figure 1. The conductivity intervals have width of six units for the pixels in the disk and width of three units for those in the square. All realizations of the discoidal inclusion are thus more conductive than the background, while the conductivity contrasts of some pixels in the square anomaly may be either positive or negative. The condition (a) of (5) is anyway satisfied: $\mathbbm{E}(\sigma_{A}^{-1})^{-1}>4.3$ for the disk, and $\mathbbm{E}(\sigma_{A}^{-1})^{-1}>1.5$ for the square. The reconstruction from noiseless data in Figure 4 accurately indicates the shapes and sizes of the two inclusions. With noisy data, the reconstruction ($\tau=12$) tends to produce one large inclusion instead of two smaller ones. The Factorization Method has been observed to exhibit similar bridging behavior also in the deterministic setting (see, e.g., [16]).

Test 4. Finally, we test the Factorization Method for three well separated disk-like inclusions, each characterized by a uniformly distributed conductivity in the interval $[0.5,3.5]$; see the bottom right image of Figure 1. Note that the condition (a) of (5) is satisfied as $\mathbbm{E}(\sigma_{A}^{-1})^{-1}\approx 1.54$, although each realization of the anomaly can have positive, negative or indefinite contrast compared to the unit background. Once again, the noiseless reconstruction in Figure 5 is very accurate, but measurement noise deteriorates the representation of the inclusion boundaries in the plot of the indicator function ($\tau=11$).

5 Concluding remarks

We have shown that an application of the Factorization Method or the Monotonicity Method to the mean value of the Neumann-to-Dirichlet boundary map characterizes stochastic anomalies, i.e., inclusions with random conductivities but fixed shapes, assuming that the contrast of the anomalies is high enough in the sense of expectation. This theoretical result was complemented by two-dimensional numerical examples demonstrating the functionality of the Factorization Method in this stochastic setting.