Abstract
We presented a 2.5D inversion algorithm with topography for frequency-domain airborne electromagnetic data. The forward modeling is based on edge finite element method and uses the irregular hexahedron to adapt the topography. The electric and magnetic fields are split into primary (background) and secondary (scattered) field to eliminate the source singularity. For the multisources of frequency-domain airborne electromagnetic method, we use the large-scale sparse matrix parallel shared memory direct solver PARDISO to solve the linear system of equations efficiently. The inversion algorithm is based on Gauss-Newton method, which has the efficient convergence rate. The Jacobian matrix is calculated by “adjoint forward modelling” efficiently. The synthetic inversion examples indicated that our proposed method is correct and effective. Furthermore, ignoring the topography effect can lead to incorrect results and interpretations.
1. Introduction
As an important geophysical prospecting method, the frequency-domain airborne electromagnetics is widely used in fields of mineral survey, geological mapping, groundwater resource exploration, environmental monitoring, and so forth. Since the prospecting is often made in mountainous areas with large undulations which can cause severe impact upon the aeroelectromagnetic response, ignorance of the landform will lead to errors in aeroelectromagnetic data interpretation. By boundary element method, Liu and Becker simulated in 1992 the topography’s impact upon the frequency-domain airborne electromagnetics [1]. Newman and Alumbaugh made 3D finite difference calculation and simulated the impact from simple topography upon the observed data results in 1995 [2]. Sasaki and Nakazato simulated the topography’s impact upon the observed result by finite difference calculation and made an attempt to find a way to eliminate this impact in 2003, concluding that no way but making inversion with the topography taken into account can eliminate the impact from the topography [3]. In a word, implementing inversion with the topography taken into account is imperative in airborne electromagnetics.
The conductivity imaging technique and 1D layered medium inversions are still the main interpretation means in frequency-domain airborne electromagnetics [4–6]. Due to the bad reconstruction capacity of 1D inversion for the high-dimensional models, lots of scholars made research on the inversion for high-dimensional models. Wilson and others realized 2.5D inversion in 2006 and applied it in theoretical and measured data [7]. Cox and others made 3D inversion based on footprint technique to inverse the measured data in 2010 [8]. Liu and Yin developed a finite difference 3D inversion program in 2013, compared the NLCG with LBFGS method, and drew the conclusion that the latter one is more applicable for the inversion in 3D frequency-domain airborne electromagnetic [9]. Yi and Sasaki put forward a joint inversion program making use of aeroelectromagnetic data and ground DC data in 2015 to restore the real model better and promote the inversion resolution ratio [10]. Nevertheless, all the above are conducted under the premise of a level surface, with no consideration of the topographical impact; thus in fact the research about inversion algorithm for the airborne electromagnetics is still quite limited.
This paper, based on vector finite element method, adopts PARDISO to solve the linear system of equations, subdivides the model in irregular hexahedral mesh, makes 3D forward modeling in topography, and employs Gauss-Newton method to carry out 2D inversion of frequency-domain airborne electromagnetics with topography, providing an effective mean to calculate, process, and interpret the aeroelectromagnetic data for mountainous areas.
2. Forward Modeling
To eliminate the source singularity, assume as the time-harmonic factor and divide the electromagnetic field into primary field (background) and secondary field (scattering). The double curl electric field equation based on the secondary field can be as follows [11]:where is the secondary electric field, is the angular frequency, is the magnetic conductivity in vacuum, is the electric conductivity, is the background field, and is the difference between the total conductivity and the background conductivity, expressed as , with as the background electric conductivity. The primary field under the uniform whole space medium or under the layered medium is calculated by analytic solution, while the secondary field (scattered field) is calculated by vector finite element method. The double curl Maxwell equation satisfies the following boundary condition:
Implement the Galerkin weighted residual method for formula (1); then
In the equation , represents a discrete element and is the number of discrete elements; thus formula (3) can be expressed in discrete form as below:where and are element stiffness matrix, expressed as
is the vector basis function [12]. 27-point Gauss integral is used for the calculations in formula (5).
The element global coordinates mapping on local coordinates is shown in Figure 1. Allocate the element stiffness matrix to the overall stiffness matrix according to certain standard, and then the large-scale linear system of equations can be drawn:
is sparse complex symmetric matrix, is the field value, and is the source item.
To guarantee a unique solution, the Dirichlet boundary condition is employed, deeming that the secondary abnormal field at the boundary which is far enough from the anomalous body has attenuated nearly to zero. In other words, at the boundary,
Once the electric field is calculated, the magnetic field can be solved by Faraday’s law:
2.1. Precision Validation of Forward Modeling
To verify the correctness of this procedure, a 2D trapezoid mountain model the same as that used by Sasaki and Nakazato [3] is adopted with the medium electric conductivity of 100 ohm-m, as shown in Figure 2. The top and bottom width of the trapezoid are, respectively, 20 m and 220 m, the height of the trapezoid is 50 m, and the gradient is 26.5°. Subdivide the modeling area in , and direction into grids, with the one in the middle sized and the largest one at the boundary sized . HCP device is employed for calculation, with the transmitter-receiver coil pair keeping 30 m above the ground in 10 m transmitter-receiver spacing in three frequencies of 1 kHz, 4 kHz, and 16 kHz. The results in Figure 3 highly coincide with Sasaki’s finite difference method results, which illustrated that the forward modeling program in this paper is highly accurate.
3. Inversion Algorithm
3.1. Objective Function Definition
The objective function can be written aswhere denotes the Euclidean norm, is the data vector obtained from forward modeling, is the measured data vector, is the model parameter vector, is the reference model vector or the prior information model vector, is the model smoothness matrix, and is the regularization parameter that balances the effect of data misfit and model regularization during minimization.
The objective function defined by formula (9) is divided into data fitting item and regularization item. The former guarantees the inversed model to fit the measured data, while the latter controls the inversion result to draw close to the known prior model. There are normally three cases for the regularization item: (1) smallest model, that is, smallest ; (2) flattest model, that is, smallest first-order difference operator for ; (3) most smooth model, that is, smallest second-order difference operator for . Due to the inversion algorithm in this paper, the third type is adopted, that is, with as the second-order difference operator, and has the following form:
Assuming as the initial model, substitute it into formula (9) and linearize it by Taylor expansion, ignoring the higher order terms to reachwhere is . Derive formula (11) and order it equal to zero so as to obtain the model updated iterative formula for Gauss-Newton method:where is the approximate Hessian matrix, is the objective function gradient, and is Jacobian matrix or sensitivity matrix, with its elements expressed as
3.2. Jacobian Matrix Calculation
From formula (12), it is known that the key to the model update iteration is to solve the linear system of equations, which in turn is decided by Jacobian matrix or sensitivity matrix ; thus the calculation for Jacobian matrix is crucial in the inversion algorithm, and its calculation efficiency is decisive to some extent for the inversion efficiency. The calculation methods for Jacobian matrix include analytic method, difference method, and adjoint equation solution method [13–15]. Since the analytic method cannot solve the matrix in high-dimensional inversion and the difference method, though simple, is inefficient, here we choose to solve the adjoint equation for the solution of Jacobian matrix. The calculation formula for the magnetic field is expressed as
The magnetic field’s model parameter derivative is expressed as
The magnetic field’s model parameter derivative can be converted into the electric field’s model parameter derivative. The forward modeling can yield a large-scale complex number linear system of equations:
Derive simultaneously from both sides of formula (15) to work out
The electric field’s model parameter derivative can be worked out by solving formula (17). Formula (17) has nearly the same coefficient matrix as the forward modeling linear system of equations, with difference at the right-hand item. The solving process is similar to the forward modeling process, named as “accompanying forward modeling.” According to interchange theorem, to obtain the model parameter sensitivity of the magnetic field at the receiving point, all that is needed is to put the unit source at the receiving point and conduct “accompanying forward modeling” once. It is equivalent to the simple weighted summation with the model parameter-correlated element field value when putting the unit source at the exact receiving point, which can be worked out from the right-hand item of formula (17). “Accompanying forward modeling” is needed to be carried out just once, which refers to solving all the measured points and frequency points, and then the Jacobian matrix can be obtained.
3.3. Solution for Model Updated Equations
As known from formula (12), to obtain the model parameter update vector quantity, the linear system of equations must be solved first, either solved directly or solved by iterative method. needs be worked out for the solution of formula (12). If selecting direct solution, Jacobian matrix is needed to be stored. Here the preconditioned conjugate gradient iteration method is adopted, only with the multiplication of the matrix and the vector required to be worked out. As for whether to store the Jacobian matrix in the iteration method, if storing, though unnecessary to solve in each time of iteration, the memory demand will be large, and Jacobian matrix must be stored entirely due to its denseness, while if not storing, the memory insufficiency can be avoided, but for each time of iteration, Jacobian matrix must be calculated once again, plus the forward modeling and “accompanying forward modeling” being conducted at each time.
3.3.1. Step Length Selection
After obtaining the model parameter update vector by working out formula (12), the next iterative model is updated by where is the step length with the range as . The suitable step length is defined by the following steepest descent formula [16]:where is a constant that in practice takes a very small value, for example, . Satisfying (19) guarantees that the our new model update adequately reduces the objective function. When performing the model update in (18), the initial value of is set to one and the new model is generated. The inequality in (19) was evaluated, and if the condition is satisfied, the new model is acceptable., However, if the condition is not satisfied, is reduced by one-half, and (19) was reevaluated. This procedure is repeated until we find a new model that meets the criterion of (19).
The iterative process is repeated until the objective function and/or the data misfit have reduced to a specified level. Alternatively, the inversion can be stopped when the magnitude of the objective functional gradient falls below a given tolerance, for example, , which indicates that the objective function has flattened out and that no further model improvements can be made. Moreover, the inversion can be stopped when the iteration number exceeds the maximum iteration number. Our experience suggests that it is usually not more than ten iterations. Otherwise, stopping the inversion once the data misfit reaches the noise level of the measure data usually provides reasonable results.
3.3.2. Regularization Factor Selection
The regularization factor can balance the data misfit and model constraints. There are lots of schemes that have been recommended in literature and used in practice for the selection of [17, 18]. To deal with , there are three most common approaches: (1) one fixed value in the whole inversion process; (2) decreasing value gradually in each time of iteration; (3) choosing the most appropriate value by L curve in each time of iteration. Normally, scheme (3) has the best effect but requires additional calculation work to seek the appropriate value each time. In this paper, schemes (1) and (2) are adopted in the inversion and are found capable of yielding satisfactory results as well. Furthermore, for a given problem one must often determine through trial and error. When using trial and error, one should seek the largest value which still permits that data to be fit down to the noise level.
3.4. Inversion Process
The inversion procedure is shown in Figure 4:(1)Set the number of iterations as , the accuracy of fitting, and the largest number of iterations, and input the initial model and inversion data.(2)Conduct the forward modeling and solve the forward modeling equation to work out the secondary magnetic field and .(3)Calculate the fitting error, exit the program if the set accuracy or the largest number of iterations is achieved, and continue the program if not.(4)Calculate Jacobian matrix by forward simulation to obtain the updated step length.(5)Update the model parameter: .
4. Synthetic Inversion Examples
The synthetic data sets from two 2D models were used to test the inversion algorithm. The first is a flat-earth model. The second example represents a 2D feature with topography. The data sets are generated for frequency-domain airborne electromagnetic with HCP or VCX using the vector finite element method.
4.1. Inversion in Flat-Earth Model
Figure 5 shows the model used to test the inversion method for flat-earth model. There are four anomalies bodies with 10 ohm-m resistivity buried in 1000 ohm-m half-space. HCP and VCX devices are employed under the frequencies of 1000 Hz, 2700 Hz, 7400 Hz, and 20000 Hz. The coil pair is keeping 20 m above the ground with the transmitter-receiver distance of 8 m, and there are totally 51 dots with spacing of 10 m. The model domain is subdivided into grids, with the one in the middle sized and the largest one at the boundary sized . The forward modeling results for HCP and VCX are shown in Figures 6 and 7, respectively, with the value normalized by the primary field in ppm.
(a) Real component
(b) Imaginary component
(a) Real component
(b) Imaginary component
The inversion model is subdivided into pieces, with the 520 electric resistance unknown, while the initial model is even half-spatial with resistance of 300 ohm-m. By 10 times of iterations, HCP and VCX are converged to 1.5% and 0.7% of the origin, respectively. Figure 8 shows the RMS variation curve along with the increasing iterations for HCP and VCX device. Their inversion results are shown in Figures 9 and 10, respectively, from which it is observed that both HCP and VCX have satisfactorily reconstructed the real model. Comparatively, HCP has better lateral resolution while VCX acts better in vertical resolution.
(a) HCP device
(b) VCX device
4.2. Inversion in the Model with Topography
The following two models are designed to analyze and check the inversion effect in undulating surface condition.
4.2.1. Low Resistance Composite Pattern
The model is an even half-spatial mountain type as shown in Figure 11, with two low resistance anomalous bodies buried therein. The electric resistance of the anomalous bodies and of the background is, respectively, 10 ohm-m and 300 ohm-m. The model top is 20 m wide and 50 m from the bottom, and the gradient is 26.6°. Here we change the coordinate of -axis in the grids to simulate the undulating surface. HCP device is employed under the frequencies of 1000 Hz, 2700 Hz, 7400 Hz, and 20000 Hz. The coil is 20 m above the ground with transmit-receive distance of 8 m, and there are totally 51 dots in spacing of 10 m. The model domain is subdivided into grids, with the one in the middle sized and the largest one at the boundary sized . The forward modeling results for HCP are shown in Figure 12, with the value normalized by the primary field in ppm.
(a) Real component
(b) Imaginary component
In the inversion, the forward modeling uses a grid of cells. The inversion model region is subdivided into pieces, with the 520 unknown electric resistivity, and the starting model is half-space with resistance of 300 ohm-m. After 10 iterations, the HCP RMS is reduced to 0.5% of the first iteration. The RMS variation with the increasing of iteration is shown in Figure 13. The inversion result is shown in Figure 14, from which it is observed that HCP has satisfactorily reconstructed the real model.
To illustrate the effectiveness of the inversion algorithm with the landform being taken into account, the inversion results for different landforms are shown in Figure 15. It is observed that ignorance of the landform will cause undesirable inversion result; for example, the low resistance anomalous body, though being reflected, is far different from the real model, especially for the lateral stretch. Besides, much falseness and abnormities arise in the inversion profile, and particularly a high resistance anomalous body arises right below the mountain due to the topographic impact.
4.2.2. High-Low Resistance Composite Pattern
The high-low resistance composite pattern is shown in Figure 16, with the anomalous bodies in 10 ohm-m low and 3000 ohm-m high resistance. The other parameters are the same as those in low resistance composite pattern. The forward modeling result is shown in Figure 17.
(a) Real component
(b) Imaginary component
The initial model for inversion is even half-spatial with resistance of 300 ohm-m. By 10 times of iterations, the fitting difference reduces to 2.2% of the first iteration misfit, referring to Figure 18 for its variation along with the iterations. From Figure 19 showing inversion result for high-low resistance composite pattern (with the landform being taken into account), it is observed that the inversion result has satisfactorily reconstructed the real geoelectric model with no redundant abnormity, whereas the high resistance anomalous body, though being located correctly, is observed with its resistance value as just 600 ohm-m which is far different from the real value of 3000 ohm-m. Oldenburg (2012) and others have asserted that this phenomenon is normal [19], because the pure secondary abnormal field collected by frequency-domain airborne electromagnetics is insensitive to high resistance body, and, consequently, the inductive source magnetic field data is essentially difficult to inverse the high resistance body. In Figure 20 showing inversion result for high-low resistance composite pattern (without taking into account the landform), it is observed that the inversion result, though capable of reflecting the outline of high or low resistance anomalous bodies, is far different from the real model, particularly with large lateral stretch arising for the low resistance anomalous body. Moreover, the same as the low resistance composite pattern, much falseness and abnormities arise, and particularly a high resistance anomalous body arises right below the mountain.
5. Conclusion
We have developed the 2.5D inversion algorithm for frequency-domain airborne electromagnetic data and applied it to synthetic data, such as flat-earth model and the model with topography. Commencing from analysis for the inversion effect of HCP and VCX device at horizontal landform, the paper makes research for multiple low resistance composite pattern of anomalous bodies and indicates that HCP has better lateral resolution while VCX acts better in vertical resolution. Therefore, data for various measuring devices shall be made use of to promote the inversion efficiency. Afterwards, through research for the mountain low resistance composite pattern and high-low resistance composite pattern, it is concluded that inversion with the topography being taken into account can better reconstruct the real model; otherwise much falseness and abnormities will arise, and particularly a high resistance anomalous body will arise right below the mountain. Thus to obtain an accurate geoelectric model, the impact from the topography must be taken into account in the measured data processing.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The authors thank the financial support of the Doctoral Fund Project of the Ministry of Education (no. 20130061110060, class tutors), the National Natural Science Foundation of China (no. 41504083), and National Basic Research Program of China (973 Program) (no. 2013CB429805).