See discussions, stats, and author profiles for this publication at: [Link]

net/publication/333667153

Statistique et Modelisation de données spatiales

Research · June 2019

CITATIONS READS
0 2,303

1 author:

Nasradine Ibrahim Ahmed

Université de Lille
1 PUBLICATION 0 CITATIONS

SEE PROFILE

All content following this page was uploaded by Nasradine Ibrahim Ahmed on 09 June 2019.

The user has requested enhancement of the downloaded file.

 Université de Lille 3

Rapport de Recherche
2018-2019

Statistique et Modelisation de
données spatiales

Realisé par : Encadré par :

Nasradine IBRAHIM Pr. Sophie DABO

9 juin 2019
 TABLE DES MATIÈRES

1 Introduction 4

2 Types de données spatiales 6

2.1 Données de type geostatistique . . . . . . . . . . . . . . . . . . . 7
2.1.1 Les précipitations de la république de Djibouti . . . . . . 7
2.2 Données de type surfacique . . . . . . . . . . . . . . . . . . . . . 8
2.2.1 Mortalité infantile en Afrique . . . . . . . . . . . . . . . . 9

3 Modélisation geostatistique 11
3.1 Rappels sur les processus stochastiques . . . . . . . . . . . . . . . 11
3.2 Processus stationnaire au second ordre . . . . . . . . . . . . . . . 12

4 Variabilité spatiale 13
4.1 Le variogramme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.1.1 Processus intrinsèques . . . . . . . . . . . . . . . . . . . . 13
4.1.2 Variogramme d’un processus intrinsèque . . . . . . . . . . 14
4.1.3 Variogrammes usuelles . . . . . . . . . . . . . . . . . . . . 15
4.2 Analyse du correlogramme . . . . . . . . . . . . . . . . . . . . . . 17
4.2.1 Estimation empirique d’un variogramme . . . . . . . . . . 17

1
 4.2.2 L’ajustement du variogramme . . . . . . . . . . . . . . . . 18

5 Interpolation par krigeage 20

5.1 Origine et principe du krigeage . . . . . . . . . . . . . . . . . . . 20
5.2 Le krigeage simple . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.2.1 Minimisation de la variance de l’erreur de prévision . . . . 23
5.3 Le krigeage ordinaire . . . . . . . . . . . . . . . . . . . . . . . . . 26

6 Tests d’autocorrelation spatiale 29

6.1 Indice de Moran . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
6.2 indice de Geary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

7 Modélisation spatiale par krigeage sur R 33

7.1 The meuse data set . . . . . . . . . . . . . . . . . . . . . . . . . . 33
7.2 Analyse variographique . . . . . . . . . . . . . . . . . . . . . . . 35
7.2.1 Le variogramme expérimental . . . . . . . . . . . . . . . . 35
7.2.2 Modélisation du variogramme . . . . . . . . . . . . . . . . 38
7.3 Le krigeage avec R . . . . . . . . . . . . . . . . . . . . . . . . . . 40

A Le krigeage sur le langage R 42

B bibliographie 46

2
 Remerciement

Tout d’abord, ce travail est le fruit d’une collaboration avec So-

phie DABO ,responsable du parcours MIASHS de l’UFR MIME à
l’Université de Lille 3 .Je remercie énormément son soutien qu’elle
m’a donné et également pour sa gentillesse .

3
CHAPITRE 1
INTRODUCTION

Ce sujet consiste à mettre en lumière la notion de statistique spatiale ou

modélisation de données spatiales .

Plusieurs domaines de recherches comme la foresterie ,la température,l’épidémiologie,la

science des cartes géographiques,...., sont à la fois localisables dans l’espace géo-
graphique et présentent des dépendances fortes ou faibles ou parfois absences
de dépendances.

Soit {Y (s), s ∈ Ω} ou Ω ⊂ <d une collection aléatoire , s est une zone ou

localisation spatiale et Y une quantité aléatoire.

Notre but est de caractériser la dépendance spatiale entres différentes obser-

vations à l’aide des outils statistiques spatiales .

Dans un premier temps nous allons aborder l’approche géostatistique qui

es une méthode spécifique pour étudier les données référencées de types géosta-
tistiques.

Après passage au second ordre des processus stationnaires, on devra analyser

empiriquement le variogramme ou fonction de covariance afin de le modéliser

4
ou de l’ajuster pour finalement se servir de l’interpolation par krigeage.

On pourra ensuite détailler l’étude d’auto-corrélation spatiale des données

spatiales en effectuant des tests spécifiques dites tests d’auto-corrélation spa-
tiales.

Dans un second temps , nous nous évertuerons la partie pratique du sujet

en effectuant la modélisation spatiales de données par krigeage sur le langage
R.

5
CHAPITRE 2

TYPES DE DONNÉES SPATIALES

Le principe ou le but de l’analyse spatiale est d’explorer le positionnement

spatial des vecteurs spatiaux,des phénomènes afin de comprendre leurs caracté-
ristiques.
Quel que soit le type de données spatiales ,elle nous pousse à comprendre les
phénomènes qu’on associent à des processus aléatoires {Y (s), s ∈ Ω} .
Avant toute analyse ,il faut prendre en mains les données brutes afin de les
adapter au processus spatial qui leurs [Link] permet non seule-
ment dans la démarche de l’analyse, de voir des problèmes particuliers dans les
données(données manquantes,problème lors de la collecte,ou valeurs aberrantes
et de laisser des discussions les hypothèses nécessaires des méthodes économé-
triques .

On distingue alors deux grandes structures de données spatiales : Données

géo-statistiques ,données surfaciques.

Pour tout type de données ,on lui apporte des méthodes statistiques spéci-
fiques.

6
2.1 Données de type geostatistique

Les données spatiales sont des informations (ou observations) qui nous donnent
à la fois les valeurs mais aussi les localisations geo-référencées .Ces objets geo-
localisés définissent un ensemble des coordonnées spatiales et constituent des
informations riches lors de l’analyse.

Les données de type geostatistique ,appelées également données continues

sont modélisées à partir des processus aléatoires {Y (s), s ∈ Ω} .A chaque mesure
de position s obtenue on l’associe à une variable aléatoire discret Ys .Ainsi il
existe une valeur pour la variable d’intérêt en tout point s du territoire étudié.

Ils sont utiles lorsqu’on analyse la composition chimique du sol ,la qualité
de l’eau ou encore la précipitation des pluies .

L’analyse de ces types de données spatiales consistent à prédire la valeur

d’une variable en un point ou elle n’a pas été échantillonné.

La méthode d’interpolations la plus utilisée pour analyser spatialement est

le krigeage .Elle s’agit de prédire Y sur l’ensemble des sites S en site s ∈ S
inconnu

2.1.1 Les précipitations de la république de Djibouti

Le 15 juin 2016,des chercheurs geostaticiens dans le domaine de la météoro-

logie étudient la dépendance spatiale entre les données obtenues dans des sites
différents grâce la technique du krigeage .Finalement ils ont abouti à la prédic-
tion des precipitations à l’interieur du pays (Djibouti), taux de précipitations
en des lieux ou le sol n’a pas été échantillonné .

7
 Figure 2.1 – Les précipitations de la République de Djibouti
La source des données nous provient de l’entreprise Own work .

2.2 Données de type surfacique

Les données de type surfacique appelées également données laticielles sont

des données dont la localisation des observations est [Link] caractérisent une
partition des régions en zones contiguë[Link] peuvent représenter des points
fixes du territoires. Lors de l’analyse du PIB par région ou le nombre de mariage
par mairie ,ces types de données nous donnent des informations potentiellement

8
riches. D’abord on s’intéresse les valeurs des observations voisines.

Pour analyser spatialement ces types de données surfaciques ,on définit la

notion de voisinage des observations puis l’influence qu’exerce une obser-
vation i sur sa voisine j enfin on évalue la significativité de cette influence
par des tests spécifiques qu’on pourra voir sur le chapitre 6.

2.2.1 Mortalité infantile en Afrique

En Afrique ,on s’est intéressé de comprendre la mortalité subite du nour-

risson . Pour cela commence par définir notre processus aléatoire Ys et on
lui associe au nombre de cas de mort subite dans chaque pays S (S regroupe
l’ensemble des pays africains) .

Après une analyse spatiale basée sur les relations entre les observations ,on
obtient cette figure suivante :

Figure 2.2 – Données surfaciques :Mortalité infantile en Afrique

9
 Une carte présentant la mortalité infantile en Afrique. Plus la couleur est
rouge foncé plus la mortalité infantile du pays est importante. Inversement plus
la couleur est clair moins la mortalité du pays est importante.

Les pays du continent Africain où la mortalité infantile sont les plus impor-
tantes sont le Niger et le Mali avec près de 110 mort pour mille naissances. Au
contraire, les pays où la mortalité infantile est moins la importante sont Mayotte,
l’Ile-Maurice et les Seychelles avec environ 11 mort pour mille naissances.

Les données proviennent du rapport annuel de la CIA en Avril 2012.

10
CHAPITRE 3

MODÉLISATION GEOSTATISTIQUE

3.1 Rappels sur les processus stochastiques

Soit (Ω ,F,P) un espace de probabilité ,S un ensemble de sites, et (E,ε)

un espace d’état mesurables .

Définition 3.1.1 Un processus stochastique à valeur dans E est une

famille {Y (s), s ∈ S} de variable aléatoire définie sur (Ω ,F,P) et à valeurs
dans (E,ε).

Définition 3.1.2 Y est processus du second ordre si pour tout s

∈ S , E(Ys2 ) < ∞.La fonction m :S → R tel que m(s)=E(Ys ) représente la
moyenne de Y. La covariance de Y est la fonction c :S × S → R tel que pour
tout s,t c(s,t)=ov(Ys , Yt ).

11
3.2 Processus stationnaire au second ordre

Dans un premier temps ,il va falloir étudier la stationnarité du processus

spatial pour voir si la variabilité du processus obéit la même loi dans les obser-
vations .Tel est le but de l’hypothèse stationnaire que nous allons définir :Une
fonction stationnaire reste invariante dans l’espace et cette répétition d’inva-
riance rend à nouveau possible l’inférence statistique à partir d’une réalisation
.

Soit Y un processus du second ordre ,la fonction m une composante déter-

ministe et la fonctionR une composante aléatoire tel que Y(s) = m(s)+R(s)
.

La stationnarité au second ordre n’a plus d’influence sur la loi de probabi-

lité,mais seulement sur la moyenne et la covariance.

Elle vérifie les conditions suivantes :

— E [Y (s)] = m(s)
L’invariance de l’espérance par translation entraîne la constance de la
composante déterministe ∀ s ,on a m(s+h) = m(s) .
— La variance est constante : E((Y (s) − m)2 ) = σ 2
— La covariance qui ne dépend que le décalage spatial h
Cov [Y (s + h), Y (s)] = (E [Y (s + h) − m)] (E [Y (s) − m)]= C(h)

12
CHAPITRE 4

VARIABILITÉ SPATIALE

4.1 Le variogramme

4.1.1 Processus intrinsèques

Lorsque la notion de stationnarité devient trop forte .Une de ces est sou-
vent quand la moyenne change sur le territoire d’intérêt et que la variance ne
peut pas toujours être borné lorsque cette région d’intérêt s’agrandit .Georges
Matheron établit alors après avoir tiré les conséquences de ces limites de la
stationnarité au second ordre en proposant une encore plus faible : la station-
narité intrinsèque .

Définition 4.1.1 Processus intrinsèque

Un processus(ou champ) Y est dit intrinsèquement stationnaire ou intrin-

sèque si ses accroissements respectent les conditions suivantes ∀s, h :

— D’espérance nulle : E [( Y(s+h)-Y(s)] = 0

— La covariance qui ne dépend que le décalage spatial h

13
 Cov [Y (s + h), Y (s)] = (E [Y (s + h) − m)] (E [Y (s) − m)]= C(h)

4.1.2 Variogramme d’un processus intrinsèque

Définition 4.1.2 Continuité en moyenne quadratique

Soit Y un processus spatial ,on dit que Y est continu en moyenne quadratique
si,pour tout s,

2
lim E [(Y (s + h) − Y (s)] = 0 , ∀s, h .
x→0

2
si on note E [(Y (s + h) − Y (s)] = V ar [Y (s + h) − Y (s)]) .Cette quan-
tité mesure alors la dispersion de l’écart entre Y(s+h) et Y(s).

Ensuite ,dire que Y est continue en moyenne quadratique revient donc à

constater que la dispersion de l’écart entre Y(s+h) et Y(s) est d’autant plus
proche de zéro que le décalage spatial h devient faible.

Dans le cas d’un champ stationnaire à l’ordre 2,

V ar [Y (s + h) − Y (s)))] = V ar(Y (s + h) + Var(Y (s) − 2Cov(Y (s + h), Y (s))=

2(C(0) − C(h)) ,ou C est la fonction d’auto-covariance du processus Y.

A partir de maintenant ,on considère la fonction Γ(h) = (C(0) − C(h))

est appelée variogramme du processus [Link] fonction est symétrique et
positive , Γ(0) = 0.

Dans le cas pratique ,souvent le variogramme est présenté comme la figure

ci-dessous :

14
 Figure 4.1 – Variogramme d’un processus

Porté a : La portée du variogramme (la portée pratique) est la distance à

partir de laquelle le variogramme atteint son palier . À cette distance, la valeur
du variogramme correspond à la variance de la variable aléatoire.

Palier C0 +C : Correspond à la valeur de la variance aléatoire du processus

.

Effet pépite C0 : Variation à très courte échelle, erreurs de localisation,

erreurs d’analyse et précision analytique.

4.1.3 Variogrammes usuelles

La littérature géostatistique propose de nombreuses fonctions qui satisfont

les propriétés du semi-variogramme tel qu’il est présenté dans la ci-dessous. Ces
fonctions paramétrées doivent permettre de décrire les différentes composantes

15
(portée, palier, pépite). Elles doivent aussi gérer le comportement de la fonction
à l’origine (tendance linéaire, tangence horizontale ou verticale).

Comme l’expliquent Matheron et les géostatisticiens , la modélisation est une

affaire de choix. Le choix du modèle théorique est un moment décisif dans la dé-
marche du géostatisticien, mais on ne peut pas associer a priori un variogramme
théorique à tel ou tel type de processus.

Figure 4.2 – Variogrammes usuelles

16
4.2 Analyse du correlogramme

Après avoir détaillé le choix du modèle probabiliste, nous nous sommes fo-
calisés sur un outil indispensable pour la modélisation spatiale des données : le
variogramme. Le variogramme reflète la structure de la régionalisation. On ap-
pelle analyse variographique l’inférence du variogramme à partir de données
expérimentales.

Elle permet de restituer des informations quant à la distribution spatiale

de la variable régionalisée, c’est-à-dire de la réalisation de la fonction aléatoire
d’intérêt. On s’attache usuellement à l’étude du variogramme plutôt que de
la covariance car, comme nous l’avons vu précédemment, celle-ci n’étant pas
nécessairement définie dans le cadre intrinsèque.

L’analyse variographique constitue une étape cruciale dans une étude géo-
statistique : elle permet postérieurement d’estimer les valeurs inconnues de la
variable régionalisée et d’assortir leur estimation d’une précision.

4.2.1 Estimation empirique d’un variogramme

On cherche maintenant à estimer le variogramme à partir d’un ensemble de

données expérimentales. Le variogramme expérimental est calculé empirique-
ment à partir des observations de la variable régionalisée. Il est donné par la
formule pour une valeur h donné :

1 1 Pn
Γ̂(h) = 2 N (h) si −sj=h (Y (si ) − Y (sj )2 )

ou N(h) est le nombre de paires de points (si ; sj ) séparés par la distance

h. En pratique, le variogramme expérimental est calculé pour des classes de
distances, et par classes de direction si nécessaire.

Preuve

17
 Pn
E[Γ̂(h)] = E[ 12 N 1(h) h (Y (si ) − Y (sj ))2 ]

Pn
= N 1(h) E[ 12 h (Y (si ) − Y (sj ))2 ]

1 Pn
= N (h)
[ h 12 V(Y (si ) − Y (sj )]

1 Pn
= N (h) h Γ(h)

= Γ(h)

Ainsi , on va bien que l’estimateur est sans biais car son espérance est égale
à l’expression générale du variogramme.

4.2.2 L’ajustement du variogramme

Il faut ensuite trouver une forme fonctionnelle adaptée au variogramme expé-

rimental. Une première étape importante est d’obtenir une représentation lissée
du variogramme, qui peut être obtenue avec la fonction variog de geoR.

En fin de compte, le choix final du modèle devrait refléter à la fois le résultat

de la procédure d’ajustement au modèle statistique et une interprétation cohé-
rente avec la compréhension scientifique du processus à l’étude". De nombreuses
méthodes sont proposées pour ajuster le variogramme. On trouve des méthodes
fondées sur les moindres carrés, ordinaires ou pondérés, des méthodes fondées
sur la vraisemblance, ainsi que des méthodes bayésiennes.

Ajustement par les moindres carrés ordinaires

On cherche le vecteur des paramètres de la fonction qui minimise une fonction

d’objectif simple, la somme des carrés des distances entre la valeur du semi-
variogramme expérimental et celle du variogramme théorique.

18
 k
ˆ i ) − Γ(hi , θ)2 )
X
θ̂M CO = argmin (Γ(h
θ∈Θ
i=1

Ajustement par les moindres carrés pondérés

Les moindres carrés ordinaires ne tiennent pas compte du nombre de couples

de points qui interviennent dans le calcul de chacun des points du variogramme
expérimental, contrairement aux moindres carrés pondérés.

k
X N (hi ) ˆ i ) − Γ(hi , θ)2 )
θ̂M CP = argmin (Γ(h
θ∈Θ
i=1
(Γ(hi , θ)2 )

Ces deux méthodes sont utilisées par deux écoles de pensées concernant
l’ajustement automatiques des modèles de variogramme. Certains prétendent
que l’ajustement visuel demeure la meilleure méthode, alors que d’autres pré-
fèrent les méthodes automatiques.

19
CHAPITRE 5
INTERPOLATION PAR KRIGEAGE

5.1 Origine et principe du krigeage

Le krigeage est une méthode d’interpolation applicable à des données spa-

tiales. Elle s’appuie sur la géostatistique linéaire, notamment le variogramme.
La théorie du krigeage a été développée par un mathématicien français (G. Ma-
theron) à partir des travaux de l’ingénieur minier sud-africain [Link]. Dans
les années 50, Krige a développé une série de méthodes statistiques empiriques
afin de déterminer la distribution de minerais à partir d’un ensemble de forages.

Le krigeage est un interpolateur exact (la valeur estimée sur un point de

mesure est égale à la valeur du point de mesure) et optimal (il est sans biais et
minimise la variance sur l’erreur d’estimation).

Contraintes : L’objectif est d’estimer le plus précisément possible la

grandeur recherchée par une combinaison linéaire pondérée des observations.
Ce critère de précision se traduit par une minimisation de l’erreur quadratique
moyenne. Ces conditions se traduisent mathématiquement en 4 contraintes.

20
— Contrainte de linéarité : L’estimé doit etre une combinaison lineaire
des données :
Pn
Yi∗ = i=1 λi Yi
On pourra écrire sous cette forme
Yi∗ = λ1 Y1 + λ2 Y2 + λ3 Y3 + ........
˙ + λn−1 Yi−1 + λn Yi
Proposition On considère Y une variable aléatoire et stationnaire C(h)
sa covariance , et Γ(h) son [Link] Y ∗ une combinaison linéaire
Pn
tel que Y ∗ = i=1 λi Y (xi ) .Les λi sont les poids et les xi représentent
des sites geolocalisé[Link] dit que la combinaison lineaire Y ∗ est admissible
si la somme des poids λi est égale à 0 .On peut voir que sa variance existe
Pn Pn
V(Y ∗ ) = i=1 j=1 λi λj C(xi − xj )
Preuve
Par définition, V(Y ∗ ) = E(Y ∗ − E(Y ∗ ))2 , comme Y est stationnaire
Pn
(E(Y ) = m) on a E(Y ∗ ) = m i=1 λi .
Alors
Pn
V(Y ∗ ) = E( i=1 λi (Y (xi ) − m))2
= λ21 C(x1 − x1 ) + λ22 C(x2 − x2 )
˙ + λ2n C(xn − xn ) + ........
+ ....... ˙ + 2λn−1 λn C(xn−1 − xn )
Notre objectif étant de calculer la variance de la combinaison lineaire,la
variance de la variable aléatoire estimée dans des sites différents et comme
n X
X n
Pn ∗
i=1 λ i = 0 en respectant la condition d’admissibilité ,obtient V(Y ) = λi λj C(xi − xj )
i=1 j=1
∗
Pn Pn
De même on peut écrire , Y = i=1 λi [Y (xi ) − Y (0)]
λi Y (xi ) = i=1
∗
Pn Pn ∗
or on a montrer que la variance de Y est égale à V(Y ) = i=1 j=1 λi λj C(xi − xj )
donc

n
X n
X n X
X n
V( λi Y (xi )) = V( λi [Y (xi ) − Y (0))] = λi λj C(xi − Y (0), xj − Y (0))
i=1 i=1 i=1 j=1

alors on peut calculer la variation ou la dispersion des variable aléatoires

entre deux sites géolocalisés distinctes :
V[Y (xi ) − Y (xj )] = V[Y (xi ) − Y (0) + Y (0) − Y (xj )]

21
 = V[Y (xi ) − Y (0)] + V[Y (0) − Y (xj )] − 2C(xi − Y (0), xj − Y (0))
On peut en déduire d’après cette expression en passant la covariance en
fonction des autres variances ,

2C(xi − Y (0), xj − Y (0)) = V[Y (xi ) − Y (0)] + V[Y (0) − Y (xj )] − V[Y (xi ) − Y (xj )]

C’est à dire en faisant apparaître le variogramme,

2C(xi − Y (0), xj − Y (0)) = 2Γ(xi ) + 2Γ(xj ) − 2Γ(xi − xj )

Ensuite , Il en résulte ,

n
X n X
X n
V λi Y (xi ) = λi λj Γ(xi ) + Γ(xj ) − Γ(xi − xj )
i=1 i=1 j=1

n
X n
X n
X n
X n X
X n
λj Γ(xi ) + λi Γ(xj ) − λi λj Γ(xi − xj )
j=1 i=1 i=1 j=1 i=1 j=1

Par conséquent , grâce au condition d’admissibilité de variance de la

combinaison aléatoire ,la somme de chacun des poids λi etλj étant nulle
Pn Pn
, j=1 λi = j=1 λj = 0,la variance de la combinaison existe et elle
est égale à

n
X n X
X n
V λi Y (xi ) = − λi λj Γ(xi − xj )
i=1 i=1 j=1

— Contrainte d’autorisation : On doit assurer l’existence de l’espérance

et de la variance de l’erreur de prévision.
— Contrainte de non biais : On doit s’assurer que l’estimateur de l’er-
reur de prévision est non biaisé :

E[Ŷ (s0 ) − Y (s0 )] = 0

22
 s0 ∈ S représente un site .
— Contrainte d’optimalité : La minimisation de la variance de l’erreur
après le krigeage V[Ŷ (s0 ) − Y (s0 )] s0 ∈ S

5.2 Le krigeage simple

Le krigeage simple est une forme d’interpolation dont sa prévision s’écrit de

la forme suivante :
n
X
Ŷi = µ + λi (Yi − µ)
i=1

avec µ une constante connue ou l’espérance de la variable aléatoire Y stationnaire-

ment intrinsèque .

Proposition :Soit Y une variable aléatoire stationnaire d’espérance µ ,s0 ∈

S un site dans l’ensemble de sites geolocalisés ,Ŷ la prédiction de Y par le
Pn
[Link] suppose par hypothèse i=1 λi = 1 .Alors l’espérance de l’erreur
de prévision est donc nulle c’est à dire E[Ŷ (s0 ) − Y (s0 )] = 0 .

Preuve

h i
On écrit , E Ŷ (s0 ) − Y (s0 )
Xn    
= µ+ λi E (Y (s0 ) − µ − E Y (s0 ) = µ − µ = 0
i=1 | {z } | {z }
=0 =µ

5.2.1 Minimisation de la variance de l’erreur de prévision

Après avoir calculé l’espérance de l’erreur de prévision ,on peut aussi expri-
mer
" sa variance : #    2
2
V Ŷ (s0 ) − Y (s0 ) = E (Ŷ (s0 ) − Y (s0 )) − E(Ŷ (s0 ) − Y (s0 ))
| {z }
=0

23
 " n #
X
2
=E ( λi (Ŷ (xi ) − µ)) + (µ − Y (s0 )))
i=1

n n " # n " # "

X X X
= λi λj E (Y (xi ) − µ)(Y (xj ) − µ) − 2 λi E (Y (xi ) − µ)(Y (x0 ) − µ) + E (Y (x0 ) − µ)2
i=1 j=1 i=1

n n " # n " # " #

X X X
2 2
= λi λj E (Y (xi )(Y (xj ) − µ −2 λi E (Y (xi )(Y (x0 ) − µ + V x0
i=1 j=1 | {z } i=1 | {z }
=C(Y (xi ),Y (xj )) =C(Y (xi ),Y (x0 ))

n n n
X X X
= λi λj C(Y (xi ), Y (xj )) − 2 λi C(Y (xi ), Y (x0 )) + C(0)
i=1 j=1 i=1

Grâce à la stationnarité de la variable aléatoire Y ,on peut en déduire cette

expression :

" # n n n
X X X
V Ŷ (s0 ) − Y (s0 ) = λi λj C(xi − xj ) − 2 λi C(xi − x0 ) + C(0)
i=1 j=1 i=1

Finalement ,on peut réécrire la variance de l’erreur de prévision sous la forme

" #
X
2
V Ŷ (s0 ) − Y (s0 ) = σE = C(0) + t λ λ − 2λt S0 ≥ 0

P
où est la matrice de variance-covariance et S0 = C(xi − x0 )

Rappel : Soit a un vecteur a ∈ Rk et M une matrice ∈ Rk∗k :

t
∂ aM a
Alors ∂a = (M + t M )a .
t
∂ aM a
En particulier,si M est symétrique alors ∂a = 2M a.

24
 Pour minimiser la variance de l’erreur de prévision , il suffit d’observer pour
quelle valeur de λ s’annule la dérivée de la variance de l’erreur de prévision ,c’est
2
∂σE
à dire ∂λ
= 0.
P
Comme la matrice est symétrique ,on peut écrire :
" #
∂ C(0) + t λ λ − 2λt S0
P
" #
2
∂σE X
t
X X
= = + λ − 2S0 = 2 λ − 2S0 = 0
∂λ ∂λ

P
Ce qui équivaut à λ = S0
P P
La dérivée seconde est égale à 2 .La matrice hessienne étant semi-
définie et positive ,ce qui assure la convexité de la variance de l’erreur de prévi-
P −1
sion .On peut en conclure que la valeur λ = S0 est le point critique qui
minimise la variance de l’erreur de prévision .

Rappelons la prévision du krigeage simple :

n
X
Ŷ = µ + λ(Y − µ) = µ + λt (Y − µ)
i=1

Ceci devient quand on remplace λ par sa valeur

X t   X  
−1 t −1
Ŷ = µ + S0 Y − µ = µ + S0 Y −µ

En tenant compte aussi la nouvelle valeur de λ dans l’expression de la va-

riance de l’erreur de prévision,on obtient une nouvelle expression,c’est ce qu’on
appelle la variance du krigeage simple .Elle se calcule de la manière sui-
vante :

Ce qu’on sait :
" #
X
2
V Ŷ (s0 ) − Y (s0 ) = σE (s0 ) = C(0) + λt λ − 2λt S0

25
Ce qui équivaut à

X X X
2 −1 −1 −1
σE (s0 ) = C(0) + S0t S0 − 2S0t S0 = σ 2 − S0t S0

Par conséquent , la variance du krigeage simple en un site s0 géolocalisé

concernant sa valeur predite par rapport à sa valeur observée est donnée par
l’expression suivante :

X
2 −1
σKrigeageSimple = σ 2 − S0t S0

5.3 Le krigeage ordinaire

Le krigeage ordinaire appelé aussi krigeage ponctuel est le plus fréquemment

utilisé car il généralise le krigeage simple c’est à dire cette fois-ci l’espérance de
la variable stationnaire semble inconnue .
Le modèle de cette méthode est une combinaison lineaire de type :

n
X
Ŷi = λi Y i
i=1

avec Y une variable aléatoire stationnaire-ment intrinsèque et les λi tel que

Pn
i=1 λi = 1 sont les poids qui nous permettent que l’espérance et la variance

de l’estimateur du krigeage existent et qui le rendent Ŷ sans biais .

L’espérance de l’erreur de prévision du krigeage ordinaire

Pn
Puisque nous avons supposé que i=1 λi = 1 ,alors on a :
" # " n n
# n
" #
X X X
E Ŷ (x0 ) − Y (x0 ) = E λi Y (xi ) − Y (x0 ) λi = λi E Y (xi ) − Y (x0 ) = 0
i=1 i=1 i=1
| {z } | {z }
=1 =0

26
 La variance de l’erreur de prévision du krigeage ordinaire

Petite astuce :
Tout abord on va débuter une petite
 astuce pour
 calculer la variance de l’erreur
de prévision du krigeage : ∀i, j = 1, ......, n ,on peut écrire :

" #2 " #2

n X
n
Y (xi ) − Y (xj ) n
Y (xi ) − Y (x0 )
X X
− λi λj +2 λi
i=1 j=1
2 i=1
2

n n n X
n n n
X (Y (xi ))2 X X X (Y (xj ))2 X
=− λi λj + λi λj Y (xi )Y (xj ) − λi
i=1
2 j=1 i=1 j=1 j=1
2 i=1
| {z } | {z }
=1 =1

n
X n
X n
X
+ λi (Y (xi )2 − 2 λi Y (xi )Y (x0 ) + (Y (x0 ))2 λi
i=1 i=1 i=1
| {z }
=1

X n
n X n
X
= λi λj Y (xi )Y (xj ) − 2 λi Y (xi )Y (x0 ) + (Y (x0 ))2
i=1 j=1 i=1

Nous aurons donc cette égalité :

n X
X n n
X
λi λj Y (xi )Y (xj ) − 2 λi Y (xi )Y (x0 ) + (Y (x0 ))2
i=1 j=1 i=1

est donc équivalent à

" #2 " #2

n X
n
Y (xi ) − Y (xj ) n
Y (xi ) − Y (x0 )
X X
− λi λj +2 λi
i=1 j=1
2 i=1
2

27
 Maintenant ,on peut commencer à calculer la variance de l’erreur de prévision
du krigeage ordinaire :

"  #   2 " n !2 #
X
2
V Ŷ (s0 ) − Y (s0 ) = E (Ŷ (s0 ) − Y (s0 )) − E(Ŷ (s0 ) − Y (s0 )) = E λi Y (xi ) − Y (x0 )
| {z } i=1
=0

Grâce à l’égalité obtenue du petite astuce détaillée au dessus,on peut donc écrire

" n
!2 # n X
n   n   
X X X
E λi Y (xi ) − Y (x0 ) = λi λj E Y (xi )Y (xj ) − 2 λi E Y (xi )Y (x0 ) + E (Y (x0 ))2
i=1 i=1 j=1 i=1

A partir de cette expression,on peut faire apparaître le variogramme de la va-

riable stationnaire Y
" !2 # " !2 #

n X
n
E Y (xi ) − Y (xj ) n
E Y (xi ) − Y (x0 )
X X
=− λi λ j +2 λi
i=1 j=1 | 2
{z } i=1 | 2
{z }
Γ(xi −xj ) Γ(xi −x0 )

Finalement,la variance de l’erreur de prévision du krigeage ordinaire est don-

née par :

n X
X n n
X
2
σE =− λi λj Γ(xi − xj ) + 2 λi Γ(xi − x0 )
i=1 j=1 i=1

28
CHAPITRE 6
TESTS D’AUTOCORRELATION SPATIALE

L’auto-corrélation mesure la corrélation d’une variable avec elle-même, lorsque

les observations sont considérées avec un décalage dans le temps (autocorréla-
tion temporelle) ou dans l’espace (autocorrélation spatiale). Si la présence d’une
qualité dans une partie d’un territoire rend sa présence dans les zones voisines
plus ou moins probable, il existe un effet de contiguïté dans la structure spatiale,
le phénomène montre une autocorrélation spatiale. Lorsque l’autocorrélation est
positive, des régions voisines tendent à avoir des propriétés identiques ou des
valeurs semblables : de telles configurations sont par exemple des régions homo-
gènes*, ou bien des gradients réguliers. Lorsque l’autocorrélation est négative,
des régions voisines ont des qualités différentes, ou bien les valeurs fortes et
faibles alternent.

6.1 Indice de Moran

L’indice de Moran (Moran 1950) permet de mesurer le niveau d’autocorré-

lation spatiale d’une variable et de tester sa significativité. Il est égal au ratio

29
de la covariance entre observations contiguës (définies par la matrice d’interac-
tions spatiales) à la variance totale de l’échantillon (Jayet, 2001). L’indice a des
valeurs comprises entre -1 (indiquant une dispersion parfaite) à 1 (corrélation
parfaite). Une valeur nulle signifie que la distribution spatiale de la variable étu-
diée est parfaitement aléatoire dans le territoire. Les valeurs négatives (positives)
de l’indice indiquent une autocorrélation spatiale négative (positive).

La formulation de l’indice de Moran est

Pn Pn
n i=1 j=1 Φij [(Yi − Ȳ ) − (Yj − Ȳ )]
IM ORAN = Pn
m i=1 [(Yi − Ȳ )]

— "x" et "y" sont réservés aux coordonnées des entités géographiques.

— Yi = valeur de la variable pour l’entité "i" , sa moyenne est Ȳ
— i = entité géographique
— j = voisins des entités "i".
— n = nombre total d’entités géographiques dans l’échantillon
— m = nombre total de paires de voisins
— Φij = matrice de pondération, dont les éléments prennent, par exemple,
la valeur "1" pour les "i,j" voisins et "0" autrement.
Indice de Moran sur R
Sur le logiciel R ,on peut calculer l’indice de Moran et faire des tests
d’autocorrelation spatiale :

30
6.2 indice de Geary
Une autre alternative pour mesurer l’autocorrélation spatiale est l’indice
de Geary qui est, à un facteur près, égal au ratio de la variance des écarts
entre observations contigües à la variance totale (Jayet, 2001). Si l’indice
de Moran est une mesure de l’autocorrélation spatiale globale, l’indice
de Geary est plus sensible à l’autocorrélation spatiale locale. L’indice de
Geary varie de 0 à l’infini et vaut 1 s’il y a indépendance spatiale. Pour
toute valeur inférieur à l’unité, il y a une autocorrélation spatiale positive,
et inversement pour des valeurs supérieures à l’unité. L’indice de Geary
varie en sens inverse de l’indice de Moran. Comme pour l’indice de Moran,
l’indice de Geary peut être testé statistiquement, par une transformation
en Z- scores et en déterminant son seuil de significativité.

31
 La formulation de l’indice de Geary est

Pn Pn
n−1 i=1 j=1 Φij [(Yi − Yj )2 ]
IGEARY = Pn
2m i=1 [(Yi − Ȳ )2 ]

— "x" et "y" sont réservés aux coordonnées des entités géographiques.

— Yi = valeur de la variable pour l’entité "i" , sa moyenne est Ȳ
— i = entité géographique
— j = voisins des entités "i".
— n = nombre total d’entités géographiques dans l’échantillon
— m = nombre total de paires de voisins
— Φij = matrice de pondération, dont les éléments prennent, par exemple,
la valeur "1" pour les "i,j" voisins et "0" autrement.
Sur le logiciel R ,on peut calculer l’indice de Geary et faire des tests d’autocor-
relation spatiale :

32
CHAPITRE 7

MODÉLISATION SPATIALE PAR KRIGEAGE SUR R

Dans cette partie , nous allons utiliser le logiciel R ,des fonctions permettant
d’effectuer une analyse variographique ainsi que le krigeage.

7.1 The meuse data set

Avant de faire une analyse variographique ,nos données géographiques s’ap-

pelle "Meuse data set" .Ces données regroupes plusieurs variables :

33
 x et y sont les coordonnées géographiques de nos points sur l’espace . La
variable zinc représente la concentration en zinc mesuré à la surface du sol dans
une plaine inondable le long de la rivière Meuse . Ces métaux polluent des
sédiments qui sont transportés par le fleuve et se déposent principalement près
de la rive du fleuve et des zones de faible altitude.

Notre objectif est d’analyser ces données spatiales afin de connaître les en-
droits en forte concentration de zinc.

Lorsqu’on fait une analyse préalable sur la concentration en zinc selon les
cinq premiers coordonnées x et y on obtient :

34
7.2 Analyse variographique

Cette étape nous permet de calculer et tracer le variogramme estimé à l’aide

de nos observations.

7.2.1 Le variogramme expérimental

La fonction R permettant de calculer le variogramme est la suivante :

>variogram(formule,data,........)
˙

cette fonction prend des arguments de types :

— formule : variabler ponse.
— data : Objet de typeataframe où se trouve nos variables utilisées.

35
— cutoff : Distance maximale entre deux points au-delà de laquelle le va-
riogramme n’est plus calculé.
— cloud : TRUE/FALSE. Si True, cette option calcule la nuée variogra-
phique.
Et elle renvoie ces valeurs suivantes :
Si le cloud = TRUE ,elle renvoie un objet de type "VariogramCloud".Sinon
elle renvoie :
— np : Nombre de paires de points utilisées pour chaque point du vario-
gramme.
— dist :Distance moyenne entre les points d’une paire pour chaque point
du variogramme.
— gamma : Donne les valeurs calculées du variogramme.
Lorsqu’on finit de calculer le variogramme expérimental , il suffit de le
représenter graphiquement. Cette fonction plot(x,...) permet de tracer
le graphique .Elle prend des arguments de types :
— type : Permet de tracer différents types de variogramme.
— identify :TRUE/FALSE : Si on a calculé la nuée variographique (c’est-
à-dire si on ajouté l’option cloud=TRUE dans la fonction variogram()),
cette option permet d’identifier pour chaque point du nuage, quelle paire
de points correspond à ce point.
La variable aléatoire zinc n’étant pas stationnaire ,sa transformation lo-
garithmique est stationnaire c’est à dire log(zinc) .

36
Et sa représentation graphique par la fonction plot est :

37
7.2.2 Modélisation du variogramme

Après avoir calculé le variogramme expérimental, l’étape suivante constitue

à modéliser ce variogramme à l’aide de modèles connus afin de pouvoir lui appli-
quer le krigeage. Pour cela, on a besoin de deux fonctions du logiciel R : vgm(...
) et [Link](...).

La fonction vgm(psill,model,..... ) permet de générer un modèle de va-

[Link] renvoie un objet de type «variogramModel». Les arguments de
cette fonction sont les suivants :

— psill : Palier (ou palier asymptotique) du modèle.

— model : Type de modèle (“Exp”,”Sph”,...).
— range : Portée du modèle.
— range : effet pépite du modèle.
La fonction [Link](object,model,...) permet d’ajuster le va-
riogramme expérimental par un modèle de variogramme .Elle prend des
arguments de types :
— object : Variogramme expérimental de type «[Link]» (issu de
la fonction variogram(...)).
— model : : Modèle de variogramme de type «variogramModel» (issu de
la fonction vgm(...)) .
— [Link] :TRUE/FALSE→détermine si le palier du modèle doit être ajusté.
— [Link] :TRUE/FALSE→détermine si le portée du modèle doit être
ajusté
Lorsque l’on termine d’ajuster le variogramme expérimental, on trace ce
dernier avec le variogramme ajusté à l’aide la commande plot(variogramme
expérimental,variogramme ajusté) .

38
On peut représenter graphiquement ensemble le variogramme expérimen-
tal et sa variogramme ajusté (ici le modèle utilisé est le modèle sphérique
).

39
7.3 Le krigeage avec R

Pour pouvoir kriger la variable régionalisée zinc avec R, il faut tout d’abord
créer une grille sur laquelle on pourra effectuer le krigeage.

Nous allons alors créer une grille sur R :

Une fois la grille créée, on peut effectuer le krigeage sur nos données à l’aide
de la fonction > krige(formula, data,...)

Les arguments de cette fonction sont :

— formula : V ariabler ponseouvariableexpliquer.
— model : Variogramme modélisé du type «variogramModel» (issu de la
fonction [Link]).
— data : Objet de type dataframe où se trouve nos variables utilisées.
— beta : Si on effectue un krigeage simple, beta correspond à espérance
de la variable régionalisée étudiée.
Cette fonction renvoie un objet de type «SpatialPixelsDataFrame» avec

40
 les données suivantes :
— les coordonnées des points qui ont été krigés.
— les valeurs prédites en ces points ([Link]).
Pour obtenir une visualisation du krigeage effectué sur nos données, on
utilise la fonction > spplot(obj, ...).
Cette fonction des arguments de type suivants :
— obj : : Objet que l’on veut représenter de type «Spatial».
— xlab :Nommer l’axe x.
— ylab : Nommer l’axe x.
Finalement ,on voici le krigeage qu’on a fait à partir de nos [Link]
figure montre la concentration en zinc sur le sol de la fleuve Meuse de la plus
faible au plus fort.

41
ANNEXE A
LE KRIGEAGE SUR LE LANGAGE R

42
43
44
45
 ANNEXE B
BIBLIOGRAPHIE

— Spatial Statistics and Geostatistics Theory and Applications

for Geographic Information Science and Technology-sagepub
(2013), Yongwan Chun, Daniel A Griffith
— Basic Linear Geostatistics-Springer-Verlag Berlin Heidelberg
(1998), Dr. Margaret Armstrong (auth.)
— Spatial and Spatio-Temporal Geostatistical Modeling and Kriging-
Wiley (2015), José-Maria Montero, Gema Fernanndez-Avilés, Jorge
Mateu
— Modélisation et statistique spatiales (Mathématiques et Applications)-
Springer (2008), Carlo Gaetan, Xavier Guyon

46

View publication stats