\documentclass{article} \begin{document} \'A partir d'\'echantillons relev\'es sur le terrain, on cherche \`a d\'eterminer la concentration par exemple d'un m\'etal dans toute la zone. On utilise pour cela des mod\`eles de distribution de la quantit\'e int\'eressante. Pour fixer les id\'ees, notons $z(s)$ la quantit\'e de m\'etal (par exemple en parties par million) en fonction de la position $s$ (en dimension 2 ou 3). Les mod\`eles \'etudi\'es ici approchent $z(s)$ par $Z(s)$ avec~: \[ Z(s) = m(s) + e(s) \] o\`u $m(s)$ est la tendance ({\em trend} en anglais) qui est d\'eterministe (par exemple $m(s)$ peut \^etre constant ou \'egal \`a la distance \`a un fleuve, etc.) et $e(s)$ repr\'esente la partie al\'eatoire de la distribution. On dispose de $n$ \'echantillons, donc de $n$ mesures $z(s_i)$ ($1\leq i \leq n$) et il s'agit de trouver $m(s)$ et $e(s)$ qui permettent de pr\'evoir $z(s)$ le plus pr\'ecis\'ement possible (krigeage) \section{Interpolation} \label{sec:inv_distance} Une m\'ethode d'estimation simple consiste \`a interpoler $z(s)$ par les $z(s_i)$ en mettant un poids proportionnel \`a l'inverse de la distance $|s-s_i|$ au carr\'e~: \[ z(s)= \frac{1}{C} \sum _{i=1}^n \frac{1}{|s-s_i|^2} z(s_i), \quad C=\sum _{i=1}^n \frac{1}{|s-s_i|^2} \] \section{Approcher la tendance par des mod\`eles lin\'eaires au sens des moindres carr\'es.} \label{sec:moindres_carres} On se donne $p$ fonctions $f_j(s)$ et on cherche des coefficients $\beta _j$ tels que~: \[ m(s)= \sum _{j=1}^p \beta _j f_j(s) \] approche $z(s)$ le mieux possible, ``mieux'' \'etant d\'efini au sens des moindres carr\'es, on veut minimiser $(z(s)-m(s))^2$ aux points d'\'echantillonage. Le plus \'equilibr\'e consiste donc \`a minimiser~: \[ g(\beta _1,...,\beta _p)=\sum _{i=1}^n (z(s_i)-m(s_i))^2 \] Pour trouver les $\beta _j$, on d\'erive $g$ par rapport \`a $\beta _1,...,\beta _p$, on doit trouver 0 puisque la fonction poss\`ede un extr\^emum en ce point. Exemple si $p=1$~: \[ g(\beta _1)=\sum _{i=1}^n (z(s_i)-\beta _1 f_1(s_i))^2 \] On d\'erive par rapport \`a $\beta _1 $~: \[ g'(\beta _1)=\sum _{i=1}^n -2 f_1(s_i) (z(s_i)-\beta _1 f_1(s_i)) \] et on cherche $\beta _1$ tel que $g'(\beta _1)=0$~: \[ \beta _1 = \frac{\sum _{i=1}^n f_1(s_i) z(s_i) }{\sum _{i=1}^n f_1(s_i)^2} \] Si $p=2$, on a~: \[ g(\beta _1,\beta _2)= \sum _{i=1}^n (z(s_i)-\beta _1 f_1(s_i)-\beta _2 f_2(s_i))^2 \] et en d\'erivant par rapport \`a $\beta _1$ et $\beta _2$, on obtient un syst\`eme lin\'eaire de deux \'equations \`a deux inconnues~: \[ \left\{ \begin{array}{rcl} \sum _{i=1}^n -2 f_1(s_i) (z(s_i)-\beta _1 f_1(s_i)-\beta _2 f_2(s_i)) & = & 0 \\ \sum _{i=1}^n -2 f_2(s_i) (z(s_i)-\beta _1 f_1(s_i)-\beta _2 f_2(s_i)) & = & 0 \end{array} \right. \] donc~: \[ \left(\begin{array}{cc} \sum_i f_1(s_i)^2 & \sum_i f_1(s_i) f_2(s_i) \\ \sum _i f_1(s_i) f_2(s_i) & \sum _i f_2(s_i)^2 \end{array}\right) \left(\begin{array}{c} \beta _1 \\ \beta _2 \end{array}\right)=\left(\begin{array}{c} \sum _i f_1(s_i) z(s_i) \\ \sum _i f_2(s_i) z(s_i) \end{array}\right) \] Dans le cas g\'en\'eral, on a int\'er\^et \`a utiliser les notations de l'alg\'ebre lin\'eaire (vecteurs et matrices). On note $F$ la matrice \`a $n$ lignes et $p$ colonnes dont la $j$-i\`eme colonne est form\'ee des valeurs de $f_j$ aux points $s_1$, ..., $s_n$ et on note $F^*$ sa transpos\'ee (c'est-\`a-dire la matrice obtenue par sym\'etrie par rapport \`a la diagonale). On note $z$ le vecteur colonne form\'e par la valeur de $z(s)$ en $s_1$,..., $s_n$ et $\beta $ le vecteur colonne form\'e par $\beta _1$,..., $\beta _p$. On montre alors que $\beta $ v\'erifie l'\'equation~: \[ F^* F \beta = F^* z \] Pour la r\'esoudre il faut donc inverser la matrice carr\'ee $F^* F$. Remarques~: \begin{itemize} \item On peut bien sur se limiter \`a faire des calculs avec les \'echantillons $s_i$ de la strate dans laquelle on veut faire des estimations si on connait l'existence de strates. \item Dans certaines situations, on met des poids diff\'erents $w_i$ aux diff\'erents sites d'observations et on minimise plutot~: \[ \sum _{i=1}^n w_i (z(s_i)-m(s_i))^2 \] \end{itemize} \section{Mod\'eliser les erreurs} \label{sec:variogrammes} Dans une approche naive, on peut supposer que la partie al\'eatoire $e(s)$ en un lieu est ind\'ependante de $e(s')$ en un autre lieu, et que la distribution de ces variables al\'eatoires est identique, on parle alors de mod\`ele IID ({\em independent identically distrbuted}). On peut aussi consid\'erer que les parties al\'eatoires sont corr\'el\'ees entre elles et ce d'autant plus que la distance entre $s$ et $s'$ est petite. Dans les mod\`eles g\'eostatistiques, on consid\`ere en fait que la corr\'elation ne d\'epend que de la distance entre $s$ et $s'$ et on estime cette corr\'elation \`a l'aide des \'echantillons $s_i$ en calculant ou en mod\'elisant le {\em variogramme\/} (ou le covariogramme) des $s_i$. Une fois ce calcul fait, on s'en sert pour pond\'erer le calcul de tendances au sens des moindres carr\'es. Covariogramme:\\ On veut estimer la covariance des variables al\'eatoires $e(s)$ et $e(s')$. On suppose donc qu'elle ne d\'epend que de $h=|s-s'|$, la distance de $s$ \`a $s'$, on la note Cov$(h)$. On fait alors la moyenne de $e(s_i) e(s_j)$ pour tous les couples $(i,j)$ de l'\'echantillonage tels que $|s_i-s_j| \in [h,h+\delta [$ o\`u $\delta $ est une constante, choisie de mani\`ere judicieuse~: $\delta $ doit \^etre petit pour que la fonction covariogramme varie peu sur l'intervalle $[h,h+\delta [$, mais il faut aussi que cet intervalle contienne suffisamment de couples de l'\'echantillon pour que la moyenne soit repr\'esentative. Une autre fonction, appel\'ee variogramme, est souvent utilis\'ee pour repr\'esenter la corr\'elation de $e(s)$ et $e(s')$, il s'agit de la moyenne de $\frac{1}{2} (e(s_i)-e(s_j))^2$ pour $|s_i-s_j| \in [h,h+\delta [$. Il existe diff\'erents mod\`eles de fonctions variogramme et covariogramme: le mod\`ele de la p\'epite par exemple qui d\'ecorr\`ele un \'echantillon du reste, des mod\`eles continus (sph\'eriques par exemple). On peut \`a partir du variogramme de l'\'echantillon d\'eterminer les coefficients d'un mod\`ele de variogramme approchant au mieux ce vraiogramme, par exemple au sens des moindres carr\'es, cf. l'appendice A.1 de la documentation du logiciel \verb|gstat| disponible sur~:\\ \verb|http://www.gstat.org| Lorsqu'on prend comme mod\`ele des erreurs variables al\'eatoires $e(s)$ telles que~: \[ \forall s, E(e(s))=0 , \quad \mbox{Cov}(e(s),e(s'))= \mbox{Cov}(h) \] ($E(e(s))$ d\'esigne l'esp\'erance, c'est-\`a-dire la moyenne de $e(s)$), on adapte la m\'ethode d'estimation par les moindres carr\'es \`a poids. On note $V$ la matrice de coefficients~: \[ v_{ij}=\mbox{Cov}(e(s_i),e(s_j))=\mbox{Cov}(|s_i-s_j|) \] et $z(s)$, $m(s)$ les vecteurs~: \[ z(s)=(z(s_1),...,z(s_n)), \quad m(s)=(m(s_1),...,m(s_n)) \] On minimise alors le produit scalaire~: \[ \] par rapport aux coefficients $\beta _j$ pour d\'eterminer ces coefficients $\beta _j$ (et donc $m$). La partie erreur en un lieu $s_0$ est estim\'ee par~: \[ \mbox{Cov}(e(s_0),e(s)) V^{-1} (Z(s)-m(s)) \] Les calculs sont un peu plus compliqu\'es, la valeur optimale de $\beta $ est alors~: \[ \beta =(F^* V^{-1} F)^{-1} F^* V^{-1} z(s) \] et on peut aussi calculer la variance de l'erreur $e(s_0)$, cf. la documentation logiciel {\tt gstat} (appendice A.2). \section{Programmation.} \label{sec:data} Il existe des formats standard pour repr\'esenter les donn\'ees des \'echantillons. Par exemple le format \verb|.eas| utilis\'e par \verb|gstat|, cf. sur \verb|pcm2| le fichier \verb|/usr/local/gstat/cmd/zinc.eas|~: \begin{verbatim} Zinc measurements on River Meuse flood plains 3 xcoord, m ycoord, m zinc, ppm 181072 333611 1022 181025 333558 1141 181165 333537 640 ... \end{verbatim} La premi\`ere ligne est un commentaire, la seconde indique le nombre de donn\'ee de chaque \'echantillon (3 ici: 2 coordonn\'ees $x$, $y$ en m\`etres, et la valeur de $z$ en ppm, ce qu'expliquent les 3 lignes suivantes), puis chaque ligne correspond aux donn\'ees d'un \'echantillon. On lit les donn\'ees dans le fichier et on les stocke dans un vecteur. L'\'el\'ement de base du vecteur sera compos\'e de $s$ et de $z(s)$, par exemple ici $s=(x,y)$~:\\ \verb|struct echantillon { double x; double y; double z; } ;|\\ (\verb|x| et \verb|y| sont les coordonn\'ees, \verb|z| repr\'esente $z(s)$). Pour lire les donn\'ees \`a partir du fichier ci-dessus, par exemple~: \begin{verbatim} ... ifstream infile("/usr/local/gstat/cmd/zinc.eas"); // on saute les 5 premieres lignes char s[256]; for (int i=0;i<5;i++) infile.getline(s,255); // on lit les donnees ligne par ligne vector v; echantillon e; for (;;){ infile >> e.x >> e.y >> e.z; if (!infile) // si le fichier est termine on arrete la boucle break; v.push_back(e); } ... \end{verbatim} Les calculs pr\'esent\'es ne posent pas de difficult\'es algorithmiques sauf le calcul de l'inverse d'une matrice pour les mod\`eles les plus sophistiqu\'es. Calcul d'inverse de matrice~: soit on dispose d'une librairie externe, soit on programme l'algorithme du pivot de Gau\ss(rappel de math\'ematiques: on colle la matrice identit\'e \`a la matrice \`a inverser, on r\'eduit, lorsqu'on a la matrice identit\'e \`a gauche, alors la matrice \`a droite est l'inverse de la matrice de d\'epart). Pour le choix du pivot \`a l'int\'erieur d'une colonne donn\'ee, on choisit le nombre de valeur absolue maximale afin de diminuer les erreurs d'arrondi. Pour le calcul du covariogramme, on pourra cr\'eer un type~:\\ \verb|struct couple { echantillon un; echantillon deux; double longueur;};|\\ On cr\'ee un vecteur \verb|v| de $n^2$ couples (\verb|vector v;|) contenant tous les couples d'\verb|echantillon| (en calculant la longueur) et on ordonne ce vecteur par ordre croissant en utilisant la fonction de tri \verb|sort| (ajoutez \verb|#include | pour pouvoir l'utiliser). On d\'efinit la fonction de tri \verb|est_plus_pres|~:\\ \verb|bool est_plus_pres(couple a,couple b){|\\ \verb| return a.longueur