Three-Dimensional Inversion of Semi-Airborne Transient Electromagnetic Data Based on a Particle Swarm Optimization-Gradient Descent Algorithm

: Semi-airborne transient electromagnetics (SATEM) is a geophysical survey tool known for its ability to perform three-dimensional (3D) observations and collect high-density data in large vol-umes. However, SATEM data processing is presently restricted to 3D model-driven inversion, which is not conducive to detailed surveys. This paper presents a new 3D model- and data-driven inversion algorithm using the particle swarm optimization (PSO) and gradient descent (GD) algorithms. PSO is used to suppress the multiplicity of solutions associated with inverse problems, and the GD algorithm is employed to accelerate the convergence of the inversion process. For the PSO-GD algorithm, a new model-updating equation is established and a cosine probability function is introduced as a weighting term for PSO and GD algorithms to ensure a smooth transition between the two algorithms in the iterative process. The α -trimmed ﬁlter function is used as a regularization constraint to smooth the model. The stability and reliability of the PSO-GD algorithm are veriﬁed through numerical simulations. Finally, the new algorithm is applied to the processing of SATEM measurements of the Qinshui coal mine in Jincheng, Shanxi Province, China.


Introduction
As an important geophysical survey tool widely applied in many fields (e.g., surveys of geological structures, prospecting of mineral resources, and hydrogeological surveys of coal fields) [1][2][3][4][5][6][7], semi-airborne transient electromagnetics (SATEM) has advantages, such as high efficiency, large data, and high-field adaptability [8][9][10]. Inversion is a key step in geophysical data processing. Currently, SATEM data are mostly inverted using 1D or quasi-2D inversion [2,3,[11][12][13][14]. The Earth is constantly suffering from faults, intrusions, and erosive processes that produce a subsurface morphology, which makes the geological structure generally complex. Therefore, it is necessary to perform a 3D inversion interpretation. Three-dimensional inversion of airborne electromagnetic data has attracted more and more attention. A flexible inversion framework was constructed, the lithology and drilling information was added to the inversion, and the multi-parameter inversion was achieved [15][16][17][18]. The above-mentioned inversion is restricted to deterministic algorithms. This paper analyzes the existing algorithms and proposes a hybrid algorithm.
Based on the driving strategy, inversion methods can be approximately categorized into two groups. One group consists of conventional, model-driven inversion methods: the deterministic inversion algorithm, which mainly includes the gradient descent (GD) algorithm [19,20]; the conjugate gradient algorithm [21,22]; and the quasi-Newton method [23,24]. All of these methods use the first-and second-order gradient directions or conjugate direction of the objective function as guidance to correct the model parameters, and these above algorithms have the characteristics of mature theories and high convergence rates. Deterministic inversion algorithms play an important role in 3D geophysical inversion [15][16][17][18]. However, model-driven methods have notable deficiencies in two areas. On the one hand, the inversion result heavily relies on the selection of an initial model [25]. Consequently, these methods are susceptible to falling into certain local optimum solutions near the initial model and face difficulties in converging with global optimum solutions for nonlinear, ill-conditioned geophysical inverse problems [26]. On the other hand, the inversion result lacks an evaluation of the model parameters. These limitations have called into question the accuracy of model-driven inversion methods [27].
The other group is composed of data-driven inversion methods based on statistics or the theory of swarm evolution. Stochastic algorithms include the particle swarm optimization (PSO) algorithm, the simulated annealing algorithm, the genetic algorithm, the bee colony algorithm, neural networks, and the Bayesian algorithm, et cetera [28][29][30][31][32][33][34][35][36][37][38][39]. Using multiple models, stochastic algorithms can not only obtain the better-fitting model, but also estimate the uncertainty of the retrieved solution. Compared to model-driven inversion methods, data-driven inversion methods are more capable of escaping from local optimum solutions, can provide an uncertainty evaluation of the global convergence and inversion result independent of the initial model, and are particularly suitable for solving geophysical inverse problems without a priori information. However, data-driven inversion is often more computationally expansive than model-driven inversion by several orders of magnitude [40]. In recent years, the development of computer science has allowed for an increasing application of data-driven inversion methods for geophysical inverse problems.
PSO is a typical data-driven inversion algorithm inspired by the optimization process involved in the search for food by bird or fish swarms in nature. Compared to other datadriven inversion algorithms, PSO requires fewer models and quickly converges during the optimization process; therefore, it has been extensively applied in many engineering fields. For the first time, Fernández-Álvarez et al. (2006) introduced PSO to the interpretation of geophysical inversion results. They employed this algorithm to interpret the inversion of electrical sounding data [41]. Later, many researchers began to study the application of PSO to geophysical inverse problems. Researchers have introduced PSO to many geophysical fields, such as magnetic techniques, seismology, and gravity, as well as demonstrated its versatility. Researchers have also systematically examined the relationships between the parameters (e.g., the inertia coefficient, as well as the individual and social acceleration coefficients) of the PSO algorithm and improved its convergence rate for geophysical inverse problems, thereby laying a foundation for higher-dimensional inversion [42][43][44][45]. By constraining the objective function with a k-means clustering model, Liu et al. (2018) identified particle swarms from multiple magnetic sources in the presence of complex magnetic anomalies [46]. Pace et al. (2019) successfully inverted two-dimensional magnetotelluric data using a hierarchical PSO algorithm with time-varying acceleration coefficients (HPSO-TVAC) [40]. However, PSO has yet to be applied in 3D electromagnetic (EM) inversion, which is mainly because the high complexity of the 3D models of the EM method restricts the convergence rate of PSO. In particular, the decrease in the modified step size and convergence rate of PSO near an optimum solution renders it susceptible to premature convergence. In addition, the increase in the number of parameters amplifies the problem of over stochasticity associated with PSO and consequently reduces the continuity of the model information in the inversion result, which necessitates the smoothing of the model [26,46].
Through a comprehensive consideration of the advantages and disadvantages of model-and data-driven inversion algorithms, a new 3D inversion algorithm, referred to as the PSO-GD algorithm, was developed for SATEM in this study. In this algorithm, PSO is first used to fit the inversion to ensure global convergence, and then the GD algorithm is employed to accelerate the process of fitting the inversion. A cosine probability function is used as a weighting function for the model-updating equation in the new algorithm to ensure a smooth transition from global to accelerated convergence. The α-trimmed filter function is selected as a regularization constraint for the new algorithm. The model is smoothed to suppress the problem of over stochasticity associated with the 3D PSObased inversion process. The rest of this paper is organized as follows: First, the basic principle and implementation of the PSO-GD algorithm are presented. Then, the stability and reliability of this algorithm are evaluated based on numerically simulated data. Finally, the PSO-GD algorithm is applied to the 3D inversion of real-world measurements, and conclusions and suggestions on further improving the algorithm are given.

Forward Simulation Theory
In the inversion process, forward simulation is a crucial factor affecting computational efficiency. In particular, 3D forward simulation of a complex model is generally extremely time consuming and requires high computational power. Therefore, through comprehensive consideration of two aspects of 3D forward simulation, namely, its adaptability in modeling complexity and computational efficiency, 3D forward simulation is performed in this study using the staggered mimetic finite volume approach [47,48].
We begin with Maxwell's equations of time-varying electromagnetic fields. In forward simulation theory, Maxwell's equations can be written as follows: where E and H are the intensities of magnetitic and electric fields, respectively, S(r, t) is the emitter term, σ is the electrical conductivity, ε is the permittivity, and µ is the magnetic conductivity (µ = µ 0 ). An electromagnetic field consists of primary and secondary fields: where E p and H p are the intensities of primary electrical and magnetic fields, respectively, and E s and H s are the secondary electrical and magnetic fields, respectively. Neglecting the displacement current allows for the transformation of Equations (1)- (5).
where σ p is the electrical conductivity of the half space. Transforming the volume integral to a surface integral based on Gauss's law gives the following discrete state-space equation: where n is the unit normal vector of ∂V.
Finally, Yee's method is employed to temporally discretize the equation. Specifically, H and E are discretized on the edges and surface, respectively.
The following equations give the initial and boundary conditions for temporal discretization: where H 0 is the initial static value of the magnetic field.

Theory of the PSO-GD Algorithm
In this section, the original PSO code developed by Miranda (2017) in Python is improved and turned into an HPSO-TVAC [40]. Then, a new model-updating equation is developed based on the correction terms in the PSO and GA models. The PSO and GA algorithms are integrated through a cosine probability function. In addition, the α-trimmed filter is selected as a regularization constraint to smooth the model.
Inspired by the search behavior of bird and fish swarms in nature, Kennedy and Eberhart proposed the PSO algorithm in 1995 [43]. The basic PSO algorithm calculates velocity and location using Equations (10) and (11) where γ 1 and γ 2 are random numbers in the range of (0, 1) α 1 and α 2 are the individual and social acceleration coefficients, respectively, i is the particle number, k is the current number of iterations, and P i are the individual optimum solutions, and G is the global optimum solution. Each particle of the particle swarm may be a possible solution (i.e., an electrical resistivity model) to a short-offset transient electromagnetic (SOTEM) inverse problem. The ability of PSO to escape from local optimum solutions is ascribed to the use of a swarm information sharing system and an individual fitness evaluation system as guidance to correct the model. A continuous search for a better-fitting model for the observation data in the solution space ultimately allows the particles to converge on the optimum solution, which is independent of the selected initial model [45]. When used to solve a 3D inverse problem, the convergence of the PSO algorithm requires many particles and multiple iterations. Therefore, the standard PSO version must be revised. Several variants of the PSO algorithm have been developed to improve the convergence rate [49][50][51][52][53]. The HPSO-TVAC algorithm considerably outperforms the standard PSO algorithm in terms of convergence and stability. The equation of the HPSO-TVAC algorithm is shown as follows: In the HPSO-TVAC algorithm, the initial value of α 1 is greater than that of α 2 . As the iterative process proceeds, α 1 and α 2 are linearly inverted. The high value of α 1 at the initial stage ensures a global optimization capability, while the high value of α 2 at the final stage ensures the convergence of the swarm towards the optimum solution. α 1 and α 2 are expressed using the following equations: where α k is the acceleration in the current (i.e., kth) iteration, α min In addition, the following conditions should be considered during the selection of accelerations to ensure stability [40]: By performing a first-order Taylor series expansion of the linear approximation equation, we obtain the following conventional GD inversion: where d obs is the observed data, e obs is the observation error, g is considered the response value of parameter m re f in the data-space model, δm true = m true − m re f , m true is the true geoelectric model, and m re f is an arbitrarily selected electrical resistance model. Simplification of Equation (18) gives where δd obs = d obs − g(m re f ). Then, the model correction term is Introducing the model-correction term in the GD algorithm to PSO yields the following equation: where α min 3 and α max 3 are the minimum and maximum values of the gradient descent acceleration parameter, respectively.
To ensure that the global PSO and GD algorithms take precedence at the early and late stages, respectively, a weight-based improvement strategy is further introduced: A cosine probability function, w k , which dynamically decays from 1 to 0 as the number of iterations increases, is used as a weighting term to realize a dynamic transition from the GD algorithm to the PSO algorithm. The calculation equation for w k is given below: where k and max(k) are the current and maximum numbers of iterations, respectively. In the PSO-GD algorithm, the GD algorithm guarantees the convergence speed in the early stage, and the global convergence of the PSO algorithm prevents it from falling into a local minimum.

α-Trimmed Filter Constraint
The missing a priori information in an inverse problem and the problem of over stochasticity associated with PSO result in abrupt changes in the information for the model parameters in the inversion result in the 3D space. To address this problem, a regularization constraint term is added to the objective function in a conventional deterministic inversion algorithm, with the goal of stabilizing the inversion process and producing a smooth and reasonable result. The equation of the objective function is as follows: where φ 0 is the observed signal, φ c is the calculated response, and λ is referred to as the Lagrangian multiplier or smoothing parameter. The right-hand side of Equation (25) consists of two terms: the first term evaluates the difference between the observed data and the response calculated by the forward simulation model, and the second term realizes the regularization constraint through the minimization of the roughness of the model.
In contrast with conventional model-driven algorithms, a mean filter is generally used to smooth the model during the inversion process of the PSO algorithm equipped with a regularization term [26]. In this study, the α-trimmed filter is used as a regularization constraint strategy for the PSO-GD algorithm. Prior to the start of each iteration, the individual optimum solutions P i in the particle swarm are successively filtered. First, a fixed window width n is selected based on the scale of the constraint. Then, 3D data centered at the location of the constraint are successively read with a sliding window. The data selected along the x-, y-, and z-axes are sorted in ascending order. Subsequently, a trimming parameter α is selected to trim off the same number of points from the frontand back-end segments of the sequence, with the goal of eliminating the interference from extreme values. Finally, the mean of the remaining points is used as a substitute for the value at the center of the window [54]. The whole process is shown in Figure 1. The 3D filtering equation is as follows: m(n x , n y , n z ) = While smoothing the model, the trim filter can also remove the impact of extreme values, thereby sharpening the boundary information and allowing the anomalies to converge.
To sum up, Figure 2 shows the flowchart of the PSO-GD algorithm developed in this study.

Experimental Function Testing
To evaluate the PSO-GD algorithm, the convergence of a potential function U on the x-y plane is first tested. The equation for U is as follows:  Table 1 summarizes the parametric information for U. Figure 3 shows the distribution pattern of the fields, which is characterized by two local minima in the upper half and a global minimum in the lower right corner. The GD, PSO, and PSO-GD algorithms are respectively employed to locate the optimum solution. The optimization capacity of each algorithm is subsequently evaluated based on its optimization result and path. Observation of the optimization results in Figure 3 reveals the following behavior: The PSO and PSO-GD algorithms converge on the global minimum, while the GD algorithm falls into a local minimum. Because the GD algorithm is designed to solve convex optimization problems, its optimization result for a nonconvex optimization problem primarily depends on the selected initial model. The GD algorithm is susceptible to converging on a local minimum near the initial model. In terms of the optimization path, the PSO algorithm searches for an optimum solution in a larger area, which involves a more complex path. When the PSO algorithm eventually converges on the global minimum, its search efficiency considerably decreases as the model approaches the optimum solution. For the PSO-GD algorithm, the addition of the GD information to the model-correction term, to some extent, improves its convergence efficiency near the optimum solution. Therefore, the new algorithm is more effective at executing a global search in the parametric space.

Numerical Simulations
The noise immunity and convergence of the PSO-GD algorithm are examined based on two inversion models and are subsequently compared with those of the conventional GD and PSO algorithms. Figure 4 shows the form of the forward simulation model and transceiver. The height of the sampling point is 50 m. The background electrical resistance in the forward simulation model is 50 Ω·m, while the electrical resistance in the high-resistance anomaly model is 250 Ω·m. The anomaly model has a dimension of 600 m × 600 m × 350 m, with its top surface located at a depth of 250 m. The fundamental frequency and length of the emitter are 12.5 Hz and 2 km, respectively. Samples are collected between 10 −3 and 20 ms. In addition, Gaussian random noise at 5%, 10%, 20%, and 30% levels is respectively added to the forward simulation data. Analysis of the inversion results in Figure 5 and the fitting error curves in Figure 6 reveals that gradually increasing the noise level increases the number of false anomalies in the inversion result and the fitting error. At a 5% noise level, the inversion result basically matches the true model. At a 10% or 20% noise level, the inversion result can still distinctly show the location and shape of the anomaly. However, this capability is no longer demonstrated at a 30% noise level. Therefore, in practice, reducing the noise level of the input signal as much as possible can help to obtain more information for the geological structure in the inversion result. The overall data noise level should not exceed 20%.

Complex Model Testing
In this section, the convergence of the PSO-GD algorithm is further examined based on a relatively complex model (see Figure 7), which has a background electrical resistance of 50 Ω·m and contains three anomalies: one shallow high-resistance (250 Ω·m) anomaly, one deep high-resistance (250 Ω·m) anomaly, and one deep low-resistance (10 Ω·m) anomaly. The shallow anomaly has a dimension of 600 m × 600 m × 250 m, with its top surface located at a depth of 20 m. Each of the two deep anomalies has a dimension of 1200 m × 800 m × 500 m, with its top surface located at a depth of 500 m. Gaussian random noise is added at a 20% level to the forward simulation data to generate an inversion-fitting objective. On this basis, the GD, PSO, and PSO-GD algorithms are respectively used to invert the data. A uniform half space with a resistance of 60 Ω·m is used as the initial model for the GD algorithm. Parameter α in the PSO and PSO-GD algorithms is set to 0.2-0.9. The HPSO-TVAC algorithm is used to update the PSO model. The inversion result yielded by the GD algorithm in Figure 8 can be used to satisfactorily reconstruct the shallow high-resistance structure, but fails to reflect the deep high-and low-resistance anomalies. The inversion result produced by the PSO algorithm in Figure 9 accurately shows the location of each deep anomaly. However, the boundary information for the anomalies requires further correction, and the model is not adequately smooth. In comparison, the inversion result given by the PSO-GD algorithm is smoother and closer to the original model ( Figure 10). In addition, note that due to the addition of noise and the relativity of global convergence, some false anomalous points appear in the deep area or at the edges where the constraint on the data decreases in the three inversion results. Figure 11 shows the fitting error curves and uncertainty analysis corresponding to the three algorithms. Compared to the conventional PSO algorithm, the PSO-GD algorithm has a markedly faster convergence rate at the late stage and ultimately yields a smaller fitting error. In an uncertainty analysis, the new algorithm has a higher capacity for model recovery and a more stable effect. This numerical test demonstrates the reliability of the new algorithm in the inversion of SATEM data and lays a foundation for its real-world exploration applications.

Case Study
Finally, the 3D PSO-GD inversion technique is employed to process and interpret SATEM measurements of the Qinshui coal mine in Jincheng, Shanxi Province, China. The primary objective of this survey was to investigate the goafs in the No. 3 coal seam of the Permian Shanxi Formation within the survey area to provide a reference for subsequent mining. Water-bearing goafs are a major safety hazard during coal-seam mining. During mining in this survey area of interest, we discovered that water-bearing karst collapse columns pose a threat to the mining project and the safety of the mine.
Based on the available geological and drilling data, as well as the distribution of the formations, a total of 12 SOTEM survey lines, each with a length of 2000 m and an offset from 100-600 m, were designed for the survey area. The emitter had a length of 1700 m and emitted a current with an intensity of 40 A at a fundamental frequency of 12.5 Hz. The flight altitude and speed were set to 50 m and 8 m/s, respectively. The verticalinduced electromotive force dBz/dt was measured. Figure 12 shows the arrangement of the survey lines. The PSO-GD algorithm was used to perform a 3D inversion. Figure 13 shows the inversion result, which satisfactorily reflects the spatial location of the goafs. The location of the abnormal body is clearly delineated. A vertical profile of the measured line S300 is shown in Figure 14. The known coal seam line is at an approximate height of 400 m, as shown in Figure 14. The results of the inversion made it possible to successfully delimit two abnormal zones of low resistivity and infer the location of the goaf, which will be verified by drilling in the future. Figure 15 shows the data fitting of 1400 points on the measured line S300 and a misfit curve of the 3D data.

Conclusions
A new algorithm is presented in this paper based on PSO and GD algorithms and used in 3D inversion of SATEM data. In the new PSO-GD algorithm, PSO ensures stable global convergence at the early stage of inversion, thereby suppressing the multiplicity of solutions associated with inverse problems, while the GD algorithm successfully improves the convergence rate at the late stage of inversion. The new algorithm capably balances the speed and accuracy of the inversion. The PSO and GD algorithms are effectively integrated through a cosine probability function, and the α-trimmed filter used in the PSO-GD algorithm successfully achieves model smoothing. Test results obtained based on simulation and field data show that the new algorithm is a feasible option for 3D inversion and interpretation of SATEM data.

Conflicts of Interest:
The authors declare no conflict of interest.