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

Authors: Sven Nordebo; Mats Gustafsson; Therese Sjöden; Francesco Soldovieri
Journal: International Journal of Geophysics
 (2011)

---

## 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.


---

## Body

## 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. Section2.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)ξmnq(1)=ψmnq(1)(θ)+wmnq(1),m=1,…,M1,n=1,…,N1,q=1,…,Q,ξmn(2)=ψmn(2)(θ)+wmn(2),m=1,…,M2,n=1,…,N2.
Here, ξmnq(1) are complex valued data (dynamic electrical fields) obtained in the frequency domain and ξmn(2) are real valued data (static electrical potentials), which are obtained from simultaneous EM and ER tomographic measurements, respectively. It is assumed that there are Ni excitations with n=1,…,Ni, and Mi measurements for each excitation with m=1,…Mi, where i=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 (mn). For the EM modality there are also Q frequency domain measurements with q=1,…,Q.The unknown, real parameter function (distributed parameter) of interest isθ(r)=[ϵr(r)σ(r)]T where ϵr(r) is the relative permittivity and σ(r) the conductivity, and where both parameter functions are defined over some spatial domain Ω with r∈Ω. Here, (·)T denotes the transpose operation.The physical forward models are denotedψ(1)(θ)={ψmnq(1)(θ)} and ψ(2)(θ)={ψmn(2)(θ)} in the two EM and ER tomography problems, respectively. The Fréchet derivative [12], or the Jacobian [13] of the forward operator ψ(i), is the bounded linear operator 𝒥(i) that yields the first order variation of ψ(i), denoted as(2)δψ(i)=J(i)δθ,
where i=1,2, and which are defined for each modality and measurement by(3)δψmnq(1)(θ)=∫ΩδθT(r)ψmnq(1)(θ)δθ(r)dv=∫Ωδϵr(r)ψmnq(1)(θ)δϵr(r)dv+∫Ωδσ(r)ψmnq(1)(θ)δσ(r)dv,δψmn(2)(θ)=∫ΩδθT(r)ψmn(2)(θ)δθ(r)dv=∫Ωδϵr(r)ψmn(2)(θ)δϵr(r)dv+∫Ωδσ(r)ψmn(2)(θ)δσ(r)dv,
where δθ(r)=[δϵr(r)δσ(r)]T denotes a perturbation of θ(r) and δθ(r)ψmnq(1)(θ) and δθ(r)ψmn(2)(θ) the corresponding gradients [12]. The following notation is employed for the Jacobian 𝒥 of the combined multiphysics problem(4)δψ=(δψ(1)δψ(2))=(J(1)J(2))δθ=(Jϵr(1)Jσ(1)0Jσ(2))(δϵrδσ)=Jδθ,
where ψ(θ)=[ψ(1)(θ)ψ(2)(θ)]T, and where it has been emphasized that the electric resistive (ER) tomography measurements depend only on the conductivity parameter σ(r).The measurement noise in the statistical model (1) is denoted wmnq(1) and wmn(2), and is assumed to be the samples of a zero mean and uncorrelated multivariate Gaussian random variable [14] with variance(5)Rmnq(1)=E{|wmnq(1)|2},Rmn(2)=E{|wmn(2)|2},
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 ℰ{(wmnq(1))2}=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

Letlnp(ξ∣θ) denote the loglikelihood function of the measurement statistics [14]. The Fisher information integral operator is then defined by the outer product(6)I(r′,r′′)=E{δθ(r′)ln⁡p(ξ∣θ)δθ(r′′)ln⁡p(ξ∣θ)},r′,r′′∈Ω,
where ℰ{·} denotes the expectation and δθ(r) 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(7)F(θ)=-ln⁡p(ξ∣θ)=∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)|ψmnq(1)(θ)-ξmnq(1)|2+∑m,n=1,1M2,N212Rmn(2)|ψmn(2)(θ)-ξmn(2)|2,(8)I(r′,r′′)=2Re{∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)δθ(r′)ψmnq(1)*(θ)δθT(r′′)ψmnq(1)(θ)}+∑m,n=1,1M2,N21Rmn(2)δθ(r′)ψmn(2)(θ)δθT(r′′)ψmn(2)(θ),
where (·)* denotes the complex conjugate, and r′,r′′∈Ω,  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=[ϵr1σ1]T and θ2=[ϵr2σ2]T, the following scalar product is employed(9)〈θ1,θ2〉=∫Ωθ1T(r)θ2(r)dv=∫Ωϵr1(r)ϵr2(r)dv+∫Ωσ1(r)σ2(r)dv.
For the mixed complex and real valued measurement data ξ=[ξ(1)ξ(2)]T and η=[η(1)η(2)]T, the following real scalar product is defined(10)〈ξ,η〉d=Re{∑m,n,q=1M1,N1,Q1Rmnq(1)ξmnq(1)*ηmnq(1)}+∑m,n=1M2,N212Rmn(2)ξmn(2)ηmn(2),
where Rmnq(1) and Rmn(2) 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 n-dimensional vector space ℂn over the complex scalar field and with scalar product xHy is isomorphic with the 2n-dimensional vector space consisting of complex vectors over the real scalar field and with real scalar product RexHy. Note that x and ix are orthogonal in the latter space.For later use, the adjoint forward operator𝒥* is now given by(11)J*ξ=〈δθ(r)ψ(θ),ξ〉d=Re{∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)δθ(r)ψmnq(1)*(θ)ξmnq(1)}+∑m,n=1,1M2,N212Rmn(2)δθ(r)ψmn(2)(θ)ξmn(2)=J(1)*ξ(1)+J(2)*ξ(2),
where 〈𝒥δθ,ξ〉d=〈δθ,𝒥*ξ〉, see also [10–12, 16].Based on the scalar product (10), the negative loglikelihood function (7) can now be written(12)F(θ)=〈ψ(θ)-ξ,ψ(θ)-ξ〉d.
By using the adjoint operator defined in (11), the first and second order variations of F(θ) are obtained as(13)δF=〈δθ,2J*(ψ-ξ)〉,δ2F=〈δθ,2J*Jδθ〉+2〈δ2ψ,ψ-ξ〉d,
where δθ is the perturbation, ψ=ψ(θ) the forward operator evaluated at some known background parameter value θ and δ2ψ the second order variation of the forward operator. It is observed that the Fisher information operator defined in (8) is given by (14)I=2J*J=2J(1)*J(1)+2J(2)*J(2)=I(1)+I(2),
where ℐ(1)=2𝒥(1)*𝒥(1) and ℐ(2)=2𝒥(2)*𝒥(2) are the Fisher information operators for each modality. The gradient δθ(r)F(θ) of the cost function F(θ) is furthermore given by(15)δθ(r)F(θ)=2J*(ψ-ξ)=2J(1)*(ψ(1)-ξ(1))+2J(2)*(ψ(2)-ξ(2)),
where(16)2J(1)*(ψ(1)-ξ(1))=2Re∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)δθ(r)ψmnq(1)*(θ)(ψmnq(1)(θ)-ξmnq(1)),2J(2)*(ψ(2)-ξ(2))=∑m,n=1,1M2,N21Rmn(2)δθ(r)ψmn(2)(θ)(ψmn(2)(θ)-ξmn(2)).A linearization is now considered where it is assumed that the second order variation of the forward operatorδ2ψ can be neglected. Consider hence the quadratic functional(17)ΔF(δθ)=δF(δθ)+12δ2F(δθ)=〈δθ,δθ(r)F(θ)〉+12〈δθ,Iδθ〉,
where the linear and quadratic terms in (13) have been used with δ2ψ=0. Note that the Hessian of F is given here by the Fisher information ℐ=2𝒥*𝒥. A regularized one-step Newton algorithm for this functional is hence given by(18)δθ+=-I+δθ(r)F(θ),
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)Jδθ=ξ-ψ,
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(20)2Jvl=σlul,2J*ul=σlvl,
where {ul} and {vl} are orthonormal systems and σl 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 μl=σl2 and eigenfunctions vl,  that is,(21)Ivl=μlvl.
The operator 𝒥 can now be represented by the singular value decomposition (SVD)(22)Jδθ=12∑lσlul〈vl,δθ〉,
and the pseudoinverse (if it exists) is given by(23)J+(ξ-ψ)=2∑l1σlvl〈ul,ξ-ψ〉d=2∑l1σlvl〈2σlJvl,ξ-ψ〉d=∑l1σl2vl〈vl,2J*(ξ-ψ)〉=-∑l1σl2vl〈vl,δθ(r)F(θ)〉=-I+δθ(r)F(θ),
which is the same as (18).The relationship between ML, SVD, and FI has been described. This will be used in Section4 to calculate the pseudoinverse and find the parameters.


## 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)ξmnq(1)=ψmnq(1)(θ)+wmnq(1),m=1,…,M1,n=1,…,N1,q=1,…,Q,ξmn(2)=ψmn(2)(θ)+wmn(2),m=1,…,M2,n=1,…,N2.
Here, ξmnq(1) are complex valued data (dynamic electrical fields) obtained in the frequency domain and ξmn(2) are real valued data (static electrical potentials), which are obtained from simultaneous EM and ER tomographic measurements, respectively. It is assumed that there are Ni excitations with n=1,…,Ni, and Mi measurements for each excitation with m=1,…Mi, where i=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 (mn). For the EM modality there are also Q frequency domain measurements with q=1,…,Q.The unknown, real parameter function (distributed parameter) of interest isθ(r)=[ϵr(r)σ(r)]T where ϵr(r) is the relative permittivity and σ(r) the conductivity, and where both parameter functions are defined over some spatial domain Ω with r∈Ω. Here, (·)T denotes the transpose operation.The physical forward models are denotedψ(1)(θ)={ψmnq(1)(θ)} and ψ(2)(θ)={ψmn(2)(θ)} in the two EM and ER tomography problems, respectively. The Fréchet derivative [12], or the Jacobian [13] of the forward operator ψ(i), is the bounded linear operator 𝒥(i) that yields the first order variation of ψ(i), denoted as(2)δψ(i)=J(i)δθ,
where i=1,2, and which are defined for each modality and measurement by(3)δψmnq(1)(θ)=∫ΩδθT(r)ψmnq(1)(θ)δθ(r)dv=∫Ωδϵr(r)ψmnq(1)(θ)δϵr(r)dv+∫Ωδσ(r)ψmnq(1)(θ)δσ(r)dv,δψmn(2)(θ)=∫ΩδθT(r)ψmn(2)(θ)δθ(r)dv=∫Ωδϵr(r)ψmn(2)(θ)δϵr(r)dv+∫Ωδσ(r)ψmn(2)(θ)δσ(r)dv,
where δθ(r)=[δϵr(r)δσ(r)]T denotes a perturbation of θ(r) and δθ(r)ψmnq(1)(θ) and δθ(r)ψmn(2)(θ) the corresponding gradients [12]. The following notation is employed for the Jacobian 𝒥 of the combined multiphysics problem(4)δψ=(δψ(1)δψ(2))=(J(1)J(2))δθ=(Jϵr(1)Jσ(1)0Jσ(2))(δϵrδσ)=Jδθ,
where ψ(θ)=[ψ(1)(θ)ψ(2)(θ)]T, and where it has been emphasized that the electric resistive (ER) tomography measurements depend only on the conductivity parameter σ(r).The measurement noise in the statistical model (1) is denoted wmnq(1) and wmn(2), and is assumed to be the samples of a zero mean and uncorrelated multivariate Gaussian random variable [14] with variance(5)Rmnq(1)=E{|wmnq(1)|2},Rmn(2)=E{|wmn(2)|2},
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 ℰ{(wmnq(1))2}=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

Letlnp(ξ∣θ) denote the loglikelihood function of the measurement statistics [14]. The Fisher information integral operator is then defined by the outer product(6)I(r′,r′′)=E{δθ(r′)ln⁡p(ξ∣θ)δθ(r′′)ln⁡p(ξ∣θ)},r′,r′′∈Ω,
where ℰ{·} denotes the expectation and δθ(r) 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(7)F(θ)=-ln⁡p(ξ∣θ)=∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)|ψmnq(1)(θ)-ξmnq(1)|2+∑m,n=1,1M2,N212Rmn(2)|ψmn(2)(θ)-ξmn(2)|2,(8)I(r′,r′′)=2Re{∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)δθ(r′)ψmnq(1)*(θ)δθT(r′′)ψmnq(1)(θ)}+∑m,n=1,1M2,N21Rmn(2)δθ(r′)ψmn(2)(θ)δθT(r′′)ψmn(2)(θ),
where (·)* denotes the complex conjugate, and r′,r′′∈Ω,  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=[ϵr1σ1]T and θ2=[ϵr2σ2]T, the following scalar product is employed(9)〈θ1,θ2〉=∫Ωθ1T(r)θ2(r)dv=∫Ωϵr1(r)ϵr2(r)dv+∫Ωσ1(r)σ2(r)dv.
For the mixed complex and real valued measurement data ξ=[ξ(1)ξ(2)]T and η=[η(1)η(2)]T, the following real scalar product is defined(10)〈ξ,η〉d=Re{∑m,n,q=1M1,N1,Q1Rmnq(1)ξmnq(1)*ηmnq(1)}+∑m,n=1M2,N212Rmn(2)ξmn(2)ηmn(2),
where Rmnq(1) and Rmn(2) 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 n-dimensional vector space ℂn over the complex scalar field and with scalar product xHy is isomorphic with the 2n-dimensional vector space consisting of complex vectors over the real scalar field and with real scalar product RexHy. Note that x and ix are orthogonal in the latter space.For later use, the adjoint forward operator𝒥* is now given by(11)J*ξ=〈δθ(r)ψ(θ),ξ〉d=Re{∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)δθ(r)ψmnq(1)*(θ)ξmnq(1)}+∑m,n=1,1M2,N212Rmn(2)δθ(r)ψmn(2)(θ)ξmn(2)=J(1)*ξ(1)+J(2)*ξ(2),
where 〈𝒥δθ,ξ〉d=〈δθ,𝒥*ξ〉, see also [10–12, 16].Based on the scalar product (10), the negative loglikelihood function (7) can now be written(12)F(θ)=〈ψ(θ)-ξ,ψ(θ)-ξ〉d.
By using the adjoint operator defined in (11), the first and second order variations of F(θ) are obtained as(13)δF=〈δθ,2J*(ψ-ξ)〉,δ2F=〈δθ,2J*Jδθ〉+2〈δ2ψ,ψ-ξ〉d,
where δθ is the perturbation, ψ=ψ(θ) the forward operator evaluated at some known background parameter value θ and δ2ψ the second order variation of the forward operator. It is observed that the Fisher information operator defined in (8) is given by (14)I=2J*J=2J(1)*J(1)+2J(2)*J(2)=I(1)+I(2),
where ℐ(1)=2𝒥(1)*𝒥(1) and ℐ(2)=2𝒥(2)*𝒥(2) are the Fisher information operators for each modality. The gradient δθ(r)F(θ) of the cost function F(θ) is furthermore given by(15)δθ(r)F(θ)=2J*(ψ-ξ)=2J(1)*(ψ(1)-ξ(1))+2J(2)*(ψ(2)-ξ(2)),
where(16)2J(1)*(ψ(1)-ξ(1))=2Re∑m,n,q=1,1,1M1,N1,Q1Rmnq(1)δθ(r)ψmnq(1)*(θ)(ψmnq(1)(θ)-ξmnq(1)),2J(2)*(ψ(2)-ξ(2))=∑m,n=1,1M2,N21Rmn(2)δθ(r)ψmn(2)(θ)(ψmn(2)(θ)-ξmn(2)).A linearization is now considered where it is assumed that the second order variation of the forward operatorδ2ψ can be neglected. Consider hence the quadratic functional(17)ΔF(δθ)=δF(δθ)+12δ2F(δθ)=〈δθ,δθ(r)F(θ)〉+12〈δθ,Iδθ〉,
where the linear and quadratic terms in (13) have been used with δ2ψ=0. Note that the Hessian of F is given here by the Fisher information ℐ=2𝒥*𝒥. A regularized one-step Newton algorithm for this functional is hence given by(18)δθ+=-I+δθ(r)F(θ),
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)Jδθ=ξ-ψ,
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(20)2Jvl=σlul,2J*ul=σlvl,
where {ul} and {vl} are orthonormal systems and σl 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 μl=σl2 and eigenfunctions vl,  that is,(21)Ivl=μlvl.
The operator 𝒥 can now be represented by the singular value decomposition (SVD)(22)Jδθ=12∑lσlul〈vl,δθ〉,
and the pseudoinverse (if it exists) is given by(23)J+(ξ-ψ)=2∑l1σlvl〈ul,ξ-ψ〉d=2∑l1σlvl〈2σlJvl,ξ-ψ〉d=∑l1σl2vl〈vl,2J*(ξ-ψ)〉=-∑l1σl2vl〈vl,δθ(r)F(θ)〉=-I+δθ(r)F(θ),
which is the same as (18).The relationship between ML, SVD, and FI has been described. This will be used in Section4 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 asE=ψ(ρ,ϕ)ẑ where the scalar field ψ(ρ,ϕ) is independent of the z-coordinate. Here, ρ and ϕ denote the radial and angular circular coordinates, respectively, and ρ=ρρ̂ the two-dimensional radius vector. Let k=ωμ0ϵ0, η0=μ0/ϵ0, μ0 and ϵ0 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 eiωt.The EM-fieldψn(1)(ρ,ϕ) resulting from excitation index n is obtained as the solution to the wave equation(24)(∇2+k2ϵ)ψn(1)=ikη0δ(ρ-ρn),n=1,…,N1,
where ϵ=ϵr-iση0/k 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 ρn.Letϵ=ϵD+δϵ where ϵD=ϵrD-iσDη0/k and δϵ=δϵr-iδση0/k denote the corresponding background parameter values and a perturbation δθ=[δϵr,δσ)]T, respectively. For the homogeneous background ϵD, the solution to (24) is given by(25)ψn(1)(ρ)=-ikη0GD(1)(ρ,ρn),
where GD(1)(ρ,ρ′) is green’s function for the homogeneous background satisfying (∇2+k2ϵD)GD(1)(ρ,ρ′)=-δ(ρ-ρ′), and which is given by(26)GD(1)(ρ,ρ′)=14i∑mJm(kϵDρ<)Hm(2)(kϵDρ>)eim(ϕ-ϕ′),
where Jm(·) and Hm(2)(·) denote the Bessel function and the Hankel function of the second kind, respectively, both of order m [18]. Here, ρ<=min{ρ,ρ′} and ρ>=max{ρ,ρ′}.A first order perturbation analysis of (24) yields the following wave equation for the perturbed field δψn(1)(27)(∇2+k2ϵ)δψn(1)=-k2δϵψn(1),
where ψn(1) 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(28)δψn(1)(ρ)=k2∫ΩGD(1)(ρ,ρ′)ψn(1)(ρ′)δϵ(ρ′)dv′=-ik3η0∫ΩGD(1)(ρ,ρ′)GD(1)(ρ′,ρn)δϵ(ρ′)dv′.
It should be noted that this is the same result, that is, obtained by employing a Born approximation [4], in which case ψn(1)(ρ) represents theincident field and δψn(1)(ρ) the scattered field.Letψmn(1) denote the measured quantity ψmn(1)=ψn(1)(ρmn) where ρmn is the measurement position at measurement m with excitation n. The corresponding gradients are then given by(29)δϵr(ρ)ψmn(1)=-ik3η0GD(1)(ρmn,ρ)GD(1)(ρn,ρ),δσ(ρ)ψmn(1)=-k2η02GD(1)(ρmn,ρ)GD(1)(ρn,ρ),
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 denotedGc(1)(ρ,ρ′), and the measured quantity is hence modeled by(30)ψmn(1)=ψn(1)(ρmn)=-ikη0Gc(1)(ρmn,ρn).
Green’s function Gc(1)(ρ,ρ′) for a homogeneous cylinder of radius a 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Φ(2)(ρ,ϕ) is independent of the z-coordinate, and where ρ and ϕ are used to denote the radial and angular circular coordinates, as above. The potential field resulting from excitation index n is obtained here as the solution to the following differential equation(31)∇⋅(σ∇Φn(2))=δ(ρ-ρn)-δ(ρ-ρN2+1),n=1,…,N2,
where σ is the conductivity and the excitation is defined by current point sources at positions ρn with a common ground connection at ρN2+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(32)Φn(2)(ρ)=GD(2)(ρ,ρn)-GD(2)(ρ,ρN2+1),
where GD(2)(ρ,ρ′) is green’s function for the homogeneous background satisfying ∇·(σD∇GD(2)(ρ,ρ′))=δ(ρ-ρ′), and which is given by(33)GD(2)(ρ,ρ′)=12πσDln⁡ρ>-14πσD∑m≠01|m|(ρ<ρ>)|m|eim(ϕ-ϕ′),
where ρ<=min{ρ,ρ′} and ρ>=max{ρ,ρ′}, see also [18].The observed voltage quantityψmn(2) corresponding to the homogeneous background is given by(34)ψmn(2)=Φn(2)(ρmn)-Φn(2)(ρ0)=GD(2)(ρmn,ρn)-GD(2)(ρmn,ρN2+1)-GD(2)(ρ0,ρn)+GD(2)(ρ0,ρN2+1),
where ρmn is the measurement position at measurement m and excitation n, and ρ0 is a position used as a voltage reference.Note that green’s functionGD(2)(ρ,ρn) satisfies the differential equation(35)∇⋅(σD∇GD(2)(ρ,ρn))=δ(ρ-ρn).
A first order perturbation analysis of (35) yields the following differential equation for the perturbed field δGD(2)(ρ,ρn)(36)∇⋅(σD∇δGD(2)(ρ,ρn))=-∇⋅(δσ∇GD(2)(ρ,ρn)).
The solution to (36) is given by(37)δGD(2)(ρ,ρn)=-∫ΩGD(2)(ρ,ρ′)∇′⋅(δσ∇′GD(2)(ρ′,ρn))dv′=∫Ω∇′GD(2)(ρ,ρ′)⋅∇′GD(2)(ρn,ρ′)δσ(ρ′)dv′,
where Green's theorem [18] has been used together with the assumption that δσ=0 on the boundary of Ω. The gradient of green’s function GD(2)(ρl,ρn) is hence given by(38)δσ(ρ)GD(2)(ρl,ρn)=∇GD(2)(ρl,ρ)⋅∇GD(2)(ρn,ρ),
where ∇GD(2)(ρn,ρ) can be obtained from (33) as(39)∇GD(2)(ρn,ρ)=ρ̂{-14πσDρn∑m≠0(ρρn)|m|-1eim(ϕn-ϕ)}+ϕ̂{i4πσDρ∑m≠0(m|m|)(ρρn)|m|eim(ϕn-ϕ)},
where (ρn,ϕn) are the circular coordinates of ρn, and ρ<ρn. The gradient of the observed voltage quantity ψmn(2) defined in (34) can finally be obtained as(40)δσ(ρ)ψmn(2)=δσ(ρ)GD(2)(ρmn,ρn)-δσ(ρ)GD(2)(ρmn,ρN2+1)-δσ(ρ)GD(2)(ρ0,ρn)+δσ(ρ)GD(2)(ρ0,ρN2+1),
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 denotedGc(2)(ρ,ρ′), and the measured quantity is hence modeled by(41)ψmn(2)=Gc(2)(ρmn,ρn)-Gc(2)(ρmn,ρN2+1)-Gc(2)(ρ0,ρn)+Gc(2)(ρ0,ρN2+1),
similar to (34). Green’s function Gc(2)(ρ,ρ′) for a homogeneous cylinder of radius a placed in a homogeneous background space is given in Appendix B, and is used to generate the synthetic data in the numerical example.


## 3.1. The Electromagnetic (EM) Tomography Problem

A two-dimensional electromagnetic tomography problem is considered where the harmonic electric field is polarized asE=ψ(ρ,ϕ)ẑ where the scalar field ψ(ρ,ϕ) is independent of the z-coordinate. Here, ρ and ϕ denote the radial and angular circular coordinates, respectively, and ρ=ρρ̂ the two-dimensional radius vector. Let k=ωμ0ϵ0, η0=μ0/ϵ0, μ0 and ϵ0 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 eiωt.The EM-fieldψn(1)(ρ,ϕ) resulting from excitation index n is obtained as the solution to the wave equation(24)(∇2+k2ϵ)ψn(1)=ikη0δ(ρ-ρn),n=1,…,N1,
where ϵ=ϵr-iση0/k 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 ρn.Letϵ=ϵD+δϵ where ϵD=ϵrD-iσDη0/k and δϵ=δϵr-iδση0/k denote the corresponding background parameter values and a perturbation δθ=[δϵr,δσ)]T, respectively. For the homogeneous background ϵD, the solution to (24) is given by(25)ψn(1)(ρ)=-ikη0GD(1)(ρ,ρn),
where GD(1)(ρ,ρ′) is green’s function for the homogeneous background satisfying (∇2+k2ϵD)GD(1)(ρ,ρ′)=-δ(ρ-ρ′), and which is given by(26)GD(1)(ρ,ρ′)=14i∑mJm(kϵDρ<)Hm(2)(kϵDρ>)eim(ϕ-ϕ′),
where Jm(·) and Hm(2)(·) denote the Bessel function and the Hankel function of the second kind, respectively, both of order m [18]. Here, ρ<=min{ρ,ρ′} and ρ>=max{ρ,ρ′}.A first order perturbation analysis of (24) yields the following wave equation for the perturbed field δψn(1)(27)(∇2+k2ϵ)δψn(1)=-k2δϵψn(1),
where ψn(1) 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(28)δψn(1)(ρ)=k2∫ΩGD(1)(ρ,ρ′)ψn(1)(ρ′)δϵ(ρ′)dv′=-ik3η0∫ΩGD(1)(ρ,ρ′)GD(1)(ρ′,ρn)δϵ(ρ′)dv′.
It should be noted that this is the same result, that is, obtained by employing a Born approximation [4], in which case ψn(1)(ρ) represents theincident field and δψn(1)(ρ) the scattered field.Letψmn(1) denote the measured quantity ψmn(1)=ψn(1)(ρmn) where ρmn is the measurement position at measurement m with excitation n. The corresponding gradients are then given by(29)δϵr(ρ)ψmn(1)=-ik3η0GD(1)(ρmn,ρ)GD(1)(ρn,ρ),δσ(ρ)ψmn(1)=-k2η02GD(1)(ρmn,ρ)GD(1)(ρn,ρ),
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 denotedGc(1)(ρ,ρ′), and the measured quantity is hence modeled by(30)ψmn(1)=ψn(1)(ρmn)=-ikη0Gc(1)(ρmn,ρn).
Green’s function Gc(1)(ρ,ρ′) for a homogeneous cylinder of radius a 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Φ(2)(ρ,ϕ) is independent of the z-coordinate, and where ρ and ϕ are used to denote the radial and angular circular coordinates, as above. The potential field resulting from excitation index n is obtained here as the solution to the following differential equation(31)∇⋅(σ∇Φn(2))=δ(ρ-ρn)-δ(ρ-ρN2+1),n=1,…,N2,
where σ is the conductivity and the excitation is defined by current point sources at positions ρn with a common ground connection at ρN2+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(32)Φn(2)(ρ)=GD(2)(ρ,ρn)-GD(2)(ρ,ρN2+1),
where GD(2)(ρ,ρ′) is green’s function for the homogeneous background satisfying ∇·(σD∇GD(2)(ρ,ρ′))=δ(ρ-ρ′), and which is given by(33)GD(2)(ρ,ρ′)=12πσDln⁡ρ>-14πσD∑m≠01|m|(ρ<ρ>)|m|eim(ϕ-ϕ′),
where ρ<=min{ρ,ρ′} and ρ>=max{ρ,ρ′}, see also [18].The observed voltage quantityψmn(2) corresponding to the homogeneous background is given by(34)ψmn(2)=Φn(2)(ρmn)-Φn(2)(ρ0)=GD(2)(ρmn,ρn)-GD(2)(ρmn,ρN2+1)-GD(2)(ρ0,ρn)+GD(2)(ρ0,ρN2+1),
where ρmn is the measurement position at measurement m and excitation n, and ρ0 is a position used as a voltage reference.Note that green’s functionGD(2)(ρ,ρn) satisfies the differential equation(35)∇⋅(σD∇GD(2)(ρ,ρn))=δ(ρ-ρn).
A first order perturbation analysis of (35) yields the following differential equation for the perturbed field δGD(2)(ρ,ρn)(36)∇⋅(σD∇δGD(2)(ρ,ρn))=-∇⋅(δσ∇GD(2)(ρ,ρn)).
The solution to (36) is given by(37)δGD(2)(ρ,ρn)=-∫ΩGD(2)(ρ,ρ′)∇′⋅(δσ∇′GD(2)(ρ′,ρn))dv′=∫Ω∇′GD(2)(ρ,ρ′)⋅∇′GD(2)(ρn,ρ′)δσ(ρ′)dv′,
where Green's theorem [18] has been used together with the assumption that δσ=0 on the boundary of Ω. The gradient of green’s function GD(2)(ρl,ρn) is hence given by(38)δσ(ρ)GD(2)(ρl,ρn)=∇GD(2)(ρl,ρ)⋅∇GD(2)(ρn,ρ),
where ∇GD(2)(ρn,ρ) can be obtained from (33) as(39)∇GD(2)(ρn,ρ)=ρ̂{-14πσDρn∑m≠0(ρρn)|m|-1eim(ϕn-ϕ)}+ϕ̂{i4πσDρ∑m≠0(m|m|)(ρρn)|m|eim(ϕn-ϕ)},
where (ρn,ϕn) are the circular coordinates of ρn, and ρ<ρn. The gradient of the observed voltage quantity ψmn(2) defined in (34) can finally be obtained as(40)δσ(ρ)ψmn(2)=δσ(ρ)GD(2)(ρmn,ρn)-δσ(ρ)GD(2)(ρmn,ρN2+1)-δσ(ρ)GD(2)(ρ0,ρn)+δσ(ρ)GD(2)(ρ0,ρN2+1),
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 denotedGc(2)(ρ,ρ′), and the measured quantity is hence modeled by(41)ψmn(2)=Gc(2)(ρmn,ρn)-Gc(2)(ρmn,ρN2+1)-Gc(2)(ρ0,ρn)+Gc(2)(ρ0,ρN2+1),
similar to (34). Green’s function Gc(2)(ρ,ρ′) for a homogeneous cylinder of radius a 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 Figure1, which also shows green’s functions Gc(1)(ρ,ρ′) and Gc(2)(ρ,ρ′) 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·106 (200 MHz) is employed with the EM modality. The large circle in the center indicates the simulated tunnel of radius a=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 ϵrD=5 and conductivity σD=1 mS/m (σD=0.001 S/m). The tunnel is modeled with relative permittivity ϵra=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.Illustration of measurement setup and Green's functionsGc(1)(ρ,ρ′) and Gc(2)(ρ,ρ′) in electromagnetic and electrical resistive tomography, respectively. The large circle in the center indicates the tunnel of radius a=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.
(a)(b)The signal-to-noise-ratio (SNR) in the measurement model (1) is defined as(42)SNR(1)=max⁡mnq⁡|ψmnq(1)|2R(1),SNR(2)=max⁡mn⁡|ψmn(2)|2R(2),
where ψmnq(1) and ψmn(2) represent the (noiseless) measured quantities, and R(1) and R(2) 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 (m,n,q), and hence Rmnq(1)=R(1) and Rmn(2)=R(2). In the simulation examples described below, the measurement noise was chosen such that SNR(1)=80 dB, SNR(2)=70 dB and R(2)/R(1)=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 R(2)/R(1) which is independent of the absolute value of the SNR.In Figures2 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 Neig=15 and Neig=30 eigenvalues, respectively.Reconstruction ofδϵr and δσ in a straightforward data fusion without statistically based weighting (the weighting is assumed to be uniform with R(2)/R(1)=1). The upper and lower rows are obtained by using Neig=15 and Neig=30 eigenvalues, respectively.
(a)(b)(c)(d)Reconstruction ofδϵr and δσ in a maximum likelihood-based data fusion where the statistically-based weighting is R(2)/R(1)=92. The upper and lower rows are obtained by using Neig=15 and Neig=30 eigenvalues, respectively.
(a)(b)(c)(d)In Figure2 is shown the reconstructions in a straightforward data fusion without statistically-based weighting. In these reconstructions, the weighting is assumed to be uniform with R(2)/R(1)=1 even though the simulation data has been generated with R(2)/R(1)=92. Hence, the correct multiphysics model has been employed, but the sensor noise statistics has not been taken into consideration.In Figure3 is shown the reconstructions in a maximum likelihood-based data fusion where the statistically-based weighting with R(2)/R(1)=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 Figures2 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, maxmn|ψmn(2)|2=9.2, maxmn|ψmn(1)|2 and SNR(2)=SNR(1)/10, yielding an imbalance in the weighting of data corresponding to the ratio R(2)/R(1)=92. As can be seen in the Figures 2 and 3, the consequence of neglecting this imbalance can be significant, especially for Neig=15.In Figure4 is shown the singular values μn(ℐ) 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μn(ℐ) 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.



---

*Source: 617089-2011-07-05.xml*