Data Fusion for Electromagnetic and Electrical Resistive Tomography Based on Maximum Likelihood

Abstract

This paper presents a maximum likelihood based approach to data fusion for electromagnetic (EM) and electrical resistive (ER) tomography. The statistical maximum likelihood criterion is closely linked to the additive Fisher information measure, and it facilitates an appropriate weighting of the measurement data which can be useful with multiphysics inverse problems. The Fisher information is particularly useful for inverse problems which can be linearized similar to the Born approximation. In this paper, a proper scalar product is defined for the measurements and a truncated Singular Value Decomposition (SVD) based algorithm is devised which combines the measurement data of the two imaging modalities in a way that is optimal in the sense of maximum likelihood. As a multiphysics problem formulation with applications in geophysics, the problem of tunnel detection based on EM and ER tomography is studied in this paper. To illustrate the connection between the Green's functions, the gradients and the Fisher information, two simple and generic forward models are described in detail regarding two-dimensional EM and ER tomography, respectively.

1. Introduction

Nondestructive monitoring systems based on information, communication, and sensor technologies will be used to provide future emergency and disasters stakeholders with high situation awareness by means of realtime and detailed information and images of the infrastructure status [1]. Important technologies include electrical resistive tomography (ERT) as well as ground penetrating radar (GPR), which have become widely applied to obtain images of the subsurface in areas of complex geology [2–4]. The development and implementation of complex ICT-based monitoring systems will also rely on new technologies, techniques, and algorithms including the integration and correlation of the electromagnetic properties corresponding to imaging modalities such as ERT and GPR [1].

Statistical maximum likelihood (ML) or Bayesian approaches towards multisensor data fusion have been considered previously in, for example, [5–8]. In this paper, the general approach to information fusion is to perform a multiphysics data fusion based on the maximum likelihood principle, and to exploit the elements of Fisher information analysis that has been developed in [9–11]. In particular, an infinite-dimensional Fisher information formulation is employed which is suitable for analytical Green's function and gradient-based approaches [11], as well as with adjoint field formulations [10].

The Fisher information is a local measure of information with respect to the state space of parameter values, and it is therefore particularly useful for inverse problems which can be linearized, such as with the Born approximation [3], and so forth. In this paper, it is shown how the principle of maximum likelihood (under the assumption of Gaussian noise) can be used to derive a proper scalar product for the measurement data. By using this scalar product, it is shown that the truncated pseudoinverse based on the singular value decomposition (SVD) of the multiphysics forward model is equivalent to a one-step Newton method based on the principle of maximum likelihood together with a truncated SVD of the Fisher information. Hence, a linearized and regularized inversion based on the proposed scalar product yields a weighting of the measurement data which is statistically optimal and which is able to deal with a diversity of multiphysics forward problems, a combination of complex valued and real valued data, and a diversity of measurement sensors and qualities, and so forth.

As a generic example concerning a multiphysics inverse problem based on geophysical sensing, this paper addresses the problem of tunnel detection [4] based on data fusion with electromagnetic (ER) and electrical resistive (ER) tomography. A two-dimensional problem is considered with radio frequency as well as electrical sensors placed horizontally above a tunnel in the ground. A Green's function approach is used to obtain the gradients for the Fisher information analysis [11] as well as for the related SVD-based inversion algorithm. The numerical examples demonstrate the potential impact of an imbalance between the singular values and the variance of the measurement noise when different imaging modalities are incorporated in the inversion. The examples furthermore illustrate the significance of taking a statistically based weighting of the measurement data into proper account.

2. Multiphysics Data Fusion Based on Maximum Likelihood

In this section multiphysics data fusion is discussed. Section 2.1 settle the statistical model and the gradients. This is then used in Section 2.2 to describe data fusion with SVD and ML approach.

2.1. Statistical Model and Gradients for the Multiphysics Problem

Consider the multiphysics inverse problem of combining electromagnetic (EM) and electrical resistive (ER) tomography based on the following physical/statistical measurement model𝜉⁽¹⁾_𝑚𝑛𝑞=𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)+𝑤⁽¹⁾_𝑚𝑛𝑞,𝑚=1,…,𝑀₁,𝑛=1,…,𝑁₁𝜉,𝑞=1,…,𝑄,_(2)𝑚𝑛=𝜓⁽²⁾_𝑚𝑛(𝜃)+𝑤⁽²⁾_𝑚𝑛,𝑚=1,…,𝑀₂,𝑛=1,…,𝑁₂. (1) Here, 𝜉⁽¹⁾_𝑚𝑛𝑞 are complex valued data (dynamic electrical fields) obtained in the frequency domain and 𝜉⁽²⁾_𝑚𝑛 are real valued data (static electrical potentials), which are obtained from simultaneous EM and ER tomographic measurements, respectively. It is assumed that there are 𝑁_𝑖 excitations with 𝑛=1,…,𝑁_𝑖, and 𝑀_𝑖 measurements for each excitation with 𝑚=1,…𝑀_𝑖, where 𝑖=1,2 corresponds to the two EM and ER imaging modalities, respectively. Hence, for each imaging modality, the measurements are uniquely identified by the double index (𝑚𝑛). For the EM modality there are also 𝑄 frequency domain measurements with 𝑞=1,…,𝑄.

The unknown, real parameter function (distributed parameter) of interest is 𝜃(𝐫)=[𝜖_r(𝐫)𝜎(𝐫)]^T where 𝜖_r(𝐫) is the relative permittivity and 𝜎(𝐫) the conductivity, and where both parameter functions are defined over some spatial domain Ω with 𝐫∈Ω. Here, (⋅)^T denotes the transpose operation.

The physical forward models are denoted 𝜓⁽¹⁾(𝜃)={𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)} and 𝜓⁽²⁾(𝜃)={𝜓⁽²⁾_𝑚𝑛(𝜃)} in the two EM and ER tomography problems, respectively. The Fréchet derivative [12], or the Jacobian [13] of the forward operator 𝜓^(𝑖), is the bounded linear operator 𝒥^(𝑖) that yields the first order variation of 𝜓^(𝑖), denoted as𝛿𝜓^(𝑖)=𝒥^(𝑖)𝛿𝜃, (2) where 𝑖=1,2, and which are defined for each modality and measurement by𝛿𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)=_Ω𝛿_𝜃^T_(𝐫)𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)(𝐫)=𝛿𝜃d𝑣_Ω𝛿_𝜖r(𝐫)𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)𝛿𝜖_r(𝐫)d𝑣+_Ω𝛿_𝜎(𝐫)𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)𝛿𝜎(𝐫)d𝑣,𝛿𝜓⁽²⁾_𝑚𝑛(𝜃)=_Ω𝛿_𝜃^𝑇_(𝐫)𝜓⁽²⁾_𝑚𝑛(𝜃)(𝐫)=𝛿𝜃d𝑣_Ω𝛿_𝜖r(𝐫)𝜓⁽²⁾_𝑚𝑛(𝜃)𝛿𝜖_r(𝐫)d𝑣+_Ω𝛿_𝜎(𝐫)𝜓⁽²⁾_𝑚𝑛(𝜃)𝛿𝜎(𝐫)d𝑣, (3) where 𝛿𝜃(𝐫)=[𝛿𝜖_r(𝐫)𝛿𝜎(𝐫)]^T denotes a perturbation of 𝜃(𝐫) and 𝛿_𝜃(𝐫)𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃) and 𝛿_𝜃(𝐫)𝜓⁽²⁾_𝑚𝑛(𝜃) the corresponding gradients [12]. The following notation is employed for the Jacobian 𝒥 of the combined multiphysics problem𝛿𝜓=𝛿𝜓⁽¹⁾𝛿𝜓⁽²⁾=𝒥⁽¹⁾𝒥⁽²⁾𝒥𝛿𝜃=^𝜖(1)r𝒥_𝜎⁽¹⁾0𝒥_𝜎⁽²⁾𝛿𝜖_r𝛿𝜎=𝒥𝛿𝜃, (4) where 𝜓(𝜃)=[𝜓⁽¹⁾(𝜃)𝜓⁽²⁾(𝜃)]^T, and where it has been emphasized that the electric resistive (ER) tomography measurements depend only on the conductivity parameter 𝜎(𝐫).

The measurement noise in the statistical model (1) is denoted 𝑤⁽¹⁾_𝑚𝑛𝑞 and 𝑤⁽²⁾_𝑚𝑛, and is assumed to be the samples of a zero mean and uncorrelated multivariate Gaussian random variable [14] with variance𝑅⁽¹⁾_𝑚𝑛𝑞=ℰ𝑤⁽¹⁾_𝑚𝑛𝑞²,𝑅⁽²⁾_𝑚𝑛=ℰ𝑤⁽²⁾_𝑚𝑛², (5) corresponding to the EM and ER imaging modalities, respectively. Here, ℰ{⋅} denotes the statistical expectation operator. For the EM-modality, the random variable is complex Gaussian with ℰ{(𝑤⁽¹⁾_𝑚𝑛𝑞)²}=0, and for the ER-modality, the random variable is real Gaussian, see, for example, [14, 15].

2.2. Maximum Likelihood and the Truncated SVD

Let ln𝑝(𝜉∣𝜃) denote the loglikelihood function of the measurement statistics [14]. The Fisher information integral operator is then defined by the outer productℐ𝐫^′,𝐫^′′=ℰ𝛿_{𝜃(𝐫′)}ln𝑝(𝜉∣𝜃)𝛿_{𝜃(𝐫′′)}(,𝐫ln𝑝𝜉∣𝜃)_′,𝐫^′′∈Ω, (6) where ℰ{⋅} denotes the expectation and 𝛿_𝜃(𝐫) the gradient operator, cf. [10, 11, 14].

The multiphysics inverse problem (1) relates to a combination of a complex valued Gaussian and a real valued Gaussian measurement statistics. The negative loglikelihood function (or cost function) and the Fisher information (FI) for this situation is given by𝐹(𝜃)()==−ln𝑝𝜉∣𝜃_{𝑀1,𝑁1,𝑄𝑚,𝑛,𝑞=1,1,1}1𝑅⁽¹⁾_𝑚𝑛𝑞𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)−𝜉⁽¹⁾_𝑚𝑛𝑞²+^𝑀2,𝑁2_{𝑚,𝑛=1,1}12𝑅⁽²⁾_𝑚𝑛𝜓⁽²⁾_𝑚𝑛(𝜃)−𝜉⁽²⁾_𝑚𝑛²,ℐ𝐫^′,𝐫^′′=2Re^{𝑀1,𝑁1,𝑄}_{𝑚,𝑛,𝑞=1,1,1}1𝑅⁽¹⁾_𝑚𝑛𝑞𝛿_{𝜃(𝐫′)}𝜓^(1)∗_𝑚𝑛𝑞(𝜃)𝛿_𝜃^T_(𝐫′′)𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)+^𝑀2,𝑁2_{𝑚,𝑛=1,1}1𝑅⁽²⁾_𝑚𝑛𝛿_{𝜃(𝐫′)}𝜓⁽²⁾_𝑚𝑛(𝜃)𝛿_𝜃^T_(𝐫′′)𝜓⁽²⁾_𝑚𝑛(𝜃), (7)(8) where (⋅)^∗ denotes the complex conjugate, and 𝐫^′,𝐫^′′∈Ω, cf. [10, 11, 14].

To obtain a useful relationship between the maximum likelihood method and the singular value decomposition (SVD) of the forward operator, the adequate scalar products must be defined, cf. [10, 11]. For the real valued parameter functions 𝜃₁=[𝜖_𝑟1𝜎₁]^T and 𝜃₂=[𝜖_r2𝜎₂]^T, the following scalar product is employed𝜃₁,𝜃₂=_Ω𝜃^T₁(𝐫)𝜃₂(𝐫)=d𝑣_Ω𝜖_r1(𝐫)𝜖_r2(𝐫)d𝑣+_Ω𝜎₁(𝐫)𝜎₂(𝐫)d𝑣. (9) For the mixed complex and real valued measurement data 𝜉=[𝜉⁽¹⁾𝜉⁽²⁾]^T and 𝜂=[𝜂⁽¹⁾𝜂⁽²⁾]^T, the following real scalar product is defined𝜉,𝜂_d=Re^{𝑀1,𝑁1,𝑄}_{𝑚,𝑛,𝑞=1}1𝑅⁽¹⁾_𝑚𝑛𝑞𝜉^(1)∗_𝑚𝑛𝑞𝜂⁽¹⁾_𝑚𝑛𝑞+^𝑀2,𝑁2_𝑚,𝑛=112𝑅⁽²⁾_𝑚𝑛𝜉⁽²⁾_𝑚𝑛𝜂⁽²⁾_𝑚𝑛, (10) where 𝑅⁽¹⁾_𝑚𝑛𝑞 and 𝑅⁽²⁾_𝑚𝑛 are the variance of the noise as defined in (5), and which are also employed in (7) and (8). It can be readily verified that (10) is a scalar product. Note, for example, that the 𝑛-dimensional vector space ℂ^𝑛 over the complex scalar field and with scalar product 𝐱^H𝐲 is isomorphic with the 2𝑛-dimensional vector space consisting of complex vectors over the real scalar field and with real scalar product Re𝐱^H𝐲. Note that 𝐱 and i𝐱 are orthogonal in the latter space.

For later use, the adjoint forward operator 𝒥^∗ is now given by𝒥^∗𝜉=𝛿_𝜃(𝐫)𝜓(𝜃),𝜉_d=Re^{𝑀1,𝑁1,𝑄}_{𝑚,𝑛,𝑞=1,1,1}1𝑅⁽¹⁾_𝑚𝑛𝑞𝛿_𝜃(𝐫)𝜓^(1)∗_𝑚𝑛𝑞(𝜃)𝜉⁽¹⁾_𝑚𝑛𝑞+^𝑀2,𝑁2_{𝑚,𝑛=1,1}12𝑅⁽²⁾_𝑚𝑛𝛿_𝜃(𝐫)𝜓⁽²⁾_𝑚𝑛(𝜃)𝜉⁽²⁾_𝑚𝑛=𝒥^(1)∗𝜉⁽¹⁾+𝒥^(2)∗𝜉⁽²⁾, (11) where ⟨𝒥𝛿𝜃,𝜉⟩_d=⟨𝛿𝜃,𝒥^∗𝜉⟩, see also [10–12, 16].

Based on the scalar product (10), the negative loglikelihood function (7) can now be written𝐹(𝜃)=𝜓(𝜃)−𝜉,𝜓(𝜃)−𝜉_d. (12) By using the adjoint operator defined in (11), the first and second order variations of 𝐹(𝜃) are obtained as𝛿𝐹=𝛿𝜃,2𝒥^∗,𝛿𝜓−𝜉₂𝐹=⟨𝛿𝜃,2𝒥^∗𝒥𝛿𝜃⟩+2⟨𝛿²𝜓,𝜓−𝜉⟩_d, (13) where 𝛿𝜃 is the perturbation, 𝜓=𝜓(𝜃) the forward operator evaluated at some known background parameter value 𝜃 and 𝛿²𝜓 the second order variation of the forward operator. It is observed that the Fisher information operator defined in (8) is given by ℐ=2𝒥^∗𝒥=2𝒥^(1)∗𝒥⁽¹⁾+2𝒥^(2)∗𝒥⁽²⁾=ℐ⁽¹⁾+ℐ⁽²⁾, (14) where ℐ⁽¹⁾=2𝒥^(1)∗𝒥⁽¹⁾ and ℐ⁽²⁾=2𝒥^(2)∗𝒥⁽²⁾ are the Fisher information operators for each modality. The gradient 𝛿_𝜃(𝐫)𝐹(𝜃) of the cost function 𝐹(𝜃) is furthermore given by𝛿_𝜃(𝐫)𝐹(𝜃)=2𝒥^∗𝜓−𝜉=2𝒥^(1)∗𝜓⁽¹⁾−𝜉⁽¹⁾+2𝒥^(2)∗𝜓⁽²⁾−𝜉⁽²⁾, (15) where2𝒥^(1)∗𝜓⁽¹⁾−𝜉⁽¹⁾=2Re^{𝑀1,𝑁1,𝑄}_{𝑚,𝑛,𝑞=1,1,1}1𝑅⁽¹⁾_𝑚𝑛𝑞𝛿_𝜃(𝐫)𝜓^(1)∗_𝑚𝑛𝑞(𝜃)𝜓⁽¹⁾_𝑚𝑛𝑞(𝜃)−𝜉⁽¹⁾_𝑚𝑛𝑞,2𝒥^(2)∗𝜓⁽²⁾−𝜉⁽²⁾=^𝑀2,𝑁2_{𝑚,𝑛=1,1}1𝑅⁽²⁾_𝑚𝑛𝛿_𝜃(𝐫)𝜓⁽²⁾_𝑚𝑛(𝜃)𝜓⁽²⁾_𝑚𝑛(𝜃)−𝜉⁽²⁾_𝑚𝑛. (16)

A linearization is now considered where it is assumed that the second order variation of the forward operator 𝛿²𝜓 can be neglected. Consider hence the quadratic functional()()+1Δ𝐹𝛿𝜃=𝛿𝐹𝛿𝜃2𝛿²𝐹()=𝛿𝜃𝛿𝜃,𝛿_𝜃(𝐫)𝐹(𝜃+1)2⟨𝛿𝜃,ℐ𝛿𝜃⟩, (17) where the linear and quadratic terms in (13) have been used with 𝛿²𝜓=0. Note that the Hessian of 𝐹 is given here by the Fisher information ℐ=2𝒥^∗𝒥. A regularized one-step Newton algorithm for this functional is hence given by𝛿𝜃⁺=−ℐ⁺𝛿_𝜃(𝐫)𝐹(𝜃), (18) where ℐ⁺ denotes a truncated pseudoinverse of the Fisher information ℐ, which is based on the singular value decomposition (SVD), see also [10, 11, 17].

It can be readily shown that the regularized maximum likelihood solution (18) is equivalent to a regularized pseudoinverse based directly on the linearized operator problem𝒥𝛿𝜃=𝜉−𝜓, (19) provided that the scalar products are properly chosen as in (9) and (10). To this end, it is assumed that the operator 𝒥 is compact and hence that the following singular system exists√2𝒥𝑣_𝑙=𝜎_𝑙𝑢_𝑙,√2𝒥^∗𝑢_𝑙=𝜎_𝑙𝑣_𝑙, (20) where {𝑢_𝑙} and {𝑣_𝑙} are orthonormal systems and 𝜎_𝑙 the singular values. (Note that 𝜎 is used here to denote singular values as well as conductivity. The correct interpretation should be evident from the context without confusion [12].) The Fisher information is the self-adjoint operator ℐ=2𝒥^∗𝒥 with eigenvalues 𝜇_𝑙=𝜎²_𝑙 and eigenfunctions 𝑣_𝑙, that is,ℐ𝑣_𝑙=𝜇_𝑙𝑣_𝑙. (21) The operator 𝒥 can now be represented by the singular value decomposition (SVD)1𝒥𝛿𝜃=√2_𝑙𝜎_𝑙𝑢_𝑙𝑣_𝑙,𝛿𝜃, (22) and the pseudoinverse (if it exists) is given by𝒥⁺𝜉−𝜓=√2_𝑙1𝜎_𝑙𝑣_𝑙𝑢_𝑙,𝜉−𝜓_d=√2_𝑙1𝜎_𝑙𝑣_𝑙√2𝜎_𝑙𝒥𝑣_𝑙,𝜉−𝜓_d=_𝑙1𝜎²_𝑙𝑣_𝑙𝑣_𝑙,2𝒥^∗𝜉−𝜓=−_𝑙1𝜎²_𝑙𝑣_𝑙𝑣_𝑙,𝛿_𝜃(𝐫)𝐹(𝜃)=−ℐ⁺𝛿_𝜃(𝐫)𝐹(𝜃), (23) which is the same as (18).

The relationship between ML, SVD, and FI has been described. This will be used in Section 4 to calculate the pseudoinverse and find the parameters.

3. Green's Functions and Gradients

Below is given a description of green’s functions and the gradients used with the EM and ER tomographic problems in this paper.

3.1. The Electromagnetic (EM) Tomography Problem

A two-dimensional electromagnetic tomography problem is considered where the harmonic electric field is polarized as 𝐄=𝜓(𝜌,𝜙)𝐳̂ where the scalar field 𝜓(𝜌,𝜙) is independent of the 𝑧-coordinate. Here, 𝜌 and 𝜙 denote the radial and angular circular coordinates, respectively, and 𝝆=𝜌𝝆 the two-dimensional radius vector. Let 𝑘=𝜔√𝜇₀𝜖₀, 𝜂₀=𝜇₀/𝜖₀, 𝜇₀ and 𝜖₀ denote the wave number, the wave impedance, the permeability, and the permittivity of free space, respectively, and where 𝜔 is the angular frequency. The time convention employed here is given by e^i𝜔𝑡.

The EM-field 𝜓_𝑛⁽¹⁾(𝜌,𝜙) resulting from excitation index 𝑛 is obtained as the solution to the wave equation∇²+𝑘²𝜖𝜓_𝑛⁽¹⁾=i𝑘𝜂₀𝛿𝝆−𝝆_𝑛,𝑛=1,…,𝑁₁, (24) where 𝜖=𝜖_r−i𝜎𝜂₀/𝑘 is the complex valued relative permittivity, 𝜖_r the relative permittivity, and 𝜎 the conductivity. Here, 𝛿(⋅) denotes the Dirac-delta function indicating that the excitation is modeled here as a point source at position 𝝆_𝑛.

Let 𝜖=𝜖_D+𝛿𝜖 where 𝜖_D=𝜖_rD−i𝜎_D𝜂₀/𝑘 and 𝛿𝜖=𝛿𝜖_r−i𝛿𝜎𝜂₀/𝑘 denote the corresponding background parameter values and a perturbation 𝛿𝜃=[𝛿𝜖_r,𝛿𝜎)]^T, respectively. For the homogeneous background 𝜖_D, the solution to (24) is given by𝜓_𝑛⁽¹⁾(𝝆)=−i𝑘𝜂₀𝐺⁽¹⁾_D𝝆,𝝆_𝑛, (25) where 𝐺⁽¹⁾_D(𝝆,𝝆^′) is green’s function for the homogeneous background satisfying (∇²+𝑘²𝜖_D)𝐺⁽¹⁾_D(𝝆,𝝆^′)=−𝛿(𝝆−𝝆^′), and which is given by𝐺_𝐷⁽¹⁾𝝆,𝝆^′=14i_𝑚J_𝑚𝑘√𝜖_D𝜌_<H_𝑚⁽²⁾𝑘√𝜖_D𝜌_>e^{i𝑚(𝜙−𝜙′)}, (26) where J_𝑚(⋅) and H_𝑚⁽²⁾(⋅) denote the Bessel function and the Hankel function of the second kind, respectively, both of order 𝑚 [18]. Here, 𝜌_<=min{𝜌,𝜌^′} and 𝜌_>=max{𝜌,𝜌^′}.

A first order perturbation analysis of (24) yields the following wave equation for the perturbed field 𝛿𝜓_𝑛⁽¹⁾∇²+𝑘²𝜖𝛿𝜓_𝑛⁽¹⁾=−𝑘²𝛿𝜖𝜓_𝑛⁽¹⁾, (27) where 𝜓_𝑛⁽¹⁾ is the solution to (24) and 𝛿𝜖 the perturbation in parameter values. For the homogeneous background with 𝜖=𝜖_D, the solution to (27) is hence given by𝛿𝜓_𝑛⁽¹⁾(𝝆)=𝑘²_Ω𝐺⁽¹⁾_D𝝆,𝝆^′𝜓_𝑛⁽¹⁾𝝆^′𝛿𝜖𝝆^′d𝑣^′=−i𝑘³𝜂_0Ω𝐺⁽¹⁾_D𝝆,𝝆^′𝐺⁽¹⁾_D𝝆^′,𝝆_𝑛𝛿𝜖𝝆^′d𝑣^′. (28) It should be noted that this is the same result, that is, obtained by employing a Born approximation [4], in which case 𝜓_𝑛⁽¹⁾(𝝆) represents the incident field and 𝛿𝜓_𝑛⁽¹⁾(𝝆) the scattered field.

Let 𝜓⁽¹⁾_𝑚𝑛 denote the measured quantity 𝜓⁽¹⁾_𝑚𝑛=𝜓_𝑛⁽¹⁾(𝝆_𝑚𝑛) where 𝝆_𝑚𝑛 is the measurement position at measurement 𝑚 with excitation 𝑛. The corresponding gradients are then given by𝛿_𝜖r(𝝆)𝜓⁽¹⁾_𝑚𝑛=−i𝑘³𝜂₀𝐺⁽¹⁾_D𝝆_𝑚𝑛,𝝆𝐺⁽¹⁾_D𝝆_𝑛,𝛿,𝝆_𝜎(𝝆)𝜓⁽¹⁾_𝑚𝑛=−𝑘²𝜂²₀𝐺⁽¹⁾_D𝝆_𝑚𝑛,𝝆𝐺⁽¹⁾_D𝝆_𝑛,𝝆, (29) where the symmetry of green’s function has been employed.

In the simulation examples given below, a homogenous cylinder is used to model a tunnel buried in the ground. Green’s function for the cylinder is denoted 𝐺⁽¹⁾_c(𝝆,𝝆^′), and the measured quantity is hence modeled by𝜓⁽¹⁾_𝑚𝑛=𝜓_𝑛⁽¹⁾𝝆_𝑚𝑛=−i𝑘𝜂₀𝐺⁽¹⁾_c𝝆_𝑚𝑛,𝝆_𝑛. (30) Green’s function 𝐺⁽¹⁾_c(𝝆,𝝆^′) for a homogeneous cylinder of radius 𝑎 placed in a homogeneous background space is given in Appendix A, and is used to generate the synthetic data in the numerical example.

3.2. The Electrical Resistive (ER) Tomography Problem

A two-dimensional electrical resistive tomography problem is considered where the electric potential field Φ⁽²⁾(𝜌,𝜙) is independent of the 𝑧-coordinate, and where 𝜌 and 𝜙 are used to denote the radial and angular circular coordinates, as above. The potential field resulting from excitation index 𝑛 is obtained here as the solution to the following differential equation∇⋅𝜎∇Φ_𝑛⁽²⁾=𝛿𝝆−𝝆_𝑛−𝛿𝝆−𝝆_𝑁2+1,𝑛=1,…,𝑁₂, (31) where 𝜎 is the conductivity and the excitation is defined by current point sources at positions 𝝆_𝑛 with a common ground connection at 𝝆_𝑁2+1, see also [17].

Let 𝜎=𝜎_D+𝛿𝜎 where 𝜎_D and 𝛿𝜎 denote the corresponding background parameter value and a perturbation, respectively. For the homogeneous background 𝜎_D, the solution to (31) is given byΦ_𝑛⁽²⁾(𝝆)=𝐺⁽²⁾_D𝝆,𝝆_𝑛−𝐺⁽²⁾_D𝝆,𝝆_𝑁2+1, (32) where 𝐺⁽²⁾_D(𝝆,𝝆^′) is green’s function for the homogeneous background satisfying ∇⋅(𝜎_D∇𝐺⁽²⁾_D(𝝆,𝝆^′))=𝛿(𝝆−𝝆^′), and which is given by𝐺⁽²⁾_D𝝆,𝝆^′=12𝜋𝜎_Dln𝜌_>−14𝜋𝜎_D𝑚≠01|𝑚|𝜌_<𝜌_>^|𝑚|e^{i𝑚(𝜙−𝜙′)}, (33) where 𝜌_<=min{𝜌,𝜌^′} and 𝜌_>=max{𝜌,𝜌^′}, see also [18].

The observed voltage quantity 𝜓⁽²⁾_𝑚𝑛 corresponding to the homogeneous background is given by𝜓⁽²⁾_𝑚𝑛=Φ_𝑛⁽²⁾𝝆_𝑚𝑛−Φ_𝑛⁽²⁾𝝆₀=𝐺⁽²⁾_D𝝆_𝑚𝑛,𝝆_𝑛−𝐺⁽²⁾_D𝝆_𝑚𝑛,𝝆_𝑁2+1−𝐺⁽²⁾_D𝝆₀,𝝆_𝑛+𝐺⁽²⁾_D𝝆₀,𝝆_𝑁2+1, (34) where 𝝆_𝑚𝑛 is the measurement position at measurement 𝑚 and excitation 𝑛, and 𝝆₀ is a position used as a voltage reference.

Note that green’s function 𝐺⁽²⁾_D(𝝆,𝝆_𝑛) satisfies the differential equation∇⋅𝜎_D∇𝐺⁽²⁾_D𝝆,𝝆_𝑛=𝛿𝝆−𝝆_𝑛. (35) A first order perturbation analysis of (35) yields the following differential equation for the perturbed field 𝛿𝐺⁽²⁾_D(𝝆,𝝆_𝑛)∇⋅𝜎_D∇𝛿𝐺⁽²⁾_D𝝆,𝝆_𝑛=−∇⋅𝛿𝜎∇𝐺⁽²⁾_D𝝆,𝝆_𝑛. (36) The solution to (36) is given by𝛿𝐺⁽²⁾_D𝝆,𝝆_𝑛=−_Ω𝐺⁽²⁾_D𝝆,𝝆^′∇^′⋅𝛿𝜎∇^′𝐺⁽²⁾_D𝝆^′,𝝆_𝑛d𝑣^′=_Ω∇^′𝐺⁽²⁾_D𝝆,𝝆^′⋅∇^′𝐺⁽²⁾_D𝝆_𝑛,𝝆^′𝛿𝜎𝝆^′d𝑣^′, (37) where Green's theorem [18] has been used together with the assumption that 𝛿𝜎=0 on the boundary of Ω. The gradient of green’s function 𝐺⁽²⁾_D(𝝆_𝑙,𝝆_𝑛) is hence given by𝛿_𝜎(𝝆)𝐺⁽²⁾_D𝝆_𝑙,𝝆_𝑛=∇𝐺⁽²⁾_D𝝆_𝑙,𝝆⋅∇𝐺⁽²⁾_D𝝆_𝑛,𝝆, (38) where ∇𝐺⁽²⁾_D(𝝆_𝑛,𝝆) can be obtained from (33) as∇𝐺⁽²⁾_D𝝆_𝑛,𝝆=𝝆−14𝜋𝜎_D𝜌_𝑛𝑚≠0𝜌𝜌_𝑛^|𝑚|−1e^{i𝑚(𝜙𝑛−𝜙)}+𝜙i4𝜋𝜎_D𝜌_𝑚≠0𝑚|𝑚|𝜌𝜌_𝑛^|𝑚|e^{i𝑚(𝜙𝑛−𝜙)}, (39) where (𝜌_𝑛,𝜙_𝑛) are the circular coordinates of 𝝆_𝑛, and 𝜌<𝜌_𝑛. The gradient of the observed voltage quantity 𝜓⁽²⁾_𝑚𝑛 defined in (34) can finally be obtained as𝛿_𝜎(𝝆)𝜓⁽²⁾_𝑚𝑛=𝛿_𝜎(𝝆)𝐺⁽²⁾_D𝝆_𝑚𝑛,𝝆_𝑛−𝛿_𝜎(𝝆)𝐺⁽²⁾_D𝝆_𝑚𝑛,𝝆_𝑁2+1−𝛿_𝜎(𝝆)𝐺⁽²⁾_D𝝆₀,𝝆_𝑛+𝛿_𝜎(𝝆)𝐺⁽²⁾_D𝝆₀,𝝆_𝑁2+1, (40) which can be computed with the aid of (38) and (39).

In the simulation examples given below, a homogenous cylinder is used to model a tunnel buried in the ground. Green’s function for the cylinder is denoted 𝐺⁽²⁾_c(𝝆,𝝆^′), and the measured quantity is hence modeled by𝜓⁽²⁾_𝑚𝑛=𝐺⁽²⁾_c𝝆_𝑚𝑛,𝝆_𝑛−𝐺⁽²⁾_c𝝆_𝑚𝑛,𝝆_𝑁2+1−𝐺⁽²⁾_c𝝆₀,𝝆_𝑛+𝐺⁽²⁾_c𝝆₀,𝝆_𝑁2+1, (41) similar to (34). Green’s function 𝐺⁽²⁾_c(𝝆,𝝆^′) for a homogeneous cylinder of radius 𝑎 placed in a homogeneous background space is given in Appendix B, and is used to generate the synthetic data in the numerical example.

4. Numerical Example

As a numerical example, we consider the problem of tunnel detection based on electromagnetic (EM) and electrical resistive (ER) tomography in geophysical sensing. A two-dimensional problem is studied where the tunnel is represented by a circular structure embedded in a homogeneous background. The measurement setup is illustrated in Figure 1, which also shows green’s functions 𝐺⁽¹⁾_c(𝝆,𝝆^′) and 𝐺⁽²⁾_c(𝝆,𝝆^′) of the circular object which have been used to generate the synthetic data. The superscripts (1) and (2) refer to the EM and ER imaging modalities, respectively. One single frequency with 𝜔=2𝜋⋅200⋅10⁶ (200 MHz) is employed with the EM modality. The large circle in the center indicates the simulated tunnel of radius 𝑎=0.5 m, and the 21 small circles the horizontal sensors which are modeled here as point sources. The measurement line is 2 m long. The surrounding medium has relative permittivity 𝜖_r𝐷=5 and conductivity 𝜎_D=1 mS/m (𝜎_D=0.001 S/m). The tunnel is modeled with relative permittivity 𝜖_r𝑎=1 and conductivity 𝜎_a=0. This formulation has been chosen to illustrate the data fusion, as it represents a simple and generic multiphysics inverse problem which is severely ill-posed due to the narrow measurement aperture (narrow spatial window) and the narrow measurement bandwidth.

(a)

(b)

Figure 1 

Illustration of measurement setup and Green's functions 𝐺(1)c(𝝆,𝝆′) and 𝐺(2)c(𝝆,𝝆′) in electromagnetic and electrical resistive tomography, respectively. The large circle in the center indicates the tunnel of radius 𝑎=0.5 m, and the 21 small circles the horizontal sensors (point sources). In both plots shown above, the source is at 𝝆′ corresponding to the fourth sensor from the left.

The signal-to-noise-ratio (SNR) in the measurement model (1) is defined asSNR⁽¹⁾=max_𝑚𝑛𝑞𝜓⁽¹⁾_𝑚𝑛𝑞²𝑅⁽¹⁾,SNR⁽²⁾=max_𝑚𝑛𝜓⁽²⁾_𝑚𝑛²𝑅⁽²⁾, (42) where 𝜓⁽¹⁾_𝑚𝑛𝑞 and 𝜓⁽²⁾_𝑚𝑛 represent the (noiseless) measured quantities, and 𝑅⁽¹⁾ and 𝑅⁽²⁾ the corresponding measurement noise variances defined in (5) corresponding to the EM and ER imaging modalities, respectively. Here, the measurement noise variance is assumed to be independent of the measurement index (𝑚,𝑛,𝑞), and hence 𝑅⁽¹⁾_𝑚𝑛𝑞=𝑅⁽¹⁾ and 𝑅⁽²⁾_𝑚𝑛=𝑅⁽²⁾. In the simulation examples described below, the measurement noise was chosen such that SNR⁽¹⁾=80 dB, SNR⁽²⁾=70 dB and 𝑅⁽²⁾/𝑅⁽¹⁾=92. It should be noted that the values of the signal-to-noise-ratios are rather high for a practical implementation, but are required here due to the absence of directivity of the point sources used in this model. However, the main point here is to illustrate the significance of the ratio 𝑅⁽²⁾/𝑅⁽¹⁾ which is independent of the absolute value of the SNR.

In Figures 2 and 3 are shown reconstruction results for 𝛿𝜖_r and 𝛿𝜎 using the truncated SVD inversion technique for data fusion that has been outlined in Section 2.2. Here, the upper and lower rows show reconstruction results obtained by using 𝑁_eig=15 and 𝑁_eig=30 eigenvalues, respectively.

(a)

(b)

(c)

(d)

Figure 2 

Reconstruction of 𝛿𝜖r and 𝛿𝜎 in a straightforward data fusion without statistically based weighting (the weighting is assumed to be uniform with 𝑅(2)/𝑅(1)=1). The upper and lower rows are obtained by using 𝑁eig=15 and 𝑁eig=30 eigenvalues, respectively.

(a)

(b)

(c)

(d)

Figure 3 

Reconstruction of 𝛿𝜖r and 𝛿𝜎 in a maximum likelihood-based data fusion where the statistically-based weighting is 𝑅(2)/𝑅(1)=92. The upper and lower rows are obtained by using 𝑁eig=15 and 𝑁eig=30 eigenvalues, respectively.

In Figure 2 is shown the reconstructions in a straightforward data fusion without statistically-based weighting. In these reconstructions, the weighting is assumed to be uniform with 𝑅⁽²⁾/𝑅⁽¹⁾=1 even though the simulation data has been generated with 𝑅⁽²⁾/𝑅⁽¹⁾=92. Hence, the correct multiphysics model has been employed, but the sensor noise statistics has not been taken into consideration.

In Figure 3 is shown the reconstructions in a maximum likelihood-based data fusion where the statistically-based weighting with 𝑅⁽²⁾/𝑅⁽¹⁾=92 has been exploited. Hence, the correct multiphysics model has been employed, and the sensor noise statistics has been taken into proper consideration based on the principle of maximum likelihood.

A comparison of Figures 2 and 3 in this example illustrates the expected behavior of the data fusion for an ill-posed multiphysics inverse problem, that is, the sensitivity of the inversion with respect to the discrepancies in the statistical assumptions. In this example, max_𝑚𝑛|𝜓⁽²⁾_𝑚𝑛|²=9.2, max_𝑚𝑛|𝜓⁽¹⁾_𝑚𝑛|² and SNR⁽²⁾=SNR⁽¹⁾/10, yielding an imbalance in the weighting of data corresponding to the ratio 𝑅⁽²⁾/𝑅⁽¹⁾=92. As can be seen in the Figures 2 and 3, the consequence of neglecting this imbalance can be significant, especially for 𝑁_eig=15.

In Figure 4 is shown the singular values 𝜇_𝑛(ℐ) for the three different measurement scenarios with electromagnetic (EM), electrical resistive (ER), and combined electromagnetic and electrical resistive (EM+ER) tomography. The figure illustrates that ER tomography contributes mainly at the lower eigenvalue indices whereas EM tomography contributes at higher indices in this particular problem. It can be anticipated that this behavior is typical since the eigenvalue distribution of ER tomography usually drops of rapidly already from the first indices whereas with EM tomography there is typically a threshold in the distribution of eigenvalues after which the eigenvalues drop off in rapid descent, see, for example, [10, 11, 13, 19, 20].

Figure 4 

Singular values 𝜇𝑛(ℐ) for the three different measurement scenarios. (a) Electromagnetic tomography with SNR(1)=80 dB. (b) Electrical resistive tomography with SNR(2)=70 dB. (c) Maximum likelihood-based data fusion with electromagnetic and electrical resistive tomography combined.

5. Summary

A maximum likelihood-based approach to data fusion for electromagnetic (EM) and electrical resistive (ER) tomography has been presented. The Fisher information is an additive measure of information which is closely linked to the principle of maximum likelihood, and which can be useful in the study of inverse problems. The Fisher information is particularly useful for inverse problems which can be linearized similar to the Born approximation. In this paper, a proper scalar product has been defined for the measurements and a truncated singular value decomposition- (SVD-) based algorithm has been devised which combines the measurement data of the two imaging modalities in a way that is optimal in the sense of maximum likelihood.

As a multiphysics problem formulation with applications in geophysics, the problem of tunnel detection based on EM and ER tomography has been studied in this paper. The connection between green’s functions, the gradients, and the Fisher information has been illustrated by a detailed description of two simple and generic forward models regarding two-dimensional EM and ER tomography, respectively. Numerical examples have been included to illustrate the potential impact of an imbalance between the singular values and the variance of the measurement noise when different imaging modalities are incorporated in the inversion, and the significance of taking a statistically-based weighting of the measurement data into proper account.

Appendices

A. EM Green's Function for a Homogeneous Cylinder

Green’s function 𝐺⁽¹⁾_c(𝝆,𝝆^′) for a homogenous cylinder of radius 𝑎 is derived below. The complex valued permittivity of the background and the cylinder are given by 𝜖_D=𝜖_rD−i𝜎_D𝜂₀/𝑘 and 𝜖_c=𝜖_rc−i𝜎_c𝜂₀/𝑘, respectively.

Green’s function 𝐺⁽¹⁾_c(𝝆,𝝆^′) satisfies the scalar wave equation1𝜌𝜕𝜌𝜕𝜕𝜌+1𝜕𝜌𝜌²𝜕²𝜕𝜙²+𝑘²𝜖𝐺⁽¹⁾_c𝜌,𝜙,𝜌^′,𝜙^′𝛿=−𝜌−𝜌^′𝜌𝛿𝜙−𝜙^′, (A.1) together with appropriate boundary conditions. By introducing the Fourier series expansion𝐺⁽¹⁾_c𝜌,𝜙,𝜌^′,𝜙^′=^∞_𝑚=−∞𝐺_𝑚𝜌,𝜌^′,𝜙^′e^i𝑚𝜙, (A.2) the wave equation (A.1) transforms to𝜕²𝜕𝑝²+1𝜌𝜕𝜕𝜌+𝑘²𝑚𝜖−²𝜌²𝐺_𝑚𝜌,𝜌^′,𝜙^′𝛿=−𝜌−𝜌^′𝜌e^{−i𝑚𝜙′}.2𝜋 (A.3)

Assume that the source is outside the cylinder, that is, 𝜌^′>𝑎. The solution to (A.3) is then given by𝐺_𝑚𝜌,𝜌^′,𝜙^′=𝐴_𝑚J_𝑚𝑘√𝜖_c𝜌𝐵0<𝜌<𝑎,_𝑚H_𝑚⁽¹⁾𝑘√𝜖_D𝜌+𝐶_𝑚H_𝑚⁽²⁾𝑘√𝜖_D𝜌𝑎<𝜌<𝜌^′,𝐷_𝑚H_𝑚⁽²⁾𝑘√𝜖_D𝜌𝜌>𝜌^′. (A.4) The appropriate boundary conditions are given by−𝐺_𝑚𝑎₋,𝜌^′,𝜙^′+𝐺_𝑚𝑎₊,𝜌^′,𝜙^′−𝜕=0,𝐺𝜕𝜌^𝑚𝑎₋,𝜌^′,𝜙^′+𝜕𝐺𝜕𝜌^𝑚𝑎₊,𝜌^′,𝜙^′−𝐺=0,_𝑚𝜌^′₋,𝜌^′,𝜙^′+𝐺_𝑚𝜌^′₊,𝜌^′,𝜙^′−𝜕=0,𝐺𝜕𝜌^𝑚𝜌^′₋,𝜌^′,𝜙^′+𝜕𝐺𝜕𝜌^𝑚𝜌^′₊,𝜌^′,𝜙^′1=−𝜌^′e^{−i𝑚𝜙′}.2𝜋 (A.5) The boundary conditions (A.5) yield the linear system of equations−J_𝑚𝑘√𝜖_c𝑎H_𝑚⁽¹⁾𝑘√𝜖_D𝑎H_𝑚⁽²⁾𝑘√𝜖_D𝑎0−𝑘√𝜖_cJ^′_𝑚𝑘√𝜖_c𝑎𝑘√𝜖_DH^(1)′_𝑚𝑘√𝜖_D𝑎𝑘√𝜖_DH^(2)′_𝑚𝑘√𝜖_D𝑎00−H_𝑚⁽¹⁾𝑘√𝜖_D𝜌^′−H_𝑚⁽²⁾𝑘√𝜖_D𝜌^′H_𝑚⁽²⁾𝑘√𝜖_D𝜌^′0−𝑘√𝜖_DH^(1)′_𝑚𝑘√𝜖_D𝜌^′−𝑘√𝜖_DH^(2)′_𝑚𝑘√𝜖_D𝜌^′𝑘√𝜖_DH^(2)′_𝑚𝑘√𝜖_D𝜌^′×𝐴_𝑚𝐵_𝑚𝐶_𝑚𝐷_𝑚=000−1𝜌^′e^{−i𝑚𝜙′},2𝜋 (A.6)which is solved for the constants 𝐴_𝑚, 𝐵_𝑚, 𝐶_𝑚, and 𝐷_𝑚 defined in (A.4).

B. ER Green's Function for a Homogeneous Cylinder

Green’s function 𝐺⁽²⁾_c(𝝆,𝝆^′) for a homogenous cylinder of radius 𝑎 is derived below. The conductivity of the background and the cylinder are denoted by 𝜎_D and 𝜎_c, respectively.

Green’s function 𝐺⁽²⁾_c(𝝆,𝝆^′) satisfies the partial differential equation∇⋅𝜎∇𝐺⁽²⁾_c𝝆,𝝆^′=𝜎∇²𝐺⁽²⁾_c𝝆,𝝆^′+∇𝜎⋅∇𝐺⁽²⁾_c𝝆,𝝆^′=𝛿𝝆−𝝆^′, (B.1) or𝜎1𝜌𝜕𝜌𝜕𝜕𝜌+1𝜕𝜌𝜌²𝜕²𝜕𝜙²+𝜕𝜎𝜕𝜕𝜌𝐺𝜕𝜌^(2)c𝜌,𝜙,𝜌^′,𝜙^′=𝛿𝜌−𝜌^′𝜌𝛿𝜙−𝜙^′, (B.2) where 𝜕𝜎/𝜕𝜌=(𝜎_D−𝜎_c)𝛿(𝜌−𝑎). By introducing the Fourier series expansion𝐺⁽²⁾_c𝜌,𝜙,𝜌^′,𝜙^′=^∞_𝑚=−∞𝐺_𝑚𝜌,𝜌^′,𝜙^′e^i𝑚𝜙, (B.3) the differential equation (B.2) transforms to𝜎𝜕²𝜕𝜌²+1𝜌𝜕−𝑚𝜕𝜌𝜌²²+𝜎_D−𝜎_c𝛿𝜕𝜌−𝑎×𝐺𝜕𝜌_𝑚𝜌,𝜌^′,𝜙^′=𝛿𝜌−𝜌^′𝜌e^{−i𝑚𝜙′}.2𝜋 (B.4)

Assume that the source is outside the cylinder, that is, 𝜌^′>𝑎. The solution to (B.4) for 𝑚≠0 is then given by𝐺_𝑚𝜌,𝜌^′,𝜙^′=𝐴_𝑚𝜌^|𝑚|𝐵,0<𝜌<𝑎,_𝑚𝜌^|𝑚|+𝐶_𝑚𝜌^−|𝑚|,𝑎<𝜌<𝜌^′,𝐷_𝑚𝜌^−|𝑚|,𝜌>𝜌^′, (B.5) and for 𝑚=0𝐺₀𝜌,𝜌^′,𝜙^′=𝐴₀𝐵0<𝜌<𝑎,₀+𝐶₀ln𝜌𝑎<𝜌<𝜌^′,𝐷₀ln𝜌𝜌>𝜌^′. (B.6)

The appropriate boundary conditions are obtained from the continuity of 𝐺_𝑚(𝜌,𝜌^′,𝜙^′) and the discontinuity of (𝜕/𝜕𝜌)𝐺_𝑚(𝜌,𝜌^′,𝜙^′) at 𝜌=𝑎 and 𝜌=𝜌^′. Hence,−𝐺_𝑚𝑎₋,𝜌^′,𝜙^′+𝐺_𝑚𝑎₊,𝜌^′,𝜙^′−𝜕=0,𝐺𝜕𝜌^𝑚𝑎₋,𝜌^′,𝜙^′𝜎_c+𝜕𝐺𝜕𝜌^𝑚𝑎₊,𝜌^′,𝜙^′𝜎_D−𝐺=0,_𝑚𝜌^′₋,𝜌^′,𝜙^′+𝐺_𝑚𝜌^′₊,𝜌^′,𝜙^′−𝜕=0,𝐺𝜕𝜌^𝑚𝜌^′₋,𝜌^′,𝜙^′+𝜕𝐺𝜕𝜌^𝑚𝜌^′₊,𝜌^′,𝜙^′=1𝜎_D𝜌^′e^{−i𝑚𝜙′}.2𝜋 (B.7) The boundary conditions (B.7) yield the following linear system of equations for 𝑚≠0−𝑎^|𝑚|𝑎^|𝑚|𝑎^−|𝑚|0−|𝑚|𝑎^|𝑚|−1𝜎_c|𝑚|𝑎^|𝑚|−1𝜎_D−|𝑚|𝑎^{−|𝑚|−1}𝜎_D00−𝜌^′|𝑚|−𝜌^′−|𝑚|𝜌^′−|𝑚|0−|𝑚|𝜌^{′|𝑚|−1}|𝑚|𝜌^{′−|𝑚|−1}−|𝑚|𝜌^{′−|𝑚|−1}×𝐴_𝑚𝐵_𝑚𝐶_𝑚𝐷_𝑚=000−1𝜎_D𝜌^′e^{−i𝑚𝜙′},2𝜋 (B.8)which is solved for the constants 𝐴_𝑚, 𝐵_𝑚, 𝐶_𝑚, and 𝐷_𝑚 defined in (B.5). For 𝑚=0, the boundary conditions (B.7) yield the following linear system of equations1−11ln𝑎000𝑎𝜎_D00−1−ln𝜌^′ln𝜌^′100−𝜌^′1𝜌^′×𝐴₀𝐵₀𝐶₀𝐷₀=0001𝜎_D𝜌^′12𝜋, (B.9) which is solved for the constants 𝐴₀, 𝐵₀, 𝐶₀, and 𝐷₀ defined in (B.6). The solution is readily obtained as 𝐴₀=𝐵₀=(1/𝜎_D)(1/2𝜋)ln𝜌^′, 𝐶₀=0 and 𝐷₀=(1/𝜎_D)(1/2𝜋).

Acknowledgment

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007–2013) under Grant Agreement no. 225663.