Global optimization of controlled source audio-frequency magnetotelluric data with an improved artificial bee colony algorithm

Journal Pre-proof Global optimization of controlled source audio-frequency magnetotelluric data with an improved artificial bee colony algorithm

Laifu Wen, Jiulong Cheng, Fei Li, Jiahong Zhao, Zhihao Shi, Hongchuan Zhang PII:

S0926-9851(18)30827-9

DOI:

https://doi.org/10.1016/j.jappgeo.2019.103845

Reference:

APPGEO 103845

To appear in:

Journal of Applied Geophysics

Received date:

19 September 2018

Revised date:

10 September 2019

Accepted date:

10 September 2019

Please cite this article as: L. Wen, J. Cheng, F. Li, et al., Global optimization of controlled source audio-frequency magnetotelluric data with an improved artificial bee colony algorithm, Journal of Applied Geophysics(2018), https://doi.org/10.1016/ j.jappgeo.2019.103845

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2018 Published by Elsevier.

Journal Pre-proof

Global optimization of controlled source audio-frequency magnetotelluric data with an improved artificial bee colony algorithm Laifu Wen1, 2 , Jiulong Cheng1* , Fei Li1, 3 , Jiahong Zhao1 , Zhihao Shi1 , Hongchuan Zhang1 1.College of Geoscience and Surveying Engineering, China University of M ining & Technology (Beijing), Beijing, 100083, China 2. School of Earth Science and Engineering, Hebei University of Engineering, Handan, 056038, China 3. Key Laboratory of M ine Disaster Prevention and Control, North China Institute of Science and Technology, Beijing 101601, China

Abstract: The inversion of controlled source audio-frequency magnetotelluric (CSAMT) data is a complex nonlinear problem. A linearization of this problem is easily trapped in local minima and the complexity of the artificial source makes CSAMT data interpretation more difficult than that of magnetotelluric (MT) data. This paper presents an improved artificial bee colony (ABC) algorithm for the 2.5D inversion of CSAMT data. New

f

initialization and generation strategies are proposed to improve the optimization achieved by the original ABC

oo

algorithm. The global optimization of CSAMT by the improved ABC algorithm is realized based on 2.5D forward modeling theory and is used in the inversion of a complex model of water-bearing anomalous bodies in sandstone.

pr

Results show that the algorithm can accurately recover the resistivity and spatial distribution of strata and anomalous bodies. The survey data for a suspected collapse column in Shandong Province are also processed

e-

using the proposed method, and inversion via the algorithm accurately shows the water abundance of the suspected collapse column. Thus, the results of the theoretical modeling and practical data indicate that the

Pr

improved ABC algorithm is effective for analyzing CSAMT data. Moreover, this algorithm improves the interpretational accuracy and resolution of CSAMT data.

Keywords: Controlled source audio-frequency magnetotelluric; Improved artificial bee colony algorithm; 2.5D;

rn

1. Introduction

al

Nonlinear inversion

Controlled source audio-frequency magnetotelluric (CSAMT, Goldstein and Strangway, 1975; Zonge and

Jo u

Hughes, 1991) survey is an electromagnetic prospecting method that works in the frequency domain. The use of an artific ial source gives the method a good anti-interference ability. CSAMT has been widely applied for oil and gas exploration, mineral, geothermal and structural mapping (Boerner et al., 1993; Savin et al., 2001; Spichak et al., 2002; Spichak and Manzella, 2009; Bastani et al., 2011; Wu et al., 2012; Aykaç et al., 2015; An et al., 2016) and aquifer exploration (Pedrera et al., 2016; Šumanovac et al., 2018). Traditional CSAMT data processing methods can be divided into two categories: 1) where the data are inverted by MT formulations because the far field satisfies the plane wave condition (Sasaki et al., 1992; Wannamaker, 1997) and 2) invert the CSAMT data with controlled source taken into consideration. There are, however, problems with both methods: direct inversion of far field data does not take the horizontal attenuation of CSAMT into account and the boundary of far field data is difficult to determine. Although the existing correction methods in the near field have been successfully applied in some cases, the correction effect is not ideal for unknown geoelectrical structures and it is probable that there will be correction errors. Routh et al. (1999) and

*

Corresponding author E-mail addresses: [email protected] (Laifu Wen), [email protected] (Jiulong Cheng) 1

Journal Pre-proof Wang et al. (2015) have conducted 1D full CSAMT inversion without near field correction. However, this method is only suitable for horizontally layered media, which is not consistent with the actual state of underground media and is not accurate enough for the inversion of geoelectrical structures with a high degree of lateral variation. For practical 3D inversion problems, although successive scholars have implemented effective CSAMT/MT 3D inversion methods including the non-linear conjugate gradient (NLCG) (Newman et al., 2000; Rodi et al., 2001; Mackie et al., 2007; Commer et al., 2008; Lin et al., 2011, 2018), the Broyden-Fletcher-Goldfarb-Shanno (BFGS) (Haber, 2007; Plessix et al., 2008; Cao et al., 2016; Cao et al., 2018; Wang et al., 2017; Wang et al., 2018), the quasi-Newton method (QN) (Schwarzbach and Haber, 2013), and Gauss-Newton (GN) (Grayver et al., 2013; Oldenburg et al., 2013) algorithms, the computational efficiency of 3D inversion is still enormously challenging, and it has not been widely used in actual data processing. Thus, a 2.5D CSAMT inversion algorithm is alternative

f

in this case (Li et al., 2016).

oo

In addition, the CSAMT inversion problem is nonlinear and traditional linearization or local linearization methods may easily trapped into local minima. The computational intelligence algorithms have been applied that

pr

can achieve global optimization, having little dependence on the initial model and being more able to avoid converging to local minima. The rise of computational intelligence algorithms has provided a new method for

e-

CSAMT data inversion. An artificial bee colony (ABC) algorithm was proposed by Karaboga (2005) for solving multidimensional and multimodal optimization problems. This algorithm is inspired by the intelligent foraging

Pr

behavior of a honeybee swarm, with the position of food sources representing the optimal solution. The food source is exploited by information exchange and collaboration between bees with three different roles (employed

its use for CSAMT inversion.

al

bees, onlookers, and scouts) (Karaboga et al. 2008, 2011; Singh, 2009; Akay, 2013). No report has been made on

rn

This study proposes an optimized inversion method for 2.5D CSAMT data using an improved ABC algorithm to simulate the 3D characteristics of the artificial source and achieve inversion of CSAMT data in

Jo u

complex media. A new initialization and generation strategy is proposed to improve the optimization performance of the standard ABC algorithm. The validity of the proposed algorithm for CSAMT inversion is verified by applying it to a complex synthetic model containing abnormal water-bearing bodies and for the water abundance detection of a suspected collapse column in a coal mine in Shandong Province, China.

2. Forward modeling 2.1 Forward theory

As shown in Fig.1, a 3D source over 2D earth is the so-called 2.5D model. The conductivity  , magnetic permeability  and permittivity  are constant in this direction and only change in the x  z plane. Assuming the time dependence is eiwt , Maxwell’s equations with an electric source can be written as (Ward and Hohmann, 1988): {

⃑ ∇ × 𝐸⃑ = −𝑖𝑤𝜇𝐻 ⃑ = 𝑖𝑤𝜀𝐸⃑ + 𝜎𝐸⃑ + ⃑⃑𝐽𝑠 ∇×𝐻

(1)

⃑ are the electric field and magnetic field vector. On the left sides of the equations are the curl of where 𝐸⃑ and 𝐻 the electric field and the magnetic field, respectively. w is the angular frequency, and ⃑⃑𝐽𝑠 is the imposed current density. 2

Journal Pre-proof

Fig. 1. Coordinate system for 2.5D modeling (The y -direction is taken as the strike direction of the structure in the modeling)

By applying the Fourier transform to equation (1) with respect to the y-direction (is the wave number) 

^

F ( x, k y , z, w)   F ( x, y, z, w)e

 ik y y



(2)

dy

we can obtain two coupled governing equations for Eˆ y and Hˆ y in the wave number k y domain (Stoyer and Greenfield, 1976; Hohmann, 1988; Mitsuhata, 2000, 2002):   1 Hˆ y  yˆ Eˆ y  yˆ Eˆ y  1 Hˆ y  ˆ ˆ y  ik y  ( 2 ( 2 ) ( 2 )  yE ) ( 2 )  x ke z x ke x z ke z z ke x  

(3)

f

  1   1  Jˆsy  ik y  ( 2 Jˆsx )  ( 2 Jˆsz )   x  z k k e e  

oo

  1 Eˆ y  zˆ Hˆ y  zˆ Hˆ y  1 Eˆ y  ˆ ˆ y  ik y   ( 2 ( 2 ) ( 2 )  zH ) ( 2 )  x ke z x ke x z ke z z ke x  

(4)

e-

pr

  1   1  zˆ  ( 2 Jˆsz )  ( 2 Jˆsx )  z ke  x ke 

derived from Eˆ y and Hˆ y (Mitsuhata, 2000).

 yˆ Eˆ y ik y Hˆ y ik y ˆ  2  2 J sz  Hˆ x  2 z ke z ke x ke  ，  Hˆ y ik y ik y Eˆ y zˆ Hˆ y zˆ ˆ ˆ  2 Jˆsx E     J sz z  z ke ke2 z ke2 x ke2  

zˆ ˆ J sx ke2

al

Hˆ y

(5)

rn

 ik y Eˆ y zˆ   Eˆ x   2 ke x ke2   yˆ Eˆ y ik y ˆ H    z  ke2 x ke2 

Pr

where zˆ  iw, yˆ  iw   , ke2  k y2  k 2 , k 2   zy ˆ ˆ  w2  iw . The other electric and magnetic fields can be

Jo u

In each element, Eˆ y and Hˆ y are represented by the bilinear interpolation (Mullen and Belytschko, 1982) and the isoparametric element is shown in Fig.2. Fig. 2. Isoparametric element ( De is the domain of an element)

The global coordinates and EM field u are represented as: 4

4

x   Nie xi ,

4

z   Nie zi ,

i 1

u   Nieui

i 1

(6)

i 1

The Galerkin weighted residual method (Zienkiewicz, 1977) was applied to equations (3) and (4), we get the following equation: (7) where Ne is the number of elements,

N ie is

(8) the interpolation function of the i-th node in the e-th element. The

integrals on the left sides are the area integrals of region De and the second integrals on the right s ides are the line integrals around the boundary of region De . Applying the equation (6) to equations (7) and (8) and assembling the elemental matrices into a global matrix, we obtain: 3

Journal Pre-proof KFˆ  B

(9) where K is the global coefficient matrix, Fˆ is the EM fields to be solved and B is the source term. Bi-conjugate gradient stabilized method (Van der Vorst, 1992) was applied to solve the equation (9). Once Eˆ y and Hˆ y are computed, other EM fields can be obtained using equation (5).Then the space-domain response E y

and H y are obtained through the inverse Fourier transform (Leppin 1992; Unsworth et al., 1993). 2.1 Verification of the 2.5D forward modeling The 2.5D FEM forward modeling is verified by the model of Farquharson and Oldenburg (2002) in Fig. 3a.

oo

f

The results were compared with those from the IE (integral equation) solution of Farquharson and Oldenburg (2002), and those from the FEM of Ansari and Farquharson (2014). As shown in Fig. 3b, the agreement is generally good. For a frequency of 3Hz, the maximum relative error between this paper’s FEM and IE is only 8.21% for the real part and 9.17% for the imaginary part. The mismatch (e.g., inobs 1

/ nobs

. Ansari and

pr

Farquharson, 2014) is only 1.75% and 2.37% respectively.

E2  E1 max( E2 , E1 )

e-

Fig. 3. Comparison of the 2.5D model results. (a) The 2D geometry for the prism. (b) Ex along a profile over the top of prism.

Pr

3. Inversion scheme 3.1 Standard ABC algorithm

al

The ABC algorithm was proposed by Karaboga (2005) for solving optimization problems. In ABC algorithm, the colony of artificial bees contains three groups of bees: employed bees, onlookers and scouts (Karaboga et al.

rn

2008, 2011). First half of the colony consists of the employed bees and the second half includes the onlookers. The bees who carried food sources information are called employed bees and who waited for choosing a food

Jo u

source are called onlookers. The scouts carry out random search for discovering new sources. When an employed bee abandons its food source, it becomes a scout. When solving optimization problems, the first step is to randomly generate FN solutions (food source positions) as the initial population in each cycle. Every solution is a D-dimens ional vector. Here, FN is the size of the food source, and D is the number of model parameters to be optimized. The position of the food source represents a possible solution to the optimization problem and the fitness calculated with equation (10) is used to evaluate the quality of the food source. After recording the optimal value of the food source at this cycle, all of the food sources are subjected to repeated cycles up to a maximum of MC by the artificial bees (Karaboga et al. 2011). fiti  1/ (1  fi )

(10)

where fiti denotes the quality (fitness) of the associated food source and fi is the objective function value of the optimization problem. After initialization, an employed bee generates a new food source by using equation (11) and evaluates its nectar quality (fitness value) by using equation (10). It compares this value with that of the old food source and replaces the old source with the new one if it is better. Otherwise, the position of the previous source is retained. 4

Journal Pre-proof This is called greedy selection (Karaboga et al. 2011). vij  xij  ij  xij  xkj 

(11)

where ij is a random number between [-1,1]. i {1, 2....., FN}, k {1, 2....., FN} , and k  i , j {1, 2....., D} . After the employed bees complete their search for all food sources, they return to the dance area and share the fitness value and position information of the memorized food sources with the onlooker bees. The onlooker bees evaluate the fitness values provided by employed bees and select a food source according to probability equation (12) (Karaboga et al. 2011). FN

Pi  fiti /  fitn

(12)

n 1

f

After the onlooker bees have chosen food sources, the employed bees generate new food sources by using

oo

equation (11) and evaluate their fitness values. Provided that the fitness value of the new food source is better than that of the previous one, the old food source is replaced with the new food source.

pr

In the ABC algorithm, a lack of improvement in the fitness value of a food source after a limited number of cycles indicates that the food source has become trapped in a local minimum and should be abandoned. The

e-

employed bee associated with the abandoned food source then becomes a scout, and a new food source will be randomly generated by equation (13) to replace the abandoned food source (Karaboga et al. 2011).

Pr

j j j xij  xmin  rand  0,1  xmax  xmin 

(13)

The food sources searched by bees represent all the possible solutions for the CSAMT inversion problem.

rn

al

That is, the search process for food sources by scout bees represents the search process for the minimum objective function. Thus, the key steps of this algorithm for the CSAMT inversion problem are described as follows. 1) The initial food sources are randomly generated. Each food source represents an inverted model, and the D-dimensional vector in each food source denotes the resistivity that must be inverted. 2) The fitness of all the

Jo u

food sources evaluated using equation (10), and the best one is memorized. 3) Global search is performed, and a new food source is generated using equation (11). 4) The probability values are calculated using equation (12), and the food source is selected. 5) When the algorithm is trapped, a new food source is randomly generated using equation (13) to jump out of the local minima. 6) The best food source, which is the resistivity obtained by inversion is memorized. 7) Cycling is performed and terminated with maximum MC. 3.2 Improved ABC algorithm 3.2.1 Improved initialization strategy In the standard ABC algorithm, food sources are randomly (uniform distribution) generated as the initial population. However, this randomness will limit the efficiency of the algorithm (Penev, 2014). As the probability theory indicates that the food sources have a 50% probability of being further from the optimal solution rather than closer to it (Rahnamayan et al., 2008). This paper has improved the initialization strategy based on the concept of Opposition-Based Learning (OBL) (Tizhoosh, 2005): the ABC generates a random initial population of food source positions, and their opposite solutions are calculated on the basis of the OBL and the fitness values are evaluated and compared. Food sources with high fitness values are retained as the initial population. The improved initialization strategy enhances the diversity of the population at the beginning of the cycle and 5

Journal Pre-proof increases the search field. The population with lower fitness values is abandoned, which is beneficial to improving the efficiency of the algorithm. 3.2.2 Improved generation strategy The standard ABC algorithm generates new food sources according to equation (11), with ij as the random perturbation. This controls the diversity of the population and avoids the solution becoming trapped in a local minimum. The value of ij directly affects the convergence rate of the algorithm: in the initial stage, the difference between each food source is relatively large so a larger ij can be used to decrease the number of iterations. As the difference between each food source decreases, the search approaches the optimum solution. Thus, a smaller ij is needed to decrease the perturbation on the food source position. Therefore, as the number of iterations increases, the perturbation ij should be adaptively modified. Gu et al. (2012) proposed an adaptive ABC

oo

f

algorithm to modify the perturbation ij that causes ij to change adaptively depending on the situation of the food source position in the iteration process. For employed bees and onlooker bees, is ( iter is the number of iterations): e



e

 rand  0.5  ， ij    other 

0.08iter max iter

0.08iter  max iter

e



4iter max iter

rand  0.5

4iter  max iter

pr

  ij    

e

(14)

other

e-

To avoid obtaining local minima, Gao et al. (2012) combined ABC and differential evolution (DE, Storn and Price, 1995) algorithm to propose an ABC algorithm based on the disturbance mechanism. That is to say, the best

Pr

solutionis selected for perturbation.

vij  xbest , j  ij  xij  xkj 

(15)

et al. (2012) and Gao et al. (2012):

al

In order to further improve efficiency, this paper proposes the following improvement strategies based on Gu

rn

vij  xbest , j n1  j n  ij  xij  xkj 

(16)

Jo u

In equation (15), the perturbation for the food source xbest , j is the J-th variable on the N-th cycle. Its fitness value is first compared with the old one. If the new food source has better fitness than the old source, the new food source should be retained, and the J-th variable cannot be changed in the (N+1)-th cycle. Conversely, if the new food source has a lower fitness, changing the J-th variable does not cause the iteration to move towards the optimal solution. Therefore, the J-th variable cannot be changed in the (N+1)-th cycle in this situation either. i.e., as shown in (16), j n 1  j n . The workflow of improved ABC algorithm is shown in Fig.4. Fig. 4. Workflow of improved ABC algorithm

3.2.3 The performance of the improved ABC algorithm In order to test the performance of the improved ABC algorithm, we compare it with the standard ABC algorithm and the particle swarm optimization (PSO, Kennedy and Eberhart, 1995; Fernández Martínez et al., 2010; Santilano et al., 2018; Godio and Santilano, 2018) algorithm. The Griewank and Rastrigin functions given by Krink et al. (2004) as shown in Table1 are used to evaluate their optimization performance. Table 1 Numerical benchmark functions 6

Journal Pre-proof Function name

Function

x 1 n 2 n cos( i )  1  xi   4000 i 1 i i 1

Griewank

f 

Rastrigin

f   ( xi2  10cos(2 xi )  10)

n

Ranges

M inimum value

xi  600

⃑ )=0 𝑓(0

xi  5.12

⃑ )=0 𝑓(0

i 1

The values of the control parameters of improved ABC algorithm, standard ABC algorithm and PSO are given in Table 2. The assigned values for PSO are the recommended ones from Krink et al. (2004). Table 2 Parameter values used in the experiments (w, inertia weight;  min ,  max , lower and upper bounds of the random velocity rule weight; dimension of the problem D=50; the maximum cycle M C=1000.) Improved ABC

Standard ABC

100

Colony size

200

w

1.0-0.7

employed bee

50% of the colony

Colony size

200

employed bee

50% of the colony

 min

0

onlooker bee

50% of the colony

onlooker bee

50% of the colony

 max

2.0

Limit

M C/2

Limit

M C/2

pr

oo

popSize

f

PSO

The convergence curves of the two functions are shown in Fig. 5a and Fig. 5b. For Griewank function, in the

e-

initial stage of evolutionary, PSO performs better convergence. With the increase of iterations, the convergence of the improved ABC is more prominent. Late in the iteration, its convergence ability is obviously better than that of

Pr

PSO and standard ABC. On Rastrigin function, the convergence ability of the three algorithms is at the same level at the beginning of iteration. However, with the succeeding iterations, the improved ABC shows a strong

al

convergence ability. Therefore, if a fixed value of fitting error is given, the improved ABC is more efficient.

rn

Fig. 5. Convergence curves of test functions. (a) Griewank function. (b) Rastrigin function.

3.3 Objective function

Jo u

The objective function of CSAMT inversion is 

obs cal 1 1 M N  ij  ij   M N j 1 i 1  ijobs

2

obs cal  1 1 M N  ij  ij     obs  M N j 1 i 1  ij

  

2

(17)

where M is the number of survey stations, N is the number of frequencies, and ijobs , ijcal , ijobs , ijcal are the measured and calculated apparent resistivity and phase, respectively. 3.4 Selection of inversion parameters Three control parameters must be set in the ABC algorithm: the number of food sources FN, the limit number of cycles L, and the maximum cycle number MC. The number of food sources can be used to balance the probability of obtaining the global optimum with the computation time (Apalak et al., 2014). The limit number is related to the ability to search globally. When the limit is too high, the ABC will have a weaker ability to search globally due to reduced randomness. Where it is too low, the randomness can be enhanced, but there will be a slower convergence rate. As one of the terminating conditions, the larger the value of MC, the more likely it is that global optimization will be reached but the longer the calculation time will be. If it is smaller, although the 7

Journal Pre-proof calculation time is shorter, the algorithm will be unable to converge (Irani and Nasimi, 2011). The optimal parameters were determined by the trial method. In other words, the third parameter is searched and optimized on the basis of other two given parameters. For a H-type layered model (Fig. 6a), Fig. 6b, c and d are the convergence of objective function with the number of food sources (FN), the convergence of objective function with the maximum number of iterations (MC) and the convergence of objective function with the limited number of times (L) respectively. Taking computing efficiency and precision into consideration，the parameters are selected as FN = D/2, MC = 80, L = 2MC/3. Fig. 6. Convergence curves of objective function for a test model (D is the number of model parameters to be optimized). (a) H-type layered model. (b) The convergence of objective function with the number of food sources. (c) The convergence of objective function

f

with the maximum number of iterations. (d) The convergence of objective function with the limited number of times.

oo

4. Inversion of synthetic data

To demonstrate the effectiveness of the improved ABC algorithm for CSAMT inversion, a three-layer model

pr

simplified from typical coalfields in North China was designed. It comprises quaternary, sandstone and limestone layers and contains two water-bearing anomalies in the sandstone, as shown in Fig. 7a. The observation system is

e-

shown in Fig. 7b. The horizontal electric dipole (HED) is located at the origin of the coordinate system; Idx refers to a 1000-m long x-direction HED. The survey stations are located between 4000 m and 5000 m in the x-direction

Pr

with spacing of 40 m. The distance in the y-direction is 10km. The frequencies were 1Hz, 2Hz, 4Hz, 8Hz, 16Hz, 32Hz, 64Hz, 128Hz, 256Hz, 512Hz, 1024Hz, and 2048Hz.

al

Fig. 7. Geoelectrical model and transmitter–receiver geometry. (a) Geoelectrical model. (b) Observation system

The initial model of the inversion uses the 1D Occam (Constable et al., 1987) inversion results (Fig. 8a),

rn

taking values according to the variation range between 50% and 150%. In addition, anomalous response data used

Jo u

for inversion are obtained by adding 5% noise to the forward results. For the 2.5D inversion algorithm, we used a desktop computer with a 3.20GHz Intel (R) Core (TM) i7-8700 CPU. The improved ABC codes was modified from the original code from Karaboga (2005). We carried out the inversion using a non-uniform grid with 93×73 elements in the x and z direction, respectively. It took about 116 hours for the inversion to convergence to the data misfit in equation 17 of 0.063 after 80 iterations. Fig. 8a, b and c are the inversion results obtained from 1D Occam, the standard ABC algorithm and the improved ABC algorithm respectively. It can be seen from the figures that both the low-resistance anomalies (black dashed boxes) and layers (black dashed lines) agree well with the true model in Fig. 8b and c. And the results in Fig. 8b and c are significantly better than those in Fig. 8a. What’s more, the location of anomalous bodies in Fig. 8c is more consistent with the true model than Fig. 8b. Fig. 8. Inverted results for a layered model containing low-resistance anomalies. (a) The result from 1D Occam. (b) The result from standard ABC. (c) The result from improved ABC.

Fig. 9a shows convergence of objective function obtained by standard ABC and the improved ABC. It can be seen that the convergence speed of the improved ABC is faster than standard ABC. Fig. 9b shows the fitting curves of the inverted result at x=4600 m for standard ABC and improved ABC respectively. The curves show that the inverted result at the survey points agrees well with the true model response over the entire frequency band. 8

Journal Pre-proof The comparison results show that the improved ABC agree well with the true model, and is obviously better than the 1D Occam and standard ABC. Fig. 9. Convergence curves and data fitting. (a) Convergence of the objective function. (b) Data fitting at x=4600 m.

5. Field data inversion To test the practicability of the improved ABC algorithm for CSAMT inversion, we applied it to the field data from a coal mine to investigate the water abundance in a suspected collapse column. 5.1 Geological overview The CSAMT data were collected from a coal mine in Shandong Province, China. The coal mine is a typical North China Permo-Carboniferous coalfield with a seam structure consisting of two mineable coal seams (Coal 3

oo

f

and Coal 16). Borehole data in Fig.10 show that the coal measure strata in the survey area have a clear sedimentary sequence and that the strata are relatively stable. Logging data in Fig.10 show that there is a gradual increase in the resistivity of each successive layer with depth.

pr

3D seismic survey data indicate a suspected collapse column in mining area No. 5. The anomaly is located in the western part of the mining area. The seismic section (Fig. 11) shows good continuity of reflection waves in

e-

coal 3, ash 3 and coal 16 outside the abnormal body. However, the reflection waves inside the anomalous body have poor continuity, and the anomalous body extends in the vertical direction. The planform shape of the

Pr

suspected collapse column is approximately elliptical in the floor of coal 3. The lengths of its long and short axes are 140 m and 80 m, respectively. CSAMT exploration was carried out to identify the water abundance in the

al

suspected collapse column.

rn

5.2 Data acquisition and setup for field work

There are 10 survey lines in the study area with a net extent of 100 m × 50 m (Fig. 12). In addition, the net

Jo u

extent in the area of the suspected collapse column is 25 m × 25 m. Frequencies between 0.125 Hz and 2048 Hz n/2

are used on the basis of the depth of the target layer, and the frequency point series is 2 . The transmitting dipole is 1300 m long, the transmitting current is 16 A, and it is 5400 m away from the first survey line. Line L1150, which passes through the suspected collapse column, was selected for analysis; it has 11 points with a spacing of 25 m. The cell grid is 99×73 for computing. Firstly, the field data were preprocessed by de-noising and filtering. And then the data were inverted by SCS2D and improved ABC algorithm respectively. SCS2D is a CSAMT inversion software developed by Zonge Company (USA) based on smooth-model inversion algorithm. Smooth-model inversion is a robust method for converting far-field CSAMT data to resistivity model cross-sections. SCS2D inverts observed apparent resistivity and impedance phase data from a line of soundings to determine resistivity in a model cross-section. For the improved ABC algorithm, the initial model is determined on the basis of the 50-150% variation range of the 1D Occam inversion results. Fig. 10. Geological section of the survey area and borehole data

Fig. 11. Suspected collapse column in the seismic section (the left is InLine168 and the right is CrossLine1136) 9

Journal Pre-proof Fig. 12. Layout of survey lines in study area (The coordinate system in the figure uses the kilometer net coordinates)

5.3 Inversion result The results of the inversion of line L1150 data are presented in Fig. 13. Fig. 13a shows the inverted results of SCS2D. The three black solid lines represent the Cenozoic floor interface and the coal floors of No. 3 and No. 16, respectively. The red solid lines show the spatial positions of the suspected collapse column based on 3D seismic interpretation. As observed in Fig. 13a, low resistivity is the key feature in the shallow area of 600 m, and the variation range is small, thereby indicating the electrical characteristics of the clay rocks in Quaternary and Neogene strata. At depths below 600 m, the resistivity values gradually increase, thereby reflecting the electrical characteristics of Permian and Carboniferous strata. At depths of 1100–1300 m, the resistivity values are generally

f

greater than 40 Ω.m, which reflect the electrical characteristics of Ordovician strata. At the measuring points of

oo

1300–1375, the resistivity contour near the coal 3 is evidently concave, and a closed low-resistance anomaly occurs in the longitudinal direction of 780–1100 m. Compared with the suspected collapse column represented by

pr

the red solid line from the 3D seismic interpretation, the spatial location of the low -resistance anomaly basically coincides with that of the collapse column. It is inferred that the suspected collapse column is rich in water.

e-

However, the shape of the low-resistance anomaly considerably differs from that of the seismic interpretation. Fig. 13b presents the inverted results of the improved ABC algorithm. The overall electrical characteristics of

Pr

the inverted resistivity profile are basically consistent with the results of SCS2D inversion. However, they provide additional details regarding resistivity to geological interpreters, primarily in the form of low -resistance anomaly. The closed low-resistance anomaly can be easily delineated from the figure at the measuring points of 1275–1375

al

in the longitudinal direction of 780–1100 m. This anomaly is conical in cross section, which is more consistent with the shape of the collapsed column from the 3D seismic interpretation represented by the red solid lines. And

rn

at the boundary of the inferred collapse column, the difference in resistivity is more pronounced. The distribution of the low-resistance anomaly infers that water is not abundant in all areas of the suspected subsidence column but

Jo u

only at the large measuring point of the CSAMT line. Data fits for the survey station 1350 is shown in Fig. 13d. According to the first term in equation (17), the result is 6.8% and 2.7% for SCS2D and improved ABC respectively. The inverted data by improved ABC fit better. This further confirms that the inversion is capable of effectively recovering the electrical characteristics in the suspected collapse column area and clearly present the underground electrical structure.

Fig. 13. Inverted results. (a) SCS2D inverted geoelectrical section. (b) Improved ABC inverted geoelectrical section. (c) Convergence of the objective function. (d) Data fitting at station 1350. (In figure 13a and 13b, the abscissa represents the survey stations, and the vertical coordinate denotes the depth)

6. Discussion The CSAMT inverse problem is nonlinear, but traditional inversion algorithms frequently linearize this problem. However, CSAMT inversion is nonunique, and many local minima are present. The ABC algorithm is a global optimization nonlinear algorithm that does not require calculating the Jacobian and Hessian matrices. Simultaneously, search is randomly performed in a given model space without restriction, and the mechanism of 10

Journal Pre-proof jumping out of local minima makes getting trapped in local minima difficult. However, the ABC algorithm also exhibits shortcomings, which focus on the selection of the initial model space, inversion parameters, and computational efficiency. For the ABC algorithm, the initial model space needs to be provided beforehand. An appropriate initial model space can effectively constrain the model parameters and enable the inversion to converge rapidly to the global minima. If the initial model space is larger than the true model, then the algorithm may easily get trapped in the local minima. Although the algorithm can theoretically search global minima after numerous iterations, this feature is inconsistent in solving practical problems. In the synthetic and field data inversion process, the initial model space is determined on the basis of 1D Occam’s inversion results in the current study. The inverted parameters of the algorithm are also critical. For example, when the food sources of the ABC algorithm are evidently less than the number of model parameters, only a small number of model parameters can be disturbed in each iteration, and thus, achieving global minima becomes difficult. If the number of iterations is too small, then the objective function cannot converge to the global minim a,

oo

f

either. In contrast with the method used by Karaboga et al. (2005), the current study adopts a trial method, and the selected inverted parameters are suitable only for the model used in this study. Although the generalization of the

pr

inverted parameters is poor for all the models, the inverted parameters selected using the trial method may be optimal relative to a specific problem. The last problem is computational time. Numerous forward modeling processes are required to search for the optimal food source of a bee colony. For example, 116 hours were spent on synthetic data inversion in this study, which may be a key factor that restricts the application of this approach

e-

to high-dimensional problems. The mechanism behind the use of supercomputers or cloud computing to accelerate the convergence of the ABC algorithm is a meaningful research direction.

Pr

To enhance the optimization ability of the ABC algorithm, this study improves the algorithm and then tests the improved version. The improved algorithm theoretically enhances optimization efficiency, but single calculation efficiency is insignificantly improved because of the increase in the number of calculation steps.

al

However, the comparison results shown in Fig. 5 indicate that the optimization ability of the improved algorithm is considerably improved. After performing the same number of iterations, the value of the objective function is

rn

smaller. That is, the improved ABC algorithm can achieve the set objective function value within a relatively small number of iterations. Thus, we use the same number of iterations in the theoretical model. The inversion

Jo u

results shown in Fig. 8 indicate that the improved algorithm is better for recovering the position and shape of low-resistance anomaly. After 80 iterations, the convergence curve of the objective function also shows that the improved algorithm is considerably superior to the standard algorithm. Given that the superiority of the improved algorithm has already been proven, only the existing commercial software and the improved algorithm inversion effect are compared in the field data processing. For the results of CSAMT inversion, low -resistance anomalies are found in the location of the collapse column during seismic interpretation. Moreover, the improved ABC algorithm inversion results provide more detailed information regarding the shape and scope of low -resistance anomaly. However, the range of low-resistance anomalies cannot completely coincide with seismic interpretation because seismic exploration can provide only structural information and cannot detect whether the collapse column is rich in water. This study interprets the results of CSAMT inversion on the basis of the results of seismic exploration. Although the accuracy of seismic exploration is high, interpretation errors may still occur. Thus , the CSAMT inversion results can be used to modify the seismic interpretation results. That is, the two methods can be interpreted on the basis of each other. However, a simple comprehensive interpretation does not make considerable sense; hence, joint inversion may perform better, which is investigated later in this study.

7. Conclusions 11

Journal Pre-proof This paper presented an ABC algorithm for 2.5D CSAMT inversion and improved the standard ABC algorithm in the following two aspects: initialization and generation strategies. The inverted parameters were determined using the trial method. The results of different test functions indicate that the improved ABC algorithm is superior to the standard ABC algorithm. Subsequently, the improved ABC algorithm was used to invert synthetic data computed from a complex model. The recovered model is highly consistent with the true model. Compared with the inverted results from 1D Occam’s method and the standard ABC algorithm, the proposed method improved inverted accuracy and provided more reasonable results that are consistent with the true model. We applied the improved ABC algorithm to field data collected in a coal mine. The obtained results exhibit features similar to those from synthetic data. This finding indicates that the propos ed method can deal well with measured data and achieve good results. In addition, although this algorithm has some problems, including the

f

selection of the initial model space, inversion parameters, and computational efficiency, the proposed method

oo

exhibits considerable application prospect and can be extended to invert other geophysical data.

pr

Acknowledgements

This research was supported by the National Key Research and Development Program of China (Grant

e-

no.2017YFC0801405), the National Natural Science Foundation of China (Grant no.51574250).

References

Pr

Akay, B., 2013. A study on particle swarm optimization and artificial bee colony algorithms for multilevel thresholding. Applied Soft Computing 13, 3066-3091. https://doi.org/ 10.1016/j.asoc.2012.03.072. An, Z.G., Di, Q.Y., 2016. Investigation of geological structures with a view to HLRW disposal, as revealed through 3D inversion of

al

aeromagnetic and gravity data and the results of CSAM T exploration. Journal of Applied Geophysics 135, 204-211. https://doi.org/ 10.1016/j.jappgeo.2016.10.013.

rn

Ansari, S., Farquharson, C. G., 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids. Geophysics 79, E149-E165.

Jo u

Apalak, M . K., Karaboga, D., and Akay, B., 2014. The artificial bee colony algorithm in layer optimization for the maximum fundamental frequency of symmetrical laminated composite plates. Engineering Optimization 46, 420-437. https://doi.org/ 10.1080/0305215X.2013.776551.

Aykaç, S., Timur, E., Sari, C., et al., 2015. CSAMT investigations of the Caferbeyli (M anisa/Turkey) geothermal area. Journal of Earth System Science 124, 149-159. https://doi.org/ 10.1007/s12040-014-0521-1. Bastani, M ., Savvaidis, A., Pedersen, L.B., Kalscheuer, T., 2011. CSRM T measurements in the frequency range of 1–250 kHz to map a normal fault in the Volvi basin, Greece. Journal of Applied Geophysics 75, 180-195. https://doi.org/10.1016/j.jappgeo. 2011.07. 001. Boerner, D. E., J. A. Wright, J. G. Thurlow, and L. E. Reed, 1993. Tensor CSAMT studies at the Buchans M ine in central Newfoundland. Geophysics 58, 12-19. https://doi.org/10.1190/1.1443342. Cao, H., Wang. K. P., Wang, T., Hua, B. G., 2018. Three-dimensional magnetotelluric axial anisotropic forward modeling and inversion. Journal of Applied Geophysics 153, 75-89. Cao, M ., Tan, H. D., Wang, K. P., 2016. 3D LBFGS inversion of controlled source extremely low frequency electromagnetic data. Applied Geophysics 13, 689-700. https://doi.org/10.1007/s11770-016-0585-6. Commer, M ., Newman, G. A., 2008. New advances in three-dimensional controlled-source electromagnetic inversion. Geophysical Journal International 172, 513-535. Constable, S. C., Parker, R., L., and Constable, C. G., 1987. Occam’s inversion: A practical algor ithm for generating smooth models 12

Journal Pre-proof from electromagnetic sounding data. Geophysics 52, 289-300. Farquharson, C. G., and D. W. Oldenburg, 2002. An integral equation solution to the geophysical electromagnetic forward-modelling problems. In Zhdanov, M . S., and Wannamaker, P. E., Eds., Three-dimensional electromagnetics. Elsevier, Inc., 3-19. Fernández M artínez, J. L., García-Gonzalo, E., Fernández Álvarez, J. P., Kuzma, H. A. and M enéndez -Pérez, C. O., 2010. PSO: A powerful algorithm to solve geophysical inverse problems. Journal of Applied Geophysics 71, 13-25. https:// doi.org/10.1016/j.jappgeo.2010.02.001. Gao, W., Liu, S., and Huang, L., 2012. A global best artificial bee colony algorithm for global optimization. Journal of Comp utational and Applied M athematics 236, 2741-2753. https://doi.org/ 10.1016/j.cam.2012.01.013. Godio, A., and Santilano, A., 2018. On the optimization of electromagnetic geophysical data: Application of the PSO algorithm. Journal of Applied Geophysics 148, 163-174. https://doi.org/10.1016/j.jappgeo.2017.11.016. Goldstein, M . A., and Strangway, D. W., 1975. Audio-frequency magnetotellurics with a grounded electric dipole source. Geophysics 40, 669-683. https://doi.org/10.1190/1.1440558.

f

Grayver, A. V., Streich, R., Ritter, O., 2013. Three-dimensional parallel distributed inversion of CSEM data using a direct forward

oo

solver. Geophysical Journal International 193, 1432-1446.

Gu, W. X., Yin, M . H., and Wang, C. Y., 2012. Self-Adaptive Artificial Bee Colony for Global Numerical Optimization: 2nd

pr

International Conference on M echanical, Industrial, and M anufacturing Engineering. Ieri Procedia 1, 59-65. Haber, E., Oldenburg, D. W., Shekhtman, R., 2007. Inversion of time domain three-dimensional electromagnetic data. Geophysical Journal International 2, 1268-1277.

e-

Hohmann, G. W., 1988. Numerical modeling for electromagnetic methods of geophysics. In Nabighian, M . N., Ed., Electromagnetic methods in applied geophysics, Vol. 1, Theory: Society of exploration geophysics.

Pr

Irani, R., and Nasimi, R., 2011. Application of artificial bee colony -based neural network in bottom hole pressure prediction in underbalanced drilling. Journal of Petroleum Science and Engineering 78, 6-12. https://doi.org/ 10.1016/j.petrol.2011.05.006. Karaboga, D., 2005. An idea based on honey bee swarm for numerical optimization. Technical Report -TR06, Erciyes. University,

al

Engineering Faculty, Computer Engineering Department. Karaboga, D., and Basturk, B., 2008. On The Performance of Artificial Bee Colony (ABC) Algorithm. Applied Soft Computing 8,

rn

687-697. https://doi.org/ 10.1016/j.asoc.2007.05.007.

Karaboga, D., and Ozturk, C., 2011. A novel clustering approach: artificial bee colony (ABC) algorithm. Applied Soft Computing 11,

Jo u

652-657. https://doi.org/ 10.1016/j.asoc.2009.12.025. Kennedy, J., and Eberhart, R., 1995. Particle Swarm Optimization. Proceedings of IEEE International Conference on Neural Networks, Extended Abstracts, 1942-1948. Krink, T., Filipic, B., Fogel, G. B., Thomsen, R., 2004. Noisy optimization problems—a particular challenge for differential evolution? In: Proceedings of 2004 Congress on Evolutionary Computation, IEEE Press, 332-339. Leppin, M ., 1992. Electromagnetic modeling of 3-D sources over 2-Dinhomogeneties in the time domain. Geophysics 57, 994-1003. Li, W. B., Zeng, F. A., Li, J., Chen, X., Wang, K., Xia, Z., 2016. 2.5D forward modeling and inversion of frequency -domain airborne electromagnetic data. Applied Geophysics 13, 37-47. https://doi.org/ 10.1007/s11770-016-0548-y. Lin, C. H., Tan, H. D., Tong, T., 2011. Three-dimensional conjugate gradient inversion of magnetotelluric impedance tensor data. Journal of Earth Science 22, 386-395. Lin, C. H., Tan, H. D., Wang, W. W., Tong, T., Peng, M ., Wang, M ., Zeng, W. H., 2018. Three-dimensional inversion of CSAM T data in the presence of topography. Exploration Geophysics 49, 253-267. http://dx.doi.org/10.1071/EG16067. M ackie, R., Watts, M . D., Rodi, W., 2007. Joint 3D inversion of marine CSEM and MT data. SEG 77th Annual International M eeting, Expanded Abstracts 574–578. M itsuhata, Y., 2000. 2-D electromagnetic modeling by finite-element method with dipole source and topography. Geophysics 65, 465-475. https://doi.org/10.1190/1.1444740. M itsuhata, Y., Uchida, T. and Amano, H., 2002. 2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source. Geophysics 67, 1753-1768. https://doi.org/10.1190/1.1527076. 13

Journal Pre-proof M ullen, R., Belytschko, T., 1982. Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation. International Journal for Numerical M ethods in Engineering 18, 11-29. Newman, G. A., Alumbaugh, D. L., 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophysical Journal International 140, 410-424. Oldenburg, D. W., Haber, E., Shekhtman, R., 2013. Three dimensional inversion of multisource time domain electromagnetic data. Geophysics 78, E47-E57. Pedrera, A., M artos-Rosillo, S., Galindo-Zaldívar, J., Rodríguez-Rodríguez, M ., Benavente, J., M artín-Rodríguez, J. F., Zúñiga-López, M . I., 2016. Unravelling aquifer-wetland interaction using CSAM T and gravity methods: the M ollina-Camorra aquifer and the Fuente de Piedra playa-lake, southern Spain. Journal of Applied Geophysics 129, 17-27. Penev, K., 2014. Free Search-comparative analysis 100. International Journal of M et heuristics 3, 118-132. https://doi.org/ 10.1504 /IJM HEUR.2014.063142.

Inverse Problems 24, 34012. https://doi.org/ 10.1088/0266-5611/24/3/034012.

f

Plessix, R. É., and M ulder, W. A., 2008. Resistivity imaging with controlled-source electromagnetic data: depth and data weighting.

oo

Rahnamayan, S., Tizhoosh, H. R., and Salama, M . M . A., 2008. Opposition-Based differential evolution. IEEE Transactions on Evolutionary Computation 12, 64-79. https://doi.org/ 10.1109/TEVC.2007.894200.

pr

Rodi, W., M ackie, R. L., 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics 66:174-187. Routh, P. S., and Oldenburg, D. W., 1999. Inversion of controlled source audio-frequency magnetoteluric data for a horizontal layered earth. Geophysics 64, 1689-1697. https://doi.org/ 10.1190/1.1826613.

e-

Santilano, A., Godio, A., and M anzella, A., 2018. Particle swarm optimization for simultaneous analysis of M agnetotelluric (M T) and Time Domain EM (TDEM ) data. Geophysics 83, 1-48. https://doi.org/10.1190/geo2017-0261.1.

Pr

Sasaki, Y., Yoneda, Y., and M atsuo, K., 1992. Resistivity imaging of controlled-source audio frequency magnetotelluric data. Geophysics 57, 952-955.

Savin, C., Ritz, M ., Join, J.L., Bachelery, P., 2001. Hydrothermal system mapped by CSAMT on Karthala volcano Grand Comore

al

Island, Indian Ocean. Journal of Applied Geophysics 48, 143-152. https://doi.org/10.1016/S0926-9851(01)00078-7.

International 193, 615-634.

rn

Schwarzbach, C., Haber, E., 2013. Finite element based inversion for time-harmonic electromagnetic problems. Geophysical Journal

Singh, A., 2009. An artificial bee colony algorithm for the leaf-constrained minimum spanning tree problem. Applied Soft

Jo u

Computing 9, 625-631. https://doi.org/10.1016/j.asoc.2008.09.001. Spichak, V., Fukuoka, K., Kobayashi, T., Mogi, T., Popova, I., Shima, H., 2002. ANN reconstruction of geoelectrical parameter s of the M inou fault zone by scalar CSAM T data. Journal of Applied Geophysics 49, 75-90. https://doi.org/ 10.1016/S0926-9851(01) 00100-8.

Spichak. V., and M anzella, A., 2009. Electromagnetic sounding of geothermal zones. Journal of Applied Geophysics 68, 459 -478. https://doi.org/10.1016/j.jappgeo.2008.05.007. Storn, R., and Price, K., 1995. Differential Evolution – A Simple and Efficient Adaptive Scheme for Global Optimization over Continuous Spaces. Technical Report TR-95-012, International Computational Science Institute. Stoyer, C. H., and Greenfield, R. J., 1976. Numerical solutions of the response of a two-dimensional earth to an oscillating magnetic dipole source. Geophysics 41, 519-530. Šumanovac, F., Orešković, J., 2018. Exploration of buried carbonate aquifers by the inverse and forward modelling of the Cont rolled Source Audio-M agnetotelluric data. Journal of Applied Geophysics 153, 47-63. https://doi.org/ 10.1016/j.jappgeo.2018.04.007. Tizhoosh, H. R., 2005. Opposition-based learning: A new scheme for machine intelligence: International Conference on Computational Intelligence for M odeling, Control and Automation and International Conference on Intelligent Agents. Web Technologies and Internet Commerce 1, 695-701. Unsworth, M . J., Travis, B. J., and Chave, A. D., 1993. Electromagnetic induction by a finite electric dipole source over a 2-D earth. Geophysics 58, 198-214. Van der Vorst, H. A., 1992. Bi-CGSTAB: A Fast and Smoothly Convergin g Variant of Bi-CG for the Solution of Non-symmetric 14

Journal Pre-proof Linear Systems. SIAM Journal on Scientific and Statistical Computing 13, 631-644. Wannamaker, P. E., 1997. Tensor CSAMT survey over the Sulphur Springs thermal area, Valles Caldera, New M exico, U.S.A.; P art 1, Implications for structure of the western caldera. Geophysics 62, 451-465. Wang, K. P., Tan, H. D., Lin, C. H., Yuan, J.L., Wang, C., Tang, J., 2018. Three-dimensional tensor controlled-source audio-frequency magnetotelluric inversion using LBFGS. Exploration Geophysics 49, 268-284. https://doi.org/ 10.1071/EG16079. Wang, R., Yin, C.C., Wang, M .Y., Di, Q.Y., 2015. Laterally constrained inversion for CSAMT data interpretation. Journal of Applied Geophysics 121, 63-70. https://doi.org/ 10.1016/j.jappgeo.2015.07.009. Wang, T., Wang, K. P., Tan H. D., 2017. Forward modeling and inversion of tensor CSAMT in 3D anisotropic media. Applied Geophysics 14, 590-605. Ward, S. H., and Hohmann, G. W., 1988. Electromagnetic theory for geophysical applications Electromagnetic M ethods in Applied Geophysical, Vol. 1, Theory, in Nabighian, M . N., Ed., Society of exploration geophysics.

CSAM T in Jiangxia, Wuhan, China. Journal of Earth Science 23, 757-767.

oo

Zienkiewicz, O. C., 1977. The finite element method, third edition: M cGraw-Hill.

f

Wu, G. J., Hu, X. Y., Huo, G. P., Zhou, X.C., 2012. Geophysical exploration for geothermal resources: An application of MT an d

Zonge, K. L., Hughes, L. J., 1991. Controlled Source Audio-Frequency M agnetotellurics. In Nabighian, M . N., Ed., Electromagnetic

Jo u

rn

al

Pr

e-

pr

M ethods in Applied Geophysics, Vol. 2, Society of Exploration Geophysicists.

15

Journal Pre-proof Dear Editors: No conflict of interest exits in the submission of this manuscript, and manuscript is approved by all authors for publication. I would like to declare on behalf of my co-authors that the work described was original research that

Jo u

rn

al

Pr

e-

pr

oo

f

has not been published previously. All the authors listed have approved the manuscript that is enclosed.

16

Journal Pre-proof This paper presents an improved artificial bee colony (ABC) algorithm for 2.5D inversion of CSAMT data. New initialization and generation strategies are proposed to improve the optimization achieved by the original ABC algorithm. Global optimization of CSAMT by the improved ABC algorithm is realized based on 2.5D forward modeling theory and is used for the inversion of a complex model of water-bearing anomalous bodies in sandstone. The results show that the algorithm can accurately restore the resistivity and spatial distribution of strata and anomalous bodies. The survey data for a suspected collapse column in Shandong Province are also processed using this method, and it is shown that inversion via the algorithm accurately reveals the water abundance of the suspected collapse column. Thus, the results for theoretical modeling and practical data show that the improved ABC algorithm is effective for analyzing CSAMT data and improves the interpretational accuracy and resolution

Jo u

rn

al

Pr

e-

pr

oo

f

of CSAMT data.

17

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13