#### 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 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 , and measurements for each excitation with , where 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 .

The unknown, real parameter function (distributed parameter) of interest is where is the relative permittivity and the conductivity, and where both parameter functions are defined over some spatial domain with . Here, 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 where , and which are defined for each modality and measurement by where denotes a perturbation of and and the corresponding gradients [12]. The following notation is employed for the Jacobian of the combined multiphysics problem where , 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 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 , 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 denote the loglikelihood function of the measurement statistics [14]. The Fisher information integral operator is then defined by the outer product 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 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 and , the following scalar product is employed For the mixed complex and real valued measurement data and , the following real scalar product is defined 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 is isomorphic with the -dimensional vector space consisting of complex vectors over the real scalar field and with real scalar product . Note that and are orthogonal in the latter space.

For later use, the adjoint forward operator is now given by where , see also [10–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 are obtained as 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 where and are the Fisher information operators for each modality. The gradient of the cost function is furthermore given by where

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 where the linear and quadratic terms in (13) have been used with . Note that the Hessian of is given here by the Fisher information . A regularized one-step Newton algorithm for this functional is hence given by 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 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 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 with eigenvalues and eigenfunctions , that is, The operator 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.

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

The EM-field resulting from excitation index is obtained as the solution to the wave equation where is the complex valued relative permittivity, 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 where and denote the corresponding background parameter values and a perturbation , respectively. For the homogeneous background , the solution to (24) is given by where is green’s function for the homogeneous background satisfying , and which is given by where and denote the Bessel function and the Hankel function of the second kind, respectively, both of order [18]. Here, and .

A first order perturbation analysis of (24) yields the following wave equation for the perturbed field
where is the solution to (24) and the perturbation in parameter values. For the homogeneous background with , 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 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 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 , and the measured quantity is hence modeled by Green’s function 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 where is the conductivity and the excitation is defined by current point sources at positions with a common ground connection at , see also [17].

Let where and denote the corresponding background parameter value and a perturbation, respectively. For the homogeneous background , the solution to (31) is given by where is green’s function for the homogeneous background satisfying , and which is given by where and , see also [18].

The observed voltage quantity corresponding to the homogeneous background is given by where is the measurement position at measurement and excitation , and is a position used as a voltage reference.

Note that green’s function satisfies the differential equation A first order perturbation analysis of (35) yields the following differential equation for the perturbed field The solution to (36) is given by where Green's theorem [18] has been used together with the assumption that on the boundary of . The gradient of green’s function is hence given by where can be obtained from (33) as where are the circular coordinates of , and . The gradient of the observed voltage quantity 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 , and the measured quantity is hence modeled by similar to (34). Green’s function 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 and 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 (200 MHz) is employed with the EM modality. The large circle in the center indicates the simulated tunnel of radius 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 and conductivity mS/m ( S/m). The tunnel is modeled with relative permittivity and conductivity . 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)**

The signal-to-noise-ratio (SNR) in the measurement model (1) is defined as 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 dB, dB and . 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 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 and eigenvalues, respectively.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

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 even though the simulation data has been generated with . 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 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, , and , yielding an imbalance in the weighting of data corresponding to the ratio . As can be seen in the Figures 2 and 3, the consequence of neglecting this imbalance can be significant, especially for .

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

#### 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 for a homogenous cylinder of radius is derived below. The complex valued permittivity of the background and the cylinder are given by and , respectively.

Green’s function satisfies the scalar wave equation together with appropriate boundary conditions. By introducing the Fourier series expansion the wave equation (A.1) transforms to

Assume that the source is outside the cylinder, that is, . The solution to (A.3) is then given by The appropriate boundary conditions are given by The boundary conditions (A.5) yield the linear system of equationswhich is solved for the constants , , , and defined in (A.4).

#### B. ER Green's Function for a Homogeneous Cylinder

Green’s function for a homogenous cylinder of radius is derived below. The conductivity of the background and the cylinder are denoted by and , respectively.

Green’s function satisfies the partial differential equation or where . By introducing the Fourier series expansion the differential equation (B.2) transforms to

Assume that the source is outside the cylinder, that is, . The solution to (B.4) for is then given by and for

The appropriate boundary conditions are obtained from the continuity of and the discontinuity of at and . Hence, The boundary conditions (B.7) yield the following linear system of equations for which is solved for the constants , , , and defined in (B.5). For , the boundary conditions (B.7) yield the following linear system of equations which is solved for the constants , , , and defined in (B.6). The solution is readily obtained as , and .

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