Abstract
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 different 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 field test data, it is further demonstrated that the QPSO-CLS joint algorithm has high optimization efficiency without manually setting the initial model. The interpretation results are consistent with the theoretical model and drilling logging results, which proves the adaptability of the proposed algorithm.
1. 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]. The 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. Therefore, a transient electromagnetic method for karst cave detection has become an efficient method [3–6]. However, the apparent resistivity profile 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 difficult to be defined 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 artificial 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]. Through 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. The 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. The 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.
The 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]. This 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. The 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.
2. 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.
2.1. 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 and the thickness of each layer need to be reversely calculated. The total number of variables to be solved is 2N – 1. Then, this problem can be transformed into a 2N – 1 dimensional optimization problem. The 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π]. Then, 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. The updated rules of particle state are as follows, including the update of the qubit angle and the probability amplitude of qubit:(a)The incremental update of the qubit angle on particles is as follows: where is a random number with inertia weight; is a self-factor and a global factor, respectively; is a random number of (0,1).(b)The probability amplitude of qubit on particles is updated as follows: Among them, ; ; and are the number of population and the dimension of unknown variables, respectively.(3)Mutation treatment. The quantum nongate is used to mutate particles. It gives each particle a random number between (0,1). When the random number is lower than the set value , the [n/2] qubits are randomly selected, for the mutation operation using equations (2) and (3).where ; [21].
2.2. 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). Therefore, equations (2)–(4) are used as the inertia weight adjustment strategy of the QPSO algorithm.
(a)
(b)
Adjustment strategy 1: equations (2)–(4) are used to realize the nonlinear dynamic adjustment of inertia weight as follows:where represents the real-time objective function value of particles; and 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 is 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:where is the maximum mutation probability, is the maximum number of iterations, and 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 adaptive strategy of is better. Therefore, this strategy is used in this study.
2.3. The QPSO-CLS Joint Algorithm
The 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. The 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. The 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. The objective function is as follows:where is the apparent resistivity value of the n-th recording time trace, is the apparent resistivity value of t iterations of the i-th recording time trace, and is the actual resistivity value of the i-th recording time trace.
3. Numerical Simulation
The 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, 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. The 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 , and the thickness as , so the magnetic field intensity response in frequency domain iswhere represents the equivalent radius of the transmitting wireframe. The recurrence formula of is as follows:
The 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 I0 is the emission current. Fourier transforms the above formula to obtain the expression of emission current in frequency domain as follows:
Therefore, the magnetic field strength can be simplified as follows:
Let the inner integral kernel function be . Hankel transform is represented by , and the cosine transform is represented by . So, formulas (3)–(5) can be abbreviated as follows:
For the uniform half space geoelectric model, the analytical expression of can be derived as follows:where is the error function. Normalize to obtain the following function:
Calculate the inverse function of and derive the region-wide apparent resistivity based on magnetic field strength as follows:
Since the function 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 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 shows the layer thickness and resistivity of the model. The square loop of 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.
The 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.
(a)
(b)
The 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.
4. Field Test
4.1. Location Overview and Data Acquisition
The 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. The specific location is shown in Figure 7 as follows. The 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. The 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. The moderately weathered carbonaceous limestone with cryptocrystalline texture and layered structure is broken by joints and fractures with calcite veins filled in the fissures. The 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. The geological conditions are complex.
The 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. Then, 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 were used, as well as the common central point detection method was adopted. The circle in Figure 7 is located in borehole Y1, and the logging information can be compared with the inversion results.
4.2. Inversion Analysis of Field Data
Figure 8 shows the comparison of inversion results and borehole logging results. The comparative analysis shows that the uppermost layer is the clay layer, which is humid and has high water content. Thus, 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. The borehole revealed that there are karst caves and moderately weathered carbonaceous limestone intercalations within 25 m–54.5 m. There are both fully filled and semifilled karst caves (the filling material is plastic clay). Thus, this section presents the characteristics of high resistivity value. The 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. The strata below 54.5 m show the characteristics of medium and high resistance values. It is inferred that the rock structure in this 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. Thus, the effectiveness of the QPSO-CLS joint algorithm and its application value in engineering are verified.
5. Conclusion
(1)The 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)The 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. The 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.(3)The QPSO-CLS joint algorithm improves the interpretation accuracy of transient electromagnetic. It provides a new inversion interpretation idea for the transient electromagnetic detection of the bored pile, reduces the engineering cost, and enriches the detection means of karst geological conditions.The 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
The [TEM data.DTE] data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The 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.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (NSFC) funded project (51969023).