Numerical determination of anomalies in multifrequency electrical impedance tomography

TL;DR

This paper introduces a spectral decomposition method for detecting anomalies in multifrequency electrical impedance tomography, enabling the reconstruction of conductivity profiles with high accuracy using gradient descent algorithms.

Contribution

It presents a novel spectral decomposition approach to identify and reconstruct anomalies in multifrequency electrical impedance tomography, improving detection accuracy.

Findings

Successful reconstruction of anomalies using spectral decomposition.

Effective use of gradient descent for numerical experiments.

Enhanced detection of piecewise constant conductivities.

Abstract

The multifrequency electrical impedance tomography consists in retrieving the conductivity distribution of a sample by injecting a finite number of currents with multiple frequencies. In this paper we consider the case where the conductivity distribution is piecewise constant, takes a constant value outside a single smooth anomaly, and a frequency dependent function inside the anomaly itself. Using an original spectral decomposition of the solution of the forward conductivity problem in terms of Poincar\'e variational eigenelements, we retrieve the Cauchy data corresponding to the extreme case of a perfect conductor, and the conductivity profile. We then reconstruct the anomaly from the Cauchy data. The numerical experiments are conducted using gradient descent optimization algorithms.

Tables

Table 1. Table 1 . Errors between u 0 ​ r ​ e ​ c ​ o ​ n ​ s ​ t ​ r ​ u ​ c ​ t subscript 𝑢 0 𝑟 𝑒 𝑐 𝑜 𝑛 𝑠 𝑡 𝑟 𝑢 𝑐 𝑡 u_{0reconstruct} and u 0 subscript 𝑢 0 u_{0} .

•	ellipse	square	near-boundary	small-central
$f=f_1$	0.04707	0.11973	0.00956	0.00502
$f=f_2$	0.01583	0.09905	0.02436	0.00893

Table 2. Table 2 . Reconstructed constants in the frequency profile.

•	real value	ellipse	square	near-boundary	small-central
${\kappa }_1$	3	2.80971	3.36482	3.00287	6.65418
${\kappa }_2$	2	1.79063	2.34197	1.96926	5.14671
${\kappa }_3$	1	1.00212	0.987247	0.999658	1.13223

Table 3. Table 3 . Relative symmetric difference.

•	ellipse	square	near-boundary	small-central
$\left|D_i△D_{target}\right|/\left|D_{target}\right|$	0.07055	0.12187	0.24299	0.19471

Equations

$$\displaystyle\left\{\begin{array}[]{lllcc}-\nabla\cdot\left(\sigma(x,\omega)\nabla u(x,\omega)\right)=0&\textrm{in}\quad\Omega,\\ \sigma(x,\omega)\partial_{\nu_{\Omega}}u(x,\omega)(x)=f(x)&\textrm{on}\quad\partial\Omega,\\ \int_{\partial\Omega}u(x,\omega)ds=0,&\end{array}\right.$$

$$\sigma (x,\omega )=k_0+(k(\omega )-k_0){\chi }_D(x)$$

$$k(\omega ):={\kappa }_1-\frac{{\kappa }_2}{{\omega }^2+i\omega {\kappa }_3},$$

$${\int }_{\Omega }\nabla Tu\cdot \nabla vdx={\int }_D\nabla u\cdot \nabla vdx.$$

$$\lambda {\int }_{\Omega }\nabla w\cdot \nabla vdx={\int }_D\nabla w\cdot \nabla vdx.$$

$$w{|}_{\partial D}^{+}=w{|}_{\partial D}^{-},{\partial }_{{\nu }_D}w{|}_{\partial D}^{+}=(1-\frac{1}{\lambda }){\partial }_{{\nu }_D}w{|}_{\partial D}^{-},{\partial }_{{\nu }_{\Omega }}w=0,$$

$$\Vert u{\Vert }_{{\mathfrak{H}}_{\diamond }}={\int }_{\Omega }|\nabla u{|}^{2}dx.$$

$$0<{\lambda }_1^{-}\le {\lambda }_2^{-}\le \cdots <\frac{1}{2},$$

$$1>{\lambda }_1^{+}\ge {\lambda }_2^{+}\ge \cdots >\frac{1}{2}.$$

$$T=\frac{1}{2}I+K,$$

$${\lambda }_1^{-}$$

$${\lambda }_n^{-}$$

$${\lambda }_1^{+}$$

$${\lambda }_n^{+}$$

$$u(x,\omega )=k_0^{-1}{u}_0(x)+\sum _{n=1}^{\infty }\frac{{\int }_{\partial \Omega }f(z)w_n^{\pm }(z)ds(z)}{k_0+{\lambda }_n^{\pm }(k(\omega )-k_0)}{w}_n^{\pm }(x),x\in \Omega ,$$

$$\displaystyle\left\{\begin{array}[]{lllcc}\Delta v=0&\textrm{in}\quad D^{\prime},\\ \nabla v=0&\textrm{in}\quad D,\\ \partial_{\nu_{\Omega}}v=f&\textrm{on}\quad\partial\Omega.\end{array}\right.$$

$$u_f=u-k_0^{-1}{u}_0,$$

$$u_f(x)=\sum _{n=1}^{\infty }{\int }_{\Omega }\nabla u_f(z)\nabla w_n^{\pm }(z)dz{w}_n^{\pm }(x),x\in \Omega .$$

$${\int }_{\Omega }\nabla u_f(z)\nabla w_n^{\pm }(z)dz$$

$${\int }_{\Omega }\nabla u(z)\nabla w_n^{\pm }(z)dz$$

$$\mathfrak{D}:=\left\{D=\left\{X_0+{\rm Y}(\theta )\left(\begin{array}{c}\cos \theta \\ \sin \theta \end{array}\right),\theta \in [0;2\pi )\right\};{\rm Y}\in \Pi \right\},$$

$$\Pi :=\left\{d_0<{\rm Y}(\theta )<d_1-\delta ;{\rm Y}(2\pi )={\rm Y}(0);\Vert {\rm Y}{\Vert }_{C^{\beta }}\le m,\beta \ge 2\right\}.$$

$$\partial D=\{x=X(t),t\in [a,b]\}.$$

$$\partial D_{\epsilon }=\{\tilde{x}=\tilde{X}(t):=X(t)+\epsilon h(X(t)){\nu }_D(X(t)),t\in [a,b]\}.$$

$$\left\{\begin{array}[]{lr}\bigtriangleup v=0&\text{in }\Omega\setminus\overline{D_{\varepsilon}},\\ \nabla v=0&\text{in }D_{\varepsilon},\\ \partial_{\widetilde{\nu}}v=f&\text{on }\partial\Omega,\\ \int_{\partial\Omega}vd\sigma=0,\end{array}\right.$$

$$u_{\epsilon }=u_0+\epsilon u_h+O({\epsilon }^2),\text{ as }\epsilon \to 0.$$

$$\left\{\begin{array}[]{lr}\bigtriangleup v=0&\text{in }\Omega\setminus\overline{D},\\ \nabla v=0&\text{in }D,\\ u|_{+}-u|_{-}=-h\partial_{\nu_{D}}u_{0}|_{+}&\text{on }\partial D,\\ \partial_{\widetilde{\nu}}v=0&\text{on }\partial\Omega,\\ \int_{\partial\Omega}vd\sigma=0.\end{array}\right.$$

$${\mathcal{S}}_D:H^{-\frac{1}{2}}(\partial D)\to H^{\frac{1}{2}}(\partial D),$$

$${\mathcal{S}}_D\phi (x)={\int }_{\partial D}G(x,z)\phi (z)ds(z),$$

$${\mathcal{S}}_{\Omega }:H^{-\frac{1}{2}}(\partial \Omega )\to H^{\frac{1}{2}}(\partial \Omega ),$$

$${\mathcal{S}}_{\Omega }\psi (x)={\int }_{\partial \Omega }G(x,z)\psi (z)ds(z),$$

Full text

1. The mfEIT Mathematical Model

Experimental research has found that the conductivity of many biological tissues varies strongly with respect to the frequency of the applied electric current within certain frequency ranges [GPG]. In [AGGJS], using homogenization techniques, the authors analytically exhibit the fundamental mechanisms underlying the fact that effective biological tissue electrical properties and their frequency dependence reflect the tissue composition and physiology. The multifrequency electrical impedance tomography (mfEIT) is a diffusive imaging modality that recovers the conductivity distribution of the tissue by using electrodes to measure the resulting voltage on its boundary, induced by two known injected currents and for many frequency values. The principal idea behind the (mfEIT) is that the dependance of the effective conductivity of the tissue with respect to the frequency of the electric current is extremely related to its state. In fact, its frequency dependence changes with its composition, membrane characteristics, intra-and extra-cellular fluids and other factors [AGGJS]. Therefore, the frequency dependence of the conductivity of the tissue can provide some information about the tissue microscopic structure and its physiological and pathological conditions. In other words, the frequency dependence of the conductivity of the tissue can help to determine if it is healthy or cancerous. The advantages of the (mfEIT) is canceling out errors due to boundary shape, the electrode positions, and other systematic errors that appear in -the more conventional imaging modality- electric impedance tomography (EIT) [Bor].

In the following we introduce the mathematical model of the (mfEIT). Let $\Omega$ be the open bounded smooth domain in $\mathbb{R}^{2}$, occupied by the sample under investigation and denote by $\partial\Omega$ its boundary. The mfEIT forward problem is to determine the potential $u(\cdot,\omega)\in H^{1}(\Omega):=\{v\in L^{2}(\Omega):\nabla v\in L^{2}(\Omega)\}$, solution to

[TABLE]

where $\omega$ denotes the frequency, $\nu_{\Omega}(x)$ is the outward normal vector to $\partial\Omega$, $\sigma(x,\omega)$ is the conductivity distribution, and $f\in H^{-\frac{1}{2}}_{\diamond}(\partial\Omega):=\{g\in H^{-\frac{1}{2}}(\partial\Omega):\int_{\partial\Omega}g\,ds=0\}$ is the input current.

In this work we are interested in the case where the frequency dependent conductivity distribution takes the form

[TABLE]

with $\chi_{D}(x)$ being the characteristic function of a $C^{2}$ domain $D$ in $\Omega$ ($\overline{D}\subset\Omega$), $k_{0}$ being a fixed strictly positive constant, and $k(\omega):\mathbb{R}_{+}\rightarrow\mathbb{C}\setminus\overline{\mathbb{R}_{-}}$, being a continuous complex-valued function.

Here $k_{0}$ represents the conductivity of the background medium, is known, and $k(\omega)$ is the conductivity of the biological tissue, given by the empirical model

[TABLE]

where $\kappa_{p}>0,\;p=1,2,3,$ are constants that only depend on the biological tissue properties (see for instance [AGGJS] ). The frequency profile $k(\omega)$ is somehow a meromorphic approximation with a single pole of the graph of experimental measurements for a given biological tissue [AGGJS]. It also appears as a homogenized model for periodically distributed biological cells in the dilute limit [AGGJS], and is similar to Drude models that describes the frequency dependence of the electric permittivity of a real metal within the visible frequency range [MFZ].

The mfEIT inverse problem is to determine the anomaly $D$ and the characteristics $\kappa_{p},\;p=1,2,3,$ of the biological tissue from measurements of the boundary voltages $u(x,\omega)$ on $\partial\Omega$, for $\omega\in(\underline{\omega},\overline{\omega})$, $0\leq\underline{\omega}<\overline{\omega}$.

There have been several numerical approaches on multifrequency electrical impedance tomography. Most of them are dealing with small-volume anomalies [ABG, ABGW, GH] and frequency-difference imaging [JS, MSHA]. In small-volume-volume imaging, only the location and the multi-frequency polarization tensor can be reconstructed from boundary measurements. In frequency-difference imaging, the main idea is to compare the images for different frequencies, and consider only the frequency dependent part. It was numerically shown that the approach can accommodate geometrical errors, including imperfectly known boundary. This approach which seems more natural since it aims to identify the anomaly by only focusing on the changes in the resulting images for different frequencies, is somehow simultaneously complementary and opposite to our analysis in this paper. Taking the difference between two images associated to different frequencies will remove the frequency independent part which is the keystone of the identification path pursued in this paper. Our strategy is based on the plasmonic spectral decomposition derived in [AT] which splits the electric potential on the boundary $\partial\Omega$ into two parts $u=k_{0}^{-1}u_{0}+u_{f}$, and separate between the frequency dependent and independent parts; see also [ADM, AMRZ]. In fact the part $k_{0}^{-1}u_{0}$ corresponds to the response of the same anomaly filled with a perfect conductor, that is $k_{0}^{-1}u_{0}$ is the limit of $u$ when $k$ tends to infinity. Precisely, in [AT], it was proven that the convergence of $u$ to $k_{0}^{-1}u_{0}$ is linear in $1/k$. We first process algebraically the data on the boundary and recover the frequency dependent part in order to acquire the Cauchy data of frequency independent part $u_{0}$. Then, based on known results and approaches on the reconstruction of zero level set of harmonic functions from Cauchy data we determine the anomaly itself. The other approach is based on perturbation techniques, suppose that the contrast $k(\omega)/k_{0}$ is close to $1$, and linearize the problem around the harmonic function in the whole domain that shares the same flux on the boundary [AAJS]. Then, the sparsity issue comes into play and help to speed up the convergence of the iterative algorithm. In our approach, the notion of sparsity appears naturally in post-processing the data on the boundary to recover the frequency independent part. Precisely, it seems that only a finite number of eigenfunctions intervene in the linear inversion, and this can be completely mastered by the shape of the anomaly and its distance to the boundary where the measurement are taken. Finally, the perturbation approach is also complementary to our analysis since the higher the frequency is the better the recovery of the frequency independent part on the boundary is.

The paper is organized as follows. In section 2 we provide the spectral decomposition derived in [AT]. The linearization of the frequency independent part with respect to the shape of the anomaly which is necessary in our identification approach is studied in section 3. Section 4 is devoted to the retrieval of the frequency independent part on the boundary. Here, we will not follow the theoretical approach developed in [AT] based on the unique continuation of meromorphic functions. We solve the problem using algebraic tools under simplification assumptions inspired by the sparsity properties of the conductivity distribution and the behavior of the eigenvalues near the unique accumulation point $1/2$. We also determine in the sequel the profile constants $\kappa_{i},\,i=1,2,3$. In section 5, one we have the Cauchy data of the frequency independent part $u_{0}$ on the measurement boundary we use a conventional optimization technique to recover the anomaly. Several numerical examples are presented in section 6. Comments on the obtained results and future directions are given in the conclusion section 7.

2. Spectral decomposition of $u(x,\omega)$

We first introduce an operator whose spectral decomposition will be later the corner stone of the identification of the anomaly $D$. Let $H^{1}_{\diamond}(\Omega)$ be the space of functions $v$ in $H^{1}(\Omega)$ satisfying $\int_{\partial\Omega}vds=0$.

For $u\in H^{1}_{\diamond}(\Omega)$, we infer from the Riesz theorem that there exists a unique function $Tu\in H^{1}_{\diamond}(\Omega)$ such that for all $v\in H^{1}_{\diamond}(\Omega)$,

[TABLE]

The variational Poincaré operator $T:H^{1}_{\diamond}(\Omega)\rightarrow H^{1}_{\diamond}(\Omega)$ is easily seen to be self-adjoint and bounded with norm $\|T\|\leq 1$.

The spectral problem for $T$ reads as: Find $(\lambda,w)\in\mathbb{R}\times H^{1}_{\diamond}(\Omega)$, $w\not=0$ such that $\forall v\in H^{1}_{\diamond}(\Omega)$,

[TABLE]

Integrating by parts, one immediately obtains that any eigenfunction $w$ is harmonic in $D$ and in $D^{\prime}=\Omega\setminus\overline{D}$, and satisfies the transmission and boundary conditions

[TABLE]

where $w|^{\pm}_{\partial D}(x)=\lim_{t\rightarrow 0}w(x\pm t\nu_{D}(x))$ for $x\in\partial D$. In other words, $w$ is a solution to (4) for $k=k_{0}(1-\frac{1}{\lambda})$ and $f=0$.

Let $\mathfrak{H}_{\diamond}$ the space of harmonic functions in $D$ and $D^{\prime}$, with zero mean $\int_{\partial\Omega}uds(x)=0$, and zero normal derivative $\partial_{\nu_{\Omega}}u=0$ on $\partial\Omega$, and with finite energy semi-norm

[TABLE]

Since the functions in $\mathfrak{H}_{\diamond}$ are harmonic in $D^{\prime}$, the $\mathfrak{H}_{\diamond}$ is a closed subspace of $H^{1}(\Omega)$. Later on, we will give a new characterization of the space $\mathfrak{H}_{\diamond}$ in terms of the single layer potential on $\partial D$ associated with the Neumann function of $\Omega$.

We remark that $Tu=0$ for all $u$ in $H^{1}_{0}(D^{\prime})$, and $Tu=u$ for all $u$ in $H^{1}_{0}(D)$ (the set of functions in $H^{1}(D)$ with trace zero).

We also remark that $T\mathfrak{H}_{\diamond}\subset\mathfrak{H}_{\diamond}$ and hence the restriction of $T$ to $\mathfrak{H}_{\diamond}$ defines a linear bounded operator. Since we are interested in harmonic functions in $D$ and $D^{\prime}=\Omega\setminus\overline{D}$, we only consider the action of $T$ on the closed space $\mathfrak{H}_{\diamond}$. We further keep the notation $T$ for the restriction of $T$ to $\mathfrak{H}_{\diamond}$. We will prove later that $T$ has only isolated eigenvalues with an accumulation point $1/2$. We denote by $\left(\lambda_{n}^{-}\right)_{n\geq 1}$ the eigenvalues of $T$ repeated according to their multiplicity, and ordered as follows

[TABLE]

in $(0,1/2]$ and, similarly,

[TABLE]

the eigenvalues in $[1/2,1)$. The eigenvalue $\lambda_{\infty}=1/2$ is the unique accumulation point of the spectrum. To ease the notation we further denote the orthogonal spectral projector on the eigenspace associated to $1/2$, by $\int_{\partial\Omega}\cdot w_{\infty}^{\pm}(z)ds(z)w_{\infty}^{\pm}(x).$ Next, we will characterize the spectrum of $T$ via the mini-max principle.

Proposition 2.1**.**

The variational Poincaré operator has the following decomposition

[TABLE]

where $K$ is a compact self-adjoint operator. Let $w_{n}^{\pm},\;n\geq 1$ be the eigenfunctions associated to the eigenvalues $\left(\lambda_{n}^{\pm}\right)_{n\geq 0}$. Then

[TABLE]

and similarly

[TABLE]

We have the following decomposition of $u(x,\omega)$ in the basis of the eigenfunctions of the variational Poincaré operator $T$.

Theorem 2.1**.**

[AT*]** Let $u(x,\omega)$ be the unique solution to the system (4).

Then, the following decomposition holds:

[TABLE]

where $u_{0}(x)\in H^{1}_{\diamond}(\Omega)$ depends only on $f$ and $D$, and is the unique solution to

[TABLE]

Proof.

In order to have a self-contained document we give the proof of the theorem.

We first observe that frequency dependent part

[TABLE]

lies in $\mathfrak{H}_{\diamond}$. Since the eigenfunctions $w^{\pm}(x)$ form an orthonormal basis of $\mathfrak{H}_{\diamond}$, the frequency part $u_{f}$ posses the following spectral decomposition:

[TABLE]

A forward computation leads to

[TABLE]

On the other hand, since $u\in H_{\diamond}^{1}(\Omega)$, we obtain

[TABLE]

Using the simple fact that

[TABLE]

we obtain the desired decomposition.

∎

In [AT], assuming that the profile $k(\omega)$ is given, the spectral decomposition (9) of the solution of the forward conductivity problem has been used to retrieve the Cauchy data corresponding to the extreme case of perfect conductor $u_{0}$. Based on unique continuation techniques, the uniqueness of the mfEIT problem, and rigorous stability estimates have been obtained from the knowledge of $u_{0}|_{\partial\Omega}$, in the case where the anomaly is within a class of star shaped domains.

Assume that $X_{0}\in\Omega$, and let $d_{1}=\textrm{dist}(X_{0},\partial\Omega)$ and let $d_{0}<d_{1}$. For $\delta>0$ small enough, and $m>0$ large enough, define the set of anomalies:

[TABLE]

where

[TABLE]

In this paper we consider the reconstruction of the profile function $k(\omega)$ defined in (5), and anomalies $D$ within the set $\mathfrak{D}$. At first glance, the numerical reconstruction of the inclusion in the mfEIT problem does not need to follow the path of the theoretical results derived in [AT], that is, to determine first $u_{0}|_{\partial\Omega}$, and then find the high conductor anomaly $D$ that produces the recovered Cauchy data of $u_{0}$. In fact preliminary numerical calculations show that a blind minimization approach that searches the anomaly $D$ and the profile $k(\omega)$ using boundary multifrequency data does not converge in most cases, and if it happens to converge the rate turns out to be very slow. These difficulties are well known in inverse conductivity problem, usually it is very hard to distinguish between the conductivity value and the size of the anomaly [AK]. On the other hand the numerical identification of a high conductor anomaly is a well known inverse problem, and many works have been done on it (see for instance [KS, LL]). We can cite for example quasi-reversibility type based methods [KS, LL]. Here, we will consider the parameterization type based methods [CK, ACLZ, Ru]. Since the problem is ill-posed we will use a cut-off regularization approach that consists on taking into account in the computation only the important Fourier modes of the parametrization function $\Upsilon(\theta)$. Then, the identification is transformed into an optimization problem with a finite number of degree of freedom. The problem is still strongly nonlinear we propose here to solve it using the gradient adjoint method.

The algorithm we propose in this paper for identifying numerically the anomaly and the frequency conductivity profile is inspired by the theoretical approach developed in[AT], and it can be summarized as follows:

(i)To recover $u_{0}(x)|_{\partial\Omega}$ and $\kappa_{p},\;p=1,2,3,$ from the knowledge of $u(x,\omega)|_{\partial\Omega},\;\omega\in(\underline{\omega},\overline{\omega})$. Here, we will use an linear algebraic approach based on the understanding of the behavior of the spectrum of Neumann-Poincaré operator near its unique accumulation point, and the sparsity of the considered conductivity distribution.

(ii) To identify the anomaly $D$ from the Cauchy data $(u_{0}(x),f(x))$ on the boundary $\partial\Omega$ using a cut-off parameterization/Fourier approach. Here to further stabilize and speed up the convergence of the iterative gradient based method we will use two linearly independent boundary currents $(f_{1},f_{2})$.

3. The linearized map

We shall use gradient methods to identify the anomaly from the Cauchy data $(u_{0}(x),f(x))$ on the boundary $\partial\Omega$. In this section, we determine $u_{h}$, the derivative of $u_{0}$ with respect to a shape perturbation in the direction $h(x)\nu_{D}(x)$, where $h(x)$ is a scalar function defined on $\partial D$. We will follow the analysis in [AKLZ] for a non-degenerate conductivity inside the anomaly, based on integral equations techniques.

Let $D\in\mathfrak{D}$ be a given anomaly. We define $X(t):[a,b]\rightarrow\mathbb{R}^{2}$ to be a smooth clockwise parametrization of $\partial D$, where $a,b\in\mathbb{R}$, $a<b$. We assume that $X\in\mathcal{C}^{\beta}([a,b])$ and $|X^{\prime}(t)|=1$ for all $t\in[a,b]$. Then

[TABLE]

Let $h\in C^{2}(\partial D)$, and define the boundary of the perturbed anomaly $D_{\varepsilon}$ by

[TABLE]

Define $u_{\varepsilon}$ to be the unique solution to the system (13) associated to the perturbed anomaly $D_{\varepsilon}$, that is

[TABLE]

where $\widetilde{\nu}$ is the outward normal vector on $\partial D_{\varepsilon}$.

The objective of this section is to derive a linear correction $u_{h}$ of $u_{\varepsilon}$, such that

[TABLE]

The main result of this section is the following.

Theorem 3.1**.**

Let $h\in C^{2}(\partial D)$ be fixed. Then, $u_{h}$ is the unique solution to the system

[TABLE]

Proof.

We first derive an integral equation representation of the field.

Let $G(x,z)=\frac{1}{2\pi}\log(|x-y|)$ be the Green function for the Laplacian in $\mathbb{R}^{2}$, and define the single layer potentials respectively on $\partial D$, and $\partial\Omega$ by

[TABLE]

and

[TABLE]

They satisfy the following jump relations through the boundary of respectively $D$ and $\Omega$ [AK]

[TABLE]

where

[TABLE]

and

[TABLE]

are compact operators.

Since $G(x,z)$ is harmonic in $\mathbb{R}^{2}\setminus\{z\}$, $S_{\Omega}\psi$ has a unique harmonic extension in $\mathbb{R}^{2}\setminus\partial\Omega$. Similarly, $S_{D}\varphi$ has a unique harmonic extension in $\mathbb{R}^{2}\setminus\partial D$. In addition, we have

[TABLE]

We also define the double layer potential $\mathcal{D}_{D}$

[TABLE]

It satisfies the following jump relations

[TABLE]

where

[TABLE]

is the $L^{2}(\partial D)$-adjoint of $\mathcal{K}_{D}^{*}$.

The solution $u_{0}$ can be written as

[TABLE]

where $\phi=\partial_{\nu_{D}}u_{0}|_{+}\in H^{-\frac{1}{2}}_{\diamond}(\partial D)$ and $\psi\in H^{-\frac{1}{2}}(\partial\Omega)$.

Similarly, we have

[TABLE]

where $\phi_{\varepsilon}=\partial_{\widetilde{\nu}}u_{\varepsilon}|_{+}\in H^{-\frac{1}{2}}_{\diamond}(\partial D_{\varepsilon})$ and $\psi_{\varepsilon}\in H^{-\frac{1}{2}}(\partial\Omega)$.

Using the jump relations and the facts that $\partial_{\nu_{D}}u_{0}|_{-}=0$ on $\partial D$, and $\partial_{\nu_{\Omega}}u_{0}|_{-}=f$ on $\partial\Omega$, the densities $\phi$ and $\psi$ satisfy the following system

[TABLE]

This system can be also represented in a matrix form

[TABLE]

The same analysis leads to the system

[TABLE]

From the parameterization of $\partial D$, we deduce that the outward unit normal vector $\nu_{D}(x)$ is given by $\nu_{D}(X(t))=R_{-\frac{\pi}{2}}T(X(t))$, where $R_{-\frac{\pi}{2}}$ is the rotation with the angle $-\frac{\pi}{2}$, and $T(X(t))=X^{\prime}(t)$ is the tangential normal vector. Let $\gamma(X(t))$ be the curvature, it satisfies

[TABLE]

Using the parameterization of $\partial D_{\varepsilon}$, we deduce the following asymptotic expansion

[TABLE]

where $h^{\prime}(t)=\frac{d}{dt}h(X(t))$ (we also use $h^{\prime}(x)$ to denote this quantity). In the same way obtain the asymptotic expansion of the length element

[TABLE]

Let $\Psi_{\varepsilon}$ be the diffeomorphism from $\partial D$ onto $\partial D_{\varepsilon}$ given by $\Psi_{\varepsilon}(x)=x+\varepsilon h(x)\nu(x)$. From [AKLZ], we deduce the asymptotic expansion of $\mathcal{K}^{*}_{D_{\varepsilon}}$

[TABLE]

where $\widetilde{\phi}=\phi\circ\Psi^{-1}_{\varepsilon}$, and the operator $\mathcal{K}^{(1)}_{D}$ is defined by

[TABLE]

Now, we calculate the asymptotic expansion of the operators $\partial_{\widetilde{\nu}}\mathcal{S}_{\Omega}$ on $\partial D_{\varepsilon}$ and $\partial_{\nu_{\Omega}}\mathcal{S}_{D}$ on $\partial\Omega$.

Let $\psi\in H^{-\frac{1}{2}}(\partial\Omega)$ be fixed. Using the relation $\widetilde{x}=x+\varepsilon h(x)\nu(x)\in\partial D_{\varepsilon}$ for $x\in\partial D$, we obtain

[TABLE]

where $\partial_{T}$ denotes the tangential derivative, and $\mathcal{S}^{(1)}_{\Omega}$ is defined by

[TABLE]

for $x\in\partial D$.

We further determine the asymptotic expansion of $\partial_{\nu_{\Omega}}\mathcal{S}_{D_{\varepsilon}}$ on $\partial\Omega$. Let $\phi\in H^{-\frac{1}{2}}(\partial D)$, and $x\in\partial\Omega$, we have

[TABLE]

Consequently

[TABLE]

where the operator $M_{h}$ on $H^{-\frac{1}{2}}_{\diamond}(\partial D)\times H^{-\frac{1}{2}}(\partial\Omega)$ is defined by

[TABLE]

So, the systems (22) and (23) imply

[TABLE]

where $\begin{pmatrix}\phi_{h}\\ \psi_{h}\end{pmatrix}$ is given by

[TABLE]

Thus, using the representation formula and following the same calculus, we determine the asymptotic expansion of the solution $u_{\varepsilon}|_{\partial\Omega}$

[TABLE]

We further denote

[TABLE]

We deduce from (33)

[TABLE]

on $\partial D$, and

[TABLE]

on $\partial\Omega$, and hence

[TABLE]

The equality (18) and the fact that $\partial_{T}u_{0}=0$ on $\partial D$, lead to

[TABLE]

which implies

[TABLE]

A similar calculus gives

[TABLE]

for $x\in\partial D$.

Thus

[TABLE]

By the continuity of the normal derivative of double layer potentials and the jump relation, we have, for $x\in\partial D$,

[TABLE]

Using (35), (3) and (3), we have,

[TABLE]

which gives the desired result. ∎

4. Reconstruction of $u_{0}(x)|_{\partial\Omega}$ and $k(\omega)$

In this section we construct $u_{0}|_{\partial\Omega}$ from the knowledge of $u(x,\omega)|_{\partial\Omega},\;\omega\in(\underline{\omega},\overline{\omega})$. Here we recall the spectral decomposition (9), also valid on the boundary $\partial\Omega$, which is the keystone of our approach.

[TABLE]

The first observation is that the functions $w_{n}^{\pm}(x)$ do not need to be orthogonal on the boundary $\partial\Omega$. Then, varying the frequency, and so the coefficients of the expansion above do not guarantee the complete separation between the frequency and the non frequency parts. The second observation is that the simultaneous determination of the plasmonic resonances $\lambda_{n}^{\pm},\,n\geq 1$, the frequency profile $k(\omega)$, and $u_{0}|_{\partial\Omega}$ is strongly nonlinear while if we assume that $k(\omega)$ and $\lambda_{n}^{\pm},\,n\geq 1$ are given the problem becomes a linear one.

We further consider $M\geq 2$ frequencies $\omega_{1},\cdots,\omega_{M}$ in $(\underline{\omega},\overline{\omega})$, and their associated solutions $u(x,\omega_{1}),\cdots,u(x,\omega_{M})$. Since $1/2$ is the unique accumulation point of the eigenvalues $(\lambda_{n}^{\pm})_{n\geq 1}$, we only consider the $N_{f}\geq 0$ first eigenvalues as unknown variables, and we approximate the others eigenvalues by the limiting value $1/2$. In fact it has been shown in [MS] that if $D$ is $C^{\beta}$ with $\beta\geq 2$ then for any $\alpha>-2\beta+3,$ we have

[TABLE]

Thus the boundary regularity is essential to the decay rate of eigenvalues. Consequently if the boundary is $C^{\infty}$ smooth, then the plasmonic eigenvalues will decay faster than any power order. Recently H. Kang and his collaborators have proved the exponential convergence of the eigenvalues in the case of analytic anomalies [AKM]. This theoretical work justifies the exponential decay behavior that has been observed numerically [PP], and checked for sample geometries like ellipses. If $\alpha>0$ is the modified maximal Grauert radius of $\partial D$, then

[TABLE]

These asymptotic properties of the spectrum of the Neumann-Poincaré operator suggest to consider only a finite number of them in the spectral decomposition (9). We further make the following approximation for $x\in\Omega$, $1\leq p\leq M:$

[TABLE]

where

[TABLE]

A simple integration by parts, leads for all $n\geq 1$

[TABLE]

where $\mathfrak{f}$ is the unique solution in $H^{1}_{\diamond}(\Omega)$ to

[TABLE]

Consequently, the function

[TABLE]

where $\mathcal{P}_{0}$ is the orthogonal projection onto the space $\mathfrak{H}_{\diamond}$. On the other hand, $u_{0}$ satisfies

[TABLE]

for all $n\geq 1$.

Since $\mathcal{P}_{0}u_{0}=0$, and $\mathfrak{f}-u_{0}\in\mathfrak{H}_{\diamond}$ the orthogonal projection of $\mathfrak{f}$ onto the space $\mathfrak{H}_{\diamond}$ is $\mathfrak{f}-u_{0}$, that is $\mathfrak{f}-u_{0}=\mathcal{P}_{0}(\mathfrak{f}-u_{0})=\mathcal{P}_{0}\mathfrak{f}.$

Therefore, the formula (40) becomes

[TABLE]

Next, we reconstruct $\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$ and $u_{0}(x)$ by an optimization algorithm. In order to do so, we need an apriori estimation of the eigenvalues $\widetilde{\lambda_{n}^{\pm}}\in[0,1]$ for $n=1,\cdots,N_{f}$. Since the eigenvalues $\lambda_{n}^{\pm}$ are within a relative narrow interval, preliminary calculations showed that the reconstruction of $u_{0}$ is indeed not very sensitive to the choice of those eigenvalues. In the rest of this section, we assume that we have we fix $\widetilde{\lambda_{n}^{\pm}}\in[0,1]$ for $n=1,\cdots,N_{f}$.

Let $(x_{j})_{1\leq j\leq N_{d}}\in\partial\Omega$ a discretization of the boundary $\partial\Omega$, and define, for $n=1,\cdots,N_{f}$ the scalar functionals

[TABLE]

where $(U_{0}^{(j)})_{1\leq j\leq N_{d}}$ and $(V_{n}^{\pm(j)})_{1\leq j\leq N_{d}}$ are vectors in $\mathbb{R}^{N_{d}}$, that approximate respectively $(u_{0}(x_{j}))_{1\leq j\leq N_{d}}$ and $(v_{n}(x_{j}))_{1\leq j\leq N_{d}}$.

The scheme consists in minimizing the scalar functional

[TABLE]

So, we can easily calculate its gradient from (4) and (6), for $i=1,2,3$, $1\leq l\leq N_{d}$ and $1\leq n\leq N_{f}$,

[TABLE]

we denote here $F_{j}(U_{0}^{(j)},V^{\pm(j)}_{1},\cdots,V^{\pm(j)}_{N_{f}},\omega_{p},\kappa_{1},\kappa_{2},\kappa_{3})$ by $F_{j}(\cdot,\omega_{p})$ in order to simplify the notations.

Then, the algorithm follows the standard gradient method for $3+N_{d}(1+2N_{f})$ variables. Once we have reconstructed the conductivity profile, i.e. the approximate values of $\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$, we can use (4) again to calculate the approximate conductivity $\widetilde{k}(\omega)$ by (6) and the approximate $u_{0}$ by the following matrix formula, let $x\in\partial\Omega$,

[TABLE]

where $\widetilde{u}(x,\omega)=u(x,\omega)-\frac{2}{\widetilde{k}(\omega)+k_{0}}\mathfrak{f}(x)$, $q_{0}(\omega)=\frac{\widetilde{k}(\omega)-k_{0}}{k_{0}(\widetilde{k}(\omega)+k_{0})}$, and $q(\widetilde{\lambda},\omega)=\frac{1}{k_{0}+\widetilde{\lambda}(\widetilde{k}(\omega)-k_{0})}-\frac{2}{\widetilde{k}(\omega)+k_{0}}$. Then, the vector $V$ can be obtained by the formula

[TABLE]

where $(L^{T}L)^{{\dagger}}$ is the pseudo-inverse of the matrix $L^{T}L$. The conditioning of the matrix $L^{T}L$ depends in fact on the distance between the sampling values $\omega_{j},\,j=1\,\cdots,M,$ and the frequency profile $\eqref{empirical}$. The approximate $u_{0}(x)$ is then recovered by taking the first coefficient of the vector $V(x)$.

Finally, the algorithm to reconstruct $u_{0}$, can be summarized in the following steps:

1. (1)

Give an apriori estimation $\widetilde{\lambda_{n}^{\pm}},\,n=1,\cdots,N_{f}$, of the eigenvalues $\lambda_{n}^{\pm},\,n=1,\cdots,N_{f}$. 2. (2)

Choose a step length $\alpha_{m}>0$ for the gradient descent. 3. (3)

Initialize the vectors $U_{0}|_{0}$, $V_{1}|_{0},\cdots,V_{n}|_{0}$ and the coefficients $\kappa_{1}|_{0}$, $\kappa_{2}|_{0}$, $\kappa_{3}|_{0}$. 4. (4)

While $|\nabla J_{m}|$ is larger then a given threshold, we do

1. (a)

Calculate the values of the functions $F_{j}$ by (4), and $\nabla J_{m}$ by (46), (47), (48). 2. (b)

Update the parameters $\kappa_{i}|_{k+1}=\kappa_{i}|_{k}-\alpha_{m}\frac{\partial J_{m}}{\partial\kappa_{i}}$, $U_{0}^{(l)}|_{k+1}=U_{0}^{(l)}|_{k}-\alpha_{m}\frac{\partial J_{m}}{\partial U_{0}^{(l)}}$, and $V_{n}^{\pm(l)}|_{k+1}=V_{n}^{\pm(l)}|_{k}-\alpha_{m}\frac{\partial J_{m}}{\partial V_{n}^{\pm(l)}}$. 5. (5)

When $|\nabla J_{m}|$ is smaller then the threshold, we stop the iterations. 6. (6)

Use (50) with the approximate coefficients $\kappa_{i}$ obtained in the previous step to calculate the approximate value of $u_{0}(x)$ for every $x\in\partial\Omega$.

5. Reconstruction of the anomaly from $u_{0}$

In this section, we propose a numerical method to identify the anomaly $D$ from a finite number of Cauchy data of $(u_{0}(f_{i}),f_{i}),\,i=1,\cdots,P,$ on $\partial\Omega$, where $P\geq 1$. We further assume that the anomaly is located within an open subdomain $\Omega_{0}\subset\Omega$ with $dist(\partial\Omega_{0},\partial\Omega)\geq\delta_{0}>0$ The scheme is based on the minimizing of a non convex functional

[TABLE]

where $u_{meas}^{(i)}$ are the measured Dirichlet data corresponding to the $i$-th Neumann data and where $u$ is the solution to (4) associated to the current domain $D\subset\Omega_{0}$. In our numerical simulations we take $P=2$ with $f_{1}=\langle e_{1},\nu_{\Omega}\rangle$ and $f_{2}=\langle e_{2},\nu_{\Omega}\rangle$, where $(e_{1},e_{2})$ is the canonical base of $\mathbb{R}^{2}$.

We further assume that $D$ is within the class $\mathfrak{D}$, that is, it is star shaped and its boundary $\partial D$ can be described by the Fourier series:

[TABLE]

where $C=\begin{pmatrix}c_{-N}\\ c_{-N+1}\\ \vdots\\ c_{N}\end{pmatrix}\in\mathbb{R}^{2N+1}$, $f_{n}(\theta)=\cos(n\theta)$ for $0\leq n\leq N$ and $f_{n}(\theta)=\sin(n\theta)$ for $-N\leq n<0$.

Using (17) in Theorem 3.1, and integration by parts, we have the expressions of the shape derivative corresponding to each Fourier coefficient $c_{n}$

[TABLE]

for $-N\leq n\leq N$, where $h(\theta)=f_{n}(\theta)\langle\begin{pmatrix}\cos\theta\\ \sin\theta\end{pmatrix},\nu_{D}\rangle,$ and $w$ is the solution of the following equation

[TABLE]

Formula (52) is also valid for the shape derivative corresponding to the displacement of $X_{0}$, in these cases, $h=\langle e_{i},\nu_{D}\rangle$, $i=1,2$.

Those expressions are the basis of the following iterative algorithm:

1. (1)

Choose an initial domain $D_{0}$. 2. (2)

For each iteration, $i>0$:

1. (a)

Calculate the solution $u_{i}$ to (4), associated to the domain $D_{i}$ for which the boundary $\partial D_{i}$ is calculated by (51). 2. (b)

Calculate the shape derivatives $\frac{\partial J}{\partial x_{1}}$, $\frac{\partial J}{\partial x_{2}}$ and $\frac{\partial J}{\partial c_{n}}$ for all $-N\leq n\leq N$. 3. (c)

Choose a step length $\alpha>0$ for the gradient descent. 4. (d)

Update the parameters of the domain $X_{i+1}=X_{i}-\alpha\nabla_{X_{0}}J(X_{i},C_{i})$ and $C_{i+1}=C_{i}-\alpha\nabla_{C}J(X_{i},C_{i})$ with $\alpha>0$. 5. (e)

If the updated domain is not entirely in $\Omega_{0}$ or if $R$ becomes negative, reduce the size of $\alpha$. 3. (3)

When $J(X_{i},C_{i})$ becomes smaller than a fixed threshold, we stop.

6. Numerical examples

The numerical tests follow the steps presented here. All the numerical experiments are done using FreeFem++ [FreeFem].

1. (1)

$\Omega$ is a centered ellipse defined by the equation: $\frac{x_{1}^{2}}{4^{2}}+\frac{x_{2}^{2}}{3^{2}}\leq 1$. 2. (2)

We use two linearly independent Neumann data: $f_{1}=\langle e_{1},\nu_{\Omega}\rangle$ and $f_{2}=\langle e_{2},\nu_{\Omega}\rangle$, where $(e_{1},e_{2})$ is the canonical base of $\mathbb{R}^{2}$. 3. (3)

The multifrequence conductivity follows the model (6) with $\kappa_{1}=3$, $\kappa_{2}=2$, $\kappa_{3}=1$ and $\omega$ are integers from $1$ to $8$. 4. (4)

Only the first two eigenvalues are taken into consideration, and they are fixed as follows $\lambda^{+}_{1}=\frac{3}{4}$, $\lambda^{-}_{1}=\frac{1}{4}$ respectively in all cases. 5. (5)

In the algorithm to reconstruct $u_{0}$ and the conductivity profile, the initial guess of $u_{0}$ is the function $\mathfrak{f}$. 6. (6)

The initial estimation of domain $D$ is a centered disk with a radius $\frac{1}{2}$. 7. (7)

We consider the first $15$ Fourier coefficients: $N=15$. 8. (8)

We use P1 finite elements for the numerical resolution of the PDEs. 9. (9)

At each iteration, we remesh the domain to adapt to the new predicted position and shape of the domain. 10. (10)

The algorithms stop if $J<10^{-5}$ or the number of iterations exceed $500$. All the tests have executed $500$ iterations.

We present here several numerical simulations of the proposed algorithm. We first give the errors in the reconstruction method of $u_{0}$ in Table (1), and the errors in the reconstructed coefficients $\kappa_{1}$, $\kappa_{2}$, $\kappa_{3}$ in Table (2). The errors are computed using the $L^{2}$-norm of the difference

[TABLE]

We show in the following figures the targets and the reconstruction result. We calculate also the relative symmetric difference $|D_{i}\bigtriangleup D_{target}|/|D_{target}|$ during the iterations, and we draw the curves of the symmetric difference with respect to $\log(u_{i})$. We finally give the relative symmetric difference of each shape in Table 3. Finally, we test a reconstruction in domain $\Omega$ that has shape different from an ellipse in Figure (5).

7. Concluding remarks

In this paper, by combining the spectral decomposition derived in [AT] and the linearization of the frequency independent part with respect to the shape of the anomaly, we have provided a new and efficient approach for reconstructing both the shape and conductivity parameter of a conductivity anomaly from multifrequency boundary voltage measurements. The approach and results of this paper can be extended in several directions: (i) to reconstruct multiple anomalies from multifrequency boundary measurements; (ii) to investigate the reconstruction of anisotropic conductivity anomalies from multifrequency boundary measurements, and (iii) to study elastography imaging of visco-elastic anomalies. These new developments will be reported in forthcoming works.

8. Acknowledgments

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX- 0025-01).