One-dimensional magnetotelluric inversion using levenberg-marquardt and particle swarm optimization algorithm

Magnetotelluric is a geophysical method which uses natural time variations of the Earth’s magnetic and electric fields. One-dimensional magnetotelluric inversion was carried out to estimate the electrical resistivity distribution of rocks which varied with depth. Non-linear relationships between model parameters and observation data lead to difficulties in the inversion process. This problem can be solved with a non-linear inversion method with a global search technique. One of the global search technique that gaining interest in the inverse problem is Particle Swarm Optimization (PSO). PSO can find acceptable solutions from a very broad set of initial parameters. This study also carried out a non-linear local search technique, namely Levenberg-Marquardt (LM) which has been used globally for the geophysical problem. Then the two inversion algorithm was compared in determining the best model solution. There are two data used in this study which are generated from a three-layer model, namely data without noise and also data that has been added by Gaussian noise with 10% standard deviation. It was found that the results of the two inversions were quite good at determining the actual model. LM algorithm is an algorithm that is truly influenced by the initial value of the damping factor ε2 while PSO algorithm depends on several parameters, namely number of particles n, inertia weight w, cognitive parameter c1 and social parameter c2. Trial experiments suggested that the global best solution could be achieved with controlling parameters wmin = 0.4, wmax = 0.9, c1 = 1.3 and c2 = 1.5 without noise data and wmin = 0.7, wmax = 0.9, c2 = 0.8 and c2 = 1.4 for data with 10% noise.


Introduction
Modeling and inversion in magnetotelluric have shown numerous successes that have been applied in various fields, such as geothermal [1,2], oil and gas [3], and also to study fault and Earth's lithosphere [4,5]. However, modeling and inversion for the 1D case are still done mainly to find the initial results before modeling on a higher dimension. The interpretation of magnetotelluric data needs to go through the modeling and inversion process. The inversion process for some algorithms sometimes still encounter problems that are ill-posed problems where small causes in large changes in the data can change in the solution [6,7].
Many techniques have been developed to invert geophysical data to estimate model parameters and most practical geophysical inversion problems are non-linear. In general, there are two types of non- 3 To whom any correspondence should be addressed. linear inversion methods, namely local search techniques and global search techniques. One of the common local search techniques that is very popular in the geophysical inverse problems is Levenberg-Marquardt (LM) algorithm which is a linearized approach and uses a damping parameter in its algorithm [8]. For global search technique, Particle Swarm Optimization (PSO) is an algorithm gaining interest lately in geophysical inverse modeling proposed by Kennedy and Eberhart [9]. PSO can find an acceptable solution from a very wide initial set of parameters [10] and commonly considered simple to apply because they are easy to understand and computationally efficient [11,12].
In this paper, we will present the application of LM and PSO algorithm for the 1D magnetotelluric case. The two inversion methods will be compared to determine the best model solution. The observation data used is synthetic data in response to a simple three-layer synthetic model and also data with 10% Gaussian noise that has been added to the synthetic model.

Forward modelling of 1D magnetotelluric
One dimensional magnetotelluric (MT) inversion, also include forward modelling, is an iterative process. The formulation starts from the layer n (halfspaces) and toward to the above layer from layer n to layer 1 ( Figure 1).

Figure 1.
A homogeneous 1D layered earth with resistivity and thickness.
The forward modelling scheme for the 1D model response is the total impedance obtained from each of the 2 layers (Z j and Z j+1 ) following to the equation (1-3) [11].
where Z 0j is characteristic impedance of the j th layer, ρ j and h j are the resistivity and the thickness of the j th layer, ω is the angular frequency, j R is reflection coefficient, and µ 0 is the magnetic permeability. At layer 1 with N-layer, MT data are expressed as apparent resistivity (ρ a ) and phase (ϕ), which can be calculated by the equation (4).

Levenberg-Marquardt
The LM algorithm iterative equation is as follows [8], where n m is iteration models, I is an identity matrix and 2  is is a damping factor whose value is chosen by try and error until it gets the appropriate value. Meanwhile, J is a Jacobian matrix written as below:    (7) where N is the total number of the data point, obs d and cal d are observed data and calculated data.

Particle swarm optimization
The PSO algorithm is based on sharing social information among members of a species that offers evolutionary advantages [9]. This is inspired by the behavior of flocks of birds that are looking for food in their habitat. In implementing the algorithm, birds are represented by "particles" and operate by moving a group of particles around the search space according to a simple mathematical formula that depends on the position and velocity of the particles. The particles are consistently looking for better solutions.
The movement of particles is influenced by their own best-known position (p best ) and the best known overall position in the population (g best ) which is updated every iteration [7]. After finding the two best values, the particles update their speed and position with the following equation (8) and (9).  As a first step, a cluster of particles is randomly selected from the predetermined parameter search space, and the velocity of each particle is initialized to zero. Furthermore, the velocity and location of the particles are updated according to equation (8) and equation (9), with the constraint that the location of each particle does not exceed the specified parameter limit. Figure 3 outlines the important steps in implementing the PSO algorithm.

Results and discussion
In this study, MT data were inverted using the LM and PSO algorithm to get rock resistivity's and layer thicknesses. The data used is synthetic data from the 3 layer model, namely data without noise and data of 10% noise. The results of inversion using the LM and PSO algorithms can be seen in Figure 4 (synthetic data without noise) and Figure 5 (synthetic data with 10% noise). In general, inversion using the LM and PSO algorithms is as good as that which uses synthetic data without noise or synthetic data that has added noise. In general, inversion using the LM algorithm is quite good at determining the actual model as well as for the PSO algorithm. Inversion with the LM algorithm tends to have a smaller error compared to the PSO method as can be seen in Table 1 and Table 2. For the inversion of apparent resistivity and phase data, LM is better than PSO (see Figure 4 and Figure 5). As can be seen in Table 1, inversion using the LM algorithm for data without noise can predict the resistivity and thickness of the layers quite well. But for the second layer, the thickness of the layer obtained by using the LM algorithm is 1500 m while the actual model has a thickness of 1000 m. Similar results are also obtained with the LM algorithm for data without noise where the result in the second layer has a difference in a layer thickness of 600 m. Whereas the inversion results in determining resistivity values almost approach the actual model. For inversion using PSO algorithm, determining layer thickness is still better than using LM inversion. But in determining the resistivity value, the LM algorithm is better than the PSO algorithm. When compared to determining the initial parameters between the two inversion algorithms, LM is simple than PSO. From equation (1), LM depends only on parameters 2  as an initial parameter which can cause the solution not to converge if the given value is not appropriate. Usually, the value of the damping factor is given with a very large at the start of computing. Then, the PSO algorithm depends on several parameters such as a number of particles n, inertia weight w, cognitive parameter c 1, and social parameter c 2 as can be seen in equation (8). Determination of parameters w, c 1 and c 2 greatly affects the inversion result. Determination of these parameters was done by try and error until the inversion result is fit with the actual model. For the data without noise, the best PSO inversion result is by using the parameter w min = 0.4, w max = 0.9, c 1 = 1.3, c 2 = 1.5 and the number of particles n = 50. As for the data with noise, the best inversion results are using parameter w min = 0.7, w max = 0.9, c 1 = 0.8 and c 2 = 1.4 and the number of particles n = 50.

Conclusion
In this study, we report the capability of LM and PSO in solving MT inverse problem. We implemented the LM and PSO algorithm on 2D synthetic MT data with 10% noise and without noise. From the inversion results, LM algorithm is better than PSO algorithm for determining apparent resistivity and phase for data without noise and with noise. In addition, the LM algorithm is better for determining rock resistivity while PSO algorithm is better for determining layer thickness. LM algorithm is an algorithm that is truly influenced by the initial value of the damping factor 2  while PSO algorithm depends on several parameters, namely number of particles n, inertia weight w, cognitive parameter c 1, and social parameter c 2 . Trial experiments suggested that the global best solution could be achieved with controlling parameters w min = 0.4, w max = 0.9, c 1 = 1.3 and c 2 = 1.5 without noise data and w min = 0.7, w max = 0.9, c 1 = 0.8 dan c 2 = 1.4 for data with 10% noise.