Simulation (géostatistique)

En géostatistique, les méthodes de simulation visent à proposer une variable régionalisée reproduisant un phénomène (ou processus) désiré. On parle de simulation conditionnelle lorsque les valeurs de la variable régionalisée en certains points sont définies.

En pratique, on veut souvent que la simulation respecte les deux premiers moments du processus, son histogramme et son variogramme.

Ces méthodes sont particulièrement employées en géostatistique non linéaire, comme souvent le seul moyen techniquement disponible pour l'estimation de grandeurs. En effet, les méthodes de modélisation géostatistique donnent la meilleure estimation de la variable régionalisée, mais lissent le résultat et donc échouent à reproduire une variabilité naturelle du phénomène.

Notations

Dans la suite, les notations suivantes seront utilisées :

• N et n respectivement le nombre de points de données et le nombre de points d'estimation
• i l'indice allant de 1 à N pour les points de données, de N+1 à N+n pour les points d'estimation
• x_i la position du ième point
• z_i la valeur de la variable régionalisée au ième point
• Z la fonction aléatoire du phénomène.
• C la covariance du phénomène simulé.

Pour un vecteur v donné, on note v_|k le vecteur des k premiers éléments de v. Par exemple, z_|N est le vecteur des valeurs connues aux points x_|N. De même, pour une matrice M, M_|k note sa sous-matrice aux k premières lignes et k premières colonnes.

Pré-calculs

Aplanissement

Cette étape est un préalable possible à la modélisation à l'intérieur d'une couche géologique. Plutôt que de travailler dans un repère de coordonnées XYZ, l'on passera à un repère XYh où h est la distance selon Z entre le point courant et une surface de référence (par exemple, la base de la couche). Ceci vise à s'approcher d'un travail en lignes de niveau qui correspondraient aux isochrones de dépôt.

Anamorphose gaussienne

L'anamorphose gaussienne consiste à appliquer une bijection à la variable pour lui donner une distribution gaussienne. En effet, plusieurs méthodes de simulation exigent ce préalable. La bijection inverse devra être appliquée sur le résultat.

Méthodes de simulations

Simulation matricielle

Cette méthode utilise la décomposition de Cholesky et se prête bien aux simulations en une dimension en voisinage globale. Soit (C_i,j) la matrice de covariance, des éléments C_i,j = C(x_j−x_j) que l'on décompose par la méthode de Cholesky en C = LL^T.

Dans le cas non-conditionnel, L est d'ordre n. On tire n valeurs aléatoirement selon une loi normale centrée réduite, x₁, …, _n. Une réalisation est alors z = Lx.

Dans le cas conditionnel, L est d'ordre N+n. Le conditionnement est fixé par y_|N = L_|N⁻¹z_|N. On tire n valeurs y_N+1,…,y_N+nselon une loi normale centrée réduite. Une réalisation est alors z = Lx.

Cette technique permet de réaliser un grand nombre de simulations sans grand allongement des calculs, une seule décomposition de Cholesky est nécessaire.

Simulation gaussienne séquentielle

Cette méthode exige que la fonction aléatoire soit gaussienne ; alors, une distribution conditionnelle est également une gaussienne, dont les espérance et variance se déduisent d'un krigeage simple.

Elle s'effectue par étapes, l'ordre de visite des n points d'estimation pouvant en pratique influer sur le résultat. À l'étape k+1un krigeage simple sur les N points de données et les k points déjà simulés donne une valeur krigée z^* et un écart-type de krigeage σ^*. On tire la valeur de z_N+k+1 selon une loi normale d'espérance z^* et de variance σ^*2.

Bandes tournantes

L'idée est de transformer la simulation sur une partie de ℝ³ ou ℝ² en composée de simulations sur des parties de ℝ. On trace dans l'espace une série de bandes (lignes) S_i, sur chacune est calculée une réalisation Y_i(x) du processus. La valeur en un point x quelconque est une somme des valeurs aux projetés sur les bandes, affectée d'un facteur. $Z\left(x\right)={\frac {1}{\sqrt {n}}}\sum _{i=1}^{n}Y_{i}\left(\langle x.S_{i}\rangle \right)$

Cette méthode requiert typiquement quelques centaines de bandes tournantes pour que s'effacent les artéfacts de calcul.

Post-conditionnement

Soit une méthode de simulation non-conditionnelle. On suppose la fonction gaussienne centrée (pas forcément réduite).

Résidu de rubanage

Le conditionnement peut s'effectuer selon l'algorithme suivant :

• tirage d'une simulation conditionnelle en tout point x_i de données ou de simulation. On obtient respectivement z_D et z_S.
• krigeage (simple ou ordinaire) aux points de simulation en utilisant les valeurs observées z_|N. On obtient n valeurs z_d*
• krigeage (idem) aux points de simulation en utilisant les valeurs simulées aux points de données z_D. On obtient z_DS.
• La simulation conditionnelle est donnée par z_d*+z_S−z_DS.

Convergence

Conditionnellement aux données, l'espérance d'un grand nombre de simulations tend en tout point vers l'estimation par krigeage ; de même, leur variance tend vers la variance de krigeage.