Pareto Joint Inversion of 2D magnetometric and gravity data- synthetic study

Pareto joint inversion for two or more data sets is an attractive and promising tool which eliminates target functions weighing and scaling, providing a set of acceptable solutions composing a Pareto front. In former author’s study MARIA (Modular Approach Robust Inversion Algorithm) was created as a flexible software based on global optimization engine (PSO) to obtain model parameters in process of Pareto joint inversion of two geophysical data sets. 2D magnetotelluric and gravity data were used for preliminary tests, but the software is ready to handle data from more than two geophysical methods. In this contribution, the authors’ magnetometric forward solver was implemented and integrated with MARIA. The gravity and magnetometry forward solver was verified on synthetic models. The tests were performed for different models of a dyke and showed, that even when the starting model is a homogeneous area without anomaly, it is possible to recover the shape of a small detail of the real model. Results showed that the group analysis of models on the Pareto front gives more information than the single best model. The final stage of interpretation is the raster map of Pareto front solutions analysis.


Introduction
Nowadays, joint inversion undoubtedly is seen as a powerful tool useful for retrieving the model from two or more sets of geophysical data, where each of these sets is represented with the same coordinates of vertices and the same number of parameters. Each of the different methods used in the joint inversion is sensitive to other physical attributes, so merging these results is especially valuable as it provides more reliable information about the predicted model. Joint inversion seems to be very promising, but it is still quite a serious challenge for today's geophysics. The complexity of the issue increases with the number of model parameters required (the more parameters, the more dimensions of solution space). Additionally, classic joint inversion suffers from the necessity of weighting and scaling target functions from incomparable methods. One of the possible strategies of solving this problem is the Pareto approach [2,4,6,10,12,15], however, it does not provide one sure answer about the final model, but rather kind of set of equally feasible solutions. What is worse, the same set of parameters can produce a few different models what raises the problem of equivalence. Taking all those facts into account, we decided to explore the problem of Pareto joint inversion effectiveness in solving the problems of scaling, weighting, and final solution interpretation using the simple dyke model.

Joint inversion
When two or more sets of geophysical data are available, using them together in an inversion procedure seems to be straightforward. Surprisingly even if the idea dates back to the beginning of the 19 th century, in geophysical application it was not proposed until 1975 [17]. At early applications data from similar methods were combined into one joint target function but it is more desirable when information comes from observation based on different physical processes. Unfortunately, simple combining misfits between observed and predicted data suffer from some serious disadvantages. First of all regular optimization procedures require one result returned by target function which in turn introduces the necessity of scaling misfits from individual methods. This operation is not obvious and usually requires arbitrary decisions. The second important problem is the correct uncertainty analysis for two combined weighted objective functions. One of the possible solution of all these problems is the Pareto inversion concept that was successfully introduced to geophysical analysis at the beginning of the 21st century [10].

Pareto approach
In the basic form of joint inversion information from different sources are combined into one target function. Separated target functions are misfits between modeled and observed data: (1) where: -measured data -model data 2 -variance.
Alternatively, multiobjective optimization theory -which has its roots in econometricsutilizes simultaneous analysis of many target functions and their reciprocal relations [4]. Joint inversion of more than one set of geophysical data entails scaling and weighting, as the results of each method differ from each other and are hard to compare. Unfortunately, such operations lead to significant loss of the information about multidimensional solution space. Moreover, it is a serious challenge to conduct uncertainty analysis for a few weighted objective functions correctly. Where more than one criterion comes into play it is not so obvious which solution is the best. In Pareto approach no parameter is more important than any other, so the inversion process implies accepting only these solutions, which produce better value for at least one target function and in the same time does not worse any other. The result of one single inversion run is a vector of values for each target function, which represents the coordinates of a point in the multi-dimensional solution space. Nonetheless, to get a reliable outcome, the inversion process is conducted many times and all results are composed into one set, where the optimal (non-dominated solutions) can be joined with a curve and make up so-called Pareto optimal front [6]. In contrast to the classic approach, there is no single best result obtained, but rather a set of all equivalent Pareto-optimal solutions, from which the best can be chosen taking into account qualitative -e.g. geological information. The multiobjective minimum can be described as: where: -target function, -space of acceptable solutions.
Solution vector x ∈ S dominates solution vector y ∈ S if and only if x provides better value in at least one dimension, not worsening any other [18]. Assuming that the function is minimized, x dominates y only if in at least one dimension x has a smaller value than y and at the same time x does not have bigger value than y in any other dimension: where ⊂ and , ∈ {1, . . . , },where nnumber of dimensions Set S* of nondominated solutions composes into Pareto front (P): In this contribution more than on Pareto front we focus on the set of all solutions fulfilling the Pareto criterion, what we call Pareto optimal set, as every obtained result is valuable and can bring different information about the retrieving model. In this case, the best misfit's values do not necessarily imply the feasible model shape. Obviously, some of these solutions could be improved by enlarging a number of allowed iterations or by additionally applying a local optimization method.

PSO -Particle Swarm Optimization
Particle Swarm Optimization (PSO) is a global optimization metaheuristic and stochastic method, coupling satisfactory precision with fast solution space searching. An algorithm was originally proposed by Kennedy and Eberhart in 1995 [9] and evolved with time, as it became very popular and gained the interest of scientists who worked over its improvement, presenting a large amount of its variations [5,13,14]. The algorithm was inspired by behavior of a flock of birds or shoal of fish looking for food. Initially, in PSO, the population of particles is randomly scattered in the n-dimensional solution space, where each dimension represents the respective geophysical parameter. Position of the particle, as a set of coordinates of all dimensions, strictly defines the proposed model, where the value of every dimension's coordinate relates to the value of a geophysical parameter. Further, the whole swarm moves in searching for the optimum of the function, iteratively changing each particle's position. In MARIA PSO was implemented in bare bones mode, where velocities were eliminated and new particle's positions were set according to the normal probability distribution of mean + 2 and standard deviation | − |, where and are the previous best positions in ℎ dimension, for an individual and for the neighbourhood (in our case of the whole swarm), respectively [8,14]. In other words, in every step particles move to the new position with some probability, which has its biggest value exactly in the middle of the way between actual particle and the particle representing the best solution found so far.

Software tools
To conduct all tests we used the software MARIA (Modular Approach Robust Inversion Algorithm), previously created as a part of former studies related to the Pareto joint inversion project [11]. MARIA was written in C language, under Fedora OS, using gsl scientific library, OpenMP for parallelization and GTK+ for Graphical User Interface. One of the priorities was to deliver a solution which would be flexible enough to handle data from two or more geophysical methods with the possibility to swap modules. Therefore, the first version [11] was combining magnetotellurics (MT) and gravity data, and as the investigations moved on, the MT module was swapped with magnetometry one. The main optimization engine was based on PSO algorithm in bare bones mode [8,9,14] and parallelized using OpenMP.

Methods
The gravity and magnetic methods are non-invasive, potential field methods which are commonly used for near-surface investigations. These methods are most useful if the physical parameters of search objects differ enough from the background regarding density and magnetic susceptibility [7]. The shape and amplitude for gravity anomaly are functions of density and depth of a body. For magnetic anomaly, the shape and amplitude are a function of magnetic susceptibility, depth of a body as well as magnetic field intensity and inclination angle in the area studied. In this paper, we focus on applying Pareto scheme for solving the problem of combined magnetic and gravity data inversion. Both methods are well known for problems related to multi-modality, solution equivalence, and low resolution. It is known that the Pareto approach can provide detailed insight into details of feasible solutions and support the decision of choosing the final one [2,6,10,12,15]. Herein, we analyze the other approach in which solutions from Pareto front are taken into consideration all at once. This can provide information about preserved model details that cannot be observed and interpreted when only one model is studied. Also, this approach can limit the adverse effects of multi-modality and equivalence. It is important to stress that obtained result cannot be interpreted as a probability distribution (like a result of MCMC methods) as it does not hold required assumptions, especially about symmetrical solution space sampling.
The forward solvers for gravity and magnetic methods use the superposition principle for gravity and magnetic effects respectively. Modeled gravity and magnetic effects are calculated as a sum from all nodes of the calculation grid. Each node is treated as a rectangle with density and magnetic susceptibility assigned. The magnetic effect for a single node is calculating using equation [16]: where: -magnetic susceptibility, -the intensity of the magnetic field [nT], -inclination angle [°], , -parameters defining the body's geometry (distance and angle).

Results and discussion
The application of Pareto joint inversion is demonstrated on two synthetic data sets. A dipping prism model which simulates magmatic intrusion (dyke) was used as the first synthetic example. In Fig. 1 and Fig. 5 models used for generating synthetic data are shown. Calculation of gravity and magnetic data was made for 200 evenly spaced points along the profile. White Gaussian noise was added to the generated data with different signal-to-noise ratio (50 dB for gravity and 10 dB for magnetics). The noise-contaminated data are shown in Fig. 3-4 and Fig.7-8. The physical properties of both "real" models are presented in Table 1. In calculated models, the inclination angle was assumed as 75° and intensity of magnetic field as 58 000 nT. These values are typical for Finland [1]. The starting models for the inversion in both examples are presented in Fig. 2 and Fig. 6. In both examples, 128 particles were used in PSO engine. All vertices of the models could change their position during the inversion process but their number is fixed. The physical properties-density and magnetic susceptibility could change their value in the range typical for real rocks. Both models represent a dyke built of igneous rocks, while the surrounding is formed by sedimentary rocks.  Fig. 2. Geometry for starting model in the first example.
The analysis of this simple model was aimed at checking the correctness of the algorithm's operation. In Fig. 3-4 fits of gravity and magnetic data for the simple dyke model are shown. The field data are represented by a blue curve, the model response before inversion by the green line, and example of well-fitted model after inversion by the red one. These results are a compromise between the two data sets. The blue curves represent data used as input for the inversion. In the beginning, the curves differ significantly. After the inversion process, it is possible to select a solution from the Pareto front that has a small misfit for the magnetometry and gravimetry data. This solution with almost perfect curve fitness produces a significantly different model from the real one.  The second example includes a more complicated model with small features (Fig. 5). The aim of inversion for this model was to check the possibility of recovering the shape of small details of a real model starting from the rectangular area. The denser mesh was used.      Many models generated during inversion return a similar fit. This is related to the equivalence for both gravity and magnetic methods. Due to the fact that it is not possible to choose the best model, another method was proposed to check the correctness of the algorithm's operation in order to recover the shape of the modeled body. Instead of examining single models, we analyzed the aggregated models in the form of one solution map. All Pareto solution models were analyzed. The number of points in the raster was counted for all 1131 runs. The results are presented as a heat map from the best 10 percent of solutions in Fig.11. In order to select the best solutions, on the basis of gravity and magnetic fitting error, the distance of points on the Pareto Optimal Set from the origin of coordinates was calculated.

Conclusions
Group analysis of models on the Pareto set gives information about the location of a dyke and some of its details. It can be concluded that even if the starting model is a rectangular area, it is possible to recover the shape of the real model if the set of solutions is analyzed. As expected the upper part of the model is well recovered whereas recognition of the deeper part is poor. This result corresponds with well-known feature of potential methods, namely the ability to recover more details in a shallow part of the analyzed medium.