Numerical methods for electromagnetic inversion

The aim of electromagnetic (EM) sounding methods in geophysics is to obtain information about the subsurface of the earth by recorded measurements taken at the surface. In particular, the goal is to determine variations in the electrical conductivity of the earth with depth by employing an inversion procedure. In this work we focus on one technique, that consists of placing a magnetic dipole above the surface, composed of a transmitter coil and different couples of adjacent receiver coils. The receiver couples are placed at different distances (offsets) from the transmitter coil. In this setting, the electromagnetic induction effect, encoded in the first-order linear Maxwell’s differential equations, produce eddy alterning currents in the soil which induce response electromagnetic fields, that can be used to determine the conductivity profile of the ground by meaning of an inversion algorithm. A typical inversion strategy consists in an iterative procedure involving the computation of the EM response of a layered model (forward modelling) and the solution of the inverse problem. Then, the algorithm attempts to minimize the mismatch between the measured data and the predicted data, by updating the model parameters at each iteration. By assuming that the local subsurface structures are composed by horizontal and homogeneous layers, general integral solutions of Maxwell equations (i.e., the EM fields) for vertical and horizontal magnetic dipoles, can be derived and represented as Hankel transforms, which contain the subsurface model parameters, i.e., the conductivity and the thickness of each layer. By a mathematical point of view, in general, these Hankel transforms are not analytically computable and therefore it is necessary to employ a numerical scheme. Anyway, the slowly decay of the oscillations determined by the Bessel function makes the problem very difficult to handle, because traditional quadrature rules typically fail to converge. In this work we consider two different approaches. The first one is based on the decomposition of the integrand function in a first function for which the corresponding Hankel transform is known exactly, and an oscillating function decays exponentially. For realistic sets of parameters, the oscillations are quite rapidly damped, and the corresponding integral can be accurately computed by a classical quadrature rule on finite intervals. The second approach consists in the application of a Gaussian quadrature formula. We develop a Gaussian rule for weight functions involving fractional powers, exponentials and Bessel functions of the first kind. Moreover, we derive an analytical approximation of these integrals that has a general validity and allows to overcome the limits of common methods based on the modelling of apparent conductivity in the low induction number (LIN) approximation. Having at disposal a reliable method for evaluating the Hankel transforms, by assuming as forward model a homogeneous layered earth, here we also consider the inverse problem of computing the model parameters (i.e., conductivity and thickness of the layers) from a set of measured field values at different offsets. We focus on the specific case of the DUALEM system. We employ two optimization algorithms. The first one is based on the BFGS line-search method and, in order to reduce as much as possible the number of integral evaluations, the analytic approximation of these integrals is used to have a first estimate of the solution. For the second approach we employ the damped Gauss-Newton method. To avoid the dependence on the initial guess of the iterative procedure, we consider a set of different initial models, and we use each one to solve the optimization problem. The numerical experiments, carried out for the study of river-levees integrity, are obtained by employing a virtual machine equipped with the NVIDIA A100 Tensor Core GPU.