Géostatistique
ENSAI - CREST
https://marieetienne.github.io/statspat
2025-01-15
$$
$$
Introduction
Sources
Ce cours est construit principalement à partir de celui de Lilane Bel et du livre de Gaetan, Guyon, et al. (2010)
Exemple illustratif
La température en France
Géostatistique
- Domaine continu : \(D \subset \mathbb{R}^2\).
- Processus stochastique : \(Z = \{Z(s)\}_{s \in D}\)
- Collection de variables aléatoires indexées par \(s \in D\).
- Données \(Z(s_1), Z(s_2), \ldots, Z(s_n)\) observés.
Objectifs :
- Prédiction \(\widehat{Z}(s_0)\) de \(Z(s_0)\) ou intégrale \(\int_B Z(s)ds\).
- Estimation de la loi de \(Z(s)\) ou d’une fonctionnelle \(\varphi(Z(s))\).
- Estimation des relations de dépendance entre \(Z(s_i)\).
Intuition
- Que veut on prendre en compte ?
- Quid des répétitions ?
- Que peut-on espérer estimer, comprendre ?
Estimation de la moyenne
Hypothèses
- \(E(Z(s)) = \mu \quad \forall s \in D\)
- \(\text{cov}(Z(s_i), Z(s_j)) = c_{ij}\)
Estimation
- Observations : \(Z(s_1), Z(s_2), \ldots, Z(s_n)\).
- Estimation de \(E(Z(s)) = \mu\) ?
Contraintes
- Estimateur linéaire : \(\widehat{\mu} = \sum_{i=1}^n \lambda_i Z(s_i)\).
- Sans biais : \(\text{E}(\widehat{\mu}) = \mu\).
- Variance minimale : \(\text{Var}(\widehat{\mu} - \mu)\).
Estimation de la moyenne
Si \(\widehat{\mu}\) est sans biais alors
Quelle est la conséquence de l’hypothèse de variance minimale ?
Solution : Best Linear Unbiased Estimator (BLUE)
- Minimisation : [ _{i,j=1}^n i j c{ij} ; {i=1}^n _i = 1 ]
Utiliser la méthode du multiplicateur de Lagrange et montrer qu’on cherche à résoudre le système linéaire
\[\begin{pmatrix} c_{11} & c_{12} & \dots & c_{1n} & 1 \\ c_{21} & c_{22} & \dots & c_{2n} & 1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ c_{n1} & c_{n2} & \dots & c_{nn} & 1 \\ 1 & 1 & \dots & 1 & 0 \end{pmatrix} \begin{pmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_n \\ m \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \end{pmatrix}\]
Conclusion si \(c_{ij}\) connus, on peut estimer la moyenne
Encore faut il pouvoir estimer la covariance
Spécificité des données spatiales
on a souvent une seule réalisation du processus spatial
Il faut avoir des hypothèses de régularité spatial pour espérer dire des choses à partir d’une réalisation
Processus du second ordre
Définitions
- \(Z\) est un processus du second ordre si \(\forall s \in D, \text{E}(Z(s)^2) < \infty\).
- La moyenne de \(Z\) est la fonction \[\begin{align} m : D & \to {\mathbb{R}}\\ s & \to m(s) = {\mathbb{E}}(Z(s))) \end{align}\]
- La covariance de \(Z\) est la fonction \[\begin{align} c : D\times D & \to {\mathbb{R}}\\ (s,t) & \to c(s,t) = {\mathbb{C}\text{ov}}(Z(s), Z(t)) \end{align}\]
Propriété Une fonction de covariance est semi définie positive i.e.
\[\forall a\in {\mathbb{R}}^n, \sum_{i,j}^n a_i a_j C(s_i, s_j)\geq 0\] Idée de preuve : Regarder \({\mathbb{V}\text{ar}}(\sum_{i=1}^n a_i Z(s_i))\)
- \(Z\) est un processus gaussien si pour toute partie finie \(S \subset D\) et toute suite réelle \(a = (a_s , s\in S)\), \(\sum_{s\in S} a_s Z(s)\) est une variable gaussienne.
Processus stationnaire
Définition
Un processus de second ordre \(Z\) est stationnaire sur \(D\subset {\mathbb{R}}^d\) si
- Moyenne constante : \(\forall s\in D, \quad m(s) = \mu\).
- Covariance invariante par translation : \[\forall (s, t)\in D^2 c(s,t) = C(t - s) \quad \text{ou} \forall h; s+h \in D, \quad c(s,s+h) = C(h) \]
Isotropie
Un processus est isotrope si \(C(h)\) dépend uniquement de \(\|h\|\) i.e \[C(h) = C(\|h\|)\]
Stationarité des accroissements
\(Z\) est un processus à accroissements stationnaires si les accroissements de \(Z\) sont stationnaires au second ordre, i.e. \[{\mathbb{E}}(Z(s + h) − Z(s)) = 0\] \[{\mathbb{V}\text{ar}}(Z(s + h) − Z(s)) = 2\gamma(h)\] La stationnarité du second ordre implique la stationnarité des accroissements.
La classe des processus à accroissement stations est donc plus ….. que la classe des processus de second ordre. ### Exemple
Que dire du mouvement Brownien sur \({\mathbb{R}}\)
Variogramme
Définition
Le semi-variogramme d’un processus à accroissements stattionnaire est défini comme : \[\gamma(h) = \frac{1}{2} {\mathbb{V}\text{ar}}(Z(s+h) - Z(s))\]
Propriétés
- \(\gamma(h) \geq 0\), \(\gamma(0) = 0\), \(\gamma(-h) = \gamma(h)\).
- Si \(Z\) est stationnaire de second ordre : \[\gamma(h) = C(0) - C(h)\]
- si \(lim_{|h|\to \infty} \gamma(h) = \ell < +\infty\) alors le processus est stationnaire du second ordre et \(\ell = C(0)\)
Variogramme
Variogrammes classiques
Pépitique : \(\gamma(h) = C\)
Exponentiel : \(\gamma(h) = C(1 - \exp(-\|h\|/\rho))\).
Sphérique : \[\gamma(h) = \begin{cases} C \left( \frac{3}{2} \frac{\|h\|}{\rho} - \frac{1}{2} \left(\frac{\|h\|}{\rho}\right)^3 \right) & \text{si } \|h\| \leq \rho \\ C & \text{si } \|h\| > \rho \end{cases}\]
Gaussien : \(\gamma(h) = C(1 - \exp(-\|h\|^2/\rho))\).
Puissance : \(\gamma(h) = C\left | h\right|^\alpha, \alpha < 2\)
Un variogramme particulier
Classe de Matèrn
\[ \gamma(h) = C\left [ 1 - \frac{1}{2^{\nu -1}\Gamma(\nu)}\left(\frac{\sqrt{2\nu }h}{\rho}\right)K_{\nu}\left(\frac{\sqrt{2\nu }h}{\rho}\right)\right]\] \(K_{\nu}\) fonction de Bessel modifiée de 3ème espèce, d’ordre \(\nu\)
\(\nu\) paramètre qui règle la régularité en 0.
- \(\nu\) = 1/2 : modèle exponentiel
- \(\nu\to \infty\)
Effet des paramètres
- palier \(C = \lim_{h \to \infty} \gamma(h)\)}
- échelle \(\rho\)
- portée \(\widetilde{\rho} \mbox{ t.q. } \|h\| \ge \widetilde{\rho} \rightarrow C - \gamma(h) \le \varepsilon\)
- pépite \(\tau = \lim_{h \to 0} \gamma(h)\)
Effet des paramètres
Le palier : S’il est fini, le processus \(Z\) est stationnaire du second ordre et \[ C = \mbox{Var}(Z(s))\]
L’échelle : Règle la vitesse à laquelle le variogramme rejoint le palier. Donne une idée de l’horizon de dépendance du processus.
La portée est la distance à partir de laquelle la corrélation est nulle ou négligeable
La pépite donne la régularité du processus
- si \(\tau\neq 0\) le processus est très irrégulier (bruit blanc, erreur de mesure)
- si \(\tau=0\), la régularité du variogramme en 0 (paramètre \(\nu\) dans le variogramme de Matern)
Estimation du variogramme
Intuition
- Quelles données pour le variogramme
- Comment peut-on s’y prendre
Nuée Variographique
$(h) = (Z(s+h) - Z(s)) = ( (Z(s+h) - Z(s))^2 )
On cherche donc à associer à chaque \(h\) \({\mathbb{E}}( (Z(s+h) - Z(s))^2 )\), on va approcher cette espérance par \(\left (Z(s_i)- Z(s_j) \right^2\)
Variogramme expérimental
\[\widetilde \gamma (d_k) = \frac{1}{2n_c}\sum_{(k-1)\delta \le \|s_i-s_j\| \le k\delta}(Z(s_i)-Z(s_j))^2 \]
Propriétés du variogramme expérimental
* sans biais,\\
* sous des conditions de m\'elange, loi limite gaussienne,
$\quad$ variance en \coul{$\frac{1}{n}$},\\
* ne n\'ecessite pas l'estimation de la moyenne, \\
* si \coul{$Z$} est gaussien, somme de $\chi^2(1)$\\Estimateur robuste
\[\widetilde \gamma (d_k) = \left\{\frac{1}{2n_c}\sum_{i,j \in C(k)}(Z(s_i)-Z(s_j))^2\right\}^{1/4}\frac{1}{0.457 + \frac{0.494}{n_c}}\]
Ajustement du variogramme
\(\gamma_{\rho,C,\nu,\tau}\)} variogramme admissible, solution du problème \[\min_{\rho,C,\nu,\tau} \sum_{k=1}^K (\gamma_{\rho,C,\nu,\tau}(d_k) - \widetilde \gamma(d_k))^2\] Autres m'ethodes :\ max de vraisemblance, moindres carrés généralisés
Prédiction : le krigeage
Objectifs
On observe aux sites \(s_1,s_2,\ldots,s_n\)
Prédire \(g(Z(s))\), \(s\in D\)
Par exemple
- \(Z(s_0)\), \(s_0\) non observ'e,
- \(\frac{1}{|B|}\int_B Z(s)ds \quad B \mbox{ région de } D\)},
- \(P(Z(s_0)>\zeta)\)
Prédicteur optimal
Prédire \(Z(s_0)\) sachant \(Z(s_1)=z(s_1),\ldots,Z(s_n)=z(s_n)\)
Prédicteur optimal \[p^{opt} = \mbox{E}\big(Z(s_0)/Z(s_1)=z(s_1),\ldots,Z(s_n)=z(s_n)\big)\] minimise le risque quadratique \[\mbox{E}\big((Z(s_0)-\widehat{Z}(s_0))^2\big)\] Prédicteur linéaire optimal (BLUP) \[ p^* = \alpha + \sum_{i=1}^n\lambda_i Z(s_i)\] Si \(Z\) gaussien, \(p^{opt}=p^*\)
Solution
\[\begin{align} \lambda &= C^{-1}c \qquad C_{i,j} = \mbox{cov}\big(Z(s_i),Z(s_j)\big) \quad c_i = \mbox{cov}\big(Z(s_0),Z(s_i)\big)\\ \alpha = & \mbox{E}(Z(s_0)) - \sum_{i=1}^n\lambda_i \mbox{E}(Z(s_i)) \end{align}\]
Krigeage simple
Supposons \(C\) et \(m(s)\) connues \[\widehat{Z(s_0)} = \; c^T C^{-1}(Z-m) + m(s_0)\] Variance de krigeage \[\sigma^2_{SK} = \sigma^2(s_0) - c^TC^{-1}c\] Si \(s_0=s_i\), alors \[\widehat{Z(s_i)}=Z(s_i)\]
Il faut estimer \(C\) et \(m(s)\).
Krigeage ordinaire
Hypothèses : \(Z(s) = \mu + \delta(s)\), à accroissements stationnaires.
On cherche \(\widehat{Z}(s_0) = \sum_{i=1}^n\lambda_iZ(s_i)\)} avec \[E(\widehat{Z}(s_0)) = E(Z(s_0)) \mbox{ et } E(\widehat{Z}(s_0) - Z(s_0))^2 \mbox{ minimum}\]
Variogramme : \(\gamma(h) = \frac{1}{2}\mbox{Var}(Z(s+h)-Z(s))\) \((\lambda_i)_{i=1,n}\)} solutions du système \[\begin{pmatrix} 0 & \gamma(s_1 - s_2) & \dots & \gamma(s_1 - s_n) & 1 \\ \gamma(s_1 - s_2) & 0 & \dots & \gamma(s_2 - s_n) & 1 \\ \cdots &\cdots &\cdots &\cdots &\cdots \\ \gamma(s_1 - s_n) & \gamma(s_1 - s_n) & \dots & 0 & 1 \\ 1 & 1 & \dots & 1 & 0\\ \end{pmatrix} \begin{pmatrix} \lambda_1\\ \lambda_2\\ \dots\\ \lambda_n\\ \alpha\\ \end{pmatrix}\] = \[\begin{pmatrix} \gamma(s_0 - s_1) \\ \gamma(s_0 - s_2) \\ \dots \\ \gamma(s_0 - s_n)\\ 1\\ \end{pmatrix}\]Krigeage Ordinaire
Variance de krigeage \[\sigma_K^2(s_0) = E\big((\widehat{Z}(s_0) - Z(s_0))^2\big) = \alpha + \sum_{i=1}^{n} \lambda_i \gamma(s_i-s_0)\]
Remarques
- On n’a pas forcément \(0 \le \lambda_i \le 1\)
- Si \(s_0 \in \{ s_1,\ldots s_n\}\) alors \(\lambda_i = 1\), \(\lambda_j=0\), \(j\neq i\)}, et \(\sigma_K^2(s_i)=0\)
- les poids de krigeage \(\lambda_i\) dépendent de la disposition des sites de données, de la localisation du site de prédiction, du nombre de données et de la fonction variogramme.