Neural network stochastic simulation applied for quantifying uncertainties

Generally the geostatistical simulation methods are used to generate several realizations of physical properties in the sub-surface, these methods are based on the variogram analysis and limited to measures correlation between variables at two locations only. In this paper, we propose a simulation of properties based on supervised Neural network training at the existing drilling data set. The major advantage is that this method does not require a preliminary geostatistical study and takes into account several points. As a result, the geological information and the diverse geophysical data can be combined easily. To do this, we used a neural network with multi-layer perceptron architecture like feed-forward, then we used the back-propagation algorithm with conjugate gradient technique to minimize the error of the network output. The learning process can create links between different variables, this relationship can be used for interpolation of the properties on the one hand, or to generate several possible distribution of physical properties on the other hand, changing at each time and a random value of the input neurons, which was kept constant until the period of learning. This method was tested on real data to simulate multiple realizations of the density and the magnetic susceptibility in three-dimensions at the mining camp of Val d’Or, Québec (Canada). keyword: Artificial Neural Networks, Stochastic, Simulation, Geophysics, Density, Magnetic susceptibility


INTRODUCTION
During the last 20 years, predictions based on statistical learning processes had been applied with great success in a wide range of applications, including character recognition, shape classification or partitioning, filtering, predictions, etc.One of the most popular methods is *Corresponding author.E-mail: nfoudilbey@gmail.comartificial neural networks which mimic the neurons behavior of human brain.This work aims at applying these methods to quantify uncertainties of results obtained during geophysical inversion.
Geophysical data inversion is a key step in interpreting geological underground structures.Inversion methods are generally classified in two groups: (i) the deterministic approaches, and (ii) the stochastic ones.Deterministic inversions based on global minimum search offer a unique solution called the best model or the most likely model [1].They are based on the iterative least squares minimization of some objective function [2] often improved by introducing a weighting function [3].Another class of methods suggest solving the inverse problem by applying a neural network algorithm to faster the search of the optimum.They have been applied in geophysics for locating ground targets from electromagnetic field data [4], estimating basement resistivity distribution [5,6,7], seismic wavelength inversion [8].At the opposite of these last methods, stochastic methods such as Monte Carlo [9], genetic algorithms [10], and simulated annealing [11] offer the possibility of generating several models that can fit to the observations and therefore to estimate the a posteriori probability model.
Geostatistical simulations such as sequential Gaussian simulation (SGS) [12] or sequential indicator simulations (SIS) [13], are based on a variogram analysis.They simulate several equiprobable realizations of a random variable.The integration of geostatistical simulation methods during the inversion process increases the performance of the stochastic inversion as it accounts for the rock spatial variability [14,15].Furthermore, co-simulations simulate several possible realizations of a random variable accounting for secondary parameters, and generally give better results than the simulation of a single variate.Cosimulations integrate more information in the process of the inversion [16,17].In theory, it would be more interesting to integrate several geophysical properties to account for the geological complexity.However in practice, it is very difficult to take into account these secondary variables.In general, it is impossible to find deterministic relationships which combine all the physical properties, and requires the use of empirical relationships to link a property to another, or physical properties to the geology or to the lithology.Multiple point statistics integrating both the spatial variability and the multivariate aspects allow simulating variables accounting for several secondary parameters [18].
As suggested by [19], neural networks can be an alternative for simulating several variables.The advantages are: (i) to define statistical relationships either between two points of the spaces, or/ and between several points of the space variables, (ii) variables can be of various nature and (iii) to integrate non-linearity relationships between variables.[19] and [20] suggest using a 2D smooth gridded shape recognition method at the learning step, and then, a multi-grid approach [21,18] is used to simulate several realizations at the prediction step.In this work, we propose a different approach based on pattern recognition neural networks [20,22] for simulating the sub-surface physical properties.In this method, the analysis of physical properties is conditioned by an existing drilling data set and the results of independent inversions.Their distributions are used for a supervised training of the neural network, then the simulation is processed by injecting the distributions available at the prediction step.

NEURAL NETWORKS
The neural network is a linear combination of non-linear functions that tries to solve complex problems.Conceptually, this approach mimic the human brain activity.Their architectures are interconnected networks similarly to the synapse connections.They request a learning phase before being able to make forecasts [23,24,25].The number of input and output neurons correspond to the number of observed or predictive variables and of variable to be predicted, respectively.There are several types of neural networks, the most popular ones being the multi-layer perceptron (MLP) [26].The choice of the neural network architecture depends on phenomena and problems to be solved.In this work, we used a multi-layer perceptron (MLP) architecture with a feed-forward approach that is well adapted for spatial modeling.This approach has been successfully tested for predicting the metal content of ore deposits [27,28,29].A back-propagation method is implemented for optimizing weights at the learning stage using a mean squares criteria to minimize the errors between the observations and the predicted output, coupled to a conjugate gradient method [30] to find the optimal values for the connections.

NETWORK ARCHITECTURE
The input and output layers are of fixed size and are conditioned by the number of observed variables and by the number of variables to be predicted, respectively.Input variables are of various nature.It includes among others: coordinates (x, y, z), the magnetic susceptibility contrast χ 1 at the drilling level and χ 2 found by the inversion, the density contrast ρ at the drilling level found by the inversion, the coded geological unit κ for the different rocks types, a standard error σ = 0.5 for accounting for uncertainties of input variables [31]; the network outputs are the predicted susceptibility χ or the density ρ. Figure 1 shows the architecture used for the neural network.

TRAINING NETWORK AND BACK PROPAGATION
The neural network training process aims at calibrating the optimum weights of the different connections between neurons by minimizing the error between the predicted output and the observed values.In this work, we used a back-propagation algorithm [26] together with a conjugate gradient method.At the beginning of the learning stage, the weights are initialized by a random number generator with uniform distribution between −0.1 to 0.1 [32].The errors corresponding to the difference between the predicted and the observed output are then estimated at each step of the learning process.Then, the weights of each layer are recurrently modified during the learning stage using back-propagation techniques in order to reduce these errors.
The convergence rate of the algorithm depends on (i) the problem to be solved, (ii) the number of connections between the various entities and, of course, (iii) on the input and output layers sizes (based on the available data).The output level of a neuron X j is given by: (1) where f is the sigmoid function X i is the i th input value; w ij the synaptic weight between neurons i and j, N the total number of neurons and b j is the input bias often fixed to 0 during the learning stage.The total error ε between the output and the observations is evaluated as the mean squared error e i between the output of neurons i and the i th observation: (2) with e i = (d i -Y i ) being the difference between the observation d i and the prediction Y i at the neuron level i. Finally the error calculated at the level of a neuron is distributed on the related synapses and neurons by the conjugate gradient algorithm [33].
Let ∆w t ij be the correction to be made on weight w ij at iteration t connecting the neurons i to j, it is given by: (3) where w kj is an intermediate path between neuron i and neuron j passing through neuron k; µ and Ω are factor ranging between 0 and 1 called the rate and momentum, respectively.Ω is an added coefficient used to increase the convergence speed of the back-propagation algorithm, it accelerates the learning stage and forces to search minima out of local minima.

APPLICATIONS
A single case study was chosen to demonstrate the ability of the proposed method to integrate different types of data when interpolating and simulating multi variables.The study area is located in the Val d'Or mining camp in the Abitibi province, Québec, Canada (Figure 2).The geological formations are composed of volcanic rocks and volcano clastic Archaean mafic to ultramafic and felsic formations [34,35].These volcanic formations are separated by a narrow band of younger sedimentary units [36].The major part of the ore deposits are of Volcanogenic Massive Sulfide (VMS) type.During the learning stage, the neural network establishes non linear complex relationships between the different rock types defined in the input and the location of ore deposits in the output, especially regarding the spatial distribution of the susceptibility in the sub-surface.The neural weights obtained after the learning stage can then be studied carefully to better understand the link between each rock type and the known mineralizations, and thus can be used to derive 3D potential maps.

INTERPOLATION
Integrating several variables impacts the interpolation results given by the neural network.
In this example, only the magnetic susceptibility and the density contrast obtained from the magnetic and gravity inversions have been considered as the primary variable together with the rock types as secondary variables.inverted density contrast and magnetization, respectively.As seen in Figures 3 e) and f), properties interpolated using neural network account for rock types.Integrating the geological information at the prediction step had improved the property estimation especially in low contrast zones.These interpolations are conditioned by boreholes data as they have been included at the training stage.Predictions respect data conditioning as shown in Figures 3. d).

SIMULATION
The position of the selected points, rock types, magnetization, density contrast and neurons with a constant value of 0.5 for the learning process are taken as input to the neural network.
In order to simulate several rock types, magnetic susceptibility and density realizations after the training stage of the neural network, the prediction is performed by injecting an input data set and by randomly varying between 0 and 1 the neuron values of the last layer at each realization.Lower figures show the total magnetic field induced by the three selected realizations.The magnetic field produced by the selected models are very similar, the latter step is essential to validate the simulations, and thus for keeping the most realistic models that have similar potential field (gravity or magnetic) considered as a reference.

COMPUTER PERFORMANCES
The study area is modeled using a centered regular 3D grids (Voxet) comprising 150 ∞ 150 ∞ 75 elementary cells.The magnetic susceptibility and the density calculated from the inversion are stored at the center of each elementary cells.Ore deposits are retrieved using a regular 2D grid.The learning stage takes about 15 mn on the Intel Core 2 Duo (2.40 GHz).This is the longest step in implementing the suggested method.The prediction of a 3D cube is relatively fast, this prediction time depends on the number of simulations and the spatial resolution of the study area.

DISCUSSION AND CONCLUSIONS
This work proposes to use a stochastic simulation method based on training neural networks.The benefits of this method are: (i) to be available to integrate several variables, (ii) it does not request a preliminary geostatistical study such as the variograms analysis to capture the spatial variability, (iii) it does not need empirical relationships between the different input and output variables.The spatial correlations are studied capture during the training stage of the neural network.These same benefits could be viewed as disadvantages as neural networks behave like black box and need a set of training images including both the subsurface geophysical responses and a careful knowledge of the underground regarding ore deposits.Posterior analysis of the predictions allow to cross validate the results.However, other advantages of the proposed methodology are the ability of neural networks to be used not only as predictors but also as simulators.This is a very interesting tool to quantify the uncertainty related to a given exploration area.Once the neural network is completely trained, changing randomly weight neurons values of the output layer impacts the predicted results and allow generating several realizations which are conditioned by the drill holes data (considered as points).In some case, up-scaling procedure might be necessary depending on the cell sizes of the grid.Neural networks are well adapted to investigate mature brownfields in which several deposits are already been discovered at the learning stage.Assuming green fields behave like brown ones, it is then possible to apply this methodology to investigate new areas containing potential undiscovered deposits and to produce potential 3D map accounting for uncertainty.

Figure 1
Figure 1 Neural network architecture.
Figures 4.a .b and .cgive three selected possible realizations for the rock types.Figures 5 .a.b and .cshow three possible realizations of the sub surface magnetic susceptibility.The most upper figures show horizontal maps at a depth of 930 m, at the center, vertical cross-sections down to a depth of 1600 m.

Figure 4
Figure 4 Horizontal cross-sections showing three possible geological maps simulated by ANN at a depth of 930 m.