A partir d’echantillons releves sur le terrain, on cherche a determiner la concentration par exemple d’un metal dans toute la zone. On utilise pour cela des modeles de distribution de la quantite interessante.

Pour fixer les idees, notons z(s) la quantite de metal (par exemple en parties par million) en fonction de la position s (en dimension 2 ou 3). Les modeles etudies ici approchent z(s) par Z(s) avec:

Z(s) = m(s) + e(s)

ou m(s) est la tendance (trend en anglais) qui est deterministe (par exemple m(s) peut etre constant ou egal a la distance a un fleuve, etc.) et e(s) represente la partie aleatoire de la distribution.

On dispose de n echantillons, donc de n mesures z(s_{i}) (1 ≤ i ≤ n) et il s’agit de trouver m(s) et e(s) qui permettent de prevoir z(s) le plus precisement possible (krigeage)

1 Interpolation Une methode 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 carre:

z(s) = \frac{1}{C} ∑_{i=1}^{n} \frac{1}{|s − s_{i}|^{2}} z(s_{i}), C = ∑_{i=1}^{n} \frac{1}{|s − s_{i}|^{2}}

2 Approcher la tendance par des modeles lineaires au sens des moindres carres. On se donne p fonctions f_{j}(s) et on cherche des coefficients β_{j} tels que:

m(s) = ∑_{j=1}^{p} β_{j}f_{j}(s)

approche z(s) le mieux possible, “mieux” etant defini au sens des moindres carres, on veut minimiser (z(s) − m(s))^{2} aux points d’echantillonage. Le plus equilibre consiste donc a minimiser:

g(β_{1}, ..., β_{p}) = ∑_{i=1}^{n}(z(s_{i}) − m(s_{i}))^{2}

Pour trouver les β_{j}, on derive g par rapport a β_{1}, ..., β_{p}, on doit trouver 0 puisque la fonction possede un extremum en ce point.

Exemple si p = 1:

g(β_{1}) = ∑_{i=1}^{n}(z(s_{i}) − β_{1}f_{1}(s_{i}))^{2}

On derive par rapport a β_{1}:

g′(β_{1}) = ∑_{i=1}^{n} −2f_{1}(s_{i})(z(s_{i}) − β_{1}f_{1}(s_{i}))

et on cherche β_{1} tel que g′(β_{1}) = 0:

β_{1} = \frac{∑_{i=1}^{n} f_{1}(s_{i})z(s_{i})}{∑_{i=1}^{n} f_{1}(s_{i})^{2}}

Si p = 2, on a:

g(β_{1}, β_{2}) = ∑_{i=1}^{n}(z(s_{i}) − β_{1}f_{1}(s_{i}) − β_{2}f_{2}(s_{i}))^{2}

et en derivant par rapport a β_{1} et β_{2}, on obtient un systeme lineaire de deux equations a deux inconnues:

\{ \begin{matrix}∑_{i=1}^{n} −2f_{1}(s_{i})(z(s_{i}) − β_{1}f_{1}(s_{i}) − β_{2}f_{2}(s_{i})) & = & 0 \\ ∑_{i=1}^{n} −2f_{2}(s_{i})(z(s_{i}) − β_{1}f_{1}(s_{i}) − β_{2}f_{2}(s_{i})) & = & 0 \end{matrix}

donc:

( \begin{matrix}∑_{i} f_{1}(s_{i})^{2} & ∑_{i} f_{1}(s_{i})f_{2}(s_{i}) \\ ∑_{i} f_{1}(s_{i})f_{2}(s_{i}) & ∑_{i} f_{2}(s_{i})^{2} \end{matrix} ) ( \begin{matrix}β_{1} \\ β_{2} \end{matrix} ) = ( \begin{matrix}∑_{i} f_{1}(s_{i})z(s_{i}) \\ ∑_{i} f_{2}(s_{i})z(s_{i}) \end{matrix} )

Dans le cas general, on a interet a utiliser les notations de l’algebre lineaire (vecteurs et matrices). On note F la matrice a n lignes et p colonnes dont la j-ieme colonne est formee des valeurs de f_{j} aux points s_{1}, ..., s_{n} et on note F^{∗} sa transposee (c’est-a-dire la matrice obtenue par symetrie par rapport a la diagonale). On note z le vecteur colonne forme par la valeur de z(s) en s_{1},..., s_{n} et β le vecteur colonne forme par β_{1},..., β_{p}. On montre alors que β verifie l’equation:

F^{∗}F β = F^{∗}z

Pour la resoudre il faut donc inverser la matrice carree F^{∗}F. Remarques:

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.

Dans certaines situations, on met des poids differents w_{i} aux differents sites d’observations et on minimise plutot:

∑_{i=1}^{n} w_{i}(z(s_{i}) − m(s_{i}))^{2}

3 Modeliser les erreurs Dans une approche naive, on peut supposer que la partie aleatoire e(s) en un lieu est independante de e(s′) en un autre lieu, et que la distribution de ces variables aleatoires est identique, on parle alors de modele IID (independent identically distrbuted).

On peut aussi considerer que les parties aleatoires sont correlees entre elles et ce d’autant plus que la distance entre s et s′ est petite. Dans les modeles geostatistiques, on considere en fait que la correlation ne depend que de la distance entre s et s′ et on estime cette correlation a l’aide des echantillons s_{i} en calculant ou en modelisant le variogramme (ou le covariogramme) des s_{i}. Une fois ce calcul fait, on s’en sert pour ponderer le calcul de tendances au sens des moindres carres.

Covariogramme: On veut estimer la covariance des variables aleatoires e(s) et e(s′). On suppose donc qu’elle ne depend 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}| ∈ [h, h + δ[ ou δ est une constante, choisie de maniere judicieuse: δ doit etre petit pour que la fonction covariogramme varie peu sur l’intervalle [h, h + δ[, mais il faut aussi que cet intervalle contienne suffisamment de couples de l’echantillon pour que la moyenne soit representative.

Une autre fonction, appelee variogramme, est souvent utilisee pour representer la correlation 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}| ∈ [h, h + δ[.

Il existe differents modeles de fonctions variogramme et covariogramme: le modele de la pepite par exemple qui decorrele un echantillon du reste, des modeles continus (spheriques par exemple). On peut a partir du variogramme de l’echantillon determiner les coefficients d’un modele de variogramme approchant au mieux ce vraiogramme, par exemple au sens des moindres carres, cf. l’appendice A.1 de la documentation du logiciel gstat disponible sur:

http://www.gstat.org

Lorsqu’on prend comme modele des erreurs variables aleatoires e(s) telles que:

∀s, E(e(s)) = 0, Cov(e(s), e(s′)) = Cov(h)

(E(e(s)) designe l’esperance, c’est-a-dire la moyenne de e(s)), on adapte la methode d’estimation par les moindres carres a poids. On note V la matrice de coefficients:

v_{ij} = Cov(e(s_{i}), e(s_{j})) = Cov(|s_{i} − s_{j}|)

et z(s), m(s) les vecteurs:

z(s) = (z(s_{1}), ..., z(s_{n})), m(s) = (m(s_{1}), ..., m(s_{n}))

On minimise alors le produit scalaire:

< z(s) − m(s)|V^{−1}(z(s) − m(s)) >

par rapport aux coefficients β_{j} pour determiner ces coefficients β_{j} (et donc m). La partie erreur en un lieu s_{0} est estimee par:

Cov(e(s_{0}), e(s))V^{−1}(Z(s) − m(s))

Les calculs sont un peu plus compliques, la valeur optimale de β est alors:

β = (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 gstat (appendice A.2).

4 Programmation. Il existe des formats standard pour representer les donnees des echantillons. Par exemple le format .eas utilise par gstat, cf. sur pcm2 le fichier

/usr/local/gstat/cmd/zinc.eas:

ZincmeasurementsonRiverMeusefloodplains
3
xcoord,m
ycoord,m
zinc,ppm
1810723336111022
1810253335581141
181165333537640
...

La premiere ligne est un commentaire, la seconde indique le nombre de donnee de chaque echantillon (3 ici: 2 coordonnees x, y en metres, et la valeur de z en ppm, ce qu’expliquent les 3 lignes suivantes), puis chaque ligne correspond aux donnees d’un echantillon.

On lit les donnees dans le fichier et on les stocke dans un vecteur. L’element de base du vecteur sera compose de s et de z(s), par exemple ici s = (x, y):

structechantillon{doublex;doubley;doublez;};

(x et y sont les coordonnees, z represente z(s)). Pour lire les donnees a partir du fichier ci-dessus, par exemple:

...
ifstreaminfile("/usr/local/gstat/cmd/zinc.eas");
//onsauteles5premiereslignes
chars[256];
for(inti=0;i<5;i++)
infile.getline(s,255);
//onlitlesdonneesligneparligne
vector<echantillon>v;
echantillone;
for(;;){
infile>>e.x>>e.y>>e.z;
if(!infile)//silefichieresttermineonarretelaboucle
break;
v.push_back(e);
}
...

Les calculs presentes ne posent pas de difficultes algorithmiques sauf le calcul de l’inverse d’une matrice pour les modeles les plus sophistiques.

Calcul d’inverse de matrice: soit on dispose d’une librairie externe, soit on programme l’algorithme du pivot de Gauß(rappel de mathematiques: on colle la matrice identite a la matrice a inverser, on reduit, lorsqu’on a la matrice identite a gauche, alors la matrice a droite est l’inverse de la matrice de depart). Pour le choix du pivot a l’interieur d’une colonne donnee, on choisit le nombre de valeur absolue maximale afin de diminuer les erreurs d’arrondi.

Pour le calcul du covariogramme, on pourra creer un type:

structcouple{echantillonun;echantillondeux;doublelongueur;};

On cree un vecteur v de n^{2} couples (vector<couple>v;) contenant tous les couples d’echantillon (en calculant la longueur) et on ordonne ce vecteur par ordre croissant en utilisant la fonction de tri sort (ajoutez #include<algorithm> pour pouvoir l’utiliser). On definit la fonction de tri est_plus_pres:

boolest_plus_pres(couplea,coupleb){
returna.longueur<b.longueur;
}

On trie le vecteur par:

sort(v.begin(),v.end(),est_plus_pres);

Il est ensuite aise de determiner pour des valeurs de h et δ donnees les couples dont la longueur appartient a [h, h + δ[.

Remarque: la methode la plus efficace pour trouver le premier couple dont la longueur est superieure a h consiste a couper en deux moities egales le vecteur de couples. Si le couple du milieu a une longueur superieure a h, on recommence avec la premiere moitie du vecteur, sinon avec la deuxieme moitie, on s’arrete lorsqu’il n’y a plus qu’un couple dans le vecteur.