Transient Electromagnetic 1-Dimensional Inversion Based on the Quantum Particle Swarms Optimization-Smooth Constrained Least Squares Joint Algorithm and Its Application in Karst Exploration

Before the construction of the bridge bored pile in the karst area, geological conditions of the excavation area should be investigated. In order to avoid the karst caves in underground space making adverse impacts on the construction, bearing capacity, and stability of pile foundation, in this paper, we use the transient electromagnetic method to detect the karst development in the bearing layer of the pile foundation, which is dierent from the traditional karst survey method. To improve the interpretation accuracy of transient electromagnetic detection for karst caves, the quantum particle swarm optimization (QPSO) algorithm was combined with the smooth constrained least squares (CLS) algorithm, and the transient electromagnetic inversion based on the QPSO-CLS joint algorithm was generated. Better inversion results were achieved by the proposed method in this study. Based on the inversion calculation results of simulation data and eld test data, it is further demonstrated that the QPSO-CLS joint algorithm has high optimization eciency without manually setting the initial model. e interpretation results are consistent with the theoretical model and drilling logging results, which proves the adaptability of the proposed algorithm.


Introduction
Karst has been widely developed in Yunnan and Guizhou areas in China and also has been frequently distributed in northern Guangdong, western Hunan, western Hubei, and eastern Sichuan of China [1]. If the development of the underground unfavorable geological body is not found out when drilling cast-in-place pile in karst area, then safety accidents such as hole collapse, ground subsidence, buried drill, and cracking of surrounding building structure can be caused [2]. e transient electromagnetic method has the advantages of low cost, simple operation, large detection depth, strong sensitivity to water, and mud bearing karst cave and is less susceptible to external interference. erefore, a transient electromagnetic method for karst cave detection has become an e cient method [3][4][5][6]. However, the apparent resistivity pro le of transient electromagnetic is a comprehensive response of underground medium, and its interpretation has low accuracy. To obtain more accurate results, the geophysical inversion method is often used to process transient electromagnetic detection data.
In fact, geophysical inversion is a highly nonlinear problem, and a reliable initial model is di cult to be de ned for geophysical inversion. To this end, researchers introduced a fully nonlinear algorithm into geophysical inversion. Somanash has made some research achievements on a simulated annealing method, a genetic algorithm, and an arti cial neural network algorithm [7]. Li et al. proposed a nonlinear programming genetic algorithm and applied it to 1D inversion of ground transient electromagnetic and achieved good results [8].
rough the integral equation numerical simulation, Chen studied the transient electromagnetic response characteristics in the full space of a mine and obtained the best response component [9]. Sun et al. introduced the simulated annealing nonlinear global optimization algorithm into transient electromagnetic inversion calculation, took L1 norm as the objective function, and achieved good inversion results [10]. However, the above algorithm has the disadvantages of slow convergence speed and low accuracy.
Particle swarm optimization (PSO) is a crowd-based algorithm [11]. Compared with the algorithms mentioned above, PSO has strong adaptability and can be based on the global optimization algorithm. Many scholars have studied the application of the PSO algorithm in the field of geophysics. Shaw and Srivastava evaluated the adaptability of the PSO algorithm in geophysical data inversion by inverting synthetic data with noise interference retrieved on a multilayer one-dimensional model [12]. Monteiro Santos used the PSO algorithm to retrieve spontaneous potential data to detect shallow anomalies [13]. Cheng et al. proposed the PSO algorithm based on the transient electromagnetic method and the direct current method. e research shows that the proposed algorithm can obtain better results and has been successfully applied in the advanced exploration of coal mine roadways [14]. Li et al. combined the particle swarm optimization algorithm with a damped least square method and realized the inversion calculation of full space transient electromagnetic data. e results show that the combined algorithm can invert the transient electromagnetic detection data of roadways with high accuracy [15]. Similar to other global optimization algorithms, the PSO algorithm is also prone to fall into local extremum and premature convergence. To solve this problem, Li and Li fused the improved QEA (quantum-inspired evolutionary algorithm) with the PSO algorithm and proposed a fast convergence and abundant algorithm of quantum particle swarm optimization (QPSO) [16]. However, in recent years, scholars have found through research that the QPSO algorithm also has the disadvantages of premature convergence and is easy to fall into the local minimum [17]. Moreover, the QPSO algorithm is rarely applied in the field of geophysics. e smooth constrained least square method (CLS) is commonly used to solve nonlinearity fittings. Based on the Newton optimized nonlinearity least square method, it can adjust the damping coefficient and the smoothing filter to keep its forward value close to the true value and measure the gap between them by the mean square deviation RMS. When the RMS tends to be stable, the result is the final inversion result [18,19]. is algorithm is faster than the conventional least square method and takes up less memory [20].
Based on the above, in this manuscript, the quantum particle swarm optimization (QPSO) algorithm was combined with the smooth constrained least squares algorithm. e QPSO algorithm was used to carry out the preliminary iterative search, and the preliminary inversion results were taken as the initial model for the following inversion calculation in the smooth constrained least squares (CLS) algorithm. Subsequently, the reliability and accuracy of the proposed inversion algorithm were further verified by a series of numerical experiments of test functions, simulation synthesis, and field measured data.

The Inversion of the QPSO-CLS Joint Algorithm
Based on the strong global search ability of the QPSO algorithm, the application of the QPSO algorithm to loop source transient electromagnetic nonlinear inversion can help the algorithm deviate from the local optimal value and provide a more reliable initial value for the next CLS algorithm.

Basic Principle.
In this study, the quantum particle swarm optimization (QPSO) algorithm is applied to the inversion of layer resistivity and layer thickness simultaneously. Assuming that the initial inversion layer number is N, the resistivity value ρ 1 , ρ 2 , ρ 3 , ..., ρ N and the thickness of each layerh 1 , h 2 , h 3 , ..., h N−1 need to be reversely calculated. e total number of variables to be solved is 2N -1. en, this problem can be transformed into a 2N -1 dimensional optimization problem. e basic calculation principle is as follows: (1) Initialization and transformation. It is different from the ordinary PSO algorithm; the qubit phase plays the role of random initial population, which is in the range of [0, 2π]. en, the qubit can be calculated by probability amplitude. After then, by solving the solution space transformation formula, the qubit can be transformed into the corresponding value in the domain of the independent variable so that the corresponding appropriate value can be calculated. (2) Update. e updated rules of particle state are as follows, including the update of the qubit angle and the probability amplitude of qubit: (a) e incremental update of the qubit angle on particles is as follows: where ϖ is a random number with inertia weight; c 1 , c 2 is a self-factor and a global factor, respectively; r 1 , r 2 is a random number of (0,1).

Advances in Civil Engineering
(b) e probability amplitude of qubit on particles is updated as follows: Among them, i � 1, 2, ..., n; j � 1, 2, ..., d; n and d are the number of population and the dimension of unknown variables, respectively.
(3) Mutation treatment. e quantum nongate is used to mutate particles. It gives each particle a random number [ran d]i between (0,1). When the random number [ran d]i is lower than the set value K − m, the [n/2] qubits are randomly selected, for the mutation operation using equations (2) and (3).

Regulation Mechanism of Inertia Weight and Mutation
Operator. In particle swarm optimization, inertia weight ω is an important parameter. If the value of ω is increased, the global search ability will be enhanced. If the value of ω is decreased, the local search ability will be enhanced. In this study, the commonly used nonlinear method is employed to adjust the inertia weight. Besides, two methods of adjusting inertia weight are comprehensively compared, as shown in Figures 1(a) and 1(b). It is found that better results can be obtained by equations (2)-(4) rather than equations (2)-(5). erefore, equations (2)-(4) are used as the inertia weight adjustment strategy of the QPSO algorithm. Adjustment strategy 1: equations (2)-(4) are used to realize the nonlinear dynamic adjustment of inertia weight as follows: where frepresents the real-time objective function value of particles; f min and f avg are the minimum and average moderate values of all particles [21]. Adjustment strategy 2: equations (2)-(5) are used to realize the nonlinear dynamic adjustment of inertia weight as follows: where k is the current evolution algebra, r is a random number of (0, 1), and ais a decimal greater than zero, ranging from [0, 0.5].
To increase the diversity of the population, equations (2)-(6) are used to adaptively adjust the mutation operator as follows:

Advances in Civil Engineering
where K m is the maximum mutation probability, G m is the maximum number of iterations, and t is the current number of iterations. As shown in Figures 1 and 2 (a) and (b), better performance can be obtained by using an adaptive mutation operator.
To verify the ability of the adjusted QPSO algorithm, two functions are set for the optimization test as follows: ① Griewank function ② Ackley function As shown in Figures 1(a) and 1(b), the performance of the QPSO algorithm under the ϖ 1 adaptive strategy of K m is better. erefore, this strategy is used in this study.

e QPSO-CLS Joint Algorithm.
e QPSO algorithm is improved by equations (2)-(4) and equations (2)-(6) and then combined with smooth constrained least squares (CLS). Finally, a QPSO-CLS joint transient electromagnetic 1D inversion method is formed. In the initial stage of inversion, the improved QPSO algorithm is used to invert the layer resistivity and thickness of the model, and reasonable parameters such as the number of particles and the number of iterations are set according to the actual needs. After the algorithm iterates to a certain extent, the QPSO algorithm is terminated. e inversion results of the QPSO algorithm are taken as the initial model of the CLS algorithm, and the CLS algorithm is started for iterative inversion until the inversion results meet the requirements. e flow chart of the QPSO-CLS algorithm is shown in Figure 2.
Since the one-dimensional layered Earth model belongs to the mutation model, the resistivity and thickness of each layer are not continuous. Considering the calculation time and accuracy, the L1 norm of observation data and model data is selected as the objective function. e objective function is as follows: where ρ N is the apparent resistivity value of the n-th recording time trace, ρ (i) t is the apparent resistivity value of t iterations of the i-th recording time trace, and ρ (i) is the actual resistivity value of the i-th recording time trace.

Numerical Simulation
e data collected by the transient electromagnetic sounding instrument used in this paper are the induced EMF (electromotive force), which can be normalised to give a late apparent resistivity by the late apparent resistivity formula. However, late apparent resistivity is assumed to be derived approximately as time tends to infinity, so in the early stages, Set the population numberand calculation parameters c1=1.2,c2=1.1,ωmax=0.7,ωmin=0.3,and r1 and r2 are two random numbers between(0,1) Inertia weight adaptive abjustment Update particle state (update qubit depression angle and qubit probability amplitude) Adaptive adjustment of mutation operator and mutation processing IS fitting error less than set value? Or is the iteration reaching the upper limit?
IS fitting error less than set value? Or is the iteration reaching the upper limit?
Output CLS algorithm to reflect the results, end. the apparent resistivity curve is severely distorted and cannot be imaged well for shallow areas, becoming one of the distractions of shallow geological interpretation, so many scholars have proposed the concept of area-time apparent resistivity. e all-time apparent resistivity calculated directly from the magnetic field strength is a more realistic reflection of the shallow geological conditions and is closer to the true apparent resistivity definition. In this paper, the area-wide apparent resistivity is adopted as the fitting parameter for the QPSO-CLS inversion algorithm.
For 1D forward of the central loop TEM, the strength of the magnetic field perpendicular to the central loop can be solved directly. Suppose a geoelectric model with N layers, the resistivity of layer j as ρ j , and the thickness as h j , so the magnetic field intensity response in frequency domain is where a represents the equivalent radius of the transmitting wireframe. e recurrence formula of u is as follows: e transient detection instrument used in this paper adopts a square wave with a duty cycle of 1 :1, so only half a period of the waveform needs to be considered in the forward simulation, which can be regarded as a step-by-step wave as follows: where I 0 is the emission current. Fourier transforms the above formula to obtain the expression of emission current in frequency domain as follows:

Advances in Civil Engineering
erefore, the magnetic field strength can be simplified as follows: Let the inner integral kernel function be func � λ 2 /λ + u (1) . Hankel transform is represented by HT * { }, and the cosine transform is represented by CT * { }. So, formulas (3)-(5) can be abbreviated as follows: For the uniform half space geoelectric model, the analytical expression of H z can be derived as follows: where erf(x) is the error function. Normalize H z to obtain the following function: Calculate the inverse function x of Z(x) and derive the region-wide apparent resistivity based on magnetic field strength as follows: Since the function Z(x) is an implicit function, its inverse function cannot be obtained; therefore, it cannot be solved analytically. However, the numerical solution can be obtained by using the dichotomous finding method or the golden section method due to Z(x) being monotonically increasing between 0 and 1.
To verify the effectiveness of the inversion algorithm, five-layer model forward data are established in this study, and the model is shown in Figure 3. Table 1  10 m × 10 m is used as transmitting coil, the emission current is 1A, and the receiving mode of common center point receiving is adopted. In the inversion algorithm, 10 layers of the initial model are set for calculation. e descent curves of iterative computation obtained by different inversion algorithms are shown in Figures 4 and 5. Figure 6 shows the results of the QPSO-CLS joint algorithm.
e results of different inversion algorithms are compared. It is found that the QPSO-CLS joint algorithm can effectively break through the local optimal extremum and presents better inversion results. By comparing Figures 4  and 5, the inversion results of the QPSO-CLS joint algorithm are more consistent with the theoretical model, and the results have a higher reduction degree and higher fitting degree for the layer thickness and resistivity of the model.

Location Overview and Data Acquisition.
e field data were collected at pile Y1 (112.70 E, 25.14 N) of the Caojiaben interchange main line bridge of Linwu-Lianzhou (Hunan-Guangdong Boundary) expressway project in Yizhang County, Hunan Province. e specific location is shown in Figure 7 as follows. e terrain conditions of the Caojiaben junction interchange are relatively simple, with relatively gentle terrain and local steep area, and the natural slope is 20°-40°. According to the boring results illustrated in Figure 8, the strata revealed in site include Holocene deluvial and eluvial clay and the underlaying bedrock: Devonian carbonaceous limestone and sandstone. e grimy highly weathered carbonaceous limestone with texture and structure partially destroyed is closely sectioned by joints and fractures and drilling core recovered as detritus and gravel. e moderately weathered carbonaceous limestone with cryptocrystalline texture and layered structure is broken by joints and fractures with calcite veins filled in the fissures. e moderately weathered rock that is widely distributed within the site and karst is discovered locally by borings.
According to the data of borehole Y1, within 24.2-33.5 m in the underground, the surrounding rock is mainly moderately weathered carbonaceous limestone with broken joints and developed karst caves. e geological conditions are complex.
e transient electromagnetic survey line in Figure 7 is named L1, with a length of 20 m and a fundamental frequency of 25 Hz. A total of 21 measuring points were set up, with an interval of 1 m. en, 25 channels of data were collected at each point, and the emission line of 1.9 m × 1.9 m and the receiving coil of 110.16 m 2 were used, as well as the common central point detection method was adopted. e circle in Figure 7 is located in borehole Y1, and the logging information can be compared with the inversion results. Figure 8 shows the comparison of inversion results and borehole logging results. e comparative analysis shows that the uppermost layer is the clay layer, which is humid and has high water content. us, the formation resistivity of about 0 m-5 m is low. Due to the existence of strongly weathered sandstone with low strength, developed fractures, and poor integrity within 5 m-25 m, the resistivity value of this layer is higher than that of the overlying clay layer, and the local low resistivity area is probably caused by the development of fracture water. e borehole revealed that there are karst caves and moderately weathered carbonaceous limestone intercalations within 25 m-54.5 m. ere are both fully filled and semifilled karst caves (the filling material is plastic clay).

Inversion Analysis of Field Data.
us, this section presents the characteristics of high resistivity value. e local high resistivity anomaly area is probably caused by the karst cavity, and the local low resistivity anomaly area is probably caused by the high-water content plastic clay-filled area in the karst cave. e strata below 54.5 m show the characteristics of medium and high resistance values. It is inferred that the rock structure in this   Advances in Civil Engineering area is complete and dense. As the known borehole depth is 61.60 m, only the borehole logging information is compared within this range. When the detection distance is greater than 61.60 m, the transient electromagnetic detection results are obtained. It can be seen that the inversion results are generally consistent with all kinds of strata exposed by drilling, indicating that the inversion results are accurate in the strata division. us, the effectiveness of the QPSO-CLS joint algorithm and its application value in engineering are verified.

Conclusion
(1) e proposed transient electromagnetic 1D inversion algorithm based on the QPSO-CLS joint algorithm can be used to effectively depict the formation boundary of the model and efficiently inverse the information of formation thickness and formation resistivity. (2) e inversion calculation of simulated data shows that the QPSO-CLS joint algorithm has better inversion results and has a higher reduction degree for the inversion of layer thickness and resistivity value of the model, which is superior to a single algorithm. e inversion results of field test data show that the inversion results of the QPSO-CLS joint algorithm are consistent with the results of borehole logging, and the proposed QPSO-CLS joint algorithm can accurately reveal the formation information. e QPSO-CLS joint inversion algorithm proposed in this manuscript is only aimed at the 1-dimensional model which assumes that the strata under the exploration area are evenly distributed in layers. However, it has certain limitations in the 2-dimensional section imaging. In the future, relevant research and testing can be carried out so that the joint algorithm can be applied to the two-dimensional inversion to increase the accuracy of the exploration of the spatial position and shape of the cave.
Data Availability e [TEM data.DTE] data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.

Authors' Contributions
Xue Liu conducted conceptualization, methodology, validation, writing, review, and editing, as well as supervision.
Chunwei Pan conducted methodology, model analysis, validation, and data curation. Fangkun Zheng conducted methodology, formal analysis, investigation, and data curation. Ying Sun conducted data curation, writing, review, and editing. Qingsong Gou conducted writing, review, and editing.