3D general-measure inversion of crosswell EM data using a direct solver


 We present a three-dimensional (3D) general-measure inversion scheme of crosswell electromagnetic (EM) data in the frequency domain with a direct forward solver. In the forward problem, we discretised the EM Helmholtz equation by the staggered-grid finite difference (SGFD) scheme and solved it using the Intel MKL PARDISO direct solver. By applying a direct solver, we simultaneously solved the multisource forward problems at a given frequency. In the inversion, we integrated a general measure of data misfit and model constraints with linearised least-squares inversion. We reconstructed a model with blocky features by selecting the appropriate measure parameters and model constraints. We used the adjoint equation method to explicitly calculate the Jacobian matrix, which facilitated the determination of an appropriate initial value for the regularisation coefficient in the objective function. We illustrated the inversion scheme using synthetic crosswell EM data with a general measure, the L2 norm, and, specifically, two mixed norms.


Introduction
The three-dimensional (3D) electromagnetic (EM) inversion problem has been a critical topic facing geophysicists over the past two decades. The impetus driving people's research comes from the application of subsurface electrical conductivity distribution. In sedimentary formations, the electrical conductivity mainly depends on the porosity, fluid conductivity and saturation, and slightly depends on the formation temperature, permeability coefficient and pore geometry (Archie 1942;Keller 1988). Therefore, important information, such as geological structure, fluid saturation, chemical composition and fracture direction, can be inferred from the conductivity image. This knowledge can be applied not only to reservoir characterisation, but also to mineral exploration, hydrogeology and chemical waste site evaluations.
At the end of the 20th century, important progress was made in dealing with the crosswell EM inversion problem. Zhou (1989) carried out pioneering work on low-frequency crosshole EM inversion by extending the acoustic inversion technique to a complex wave number. On the basis of Zhou's work, Alumbaugh and Morrison developed an EM conductivity imaging technique using the second Born approximation to improve the accuracy of inverted conductivity images (Alumbaugh & Morrison 1993). To calculate models that contain large-scale anomalies, Newman (1995) applied the finite difference method to solve the differentialequation (DE) formulation of Maxwell's equations instead of using the integral-equation (IE) approach. At the same time, the first field experiment was carried out at the Devine test site that is located southwest of San Antonio, Texas, USA (Wilt et al. 1995). The reconstructed resistivity from the one-dimensional (1D) layered model inversion was consistent with the values from the induction log at the Devine test site. By separating and storing the data of different calculation domains into hundreds of processors, Newman & Alumbaugh (1997) took advantage of the message passing interface technique and the iterative solver of conjugate gradients (CG) to perform a large-scale 3D EM inversion. Using the regularisation method of linearised data and the modified Gram-Schmidt method, Shen et al. (2008) developed a 2.5D crosshole EM inversion scheme with the forward responses calculated by the finite element method and the iterative solver of biconjugate gradient (Bi-CG). Based on an IE approach and Tikhonov regularisation, MacLennan et al. (2014) investigated the effects of complex conductivity on crosswell EM data. Colombo et al. (2020) explored the use of the U-Net deep learning network for monitoring an EM-based reservoir. Zhang et al. (2020) applied minimum support regularisation and the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm to solve the inverse problem and developed an ensemble-based history-matching framework to interpret crosswell EM data. Fang et al. (2020) investigated a radial basis function neural network algorithm and compared it with other five neural networks in a crosswell configuration.
Almost all of these inversion schemes used L 2 -norm as their measure of data misfit and model structure. With the help of structure constraints, however, their reconstructed images had a smooth and minimum structure, and L 2 -norm inversion could not recover images with blocky features. To solve this problem, researchers have published a number of articles about non-L 2 norm inversion. Egbert & Booker (1986) presented an iteratively reweighted least square (IRLS) inversion for geomagnetic induction data with a hybrid measure of L 1 -norm and L 2 -norm. Farquharson & Oldenburg (1998) investigated general measures in the linear 1D magnetic data inversion and the nonlinear 1D transient EM data inversion. Following the Ekblom norm (Ekblom 1973), Zhang et al. (2000) applied 1D IRLS EM inversion to determine the bed boundary positions, which were used as a priori information of two-dimensional (2D) inversion. By using the smoothness change information, Sun & Li (2014) applied the dynamic Ekblom norm to different parts of the model structure for automatic reconstruction of both blocky and smooth models.
Over the past decade, direct solvers have gradually become popular tools for solving forward equations. There are two main reasons: (1) direct solvers only need to factorise the system matrix once to solve solutions for multiple righthand sides (Newman 2013), which greatly reduces the timeconsuming multisource forward modelling; and (2) compared with iterative techniques, direct solvers are more stable and accurate, especially when dealing with models with low-frequency, nonuniform grids and large resistivity contrast (Gould et al. 2007;Grayver et al. 2013;Long et al. 2020). Although direct solvers have many advantages, their disadvantage is obvious: huge memory consumption. With the development of computing technology, however, the out-ofcore mode can store matrices on disk instead of in RAM (Xiong et al. 2018), so that personal computers can also run a large-scale EM inversion program with a sacrifice of time for data exchange.
The main objective of this paper was to develop a general-measure inversion scheme for 3D crosswell EM data based on a direct solver. The scheme had three key components: forward modelling, general measures and the inversion algorithm. First, we used the staggered-grid finite difference (SGFD) method to discretise the DE formulation of Maxwell's equations and apply it to the Intel MKL PARDISO direct solver to calculate large-scale sparse linear system equations. Second, we presented one of the general measures, the Ekblom norm, with its particular property. Next, we showed how to solve the nonlinear inverse problem by linearised forward responses with the inversion scheme of the Ekblom norm and the minimum support (MS) constraint. Finally, we illustrated and discussed the generalmeasure inversion scheme through synthetic data of a crosswell EM configuration.

Forward modelling
We obtained model responses of crosswell configurations by applying the approach of primary and secondary field decomposition to the vector Helmholtz equations in the frequency domain. Assuming a time-harmonic dependence of e j t and displacement currents are ignored, the equation for the secondary electric field E s in the frequency domain is written as follows: where j 2 = −1, denotes the angular frequency, represents the magnetic permeability of free space, is the conductivity, p represents the background conductivity and E p is the primary or background electric field. The total electric field E is given by Following Shen (2003), we computed the primary field E p analytically, and we approximated the secondary field E s with the SGFD discretisation method. We imposed Dirichlet boundary conditions on the approximated calculation. This resulted in a linear system of equations for a given frequency, where A is a N e × N e complex, sparse coefficient matrix, with N e = 3N x N y N z (where N i is the number of cells in each direction). N e -vector b is equal to the right-hand side of equation (1) that contains the source information of the secondary field and boundary conditions.
For the 3D crosswell EM problem, we had to solve the multisource multifrequency forward equation (3) in the forward modelling and sensitivity calculation. We used the Intel MKL PARDISO direct solver to simultaneously calculate multisource large sparse linear equations at a given frequency (Schenk & Gärtner 2002;Schenk & Gärtner 2004).
Once we obtained the secondary electric field E s , we substituted it in equation (2) to yield the total electric field E. We calculated the total magnetic field H from the total electric field E as follows:

General measures
Considering a N m -vector m = (m 1 , m 2 , ..., m N m ) T , a general measure of m is given by (Huber 1964) When (m) = m 2 or (m) = |m|, the measure becomes the L 2 or L 1 norm. When (m) = |m| p , p ≥ 1, this measure becomes the L p norm of the vector. Thus, the weight on each element of a vector m is controlled by parameter p. The larger the parameter p, the more weight is imposed on large elements of a vector (Zhang et al. 2000). Compared with the L 1 norm, the L 2 norm imposed more weight on large m i during the process of minimisation. The image reconstructed from the L 2 -norm inversion was rather smooth as a result, and limited the occurrence of large elements in the inverted model (Sun & Li 2014). For this reason, we were able to recover features of sharp boundaries by applying the L 1 measure to the derivative of this model, and we smeared out features only when using L 2 norm. In this paper, we adopted a damped L p norm proposed by Ekblom (1973): where is a positive number. For brevity, it is referred to as the Ekblom norm. Figure 1 shows the Ekblom norm with different p values (i.e. 0, 0.5, 1, 1.5 and 2) at = 10 −5 . If m ≫ , the Ekblom measure puts more weight on a large m i and less weight on a small m i when the parameter p increased, which led to a rather smooth image. As became larger, the measure behaved like a sum-of-squares measure scaled by (Farquharson & Oldenburg 1998). Compared with the L p norm, the Ekblom norm removed the singularity at m = 0 when p = 1. We applied this general measure to perform crosswell EM inversion.

Linearised inversion
Considering that geophysical inversion is nonunique and the recovered model should be geologically interpretable, the objective function always consists of a data misfit function and model structure stabilising functions (Wang 2016): where Φ g d and Φ g m are, respectively, symbols of the data misfit function and model structure function defined by the Ekblom norm. The complex N d -vector d°b s is the observed the data, f(m) is a forward response, m ref denotes the reference model, and W d is a N d × N d weighting diagonal matrix that consists of the reciprocal of the data standard deviations. We used the N m × N m weighting matrices W c and W m for model constraint. The second term in equation (7) is also known as the stabiliser, which keeps the inversion stable and constrains the reconstructed model structure (Tikhonov & Arsenin 1977;Wang 2016). is the Lagrange multiplier, and parameter determines the relative contributions of the smoothness term (i.e. W c m) and the similarity term (i.e. W m (mm ref )) to the objective function (Oldenburg et al. 1993). We set = 0.1 (we discuss the selection of later).
Once we defined the objective function, we could solve the inverse problem by minimising the nonlinear objective function with the linearised data. Assuming m k+1 is the model perturbation of m k , the Taylor series of the forward responses f(m) at m k+1 without higher-order terms is given by the following: or where J k is the Jacobian matrix at kth iteration and Δm k = m k+1 − m k . Then, the objective function of Δm k is written as follows: The minimisation of the objective function (10), that is, Φ( m k )∕ Δm = 0, yields the system of the normal equations, where r d , r c and r m are diagonal matrices (Sun & Li 2014): If we set p = 2, all three matrices would become identity matrices I. So, the system of the normal equations defined by the Ekblom norm included the one defined by the conventional L 2 norm as a special case. When p = 1 and < m, the reconstructed image obtained from equation (11) was similar to the one from the L 1 -norm inversion. Matrices r d , r c and r m acted as weighting coefficients for each element of the data-weighting matrix (W d ) and the model constraints matrices (W c and W m ).
Equation (11) can be simplified to the linear least-squares system, as follows: which results in a more accurate solution (Sasaki 2001). We used the gelsy routine of the Intel MKL LAPACK software package to calculate the updated model Δm k at the kth iteration from the least-squares system (13). The gelsy routine used the complete orthogonal decomposition of the system matrix to calculate the solution of the linear least-squares problem (Anderson et al. 1999). Then, we obtain a new model as follows: Until the misfit was less than the noise level or a given number of iterations had been reached, the inversion stopped. Other than the model constraint matrices discussed in the next section, there are two matrices/parameters that we have not interpreted: the Jacobian matrix J and the Lagrange multiplier . We used the adjoint equation method to calculate the Jacobian matrix J (McGillivray et al. 1994) and, next, we discuss the selection of .
Based on the explicit calculation of the Jacobian matrix, we developed a new cooling strategy to select the Lagrange multiplier . In the traditional cool strategy, the initial value of was so large that the model constraints dominated data misfit (Oldenburg & Li 2005). More iterations needed to be carried out until the value of decreased to a proper level. To solve this problem, we determined the initial value of through the balance equation: For the least-squares system in equation (13) at p = 1, the initial is given by the following: We took the infinite norm of the matrix because elements in the Jacobian matrix were arranged in rows. So, the Lagrange multiplier at kth iteration is shown as follows: where is an empirical constant, n is a non-negative constant and k is the number of iterations.

Model constraints
For 3D inversion, the model smoothness matrix W c is composed of the following three spatial components (Li & Oldenburg 1996): where ( x , y , z ) are weighting coefficients in three directions, (W x c , W y c , W z c ) are the first-order difference operators, for the model x = 1, … , N x , y = 1, z = 1 and the first-order forward difference coefficient matrices are as follows: For the case of y = 2, … , N y , z = 2, … , N z , the coefficient matrices (W x c , W y c , W z c ) were consistent with equation (19). The smoothness matrix of the full model W c could be easily constructed using a triple-loop algorithm. This constraint is known as the flattest model (FM) constraint.
We applied the MS constraint to the last term of the objective function. The MS stabilising function was proposed by Last and Kubik to find the model with the smallest volume change compared with the reference model (Last & Kubik 1983). Here, W m is the diagonal matrix where is the stability factor when (m i − m ref i ) → 0. Such a diagonal weighting matrix has been used extensively in geophysics (Wang 1999). We applied an L-curve method that was similar to the optimisation of regularisation coefficients to find the optimal parameter (Portniaguine & Zhdanov 1999).

Synthetic data example
Next, we applied the inversion scheme to a 3D crosswell EM problem-that is, a synthetic crosswell EM survey over an area of 400 × 400 × 400 m. Two 20 Ω ⋅ m -1 anomalous objects were embedded in a homogeneous half-space of 2 Ω ⋅ m -1 (figure 2). A cuboid of size 120 × 40 × 40 m with its top at a depth of 140 m and another identical cuboid was placed 40 m below it. The two 200-m-deep source wells and two 200-m-deep receiver wells were placed around anomalies. We deployed 21 vertical magnetic dipole transmitters at each well with 10 m intervals, and the receivers for each of the receiver wells had the same layout, which yielded 1764 source-receiver combinations. We inverted the vertical magnetic field component E z of 370 Hz, which contained 2% random noise.
The dotted rectangular area in figure 2 is the inverted area used in the inversion algorithm, which was smaller than the area used in the forward problem. For forward modelling, the 128 Figure 3. Inversion results obtain by applying (a-c) Ekblom norm and (d-f) L 2 norm. model was discretised into the grid of 50 × 50 × 50 cells, with a minimum cell size of 10 × 10 × 10 m and a maximum size of 200 × 200 × 200 m, whereas the inverted region contained 8000 (20 × 20 × 20) cells. It was critically important to exclude receiver points from the inverted region in the inversion of the explicitly calculated Jacobian matrix. When we used the adjoint equation method to explicitly calculate sensitivities, the receiver points had to be excluded from the inverted region. As a result, a number of singularities from the solution of adjoint equations were introduced into Jacobian matrix, which made it difficult to solve equation (13).
The initial model for inversion was a homogeneous halfspace of 2 Ω ⋅ m -1 . We set the parameters of the Ekblom norm as p = 1 and = 10 −5 . We chose = 0.1, n = 2 and = 5000 from our experiences. Several trial-and-error tests were required to determine these values, especially beta. For these model constraints, we set x = 1, y = 0.35, z = 0.35 and = 0.15. The misfit is calculated by the mean percentage error: Figure 3 shows the models obtained by the inversion with the Ekblom norm and L 2 norm. The Ekblom-norm inversion achieved relatively blocky anomalies with the help of the FM and MS constraints, whereas the L 2 norm inversion smoothed out the boundaries of anomalies under the same model constraints. The reason for these results was that the Ekblom measure (p = 1 and = 10 −5 ) put less weight on large m i and more weight on small m i , which led to an image with relatively sharp boundaries (Zhang et al. 2000;Sun & Li 2014). Although both inversions could correctly determine the position of the anomalies, they somewhat underestimated the resistivity of the objects. Note that we set the upper limit of the resistivity in figure 3 to 10 Ω ⋅ m -1 , that 129 Figure 5. Four source-receiver combinations (first column), amplitudes of the observed data (second column), mean percentage error of the predicted data and the observed data before and after Ekblom-norm inversion, respectively (third and last column).
is, a half of the true model. The maximum resistivity value of the anomalous bodies from the Ekblom-norm inversion was 13.2 Ω ⋅ m -1 , whereas that of the L 2 -norm inversion was 11.7 Ω ⋅ m -1 . Figure 4 shows the convergence rates of the inversions with two measures. Both inversions converged to a 2% data noise floor within ten iterations. The misfit of the Ekblom-norm inversion at the first iteration showed few decreases compared with the initial misfit, which suggested that 5000 was not the optimal value of the Lagrange multiplier . For the purpose of comparison between two mea-sures, however, we had to control the variable. To fully illustrate the inversion results, we displayed amplitudes of the observed data for all observation systems and the data misfit for all source-receiver combinations before and after the Ekblom-norm inversion shown in figure 5. The first column of figure 5 shows four observation systems, and the second column shows the amplitudes of the observed data. The third and last columns in figure 5 show data misfit at the initial iteration and the sixth iteration. As shown in the last column of figure 5, misfits of almost all source-receiver combinations reached the noise level of 2%. Figure 6 shows the Figure 6. Histogram of normalised probability for data misfit between predicted and observed data before (blue) and after (red) inversion. Dashed line denotes noise level of 2%. data misfit histogram before and after the inversion. After ten iterations, 99.3% of the predicted-data misfit of the Ekblomnorm inversion was less than the noise level of 2% and 96.7% of the misfit was less than 1%. The predicted-data graphs of the L 2 norm were similar to figures 5 and 6, which are not shown because of the length of the paper.
To investigate the influence of the data misfit term and model term with different measures on the inversion, we tested the inversion of two mixed norms. The objective functions of the two mixed-norm inversions were Φ = Φ Ekblom , respectively. This was easy to implement. For the former, we set p = 2 in r c and r m , whereas for the latter, we set p = 2 in r d . Figure 7 shows the inversion results of the two mixed norms. The inversion results of the Ekblom norm and the L 2 norm also are shown for reference. The second and third rows of figure 7 show that applying the Ekblom norm to the model term alone could not reconstruct blocky anomalies, whereas the inversion results of applying the Ekblom norm to the data misfit term alone were similar to the L 2 -norm inversion results. Figure 8 shows the convergence curves of the mixed norm, Ekblom norm and L 2 norm. The convergence rates of applying the Ekblom norm to the data-fitting term (blue and yellow lines) were less than the case of L 2 norm (purple line). The convergence rate of the Ekblom norm for the model term alone (orange line) was faster than the L 2 -norm inversion. These two experiments showed that the structural characteristics of the anomalies reconstructed by the inversion completely depended on the model constraints, whereas the application of the Ekblom norm to the data misfit term changed the dataweighting matrix, thus affecting the choice of the Lagrange multiplier.

Computational performance analysis
We ran all inversion programs on a personal computer equipped with an Intel Core 4.0 GHz CPU. As a control group of iterative solvers, the forward modelling used the generalised product-type biconjugate gradient (GPBiCG) solver (Dehghan & Mohammadi-Arani 2016). We employed a preconditioner of the incomplete Cholesky decomposition and a static correction (Smith 1996;Sasaki & Meju 2009). Meanwhile, we solved the inversion equations system (equation 13) using the modified Gram-Schmidt method (Bjõrck 1996).
In the direct solver examples, one iteration of the inversion took about 10 minutes. It was composed of three parts: forward calculation of observed data for updated model, forward calculation of sensitivity matrix and solution of the least-squares system (equation 13). The maximum memory usage of the inversion was about 8569 MB. More detailed information is shown in Table 1.

Conclusions
We presented a general-measure inversion scheme for 3D crosswell EM data based on a direct solver. In the forward modelling, we applied the SGFD method and the Intel MKL PARDISO direct solver to discretise and solve the vector EM Helmholtz equation, respectively. In the inversion, we integrated the Ekblom norm into the objective function. The FM constraint and MS constraint were imposed simultaneously on the model. By testing the synthetic data, we obtained the following understanding of the 3D resistivity inversion scheme proposed in the study: (1) The Ekblom norm used in the objective function served as the weight matrices of data misfit and model structure, and appropriate norm parameters clarified the boundary of the inverted target.
(2) The direct solver not only accelerated the forward modelling and sensitivity calculation, but also had a very good acceleration effect on solving the inversion leastsquares system. (3) The model item in the objective function determined the structural characteristics of the reconstructed model, whereas the data item determined the position and conductance (i.e. volume multiplied by conductivity) of the abnormal body.