Research Progress on Geophysical Three-Dimensional Electromagnetic Inversion Methods

1. Current Research Status

3D inversion of geophysical electromagnetic data is characterized by nonlinearity, underdeterminacy, and large scale, making it a very difficult optimization problem to solve. Over the past two decades, with the continuous maturation of numerical methods and the rapid development of computing equipment, geophysical electromagnetic 3D inversion technology has made great progress. Several 3D electromagnetic inversion codes based on regularization have been successfully developed and widely applied in production and scientific research. The development of electromagnetic data 3D inversion mainly revolves around the following aspects:

(1) High-precision 3D forward modeling technology

The rapid development of electromagnetic data 3D inversion is mainly due to the continuous progress of forward modeling technology. The earliest electromagnetic 3D forward modeling was based on the integral equation method, but due to the simplicity of the computational model, it could not solve complex geological problems. In the 1990s, finite difference technology based on staggered grids was introduced into 3D electromagnetic simulation, solving the forward modeling problem of complex models. Forward modeling acceleration technology based on current density divergence correction further solved the problem of slow convergence speed in low-frequency geomagnetic forward modeling. Finite-difference-based three-dimensional magnetotelluric inversion developed in the 1990s, with the algorithm developed by Newman and Alumbaugh (2000) based on finite-difference and nonlinear conjugate gradient methods being a representative example. Subsequently, the ModEM software developed by Egbert and Kelbert (2012) further propelled this technology toward practical application. ModEM has a well-modular structure and high computational efficiency, and after being open-sourced, it quickly became a standardized three-dimensional magnetotelluric inversion tool.

However, traditional 3D forward modeling algorithms based on regular hexahedral meshes cannot meet the needs of refined geophysical exploration. To accurately fit undulating terrain and complex geological structures, 3D electromagnetic forward modeling based on unstructured finite element methods has gradually developed and become a research hotspot in the electromagnetic field. Finite element methods based on unstructured element subdivision are mainly divided into two categories: deformable hexahedral and unstructured tetrahedral. The deformable hexahedral vector finite element method is simpler to implement than the unstructured tetrahedral finite element method, but it requires handling the hanging point problem during local refinement. Rücker et al. (2006) first applied the finite element method based on unstructured tetrahedral meshes to DC electrical method 3D forward modeling simulation, demonstrating its technical advantages in solving the problem of electrical method forward modeling of undulating terrain. Schwarzbach and Haber (2013) first realized 3D inversion of marine controlled-source electromagnetic methods based on unstructured tetrahedral meshes and explored the application effects of different regularization methods. Subsequently, this technology was successfully applied to magnetotelluric, marine controlled-source electromagnetic, and time-domain electromagnetic 3D inversion. Unstructured tetrahedral elements offer greater flexibility than deformable hexahedral elements, but their implementation is more complex, and the mesh size is difficult to control. While various adaptive mesh optimization strategies exist, they all only include mesh refinement algorithms and cannot remove redundant elements to generate reasonably spaced unstructured tetrahedral meshes. Besides deformable and unstructured tetrahedral elements, the finite volume method based on octree mesh partitioning has also seen varying degrees of application in electromagnetic 3D forward and inverse modeling.

Compared to differential equation methods, 3D forward and inverse modeling techniques based on integral equations are less commonly used in the electromagnetic field. The electromagnetic research team of M. Zhadnov at the University of Utah proposed an integral equation method based on quasi-linear and quasi-analytical approximations, which can quickly solve large-scale forward and inverse problems. Other newly developed forward modeling techniques, such as the spectral element method based on high-order polynomial interpolation, meshless methods with flexible model partitioning, and hybrid methods, have made some progress, but corresponding 3D inverse algorithms have not yet been developed.

(2) Efficient Equation Solving Techniques

Traditional three-dimensional electromagnetic forward modeling equations are mainly solved based on the Krylov subspace iteration method, such as the conjugate gradient method (CG), the quasi-minimum residual method (QMR), and the stable biconjugate gradient method (BICG-STAB). These methods require a small condition number for the forward modeling equations and necessitate preprocessors and divergence correction techniques to accelerate the solution process. In recent years, direct solution techniques have strongly promoted the development of electromagnetic three-dimensional simulation technology, with mainstream tools including MUMPS, Paradiso, and SuperLU. The advantages of applying direct solution techniques to three-dimensional electromagnetic problems are: 1. This technique can overcome the influence of an excessively large condition number in the equations, achieving stable solutions; 2. In forward modeling equations for multi-source and multi-time-channel electromagnetic systems (such as mobile platform aerospace and marine electromagnetic systems), the left-hand forward modeling matrix is ​​the same, and changes in the source and time channel only affect the right-hand side terms. Therefore, acceleration can be achieved by saving the matrix decomposition results and replacing the right-hand side terms for back substitution. Since the condition number of algebraic equations formed by the finite element method based on unstructured tetrahedral meshes is very large, direct solution techniques are the mainstream technique for solving this type of equation. The drawback of direct solution techniques is their high memory consumption, making them unsuitable for solving 3D problems on ultra-large-scale models. They require dimensionality reduction to achieve large-scale computation. Iterative solutions, on the other hand, require less memory. Therefore, iterative solutions with efficient preprocessing techniques are a research hotspot and important development direction for future electromagnetic large-scale 3D forward and inverse modeling.

Besides improving equation-solving techniques to increase 3D inversion efficiency, dividing the large-scale model into several smaller-scale models is also an effective speed-up method. A typical application is airborne electromagnetics. Due to the compactness (small impact range) of airborne electromagnetic systems, a separate sub-model can be applied to each data point for electromagnetic response and sensitivity calculations. Cox et al. (2010) proposed a fast airborne electromagnetic 3D inversion technique based on "moving footprint." Since the integral equation method does not require edge expansion, the frequency domain airborne electromagnetic calculation grid for a single measurement point requires only 8×8×10 elements, with an average inversion time of 1.75 seconds per measurement point, making large-scale 3D inversion possible. Furthermore, recently proposed local mesh forward modeling algorithms have solved the problem of fast time-domain airborne electromagnetic 3D inversion. The computational complexity of this type of model decomposition method increases linearly with the number of emission sources, and it is expected to become the mainstream acceleration technology for three-dimensional electromagnetic inversion.

(3) Regularization Methods

Tikhonov regularization is the theoretical foundation of modern geophysical inversion. Its basic idea is to impose prior constraints on the parameters to be solved in the inversion objective function, thereby minimizing the objective function under the constraints. Regularization techniques have solved the ill-posedness problem in inversion to a certain extent, reducing the influence of noise on ill-posedness and ensuring that the solution does not overfit. The key to regularized inversion lies in determining the optimal regularization factor l. The L-curve method is a reliable method for selecting l, but it requires multiple forward calculations, so it is not widely used in three-dimensional inversion. The adaptive descent method has comparable application results to the L-curve method and has become the main method for selecting the regularization factor in three-dimensional electromagnetic inversion. fm has a significant impact on the structural morphology of the inversion results; different definitions produce inversion results with different characteristics. Currently, in mainstream inversion algorithms, fm is defined as the L2 norm of the first or second difference of the model parameter space. Although this definition allows the objective function to maintain good convex function characteristics, resulting in stable inversion and fast convergence, the obtained inversion model is relatively smooth and has poor ability to identify steep boundaries. Focused inversion and L1-norm-based regularized inversion can effectively improve the boundary resolution of inversion by approximating the L0 norm. Another method to improve inversion resolution by utilizing the sparsity of the solution is sparse regularized inversion based on multi-scale analysis, such as wavelet transform and curvelet transform. This type of method transforms the inversion model to the sparse domain and then inverts the coefficients in the sparse domain, finally obtaining the spatial domain model parameters through inverse transform. Since the L1 norm is used to construct the regularization term in the inversion objective function, the sparsity of the coefficient solution can be well controlled during the inversion process, realizing multi-scale inversion from coarse to fine, and balancing the stability and resolution of the inversion solution. Liu et al. (2018) realized frequency domain aero-electromagnetic three-dimensional sparse regularized inversion based on wavelet transform and compared the differences between this method and traditional methods in terms of computational speed and resolution. At present, multi-scale image processing tools are developing rapidly, such as curvelet transform, shear wave, dictionary learning and deep learning neural networks, which can realize feature extraction of images at different scales, providing development space for related inversion methods.

(4) Optimization Methods

The main optimization methods used in 3D electromagnetic data inversion include the nonlinear conjugate gradient method (NLCG), the Gauss-Newton method (GN), and the quasi-Newton method (QN). Advances in optimization methods have primarily stemmed from the mathematical field, but no breakthroughs have been achieved in 3D electromagnetic inversion. Egbert and Kelbert (2012) optimized the NLCG in ModEM, improving the stability and computational efficiency of the inversion. The NLCG and QN methods only require the gradient information of the objective function and do not require explicit calculation of the sensitivity matrix, thus requiring fewer forward iterations in each inversion iteration. The Gauss-Newton method, due to its second-order sensitivity information, has a significantly faster convergence speed than the NLCG and QN methods. However, the Gauss-Newton inversion requires a large number of forward calculations to solve the equations, so its overall computational speed is often lower than the other two methods. For electromagnetic exploration methods (such as frequency-domain airborne electromagnetic methods) where the time required to explicitly calculate the sensitivity matrix is ​​close to that required to calculate the gradient of the objective function, explicitly calculating the sensitivity matrix and constructing the Gauss-Newton inversion equation can save forward modeling calculations in the solution process, fully leveraging the speed advantage of the Gauss-Newton method. Quasi-Newton methods (QN) approximate the inverse of the Hessian matrix during iteration, offering higher efficiency in step-size search than nonlinear conjugate gradients and slightly faster inversion calculations, and have become the mainstream technique for large-scale 3D electromagnetic inversion. Recent major advancements in optimization theory have been the development of stochastic optimization algorithms, including stochastic L-BFGS and stochastic Gauss-Newton methods. These methods have the advantage of using less known information during optimization, enabling large-scale 3D inversions with limited computational resources; however, accuracy and computational speed remain key issues that need to be addressed.

(5) Joint Inversion

Different geophysical methods have different detection capabilities. Joint inversion can achieve complementary advantages of multiple data and improve inversion resolution. Joint inversion is divided into three types: 1. Joint inversion of different data with the same physical property, such as joint inversion of subsurface electrical structures using magnetotellurics and controlled-source electromagnetic methods, or time-domain and frequency-domain airborne electromagnetic methods. This type of method can combine the technical advantages of one method's high resolution and another method's large detection depth, significantly improving the imaging resolution of subsurface electrical structures. A typical example is the three-dimensional joint inversion of marine magnetotellurics and marine controlled-source electromagnetic methods. By using magnetotelluric data to control the regional tectonic background, and then using marine controlled-source electromagnetic methods to achieve accurate imaging of high-resistivity reservoirs. Currently, the most widely used method in industry is the two-dimensional joint inversion of marine magnetotellurics and marine controlled-source electromagnetic methods. Efficient and practical three-dimensional joint inversion methods are yet to be developed; 2. Based on the rock physical properties of the inversion parameters, establish empirical functional relationships between different physical properties of subsurface media, and then use these functional relationships to treat different subsurface media as the same medium for joint inversion. Establishing such functional relationships requires extensive rock physics experiments, and their accuracy determines the reliability of the joint inversion; 3. Structural similarity constraint inversion between different physical properties, such as joint inversion of gravity, magnetic, and electrical data and joint inversion of seismoelectric data. This type of method uses the structural similarity of different physical properties of anomalous bodies to constrain the inversion, achieving high-precision imaging. Currently, the most widely used technique is the cross-gradient method (Gallardo and Meju, 2003), and there are relatively mature two-dimensional joint inversion software for gravity, magnetic, and electrical data and seismoelectric data, as well as three-dimensional joint inversion of magnetotelluric data and seismic travel time. Yin Changchun et al. (2018) proposed a three-dimensional magnetotelluric and gravity data joint inversion method based on local correlation constraints, which achieves structural similarity constraints and inversion by setting Pearson correlation coefficients at different scales. Another type of structural similarity constraint is to first use high-precision methods to determine the main underground structural features (e.g., seismic or well logging methods to locate underground interfaces), and then achieve the inversion of other physical property parameters under structural constraints.

(6) Anisotropic Inversion

Electromagnetic methods have limited resolution, making it difficult to distinguish minute electrical structures or strata, and can usually only be studied from a macroscopic perspective. Stratigraphic bedding or isotropic structures and textures with directional arrangement can be described macroscopically using anisotropy. Anisotropy has a significant impact on electromagnetic observation data, and ignoring its influence will lead to erroneous interpretations. Anisotropic information extraction and inversion from electromagnetic data has become a current research hotspot in electromagnetic exploration. The successful application of marine controlled-source electromagnetic methods in seabed resource exploration has promoted the advancement of three-dimensional electromagnetic anisotropic inversion technology. Since seabed sedimentary environments have typical anisotropic characteristics, the ability to effectively invert the anisotropy of strata is key to the success of marine electromagnetic exploration. Newman et al. (2010) first achieved anisotropic inversion of marine controlled-source electromagnetic data. Comparison of isotropic and anisotropic three-dimensional inversion results revealed that using a simple isotropic model to invert electromagnetic data affected by anisotropy will produce false anomalies. Wang et al. (2018) integrated the technical advantages of unstructured finite element method, primary field separation, and direct solution method to achieve high-precision and rapid three-dimensional anisotropic inversion of marine controlled-source electromagnetic data. Based on flexible unstructured tetrahedral mesh generation, the algorithm can invert marine controlled-source electromagnetic data under complex undulating seabed conditions. Studies have shown that the deep Earth is anisotropic, therefore, the anisotropic effect often needs to be considered in the interpretation of magnetotelluric data from deep Earth exploration. However, there are currently no practical magnetotelluric anisotropic inversion algorithms and software in the electromagnetic field. Due to the complexity of anisotropy in nature and the strong underdeterminacy of anisotropic inversion, current related algorithms for anisotropic electromagnetic inversion are mainly for simple models such as TI or triaxial anisotropy. Three-dimensional inversion considering arbitrary anisotropic parameters is one of the urgent problems to be solved in the field of electromagnetic exploration.

2. Research Prospects

Electromagnetic 3D inversion has now achieved a leap from non-existence to existence, with various developed software possessing basic inversion functions. However, due to factors such as computational limitations, software complexity, and the difficulty of acquiring area-based 3D data, the practical application of 3D inversion remains low, and it has not yet become a routine data processing technique. Future development directions for electromagnetic 3D inversion mainly include:

Solving the problem of efficient forward modeling for large-scale models. Efficient and high-precision 3D forward modeling is the foundation and prerequisite for inversion. Current 3D inversion model cell partitioning is generally below the million level, making it difficult to finely characterize the electrical distribution of underground media. Therefore, developing efficient forward modeling simulation techniques, mesh optimization strategies, and fast solution methods suitable for the tens of millions and above is the primary task in the current electromagnetic field.

Fast 3D inversion based on artificial intelligence technology. Deep learning and neural networks, among other artificial intelligence methods, can construct nonlinear relationships between data, completing a rapid mapping from data to the model, thereby achieving rapid imaging of electromagnetic data and serving fields such as oil well extraction, water injection fracturing monitoring, underground engineering stability, and disaster monitoring. Although such imaging results exhibit varying degrees of approximation, they can be used to construct more accurate initial models, improving the convergence and accuracy of 3D electromagnetic inversion. The challenge lies in constructing a large training set and ensuring the universality of the trained network model.

Random sampling forward and inversion methods based on compressed sensing theory. Electromagnetic exploration exhibits strong volume effects, and the data shows smooth, continuous variations in time and space. Therefore, compressed sensing theory allows for the reconstruction of the overall data using fewer random sampling points. Based on this idea, random undersampling can be performed on measured electromagnetic data to complete response and sensitivity calculations, thereby reconstructing all measurement point information, reducing the number of forward steps in the inversion, and improving inversion efficiency. Recent developments in stochastic optimization methods have further enriched this research. However, designing suitable random sampling functions and methods to improve inversion accuracy and accelerate computation speed remains a key technical challenge.

3D anisotropic inversion of electromagnetic data. Anisotropy has always been a research hotspot in the electromagnetic field, but progress has been slow. Conventional algorithms can only perform inversion of three-axis anisotropic models, which is insufficient to meet the needs of actual deep Earth structural research. Therefore, there is an urgent need for research on the identification of response characteristics of arbitrary anisotropic models and three-dimensional electromagnetic inversion.

Three-dimensional inversion technology for transient electromagnetic induced polarization data. During transient electromagnetic measurements, induced polarization effects cause rapid decay and even sign changes in the response. Since the induced polarization effect is coupled to multiple time-channel electromagnetic signals, it cannot be effectively separated using mathematical tools. Therefore, conventional inversion techniques cannot achieve inversion of this type of electromagnetic data, necessitating the development of corresponding induced polarization multi-parameter inversion techniques. Furthermore, the parameters of the induced polarization model exhibit strong nonlinearity and non-uniqueness; obtaining the true induced polarization parameters through inversion is a challenging optimization problem.

Three-dimensional inversion technology based on multi-information fusion. Electromagnetic methods themselves have limited resolution. How to utilize regularization and related techniques to introduce prior information to achieve high-precision imaging of subsurface structures is one of the important future development directions of electromagnetic methods. Currently, the mainstream method is joint inversion technology, which can combine the advantages of multiple detection data to achieve accurate imaging of subsurface structures. The technical challenge of joint inversion of electromagnetic data with other geophysical data lies in how to couple the computational grids of different methods and how to implement reasonable structural similarity constraints.

Research on three-dimensional inversion of electromagnetic data is in a period of rapid development. Advances in various numerical simulation, image, and signal processing technologies have provided multiple possible approaches for three-dimensional electromagnetic inversion. It can be expected that in the near future, three-dimensional electromagnetic inversion will rapidly develop into a practical tool for interpreting electromagnetic exploration data.