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

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.


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][3][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][6][7][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][10][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 2 International Journal of Geophysics 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.

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.

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 Here, ξ (1) mnq are complex valued data (dynamic electrical fields) obtained in the frequency domain and ξ (2) mn are real valued data (static electrical potentials), which are obtained from simultaneous EM and ER tomographic measurements, respectively. It is assumed that there are N i excitations with n = 1, . . . , N i , and M i measurements for each excitation with m = 1, . . . M i , 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 measurement noise in the statistical model (1) is denoted w (1) mnq and w (2) mn, and is assumed to be the samples of a zero mean and uncorrelated multivariate Gaussian random variable [14] with variance 2 } = 0, and for the ER-modality, the random variable is real Gaussian, see, for example, [14,15].

Maximum Likelihood and the Truncated SVD.
Let ln p(ξ | θ) denote the loglikelihood function of the measurement statistics [14]. The Fisher information integral operator is then defined by the outer product where E {·} 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 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 For the mixed complex and real valued measurement data ξ = [ξ (1) ξ (2) ] T and η = [η (1) η (2) ] T , the following real scalar product is defined where R (1) mnq and R (2) mn 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 C n over the complex scalar field and with scalar product x H y is isomorphic with the 2n-dimensional vector space consisting of complex vectors over the real scalar field and with real scalar product Re x H y. Note that x and ix are orthogonal in the latter space.
For later use, the adjoint forward operator J * is now given by where Jδθ, ξ d = δθ, J * ξ , see also [10][11][12]16]. Based on the scalar product (10), the negative loglikelihood function (7) can now be written By using the adjoint operator defined in (11), the first and second order variations of F(θ) are obtained as 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 where I (1) = 2J (1) * J (1) and I (2) = 2J (2) * J (2) are the Fisher information operators for each modality. The gradient δ θ(r) F(θ) of the cost function F(θ) is furthermore given by where 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 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 I = 2J * J. A regularized one-step Newton algorithm for this functional is hence given by where I + denotes a truncated pseudoinverse of the Fisher information I, 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 provided that the scalar products are properly chosen as in (9) and (10). To this end, it is assumed that the operator J is compact and hence that the following singular system exists where {u l } and {v l } 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 I = 2J * J with eigenvalues μ l = σ 2 l and eigenfunctions v l , that is, The operator J can now be represented by the singular value decomposition (SVD) and the pseudoinverse (if it exists) is given by 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.

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.

The Electromagnetic (EM) Tomography Problem.
A twodimensional electromagnetic tomography problem is considered where the harmonic electric field is polarized as E = ψ(ρ, φ) z where the scalar field ψ(ρ, φ) is independent of the z-coordinate. Here, ρ and φ denote the radial and angular circular coordinates, respectively, and ρ = ρ ρ the twodimensional 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 e iωt .
The EM-field ψ (1) n (ρ, φ) resulting from excitation index n is obtained as the solution to the wave equation where = r − iση 0 /k is the complex valued relative permittivity, r the relative permittivity, and σ the conductivity.
A first order perturbation analysis of (24) yields the following wave equation for the perturbed field δψ (1) where ψ (1) n 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 It should be noted that this is the same result, that is, obtained by employing a Born approximation [4], in which case ψ (1) n (ρ) represents the incident field and δψ (1) n (ρ) the scattered field.
Let ψ (1) mn denote the measured quantity ψ (1) mn = ψ (1) n (ρ mn ) where ρ mn is the measurement position at measurement m with excitation n. The corresponding gradients are then given by 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 G (1) c (ρ, ρ ), and the measured quantity is hence modeled by Green's function G (1) c (ρ, ρ ) 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.

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 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].
The observed voltage quantity ψ (2) mn corresponding to the homogeneous background is given by 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 function G (2) D (ρ, ρ n ) satisfies the differential equation A first order perturbation analysis of (35) yields the following differential equation for the perturbed field δG (2) D (ρ, ρ n ) The solution to (36) is given by 6 International Journal of Geophysics where Green's theorem [18] has been used together with the assumption that δσ = 0 on the boundary of Ω. The gradient of green's function G (2) D (ρ l , ρ n ) is hence given by where ∇G (2) D (ρ n , ρ) can be obtained from (33) as where (ρ n , φ n ) are the circular coordinates of ρ n , and ρ < ρ n . The gradient of the observed voltage quantity ψ (2) mn defined in (34) can finally be obtained as 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 G (2) c (ρ, ρ ), and the measured quantity is hence modeled by similar to (34). Green's function G (2) c (ρ, ρ ) 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.

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 twodimensional 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 G (1) c (ρ, ρ ) and G (2) 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 6 (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 Green's function G (2) Figure 1: Illustration of measurement setup and Green's functions G (1) c (ρ, ρ ) and G (2) c (ρ, ρ ) 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.
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.
The signal-to-noise-ratio (SNR) in the measurement model (1) is defined as where ψ (1) mnq and ψ (2) mn 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 R (1) mnq = R (1) and R (2) mn = R (2) . In the simulation examples described below, the measurement noise was chosen such that SNR (1)  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 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 N eig = 15 and N 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 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 Figure 3 is shown the reconstructions in a maximum likelihood-based data fusion where the statisticallybased 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 Figures 2 and 3 in this example illustrates the expected behavior of the data fusion for an illposed multiphysics inverse problem, that is, the sensitivity of the inversion with respect to the discrepancies in the statistical assumptions. In this example, max mn |ψ (2) mn| 2 = 9.2, max mn |ψ (1) mn| 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 N eig = 15.
In Figure 4 is shown the singular values μ n (I) 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].

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. The appropriate boundary conditions are obtained from the continuity of G m (ρ, ρ , φ ) and the discontinuity of (∂/∂ρ) G m (ρ, ρ , φ ) at ρ = a and ρ = ρ . Hence, 2π .