INTRODUCTION À LA

GÉOSTATISTIQUE APPLIQUÉE

Enseignant : Dr. BOUKA Stéphane, Maître Assistant CAMES

Département de Mathématiques et Informatique de la Faculté des Sciences

Université des Sciences et Techniques de Masuku, Franceville, Gabon

1
Ob jectif

La géostatistique est une branche de la statistique adaptée à l'estimation spatiale de

propriétés du milieu physique. Elle traite les propriétés observées de façon discontinue dans
l'espace géographique (en un point, sur une petite surface).

Nous allons essayer au cours de ce cours introductif de présenter les principes de mise en
÷uvre d'une analyse géostatistique. Le TP qui prolonge ce cours, a deux objectifs : d'une
part, montrer la mise en ÷uvre pratique de la géostatistique, d'autre part, montrer que la
géostatistique se fonde sur les bases des statistiques classiques mieux connue et maitrisées
par les étudiants.

Pour un étudiant confronté à un problème de variabilité spatiale, le premier choix qu'il aura
à faire concerne le type d'approche qu'il met en ÷uvre. Deux grandes voies lui sont ouvertes :
 Employer une démarche d'interpolation déterministe. C'est le type d'approche utilisée
depuis bien longtemps et même encore de nos jours. Les observations sont implantées
dès lors que les caractéristiques du paysage changent. Les limites sont tracées en
s'appuyant sur les modications du paysage. Cette technique est souvent économe
en moyens et présente des résultats très parlants. Elle ne permet par contre jamais
d'obtenir une carte dont on connaît la précision.
 Employer une démarche statistique, basée sur les statistiques classiques (recherche de
moyennes, de variances au sein d'une population ou de strates) ou sur la géostatistique
(obtention de cartes). On reproche souvent à ces techniques leurs exigences élevées en
matière d'échantillonnage. Il demeure que ces techniques sont incontournables dès lors
que l'on désire obtenir des estimations dont on connaît la précision. Ces techniques sont
également les seules que l'on puisse mettre en ÷uvre dans certains cas : phénomène
naturelle, propriété dont la variabilité ne dépend que de l'action de l'homme (pollution)
...

Les querelles d'école entre ces deux types d'approches restent nombreuses. On peut proposer
une approche pragmatique pour faire le choix. On envisagera l'approche déterministe quand
le paysage est très contrasté et que l'on sait que ces contrastes correspondent à des états
diérents de la propriété étudiée. On préférera l'approche statistique quand il est utile d'avoir
des estimations de précision connue ou que le paysage varie peu.

C'est dans cet objectif que ce cours est inscrit. Il vise à introduire les concepts de la géosta-
tistique.

2
Chapitre 1

LA VARIOGRAPHIE

1.1 Introduction à la notion de variogramme

Considérons une propriété notée  Y  connue en  n  points de l'espace géographique,

chacun de ces points étant repérés par le vecteur  x  de ses coordonnées géographiques
(longitude et latitude). De la sorte, la notation "Y (xi )" représente la valeur observée de la
propriété Y au i ème point d'échantillonnage de coordonnées  xi .
Pour simplier, prenons deux points pour lesquels on connaît des valeurs Y (x1 ) et Y (x2 )
de la propriété Y dans un espace géographique tel que le montre la gure ci-dessous.

Pour comparer ces deux valeurs, la façon la plus simple est d'utiliser la variance entre les
observations de ces deux sites, notée  S 2 . Elle est par dénition égale à :

2  2
S 2 = Y (x1 ) − Y + Y (x2 ) − Y , (1.1)


où : Y est la moyenne entre ces deux observations.

Cette variance  S 2 , qui traduit l'importance des écarts à la moyenne, est d'autant plus
grande que les observations sont diérentes et, au contraire, si elle est faible les observations
sont de plus en plus identiques. L'équation 1.1 peut être développée pour obtenir une autre
expression de la valeur S 2 :
1
S2 = [Y (x1 ) − Y (x2 )]2 , (1.2)
2
Cette nouvelle équation pour déterminer la variance (eq.1.2) peut être écrite pour tout couple
de sites. Pour cela, considérons deux sites xi et xi + h où xi représente les coordonnées
géographiques d'un des sites et  h  est un vecteur caractérisant la distance entre les sites.
L'équation (1.2) s'écrit alors :
1
S2 = [Y (xi ) − Y (xi + h)]2 , (1.3)
2

3
Calculons à présent la distance géographique séparant y(x1 ), y(x2 ), y(x3 ), y(x4 ) chacun des
points d'observation et considérons les  m  couples de point séparés par une même distance
géographique h.

On peut comme précédemment, calculer la variance des observations pour les sites pris deux
à deux. La moyenne S 2 de ces m variances s'écrit en employant (1.3) :
m
1 X
2
S = [Y (xi ) − Y (xi + h)]2 , (1.4)
2m i=1
Où : m est le nombre de couple.
Pour une distance h séparant deux points d'observation, S 2 rend compte de la ressemblance
et/ou la dissemblance des observations faites en ces deux points : il sera d'autant plus grand
que ces observations sont diérentes et le contraire signie une grande ressemblance entre les
observations. S 2 est qualiée de "semi-variance".
De manière intuitive, on conçoit que deux observations soient en général d'autant plus sem-
blables qu'elles sont proches géographiquement l'une de l'autre. Le calcul de S 2 pour dié-
rentes distances h, va permettre de quantier cette idée : il permet de suivre l'évolution des
écarts entre des observations en fonction de la distance qui les sépare.
Mathéron (1965) a montré l'intérêt de cette notion simple et les conditions de généralisation
ont été dénies par la théorie qu'il a appelée théorie des variables régionalisées. Cette
théorie montre que la généralisation de l'équation (1.4) suppose deux conditions, regroupées
sous le terme d'hypothèse intrinsèque et qui sont :
 L'espérance de Y est constante quelle que soit la position géographique x :
E[Y (x)] = m (constante) (1.5)
 Pour toute distance h, la diérence  [Y (x) − Y (x + h)]  a une variance nie, qui ne
dépend que de la distance h séparant les points.
V AR[Y (x + h) − Y (x)] = 2λ(h) (1.6)
= E[Y (x + h) − Y (x)]2
Quand ces deux conditions sont vériées, la valeur S 2 dénie dans l'équation (1.4) consti-
tue un estimateur non biaisé de la fonction λ(h) dénie en équation (1.6). Cette fonction
λ(h) est nommée  variogramme .
m(h)
1 X
γ(h) = [Y (xi ) − Y (xi + h)]2 , (1.7)
2m(h) i=2

4
Où : m est le nombre de couple.

Intérêt du variogramme : En étudiant l'évolution du variogramme λ(h) en fonction de

la distance h séparant des couples d'observation, on va analyser la façon dont se détériore
l'information acquise en un point au fur et à mesure que l'on s'éloigne de ce point.

1.2 Le calcul du variogramme

On cherche à construire un graphique représentant en abscisse les distances h séparant

les points et en ordonnée les semi-variances [γ(h)].
La construction du variogramme est illustrée ci-dessous par des schémas établis à partir
de 8 points d'observation répartis à distance égale de 1 mètre le long d'un transect.

Figure 1.1  Illustration du calcul du varriogramme sur un transect de 8 points séparés par
une distance h = 1m
Le schéma ci-dessus (Figure 1.1) montre que le nombre de points participant au cal-
cul du variogramme diminue au fur et à mesure que la distance augmente. Les valeurs de
semivariance risquent donc d'être moins précises pour les grandes valeurs de h.

Exemple de calcul
Soit deux exemples (série A et série B) ctifs correspondant à des observations disposées le
long d'un transect à des intervalles réguliers de 1 mètre.
 Calculez pour chacun des exemples : la moyenne, la variance, γ(1), γ(2), γ(3) et γ(4).
 Que peut-on conclure ?

Solution
- La moyenne A
1X 1
XA = Xi = [4 + 3 + 2 + 1 + 0 + 1 + 2 + 3 + 4] = 2.22
n i 9
- La moyenne B

5
 1X 1
XB = Xi = [4 + 2 + 1 + 0 + 3 + 1 + 2 + 4 + 3] = 2.22
n i 9
- La variance A
1 X
SA2 = (Xi − X A )2
n−1 i
1
= [(4 − 2.22)2 + (3 − 2.22)2 + (2 − 2.22)2 + (1 − 2.22)2 + (0 − 2.22)2
8
+(1 − 2.22)2 + (2 − 2.22)2 + (3 − 2.22)2 + (4 − 2.22)2 ]
= 1.94

- La variance B
1 X
SB2 = (Xi − X B )2
n−1 i
1
= [(4 − 2.22)2 + (2 − 2.22)2 + (1 − 2.22)2 + (0 − 2.22)2 + (3 − 2.22)2
8
+(1 − 2.22)2 + (2 − 2.22)2 + (4 − 2.22)2 + (3 − 2.22)2 ]
= 1.94

La gure ci-dessus illustre la méthodologie de calcul du variogramme de l'exemple des série

A et B.

Figure 1.2  Illustration de calcul des valeurs de γ(h) de la série A et B.

- Le calcul des valeurs de λ(h) de la séries A

1
γ(1) = [(4−3)2 +(3−2)2 +(2−1)2 +(1−0)2 +(0−1)2 +(1−2)2 +(2−3)2 +(3−4)2 ] = 0.5,
2∗8
1
γ(2) = [(4 − 2)2 + (3 − 1)2 + (2 − 0)2 + (1 − 1)2 + (0 − 2)2 + (1 − 3)2 + (2 − 4)2 ] = 1.71,
2∗7

6
 γ(3) = [(4 − 1)2 + (3 − 0)2 + (2 − 1)2 + (1 − 2)2 + (0 − 3)2 + (1 − 4)2 ]/(2 ∗ 6) = 3.16,
γ(4) = [(4 − 0)2 + (3 − 1)2 + (2 − 2)2 + (1 − 3)2 + (0 − 4)2 ]/(2 ∗ 5) = 4.

- Le calcul des valeurs de λ(h) de la séries B

1
γ(1) = [(4−2)2 +(2−1)2 +(1−0)2 +(0−3)2 +(3−1)2 +(1−2)2 +(2−4)2 +(4−3)2 ] = 1.56,
2∗8
1
γ(2) = [(4 − 2)2 + (3 − 1)2 + (2 − 0)2 + (1 − 1)2 + (0 − 2)2 + (1 − 3)2 + (2 − 4)2 ] = 1.92,
2∗7
γ(3) = [(4 − 0)2 + (2 − 3)2 + (1 − 1)2 + (0 − 2)2 + (3 − 4)2 + (1 − 3)2 ]/(2 ∗ 6) = 3.08,
γ(4) = [(4 − 3)2 + (2 − 1)2 + (1 − 2)2 + (0 − 4)2 + (3 − 3)2 ]/(2 ∗ 5) = 1.9.

- Conclusions ?
ˆ Les deux séries ont même moyenne et même variance, toutefois on constate clairement
qu'elles n'ont pas le même degré de continuité spatiale, la première série (série A) étant
nettement plus continue que la seconde (série B).

Remarque
On peut aussi calculer le variogramme selon certaines directions spéciques ; pour cela le
variogramme est dénie par son pas de calcul  h  et sa direction déterminée par un angle
(θ). On parle dans ce cas du variogramme directionnel et l'équation (1.7) s'écrit :
m(h,θ)
1
(1.8)
X
γ(h, θ) = [Y (xi ) − Y (xi + h)]2 .
2m(h, θ) i=2

Exemple de calcul de variogramme directionel

Soit une matrice de données 3 x 3 ayant les valeurs suivantes : la distance horizontale et
verticale entre 2 éléments consécutifs est de 1 m et  N  indique une donnée manquante
(voir schéma ci-dessous).
- Le calcul du variogramme dans la direction horizontale : θ = 0°

7
 γ(1, 0°) = [(3 − 6)2 + (6 − 5)2 + (7 − 2)2 + (2 − 2)2 ]/(2 ∗ 4) = 4.4

γ(2, 0°) = [(3 − 5)2 + (7 − 2)2 + (4 − 0)2 ]/(2 ∗ 3) = 7.5

-Le calcul du variogramme dans la direction verticale : θ = 90°

γ(1, 90°) = [(4 − 7)2 + (7 − 3)2 + (2 − 6)2 + (0 − 2)2 + (2 − 5)2 ]/(2 ∗ 5) = 5.4,

γ(2, 90°) = [(4 − 3)2 + (0 − 5)2 ]/(2 ∗ 2) = 6.5.

1.3 Variogramme expérimental

Figure 1.3  Un exemple de variogramme expérimental : variogramme moyen du pH du sol

d'une parcelle expérimentale de 1,5 hectare (Douaoui, 1993).

La gure ci-dessus représente, à titre d'exemple, un variogramme sur des données de mesure
du pH du sol d'une parcelle expérimentale de 1,5 hectares , 150 mesures ont été eectuées

8
suivant un échantillonnage régulier (10m x 10m). Sur cette gure, on représente en abscisse
diérentes distances séparant des couples de points expérimentaux : ces distances sont nom-
mées "pas" (lag en anglais). En ordonnée, on représente les valeurs des semivariances ou γ(h)
calculées suivant l'équation (1.7).
Ce qu'on peut constater sur cette gure :
(i) Jusqu'à un pas de 20 mètres, le variogramme est croissant. Les écarts moyens entre
les observations augmentent donc quand la distance séparant ces observations augmente. Les
observations "se ressemblent donc de moins en moins", ce qui est conforme à l'intuition. On
peut dire également que les observations sont spatialement dépendantes ou liées sur cette
distance de 20m.
(ii) Au-delà de vingt mètres, le variogramme reste quasi-constant. Quelle que soit la
distance, les écarts moyens entre les observations sont identiques. On parlera pour cette
gamme de distance d'indépendance spatiale entre les observations.
(iii) La projection du variogramme à l'origine conduit à une valeur de semi-variance non
nulle bien que la distance est nulle (h = 0).

Des dénitions.
Un certain nombre de termes sont utilisés pour décrire un variogramme de la gure 1.3 et
qui sont :
a. L'eet de pépite (nugget eect) : il s'agit da la valeur de la semi-variance pour une
distance nulle. En théorie, on devrait avoir un λ(h) = 0 pour un h = 0, mais fréquemment,
le variogramme présente une ordonnée à l'origine non nulle (g. 1.3). Cet écart est qualié
"d'eet de pépite" (nugget eect en anglais). Il est interprété comme le résultat d'erreurs
de mesure de la variable étudiée, ou erreur de positionnement ou d'une variabilité spatiale
présente à une distance inférieure au pas d'échantillonnage.
b. Le palier (Sill) : valeur de la semi-variance à partir de laquelle le variogramme ne croît
plus (g. 1.3).
c. La portée (Range) : distance à partir de laquelle le palier est atteint (g. 1.3). La
portée est la distance à partir de laquelle les valeurs de la variable entre deux points sont
indépendantes (non corrélées) (g. 1.3).

Figure 1.4  Le variogramme expérimental

9
1.4 Les modèles du variogramme

L'analyse du variogramme cherche une fonction caractéristique de la structure de la va-

riable étudiée. En premier lieu, on étudie quelques caractéristiques du variogramme :

1.4.1 Comportement au voisinage de l'origine

La continuité et la régularité dans l'espace de la fonction aléatoire et donc la variable

régionalisée qu'elle présente sont liées au comportement à l'origine du variogramme
Delhomme (1976) distingue 04 types
1. Allure parabolique : comportement dérivable à l'origine, ceci est la caractéristique d'une
variabilité spatiale hautement régulière (g.1.5.a)
2. Allure linéaire : γ(h) reste continue à l'origine mais n'est plus dérivable, donc moins
régulière (g.1.5.b).
3. Discontinuité à l'origine : γ(h) ne tend pas vers (0) lorsque h tend vers (0), cette
discontinuité en h = 0 du variogramme est appelée eet de pépite (g.1.5.c) qui est dû : soit
à la présence d'une structure dont l'échelle est très inférieure à l'espacement des données et
on parle de micro régionalisation des données, soit à la présence d'erreurs de mesures, soit
au nombre insusant de couples de mesures à faible distance induisant éventuellement une
incertitude sur la détermination de l'eet de pépite.
4. Eet de pépite pur (Aléatoire pure) : c'est le cas limite du cas précèdent quand γ(h)
ne traduit plus que la seule discontinuité à l'origine (g.1.5.d) γ(0) = 0 et γ(h) = C0 dès que
h > 0. Cela indique que Y (x) et Y (x+h) sont sans corrélation quelle que soit leur distance (h)
non nulle, ce type de modèle s'explique généralement par l'absence d'une structure spatiale,
plus fréquemment, par l'existence d'une structure marquée par des erreurs expérimentales ou
inférieures au plus petit intervalle d'observation.

Figure 1.5  Comportement à l'origine des diérents variogrammes (Delhomme, 1976)

1.5 Modélisation du variogramme

Pour tenir compte des caractéristiques du variogramme dans la démarche géostatistique,

il est indispensable d'ajuster une fonction au variogramme expérimental, ce qui permet d'en
résumer les principales caractéristiques. Ces fonctions doivent présenter deux qualités :

10
  Rendre compte le mieux possible de l'information du variogramme expérimental.
 Satisfaire les conditions théoriques : elles doivent être "Semi-Positives".
L'ajustement se fait par l'emploi d'un certain nombre de modèles autorisant essentiellement
deux types :

1.5.1 Modèles croissants non bornés

Les modèles non bornés sont dénis comme suit :

* Modèle linéaire :

γ(h) = C0 + bh (1.9)

Avec : C0 : L'ordre à l'origine, b : la pente de la droite, h : distance séparant les points

- cas particulier du linéaire :

γ(h) = C0 (1.10)

C'est le cas d'un variogramme plat appelé pépidique (eet de pépite pûre)

* fonction puissance :

γ(h) = C0 + bhα (1.11)

0 < α < 2, avec : b : la pente de la droite, α : un coecient xant la forme de la courbe.

1.5.2 Modèles croissants bornés

Les modèles croissants bornés sont dénis comme suit :

* Modèle Sphérique :
(
C0 + C[3h/2a − 1/2(h/a)3 ] si h < a
γ(h) = , (1.12)
C0 + C si h > a

avec : C0 : l'ordre à l'origine, C : est le palier moins l'ordonnée à l'origine, a : la portée

* Modèle Exponentiel :

γ(h) = C0 + C[1 − exp(−h/r)], (1.13)

avec : C0 : l'ordre à l'origine, r : paramètre de la distance égale environ le tier (1/3) de la

portée.

11
1.6 Stratégie pour le calcul de variogrammes et l'a juste-

ment des modèles

Pour le calcul et l'ajustement des variogrammes il faut tenir compte des points suivants :
 On accorde plus de poids aux points du variogramme expérimental calculés avec beau-
coup de paires.
 On essaie d'avoir un nombre de couple supérieur à 30 minimum (idéal 50 couples) pour
chaque point expérimental du variogramme. Si ce n'est pas possible pour certaines
classes, on accorde moins d'importance à ces points. Si le nombre de paires est très
faible, on ne considère plus du tout le point.
 On accorde plus de poids aux premiers points du variogramme (h petit) car ce sont
ces valeurs qui ont le plus d'impact dans les calculs géostatistiques.
 Lorsque  h  dépasse environ dmax/2, on ne tient pas compte des valeurs du vario-
gramme. (dmax est la taille du phénomène étudié dans la direction considérée).
 On cherche à obtenir des modèles les plus simples possible qui rendent bien compte
des valeurs expérimentales.

12
Chapitre 2

ESTIMATION D'UNE TENEUR

PONCTUELLE (KRIGEAGE)

2.1 Introduction

Le krigeage, le deuxième outil de la géostatistique, est une méthode d'interpolation ap-

plicable à des données spatiales et qui consiste à estimer la teneur ponctuelle de la variable
étudiée en des sites non échantillonnés. La théorie du krigeage a été développée par un ma-
thématicien français G. Matheron, école des Mines de Paris au début des années 1960, à
partir des travaux de l'ingénieur minier sud-africain D. G. Krige. En eet, durant les an-
nées 50, Krige a développé une série de méthodes statistiques empiriques an de déterminer
la distribution de minerais à partir d'un ensemble de forages. Pour rendre hommage à D.
Krige, Mathéron a nommé ce deuxième outil de la géostatistique par le  Krigeage . Il existe
plusieurs types de krigeages, nous allons exposer, dans le cadre de ce chapitre, le krigeage
ordinaire (KO) et le krigeage simple (KS).

2.2 Le krigeage ordinaire

Le krigeage est déni comme un estimateur Y ∗ d'une propriété Y en un point quelconque

de l'espace géographique. Cette estimation est faite à partir des observations eectives y(x)
de cette propriété. Il fournit ensuite un indicateur de la précision de l'estimation faite à
travers une variance d'estimation.

2.3 Aspects théoriques

L'objectif du krigeage est d'estimer la valeur de la variable régionalisée à interpoler Y (.)

en un site non échantillonné noté x0 . La première étape pour atteindre ce but consiste à
déterminer le  voisinage de krigeage . Ce voisinage se dénie par le domaine du champ D
contenant x0 ainsi que les sites x1 à xn associés aux observations utilisées dans la prévision
de Y (x0 ). Ces sites doivent former un sous-ensemble de l'ensemble du site d'observation.
Le choix du voisinage de krigeage se base sur une certaine connaissance de la structure de

13
dépendance spatiale entre les observations. La taille n de ce voisinage doit cependant être
assez grande pour mener à une estimation précise.
La forme la plus simple et la plus employée de cette technique est celle du krigeage linéaire.
L'estimation y ∗ (x0 ) faite en un point x0 par le krigeage linéaire est telle que :
n
(2.1)
X
∗
Y (x0 ) = λi Y (xi ),
i=1

Où : n est le nombre expérimentaux pris en compte dans l'estimation, λi est le poids aecté
au point expérimental xi .
Plus simple, pour résoudre le système d'équation induit par la recherche des poids λi , il
faut introduire les conditions d'optimisation. Ces conditions sont les suivantes :
- non biais

E[Y ∗ (x) − Y (x)] = 0 (2.2)

- variance d'estimation minimale

V AR[Y ∗ (x) − Y (x)] minimale (2.3)

Le problème à résoudre pour estimer la valeur d'une propriété Y consiste donc à calculer le
poids λi aecté à chaque point observé.
Quand l'hypothèse intrinsèque est vériée, l'ajustement d'une fonction autorisée au vario-
gramme expérimental permet de résoudre le système déni par les équations 2.1, 2.2 et 2.3.
De la sorte, on peut calculer les poids λi de l'équation 2.1 et donc la valeur de l'estimation
Y ∗ (x0 ).
Le calcul du poids aecté à un point observé ne dépend pas du tout de la valeur de la
variable étudiée en ce point. Il dépend uniquement :
- de la structure spatiale de la variable révélée par le variogramme
- de la distance géographique du point observé au point à estimer.
On a
n
!
X
V ar λi Y (xi ) − Y (x0 )
i=1
n
! n
!
X X
= V ar λi Y (xi ) + V ar (Y (x0 )) − 2Cov λi Y (xi ), Y (x0 )
i=1 i=1
n X
X n Xn
= λi λj Cov (Y (xi ), Y (xj )) + V ar (Y ) − 2 λi Cov (Y (xi ), Y (x0 ))
i=1 j=1 i=1
n X
X n n
X
= λi λj [V ar(Y ) − γ(hij )] + V ar (Y ) − 2 λi [V ar(Y ) − γ(hi0 )]
i=1 j=1 i=1

où hij désigne la distance entre les sites d'observation xi et xj . De même hi0 désigne la
distance entre le site d'observation xi et le site x sur lequel on souhaite calculer l'estimateur.

14
En remarquant que ni=1 nj=1 λi λj V ar(Y ) = V ar(Y ), puisque les sommes de poids valent
P P
1 et en abrégeant γij = γ(hij ), on obtient :
n
! n n X
n
X X X
V ar λi Y (xi ) − Y (x0 ) =2 λi γi0 − λi λj γij
i=1 i=1 i=1 j=1

En ajoutant le terme de contrainte, nous obtenons une expression plus explicite du lagran-
gien :
n n X
n n
!
X X X
L(W, λ) = 2 λi γi0 − λi λj γij − 2λ λi − 1
i=1 i=1 j=1 i=1

Il ne reste alors plus qu'à calculer les dérivées partielles (en λi et λ, soit n + 1 dérivées au
total). La dérivée partielle par rapport à un poids λi arbitraire s'écrit :
n
∂L X
= 2γi0 − 2 λj γij − 2λ
∂λi j=1

En divisant cette équation par 2, puis ajoutant l'équation sur la dérivée partielle en λ on
obtient le système d'équations de krigeage de n + 1 équations à n + 1 inconnues :

 ∂L = γ − Pn λ γ − λ = 0
(2.4)
i0 j=1 j ij
∂λi
 n λ −1=0
P
i=1 i

Remarquons que la dernière équation n'est rien d'autre que la contrainte de somme des poids
unitaire (contrainte d'autorisation et/ou d'universalité). La résolution de 2.4 s'eectue nale-
ment à l'aide de l'algèbre linéaire : on note Γ ∈ Rn×n la matrice des variances entre sites ob-
servés : Γij = γij = γ(hij ) = γ(∥xi − xj ∥) ∀i, j ∈ [1, · · · , n]. Grâce à la modélisation variogra-
phique, cette matrice Γ possède nécessairement les propriétés algébriques adéquates, et nous
n'aurons pas de mauvaise surprise par la suite. De même on note Γ0 = [γ10 γ20 · · · γn0 ]T ∈ Rn
le vecteur des variances entre les sites observés et le site à estimer x0 . On vérie alors facile-
ment que 2.4 s'écrit matriciellement sous la forme :
        −1  
Γ 1 W Γ0 W Γ 1 Γ0
= ⇒ =
1T 0 λ 1 λ 1T 0 1

où 1 = [1, · · · , 1]T est un vecteur colonne contenant n fois la valeur 1, W = [λ1 , · · · , λn ]T .

2.4 Mise en oeuvre du krigeage

Dans cette section, nous voyons des exemples pratiques d'interpolation par krigeage à
partir d'un semis de points d'observation. Nous verrons un exemple plus opérationnel par la
suite, dans lequel nous résolverons le problème modèle et construirons notre premier modèle
numérique de terrain (MNT) par krigeage.
La méthodologie est la suivante (en supposant que γ est un modèle de variogramme licite,
dont les paramètres ont été estimés éventuellement à partir d'un variogramme expérimental) :

15
 ◦ Calcul de la matrice H contenant toutes les distances hij = ∥xi − xj ∥ entre les sites
observés (en général, cette matrice a déjà été évaluée précédemment lors de l'estimation
du variogramme expérimental), et du vecteur H0 contenant les distances hi0 = ∥xi −
x0 ∥ entre les sites observés et le site à interpoler.
◦ Transformation de H et H0 par la fonction γ : Γ = γ(H) et Γ0 = γ(H0 ).
◦ Formation de la matrice A en complétant H par une colonne et une ligne de 1 et en
ajoutant un 0 dans le coin inférieur droit.
◦ Formation de la matrice B en complétant le vecteur H0 par un 1.
◦ Calcul de la solution X du système : AX = B et récupération de ses n premiers
éléments dans un vecteur colonne de poids W.
◦ Calcul de l'estimation par le produit scalaire : Y ∗ (x0 ) = WT Y
◦ Calcul de la variance d'estimation associée : σ 2 = WT Γ0 − λ, où λ est le dernier terme
de X

2.5 Application à l'évaluation du risque de crue

Dans cette activité, qui constitue en quelques sortes un travail de synthèse des notions
vues dans ce cours, nous allons évaluer la surface inondable d'un terrain donné, dont nous
disposons d'un semis de relevés d'altitude heights.txt (disponible dans le répertoire associé à
ce cours, cf en deuxième de couverture). Pour contrôler la qualité de notre estimation, nous
utiliserons également le MNT (de résolution r = 250 m) de la zone : mnt.asc. Ces deux chiers
peuvent être inspectés à l'aide d'un éditeur de text type bloc-notes.
Jusqu'à présent, dans un but purement pédagogique, nous n'avons utilisé que les fonc-
tions de base du langage R. En pratique, le traitement numérique des problèmes concrets
de Géostatistique présentent une multitude d'écueils (erreurs d'arrondi, problème de condi-
tionnement des matrices, non-convergence des régressions paramétriques...) rendant ainsi
quasi-indispensable le recours `a des librairies dédiées qui prennent en charge ces problèmes
de manière transparente pour l'utilisateur.
Nous allons utiliser la librairie R gstat(Pebesma, 2020), qui pourra être installée à l'aide de
l'instruction install.packages("gstat"), puis en sélectionnant un serveur dans la liste proposée.
Pour pouvoir facilement manipuler des données géographiques, nous utiliserons également le
package sp (Pebesma et Bivand, 2005) : install.packages("sp"). Pour activer ces librairies, on
placera les deux lignes d'instructions suivantes en entête du code :
library("gstat") : Librairie de Geostatistique
library("sp") : Librairie de gestion de données spatiales
Malgré le soin porté à la mise-à-jour de ce document, il est n'est pas impossible que
certaines fonctionnalités de gstat présentées ci-dessous ne soient plus disponibles exactement
sous le même formalisme syntaxique. On pourra en général facilement résoudre le problème
en se référant au manuel d'utilisation (Pebesma, 2001) le plus récent.
Enn, précisons que cette activité n'est en principe pas trop sensible au caractère aléatoire
des réalisations. Malgré tout, dans un souci de reproductibilité et de débuggage éventuel des
problèmes, nous mentionnons que les résultats ci-dessous ont été obtenus avec la graine :
set.seed(1).
On considère une région côtière, d'altitude comprise entre 0 et 154 m, et soumise à un

16
risque de montée des eaux, dont une analyse préalable a permis de montrer que tous les ter-
rains situés à une altitude inférieure à 5 m pouvaient potentiellement être impactés. L'objectif
de l'étude est de déterminer la surface totale (en km2) à évacuer.
Q1. Dans un premier temps, nous allons évaluer cette surface dans l'hypothèse où le MNT
de la zone est connu. La valeur trouvée à l'issue de cette étape préliminaire constituera une
vérité terrain pour pouvoir comparer par la suite les résultats obtenus respectivement par
krigeage et par simulations.
On commence par xer deux paramètres : la résolution du MNT à disposition (r) et le seuil
d'altitude en dessous duquel un terrain est situé en zone inondable (threshold). On dénit
également une fonction estimate permettant de calculer la surface inondable totale d'un MNT.
Notons que threshold est une variable globale. Sa modication entraîne automatiquement
celle de la fonction estimate. Par ailleurs, soulignons la multiplication par r2 (qui permet
d'exprimer un nombre de cellules en une surface) puis la division par 106 (qui assure la
conversion en km2).
Charger le MNT du chier mnt.asc et calculer la surface inondable sur la zone d'étude.
Q2. On suppose à présent ne pas avoir de MNT à disposition. L'ingénieur chargé de
l'étude décide alors de relever à l'aide d'un GPS professionnel de précision centimétrique
un total de 150 points d'altitude répartis aléatoirement et uniformément sur la zone. Le
résultat de la campagne est con- signé dans le chier heights.txt, dans lequel chaque ligne
représente les coordonnées géographiques (X, Y ) d'un site, exprimées dans une projection
plane quelconque, et l'altitude Z qui y a été mesurée. L'objectif consiste à essayer d'estimer
avec la meilleure précision possible, la surface de la zone inondable à partir de ces données
réduites.
Q3. Proposer un modèle de variogramme et estimer ses paramètres.
Après inspection graphique du résultat, il paraît raisonnable de choisir un modèle de
variogramme exponentiel (notons qu'un variogramme sphérique, ou même linéaire avec pallier
semblent être également des solutions convenables). Dans gstat, chaque modèle est déni par
un code de 3 lettres. On pourra consulter la liste des modèles disponibles avec l'instruction
vgm(). La régression paramétrique se fait avec la commande t.variogram.
Q4. Calculer par krigeage le relief du terrain à partir du semis de points observés et en
déduire une estimation de la surface inondable.
Pour calculer le champ par krigeage, on doit commencer par dénir la grille d'interpo-
lation. An de travailler sur une zone identique à la vérité terrain, on lui donne les mêmes
paramètres (résolution et taille) que le MNT.
Q5. Calculer 500 simulations conditionnées au semis de points observés. Pour chaque
simulation générée, on évaluera la surface de la zone inondable et on stockera la valeur
obtenue. Calculer la moyenne des estimations obtenues (on pourra également dériver un écart-
type et des bandes de conance). Cette nouvelle valeur est-elle plus précise que l'estimateur
du krigeage ?
Pour eectuer des simulations conditionnelles avec la bibliothèque gstat, on utilise exactement
la même fonction que pour le krigeage à laquelle on ajoute une entrée nsim permettant de
spécier le nombre de simulations `a calculer

17
Code R :
library("gstat")
library("sp")
set.seed(1)
par(mfrow=c(2,4))
r = 250
threshold = 5
map = c(topo.colors(255)[50:60], terrain.colors(255))

# -------------------------------------------------
# Utilitaires
# -------------------------------------------------

estimate = function(Z){
return(length(which(Z < threshold))*r^2/10^6)
}

# -------------------------------------------------
# Préparation des données
# -------------------------------------------------

S = as.matrix(read.table("C:/Users/Desktop/cours de geostat/mnt.asc"))
PTS = read.csv("C:/Users/Desktop/cours de geostat/heights.txt",
header=1, sep=",", na.strings="NA", dec=".", strip.white=TRUE)

N = nrow(PTS)
gx = (1:nrow(S))*r-r/2
gy = (1:ncol(S))*r-r/2
image(gx, gy, S, col=terrain.colors(255))
points(PTS$X, PTS$Y, pch=16, cex=.5)
S[which(S<threshold)] = 0
image(gx, gy, S, col=map)

# -------------------------------------------------
# Calcul du variogramme
# -------------------------------------------------

coordinates(PTS)=~X+Y
Z = PTS$Z
DX = max(PTS$X)-min(PTS$X)
DY = max(PTS$Y)-min(PTS$Y)
D = sqrt(DX^2+DY^2)/2 # Choix automatique
D = 13000 # Choix manuel

18
vario = variogram(Z~1, data=PTS, cutoff=D, width=1000)
vmod = fit.variogram(vario, vgm("Exp"))
h = 0:D
g = vmod$psill[1] + vmod$psill[2]*(1-exp(-h/vmod$range[2]))
plot(vario$dist, vario$gamma, pch=3, ylim=c(0,max(vario$gamma)),
xlim=c(0,max(vario$dist)))
lines(h, g, col='blue')

# -------------------------------------------------
# Krigeage
# -------------------------------------------------

nx = length(gx)
ny = length(gy)
GRID = expand.grid(x=gx, y=gy)
gridded(GRID) = ~x+y
krigeage = krige(Z~1, PTS, GRID, model = vmod, nmax=50, debug.level=-1)
K = matrix(krigeage$var1.pred, nx, ny)
image(gx, gy, K, col=terrain.colors(255))
K[which(K<threshold)] = 0
image(gx, gy, K, col=map)

# -------------------------------------------------
# Simulation
# -------------------------------------------------

SURFACES = rep(0,500)
simulation = krige(Z~1, PTS, GRID, model = vmod, nmax = 500,
nsim = length(SURFACES), debug.level=-1)

for (repetition in 1:length(SURFACES)){

SURFACES[repetition] = estimate(simulation[[names(simulation)[repetition]]])
}

SIM = matrix(simulation$sim1, nx, ny)

image(gx, gy, SIM, col=terrain.colors(255))
SIM[which(SIM<threshold)] = 0
image(gx, gy, SIM, col=map)
d = density(SURFACES)
surface_real = estimate(S)
surface_krig = estimate(K)
m = mean(SURFACES)
b_inf = as.vector(quantile(SURFACES, 0.05))
b_sup = as.vector(quantile(SURFACES, 0.95))

19
plot(d$x, d$y, type='l', xlab="surface (km2)", ylab="freq")
abline(v=m)
abline(v=b_inf, lty=2)
abline(v=b_sup, lty=2)
abline(v=surface_real, col='blue')
abline(v=surface_krig, col='red')
Résultat du code R

1000
●
●
● ● ●
● ●
● ●
25000

25000

25000
● ● ●
●
●
● ●
● ●
● ●
● ●
● ●

800
●

● ●
●
●
● ● ●
●
20000

20000

20000
●
● ● ●
●
●
● ●
●

● ●
● ●

600
● ●

vario$gamma
● ●
●
●
15000

15000

15000
● ●
●
●
● ● ●
● ● ●
gy

gy

gy
● ●
● ● ●
●
●
●
●
● ●
●

400
● ●
● ●
●
10000

10000

10000
● ● ●
● ●
● ● ●

●
● ● ●
● ●
● ● ●
● ●
● ●
● ●
● ●
● ●
● ● ●
●
●

200
● ●
● ●
●
5000

5000

5000
● ● ●
● ● ●
● ● ● ●
●
● ●
● ● ● ● ●
●
● ● ●
● ● ●

●
●
●
●
● ●
● ● ●

0
●
0

0
0 5000 10000 15000 0 5000 10000 15000 0 5000 10000 15000 0 5000 10000 15000

gx gx vario$dist gx
25000

25000

25000

0.06
20000

20000

20000
15000

15000

15000

0.04
freq
gy

gy

gy
10000

10000

10000

0.02
5000

5000

5000

0.00
0

0 5000 10000 15000 0 5000 10000 15000 0 5000 10000 15000 15 20 25 30 35 40 45

gx gx gx surface (km2)

Figure 2.1  1 : MNT ; 2 : surface réelle inondable ; 3 : Variogramme ; 4 : MNT par krigeage ;
5 : surface prédite par Krigeage ; 6 : MNT par simulation ; 7 : surface prédite par simulation ;
8 : densité de la surface prédite par simulation (rouge : S. krigeage ; noir : S. simulation ;
bleu : S. réelle).
La vraie surface inondable : 30.8 Km2 ; La surface inondable prédite par Krigeage : 16.9 Km2 ;
La surface inondable prédite par 50 simulation du Krigeage : 28.8 Km2 .

2.6 APPLICATION TO OZONE POLLUTION FORE-

CASTING AT THE NON-VISITED SITES

In this section, we apply our methodology to predict the level of ozone pollution at the
non-visited sites of California state. For that, we consider the available data on internet site
https ://www.epa.gov/outdoor-air-quality-data. The explanatory functional variables
{Xsi (t), t = 1, · · · , 100, si = (Latitude, Longitude)i , i = 1, · · · , 51}
correspond to the measurements of ozone concentration obtained for the p = 100 rsts days,
from January 1st, 2021 to April 12th, 2021 on each of n2 = 51 sites. The response variables
{Ysi , si = (Latitude, Longitude)i , i = 1, · · · , 35}

20
correspond to the measurements of ozone concentration obtained the April 13th, 2021 on each
of 35 rsts stations. For evaluating the performances of method described from the spatial
functional linear regression model (SFLR) studied in Bouka et al (2023), we yet compare it
to that of SFRD dened from the spatial functional linear regression model with derivatives
(SFLRD) studied in bouka et al. (2018). Specially, we write the SFLRD as
Z 1 Z 1
Yiℓ = γ1 (t)Xiℓ (t)dt + γ2 (t)Xi′ℓ (t)dt + ϵiℓ ,
0 0

where γ1 and γ2 are two functions to estimate, Xi′ℓ stands for the rst derivative of Xiℓ .
Recall that the SFLR model is dened when γ2 (t) = 0 ∀t ∈ [0, 1]. So, we predict from both
methodologies
{Ysi , si = (Latitude, Longitude)i , i = 36, · · · , 51},
which correspond to the measurements ozone concentration obtained the April 13th, 2021
on these 16 others sites assumed non-visited at this date. In what follows, we illustrate our
data by graphics (see Figure 2.2).

Ozone concentration to the 35 sites Ozone concentration at the 1st site

0.07

●
0.06
Ozone concentration

Ozone concentration

● ●
0.05

● ●● ● ●
● ●● ● ● ●● ●
● ● ● ●
● ●●● ●●● ● ● ● ● ● ●●
0.04

● ● ● ● ●
● ● ● ● ● ● ● ● ●● ● ● ●
● ● ● ● ● ●
●●●● ● ● ● ●●● ●● ● ● ●
0.03

● ● ● ●●
● ●● ● ●
●
● ● ●
● ●
● ● ● ●●
0.02

●
●
0.01

● ●
● ●

0 20 40 60 80 100 0 20 40 60 80 100

Days Days

Smoothed vs. raw data at the 1st site

spatial response variable
●
0.06

0.065
−122.0 −121.6 −121.2
Ozone concnetration

● ● 0.060
● ●●
●● ● ● 0.055
● ● ● ● ● ●● ●
● ●
●●● ●●● ●
● ● ● ● ● ●● 0.050
0.04

● ● ●
Lon

● ● ● ● ● ● ● ●
● ● ●● ● ● ●
0.045
● ● ● ● ● ●
● ● ●●●● ● ● ● ●●● ●● ● ● ●
● ●●
● ●● ● ● 0.040
●
● ● ●
●
0.035
● ● ● ● ●●
0.02

● 0.030
●
● ●
● ●

37.65 37.70 37.75 37.80 37.85

0 20 40 60 80 100
Lat
Days
RMS residual = 0.004989

Figure 2.2  Ozone pollution data : daily ozone concentration at the 35 stations (left and
top), daily ozone concentration at the rst station (right and top), smoothed and raw of daily
ozone concentration data at the rst station (left and down), ozone pollution level (right and
down).

In these previous graphics, we see that these data correspond to our study. In what follows,
we present results of our study.

SFLR (nbasis=7, ρ = 0.09) SFLRD (w = 0.01, ϕ = 0.11)

Prediction error (PE) 0.0282 0.038

Table 2.1  Prediction error computed from both methods.

21
 SFLR SFLRD

0.08

0.08
0.07

0.07
Predicted ozone concentration

Predicted ozone concentration

0.06

0.06
●
●
● ●●
0.05

0.05
● ●
●
● ● ●●
● ●● ● ●●
●
● ●● ●
●
●
●
0.04

0.04
●
●
0.03

0.03
0.02

0.02
0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.02 0.03 0.04 0.05 0.06 0.07 0.08

Measured ozone concentration Measured ozone concentration

Figure 2.3  Predicted values of ozone concentration, from SLFR (left) and SLFRD (right),
versus the measured values.

Both graphics in Figure 2.3 present a small dierence and it is conrmed by computation
(see table 2.1) of prediction error (PE) given by
v
u 51  2
uX
PE = t Ysi − Ybsi .
i=36

We see a very minor advantage for prediction obtained from model SFLR studied in this
paper.

Bibliography

Bouka S, Dabo-Niang S, Nkiet GM (2018). On estimation in a spatial functional linear

regression model with derivatives. C. R. Acad. Sci. Paris Ser. I, 356 : 558562.
Bouka S, Dabo-Niang S, Nkiet GM (2023). On estimation and prediction in spatial functional
linear regression model. Lithuanian Mathematical Journal, 63(1), 1330.
Delhomme J.P. (1976). Application de la théorie de la variable régionalisée dans la science
de l'eau. Thèse Doc. Ing.. ENSM de Paris. 130P
Matheron G. (1965) - Les variables régionalisées et leur estimation. Paris. MASSON 305P.
Pebesma E. J. et Bivand R. S. (2005). Classes and methods for spatial data in R. R News,
5(2) : 913.
Pebesma E. J. (2001). Gstat user's manual. Dept. of Physical Geography, Utrecht University,
Utrecht, The Netherlands.
Pebesma E. (2020). The meuse data set : a brief tutorial for the gstat r package.

22