# Ant Colony Optimization Inversion Using the L1 Norm in Advanced Tunnel Detection

Authors: Zhao Ma; Lichao Nie; Pengfei Zhou; Zhaoyang Deng; Lei Guo
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

During the construction of the tunnel, there may be water-bearing anomalous structures such as fault fracture zone. In order to ensure the safety of the tunnel, it is necessary to carry out advanced tunnel detection. The traditional linear inversion method is highly dependent on the initial model in the tunnel resistivity inversion, which makes the inversion results falling into the local optimal optimum rather than the global one. Therefore, an inversion method for tunnel resistivity advanced detection based on ant colony algorithm is proposed in this paper. In order to improve the accuracy of tunnel advanced detection of deep anomalous bodies, an ant colony optimization (ACO) inversion is used by integrating depth weighting into the inversion function. At the same time, in view of the high efficiency and low cost of one-dimension inversion and the advantages of L1 norm in boundary characterization, a one-dimensional ant colony algorithm is adopted in this paper. In order to evaluate the performance of the algorithm, two sets of numerical simulations were carried out. Finally, the application of the actual tunnel water-bearing anomalous structure was carried out in a real example to evaluate the application effect, and it was verified by excavation exposure.


---

## Body

## 1. Introduction

Geophysical inversion is an effective method to find a model with a response similar to the actual measured values. All inversion methods are essential to determine the underground model whose response is consistent with the measured data within certain restrictions and acceptable limits [1]. For example, Liu et al. [2] performed linear inversion by using the least-square 4D resistivity inversion method, which can quickly locate and depict the change in resistivity with high accuracy. In addition, one commonly used version of the geophysical inversion is nonlinear inversion which is widely used in various fields. For example, deep leaning algorithm of the nonlinear inversion method has been applied to the resistivity and seismic fields and achieved good results [3, 4]. Song et al. [5] proposed a novel multipopulation parallel co-evolutionary differential evolution, which can achieve superior performance on the accuracy, stability, and robustness. Deng et al. [6] proposed a novel improved differential evolution algorithm and a new optimal mutation strategy based on the wavelet basis function. It can improve the search quality, accelerate convergence, and avoid falling into local optimum and stagnation. Song et al. [7] proposed a parameter optimization method based on enhanced success history adaptive differential evolution algorithm with greedy mutation strategy. Deng et al. [8] presented an improved quantum evolutionary algorithm based on the niche co-evolution strategy and enhanced particle swarm optimization, which can achieve better results on optimization performance, robustness, and stability. An enhanced improved quantum-inspired differential evolution with multistrategies was proposed by Deng et al. [9], which can improve search accuracy, accelerate convergence, and achieve the optimal solution.In particular, great progress has been made in non linear inversion methods in the field of tunnel detection. For example, Nguyen et al. [10] proposed a new hybridized global optimization method combining simulated annealing global search with the minimization of trackless Kalman filtering, which was used to solve the waveform inversion problem of advanced prediction of the underground tunnel face. The particle swarm optimization and homotopy optimization were applied to the displacement back analysis in tunnel engineering by Hu et al.[11]. Park et al. [12, 13] used the harmony search method to predict the location, permittivity ratio, and conductivity of the anomalous zone of a tunnel face by utilizing the electrical resistivity of the ground. Nie et al. [14] provide a joint inversion method based on an ant colony algorithm and least-square inversion method, which could present a better identification of the position and spatial shape of the water-bearing structure in the tunnel ahead prospecting.Ant colony optimization was introduced in the early 1990s by Marco Dorigo and colleagues [15, 16]. The ACO algorithm is inspired by the ants’ foraging behavior. The core of this behavior is the indirect communication between ants via the chemical pheromone trails, which allow them to find shortcuts between the nest and a food source through their cooperation. ACO algorithm has been widely used in optimization problems, which has advantages of positive feedback, parallelism, and robustness [17]. Since the observed data generated by different electrode intervals is related to the buried depth [18], the data weighting factor is introduced into the ACO to optimize the algorithm. In this case, the ACO algorithm based on data weighting is proposed in this paper.The common approach in the regularization optimization method is the smoothness-constrained or L2 norm method [19]. Under the condition of stable underground geological changes, the model recovered by the L2 norm method has typical smooth characteristics and will give positive results [20]. However, in regions of sharp transition in the subsurface resistivity, this method usually cannot reconstruct these sharp geological changes well [21]. In contrast, the L1 norm tends to generate piecewise continuous models and produce models separated by sharp boundaries [22]. This might be more consistent with the complex environment of the tunnel, such as the presence of fault fracture zones ahead of the tunnel face utilizing electrical resistivity.The rest of the article is arranged as follows. We first briefly introduce the ACO inversion with L1 norm or L2 norm. We also describe the procedures for ACO inversion in advanced tunnel detection. We illustrate the effectiveness of this inversion method by applying it to synthetic examples and field survey in a tunnel environment. The field data were collected from the water conveyance project from Songhua River in the middle of Jilin province.


## 2. Methodology

### 2.1. The Ant Colony Optimization

The ant colony algorithm is a method that ants find the shortest path from their nest to a food source. The core of this behavior is the indirect communication between ants through chemical pheromone trails. When the ants are looking for food, they will initially explore the area around the nest in a random manner. During exploring, ants will leave chemical pheromone on the ground. Once an ant finds a food source, it evaluates the quantity and quality of the food and brings some back to the nest. The amount of pheromones ants left on the ground during the return depends on the quantity and quality of food. Pheromone will guide other ants to find a food source [23]. The main parameters are described below.The objective function of resistivity inversion based on theLp norm ACO algorithm is as follows:(1)Φmn=∑i=1N|ρdi−ρsiρdi|p,p=1,2,where Φmn is the data constraints, N represents the number of observed data, and ρdi and ρsi represent the ith observed data and predicted data, respectively. The objective function expression is expressed as the L1 norm when p = 1, otherwise it is L2 norm. In order to better use the resistivity method to detect the water-bearing fault zone ahead of the tunnel face, the L1 norm is adopted, i.e., p = 1. When Φmn reaches the global minimum, the foraging path chosen by ants is the shortest. The model parameters corresponding to the shortest path are the inversion results.As the electrode spacing increases, the sensitivity of the observed data to the deep also increases. The data weighting factor is introduced into the inversion problem to increase the weight of effective observed data. Therefore, the objective function of the ACO algorithm based on the L1 norm is expressed as follows [18, 24, 25]:(2)Φmn=∑i=1N|wiρdi−ρsiρdi|1,wi=k0zi+z0,where wi represents the data weighting factor, k0 and z0 are constants, and zi is the distance between tunnel face and current electrode of the ith circle along z direction.The probability of ants transferring from one cell to another is determined by pheromone intensity and visibility. The transition probability is expressed as follows [26]:(3)Pijk=τijα⋅ηijβ∑i=1numτijα⋅ηijβ,where Pijk represents stochastically tour from node j is selected by ant k to be visited after node i, τij is the pheromone intensity on the jth node of the ith cell, ηij is visibility about the path that ant transited from the ith node to jth node, and α and β are weight factors of pheromone intensity and visibility, respectively. In order to calculate the visibility of the path, it can be expressed as(4)ηij=1|ρi−1j−1−ρij|+ρ0,where ρij represents the resistivity value of the jth node selected on the ith cell and ρ0 is half of the resistivity range of the ith cell. It can be shown that the smaller the difference between two cells, the larger the transition probability.At the beginning,k ants were randomly placed. Then, each ant decides the next node to visit according to the Pijk given by Equation (3). After iterating this process t times, each ant completes a tour. When ants are passing through the path, they generate new pheromone, and the intensity of the original pheromone would evaporate over time. Obviously, ants with the shorter tours will leave more pheromones than ants with the longer tours. Furthermore, when a group of ants passes this path, the pheromone will be updated. Therefore, the pheromone is updated as follows [23, 27]:(5)τijnew=λτijold+Δτijk,Δτijk=Qϕmn,if ant k pass through the nodei,j,0,otherwise,where τijnew represents the intensity after pheromone update, τijold represents the intensity of the original pheromone, λ∈0,1 is the volatility coefficient of pheromone trails, which avoids the infinite accumulation of pheromones, and Q is the total quantity of pheromone.To solve the problem of normal L1 norm inversion, the main methods include the iterative shrinkage thresholding/shaping method [28–30], the preconditioned conjugate gradient (PCG) method [31, 32], and split-Bregman method [33]. Among them, the PCG method has the advantages of solving the convenient inverse matrix and accelerating the convergence rate [34]. Therefore, the PCG method is used to solve the inversion problem in this paper.


### 2.2. Process of Ant Colony Optimization Inversion for Tunnel Detection

First, the 3D induced polarization method was used for forward modeling of advanced detection in tunnel [35]. The potential electrode M placed on the tunnel face and the current electrodes A1∼A4 are arranged on the tunnel outline, as shown in Figure 1. During the tunnel detection process, the four current electrodes on the tunnel surface are injected with the current at first. The potential electrodes on the tunnel face collect observed data, and then, the current electrode ring is moved away from the tunnel face and the current is injected again. Repeating this operation, the distance between the current electrode A and the potential electrode M gradually increases. The current electrode B and potential electrode N are arranged at a remote distance from the tunnel face.Figure 1
Schematic diagram for induced polarization advanced prediction in the tunnel. The position of electrodes B and N are at a remote distance from tunnel face.For 3D resistivity forward tunnel modeling problems, the finite element method is used in this paper. In particular, as for the model grid division, the inversion area in front of the tunnel is divided into 10 layers, and each layer is 3 m.Finally, the ACO algorithm for inversion is used, and the process is as follows:(1)
Initializing the parameters: input number of a group ant, upper limits of iteration, and other parameters.(2)
Ants choose paths based on the transition probabilityPijk.(3)
Calculate objective function valueΦmn for each ant, and select a global minimum model parameter up to now.(4)
Update pheromone intensityΔτijk and execute the next iteration. If the objective function Φmn bellows the tolerance ψ or the iteration times exceed the upper limits Ni, the final result will be outputted.


### 2.3. Numerical Simulation

To illustrate the efficiency of the ACO algorithm based on the L1 norm in the 1D layered stratum inversion of the tunnel, the algorithm was tested using two synthetic examples to simulate the water-bearing fault.The first model example is shown in Figure2(a). The size of the tunnel face was used is 8 m × 8 m. According to the calculation experience, the value of k0 and z0 are both 1 in the L1 norm and 1.8 and 1 in the L2 norm, respectively. An anomalous layer is located 12 m ahead of the tunnel face, and its thickness is 3 m. The background resistivity of the model is 1000 Ωm, and the anomalous layer with a resistivity of 400 Ωm.Figure 2
Inversion results for two models with low resistivity. (a) True resistivity of model 1. The results obtained after inversion of model 1 are shown in (b) using L1 norm and (c) using L2 norm. (d) Sketch of model 2. The results obtained after inversion of model 2 are shown in (e) using L1 norm and (f) using L2 norm. The black dotted frame represents the actual location of the low resistivity layer set by the model.
(a)(b)(c)(d)(e)(f)Figures2(b) and 2(c), respectively, show the inverted 1D geoelectrical models. Results are obtained using the L1 norm and L2 norm ACO. One can see that the location of the anomalous layer was reconstructed very well, and the depth from the tunnel face is about 12 m, which is consistent with the real distance. The resistivity of the anomalous in Figure 2(b) is about 530 Ωm, which is not much different from the real resistivity value. Comparing with the result of the L1 norm, the location of the low resistivity area is about 12 m and the resistivity is about 830 Ωm in Figure 2(c), but it is not reflecting the resistivity of the anomalous layer well.The second model consists of one anomalous layer with a low resistivity value (400 Ωm), which is located at 21 m alongz-axis. Other parameters’ settings are shown in model 1. As a result, the location of the low resistivity layer is about 21 m along z-axis, and the resistivity is about 420 Ωm in Figure 2(e); the position and region of the anomalous body are consistent with the actual one. In contrast, the result of the L2 norm deviates from the actual mode more than the result of the L1 norm. The position of the low resistivity anomalous is about 22 m along the z-axis, and the resistivity value is about 590 Ωm.From the two numerical simulations, the inversion results show that the L1 norm is more effective than the L2 norm in reflecting the position and resistivity of the low resistivity layer ahead of the tunnel face. These results will be verified in the field data.


### 2.4. Application Example

To examine the efficiency of the inversion algorithm on the tunnel detection, we inverted 1D resistivity data from the water conveyance project from Songhua River in the middle of Jilin province. The No.4 contract section of the main line construction is located between the Chalu River and Yinma River sections in Jilin City, the survey area at the tunnel mileage of 64 + 728 – 64 + 698 m. According to the geological reconnaissance, the surrounding rock in this section is mainly limestone, and the integrity of the rock mass is poor. Moreover, this section is with the scope of the fault.We adopted the induced polarization method proposed by Li et al. [16] to collect and process data. The inversion results of the L1 norm and L2 norm are illustrated by using the 1D ACO and are shown in Figure 3. The inversion grid division is consistent with numerical simulation. The number of iterations reached 20 when calculating the cut-off conditions. From Figure 3(a), it can be seen that the region of the low resistivity anomalous body is approximately from 0 m to 10 m, and the value of resistivity is approximately below 500 Ωm. Figure 3(b) shows the result of the observed data using the L2 norm. It can be observed that the low resistivity area is approximately 2 m to 14 m along the z-axis. The color of the low resistivity area is lighter than that of Figure 3(a), which means that its value is relatively high, around 500–600 Ωm.Figure 3
Inversion results for field data. (a) L1 norm inversion imaging. (b) L2 norm inversion imaging. The black dotted frame represents the low resistivity anomalous body.
(a)(b)According to the tunneling exposures shown in Figure4, tunnel exaction encountered water gushing at 64 + 728 m to 64 + 718 m. The influence range of the water body is 0 m to 10 m ahead tunnel face. It was revealed by field excavation, as shown in Figure 5, that the water inflow is 1500 m3 per hour at 64 + 727. The excavation result is consistent with the result obtained by L1 norm inversion. On the whole, the ACO inversion algorithm based on the L1 norm could better estimate the low resistivity anomalous bodies than the L2 norm.Figure 4
The complete section projection drawing of the geologic sketch and the excavation mileage is from 64 + 728 m to 64 + 698 m.Figure 5
The excavation results at 64 + 727.


## 2.1. The Ant Colony Optimization

The ant colony algorithm is a method that ants find the shortest path from their nest to a food source. The core of this behavior is the indirect communication between ants through chemical pheromone trails. When the ants are looking for food, they will initially explore the area around the nest in a random manner. During exploring, ants will leave chemical pheromone on the ground. Once an ant finds a food source, it evaluates the quantity and quality of the food and brings some back to the nest. The amount of pheromones ants left on the ground during the return depends on the quantity and quality of food. Pheromone will guide other ants to find a food source [23]. The main parameters are described below.The objective function of resistivity inversion based on theLp norm ACO algorithm is as follows:(1)Φmn=∑i=1N|ρdi−ρsiρdi|p,p=1,2,where Φmn is the data constraints, N represents the number of observed data, and ρdi and ρsi represent the ith observed data and predicted data, respectively. The objective function expression is expressed as the L1 norm when p = 1, otherwise it is L2 norm. In order to better use the resistivity method to detect the water-bearing fault zone ahead of the tunnel face, the L1 norm is adopted, i.e., p = 1. When Φmn reaches the global minimum, the foraging path chosen by ants is the shortest. The model parameters corresponding to the shortest path are the inversion results.As the electrode spacing increases, the sensitivity of the observed data to the deep also increases. The data weighting factor is introduced into the inversion problem to increase the weight of effective observed data. Therefore, the objective function of the ACO algorithm based on the L1 norm is expressed as follows [18, 24, 25]:(2)Φmn=∑i=1N|wiρdi−ρsiρdi|1,wi=k0zi+z0,where wi represents the data weighting factor, k0 and z0 are constants, and zi is the distance between tunnel face and current electrode of the ith circle along z direction.The probability of ants transferring from one cell to another is determined by pheromone intensity and visibility. The transition probability is expressed as follows [26]:(3)Pijk=τijα⋅ηijβ∑i=1numτijα⋅ηijβ,where Pijk represents stochastically tour from node j is selected by ant k to be visited after node i, τij is the pheromone intensity on the jth node of the ith cell, ηij is visibility about the path that ant transited from the ith node to jth node, and α and β are weight factors of pheromone intensity and visibility, respectively. In order to calculate the visibility of the path, it can be expressed as(4)ηij=1|ρi−1j−1−ρij|+ρ0,where ρij represents the resistivity value of the jth node selected on the ith cell and ρ0 is half of the resistivity range of the ith cell. It can be shown that the smaller the difference between two cells, the larger the transition probability.At the beginning,k ants were randomly placed. Then, each ant decides the next node to visit according to the Pijk given by Equation (3). After iterating this process t times, each ant completes a tour. When ants are passing through the path, they generate new pheromone, and the intensity of the original pheromone would evaporate over time. Obviously, ants with the shorter tours will leave more pheromones than ants with the longer tours. Furthermore, when a group of ants passes this path, the pheromone will be updated. Therefore, the pheromone is updated as follows [23, 27]:(5)τijnew=λτijold+Δτijk,Δτijk=Qϕmn,if ant k pass through the nodei,j,0,otherwise,where τijnew represents the intensity after pheromone update, τijold represents the intensity of the original pheromone, λ∈0,1 is the volatility coefficient of pheromone trails, which avoids the infinite accumulation of pheromones, and Q is the total quantity of pheromone.To solve the problem of normal L1 norm inversion, the main methods include the iterative shrinkage thresholding/shaping method [28–30], the preconditioned conjugate gradient (PCG) method [31, 32], and split-Bregman method [33]. Among them, the PCG method has the advantages of solving the convenient inverse matrix and accelerating the convergence rate [34]. Therefore, the PCG method is used to solve the inversion problem in this paper.


## 2.2. Process of Ant Colony Optimization Inversion for Tunnel Detection

First, the 3D induced polarization method was used for forward modeling of advanced detection in tunnel [35]. The potential electrode M placed on the tunnel face and the current electrodes A1∼A4 are arranged on the tunnel outline, as shown in Figure 1. During the tunnel detection process, the four current electrodes on the tunnel surface are injected with the current at first. The potential electrodes on the tunnel face collect observed data, and then, the current electrode ring is moved away from the tunnel face and the current is injected again. Repeating this operation, the distance between the current electrode A and the potential electrode M gradually increases. The current electrode B and potential electrode N are arranged at a remote distance from the tunnel face.Figure 1
Schematic diagram for induced polarization advanced prediction in the tunnel. The position of electrodes B and N are at a remote distance from tunnel face.For 3D resistivity forward tunnel modeling problems, the finite element method is used in this paper. In particular, as for the model grid division, the inversion area in front of the tunnel is divided into 10 layers, and each layer is 3 m.Finally, the ACO algorithm for inversion is used, and the process is as follows:(1)
Initializing the parameters: input number of a group ant, upper limits of iteration, and other parameters.(2)
Ants choose paths based on the transition probabilityPijk.(3)
Calculate objective function valueΦmn for each ant, and select a global minimum model parameter up to now.(4)
Update pheromone intensityΔτijk and execute the next iteration. If the objective function Φmn bellows the tolerance ψ or the iteration times exceed the upper limits Ni, the final result will be outputted.


## 2.3. Numerical Simulation

To illustrate the efficiency of the ACO algorithm based on the L1 norm in the 1D layered stratum inversion of the tunnel, the algorithm was tested using two synthetic examples to simulate the water-bearing fault.The first model example is shown in Figure2(a). The size of the tunnel face was used is 8 m × 8 m. According to the calculation experience, the value of k0 and z0 are both 1 in the L1 norm and 1.8 and 1 in the L2 norm, respectively. An anomalous layer is located 12 m ahead of the tunnel face, and its thickness is 3 m. The background resistivity of the model is 1000 Ωm, and the anomalous layer with a resistivity of 400 Ωm.Figure 2
Inversion results for two models with low resistivity. (a) True resistivity of model 1. The results obtained after inversion of model 1 are shown in (b) using L1 norm and (c) using L2 norm. (d) Sketch of model 2. The results obtained after inversion of model 2 are shown in (e) using L1 norm and (f) using L2 norm. The black dotted frame represents the actual location of the low resistivity layer set by the model.
(a)(b)(c)(d)(e)(f)Figures2(b) and 2(c), respectively, show the inverted 1D geoelectrical models. Results are obtained using the L1 norm and L2 norm ACO. One can see that the location of the anomalous layer was reconstructed very well, and the depth from the tunnel face is about 12 m, which is consistent with the real distance. The resistivity of the anomalous in Figure 2(b) is about 530 Ωm, which is not much different from the real resistivity value. Comparing with the result of the L1 norm, the location of the low resistivity area is about 12 m and the resistivity is about 830 Ωm in Figure 2(c), but it is not reflecting the resistivity of the anomalous layer well.The second model consists of one anomalous layer with a low resistivity value (400 Ωm), which is located at 21 m alongz-axis. Other parameters’ settings are shown in model 1. As a result, the location of the low resistivity layer is about 21 m along z-axis, and the resistivity is about 420 Ωm in Figure 2(e); the position and region of the anomalous body are consistent with the actual one. In contrast, the result of the L2 norm deviates from the actual mode more than the result of the L1 norm. The position of the low resistivity anomalous is about 22 m along the z-axis, and the resistivity value is about 590 Ωm.From the two numerical simulations, the inversion results show that the L1 norm is more effective than the L2 norm in reflecting the position and resistivity of the low resistivity layer ahead of the tunnel face. These results will be verified in the field data.


## 2.4. Application Example

To examine the efficiency of the inversion algorithm on the tunnel detection, we inverted 1D resistivity data from the water conveyance project from Songhua River in the middle of Jilin province. The No.4 contract section of the main line construction is located between the Chalu River and Yinma River sections in Jilin City, the survey area at the tunnel mileage of 64 + 728 – 64 + 698 m. According to the geological reconnaissance, the surrounding rock in this section is mainly limestone, and the integrity of the rock mass is poor. Moreover, this section is with the scope of the fault.We adopted the induced polarization method proposed by Li et al. [16] to collect and process data. The inversion results of the L1 norm and L2 norm are illustrated by using the 1D ACO and are shown in Figure 3. The inversion grid division is consistent with numerical simulation. The number of iterations reached 20 when calculating the cut-off conditions. From Figure 3(a), it can be seen that the region of the low resistivity anomalous body is approximately from 0 m to 10 m, and the value of resistivity is approximately below 500 Ωm. Figure 3(b) shows the result of the observed data using the L2 norm. It can be observed that the low resistivity area is approximately 2 m to 14 m along the z-axis. The color of the low resistivity area is lighter than that of Figure 3(a), which means that its value is relatively high, around 500–600 Ωm.Figure 3
Inversion results for field data. (a) L1 norm inversion imaging. (b) L2 norm inversion imaging. The black dotted frame represents the low resistivity anomalous body.
(a)(b)According to the tunneling exposures shown in Figure4, tunnel exaction encountered water gushing at 64 + 728 m to 64 + 718 m. The influence range of the water body is 0 m to 10 m ahead tunnel face. It was revealed by field excavation, as shown in Figure 5, that the water inflow is 1500 m3 per hour at 64 + 727. The excavation result is consistent with the result obtained by L1 norm inversion. On the whole, the ACO inversion algorithm based on the L1 norm could better estimate the low resistivity anomalous bodies than the L2 norm.Figure 4
The complete section projection drawing of the geologic sketch and the excavation mileage is from 64 + 728 m to 64 + 698 m.Figure 5
The excavation results at 64 + 727.


## 3. Conclusions

In this paper, an ACO algorithm based on the L1 norm is proposed for advanced tunnel detection. The key component of this approach is using the 3D induced polarization method to collect and process the data in the forward modeling. The inversion adopts the data weighting L1 norm ACO, which can increase the weight of effective data in the deep, so as to locate the deep anomalous bodies.The proposed method provides an ACO algorithm using the L1 norm for tunnel advance detection. Compared to the L2 norm method, the developed method produces a better effect in identifying the sharp boundary of low resistivity anomalous bodies.The developed method was tested using the two low resistivity models to simulate the water-bearing fault, and the 1D distribution of low resistivity body can be recovered by this inversion method. The method was applied to the field detection of the water conveyance project from Songhua River in the middle of Jilin province, which was verified through tunnel excavation.



---

*Source: 6649454-2021-04-13.xml*