Géostatistique
ENSAI - CREST
https://marieetienne.github.io/statspat
2025-03-11
$$
$$
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 : \[\min \sum_{i,j=1}^n \lambda_i \lambda_j c_{ij} \quad ; \quad \sum_{i=1}^n \lambda_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)\)
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
\[\gamma(h) =\frac{1}{2} {\mathbb{V}\text{ar}}(Z(s+h) - Z(s)) = \frac{1}{2} {\mathbb{E}}( (Z(s+h) - Z(s))^2 )\]
On cherche donc à associer à chaque \(h\), \(0.5 {\mathbb{E}}( (Z(s+h) - Z(s))^2 )\), on va approcher cette espérance par \(0.5 \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élange, loi limite gaussienne, $\quad$ variance en $\frac{1}{n}$,
* ne nécessite pas l'estimation de la moyenne,
* si $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, \((\rho,\nu, C, \tau)\) 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éthodes possibles : max de vraisemblance, moindres carrés généralisés
Prédiction : le krigeage
Objectif
On observe \(z(s_1),z(s_2),\ldots,z(s_n)\)} aux sites \(s_1,s_2,\ldots,s_n\) et on veut prédire \(g(Z(s))\), \(s\in D\)
Par exemple
- \(Z(s_0)\), \(s_0\) non observé,
- \(\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
Par la propriété de l’espérance conditionnelle, le prédicteur optimal est
\[p^{opt}_{s_0}(Z(s_1),\ldots,Z(s_n)) = \mbox{E}\big(Z(s_0)/Z(s_1),\ldots,Z(s_n)\big) \] puisqu’il 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 universel
Hypothèses : \(Z(s) = X\beta + \epsilon(s)\),
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}\]
Covariance : \(C(h) = Cov(Z(s+h),Z(s))\), \(\boldsymbol{\Sigma}\) the covariance matrix (\(\Sigma_{ij}= C(s_j-s_i)\)), \(\boldsymbol{c_0}\) the covariance vector (\(c_j =C(s_o-s_j)\))
Proposition :
\[\widehat{Z}(s_0) = \left \lbrace \boldsymbol{c_0}^T \boldsymbol{\Sigma}^{-1} + (X_0 - X^T \boldsymbol{\Sigma}^{-1} \boldsymbol{c_0})^T (X^T \boldsymbol{\Sigma}^{-1} X)^{-1} X^T \boldsymbol{\Sigma}^{-1} \right\rbrace Z\] \[ \tau^2(s_0) = \sigma_0^2 - \mathbf{c_0}^T \mathbf{\Sigma}^{-1} \mathbf{c_0} + \left( X_0 - \mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{c_0} \right)^T \left( \mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{X} \right)^{-1} \left( X_0 - \mathbf{X}^T \mathbf{\Sigma}^{-1} \mathbf{c_0} \right)\]
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.
En pratique
Voir le TD