Le Krigeage PDF
Le Krigeage PDF
Le krigeage :
revue de la théorie et application à l’interpolation spatiale de
données de précipitations
Mémoire présenté
à la Faculté des études supérieures de l’Université Laval
dans le cadre du programme de maı̂trise en statistique
pour l’obtention du grade de Maı̂tre ès sciences (M.Sc.)
Avril 2005
c
°Sophie Baillargeon, 2005
Résumé
Je tiens tout d’abord à remercier les personnes qui ont contribué, de près ou de loin,
à la réalisation de ce mémoire. Premièrement, merci à mon directeur de recherche, mon-
sieur Louis-Paul Rivest, professeur au Département de mathématiques et statistique de
l’Université Laval. Il s’est toujours montré disponible et intéressé par mon projet. Il
a su me diriger tout en me laissant une grande liberté dans mes travaux, ce que j’ai
beaucoup apprécié. Merci aussi à ma co-directrice de recherche, madame Jacynthe Pou-
liot, professeure au Département des sciences géomatiques de l’Université Laval. Elle
a suscité mon intérêt pour la géostatistique et m’a proposé ce projet. Ses précieux
conseils de rédaction m’ont également permis d’écrire un document qui, je l’espère,
pourrait intéresser autant un géomaticien qu’un statisticien. De plus, je remercie mon-
sieur Vincent Fortin, chercheur à l’Institut de recherche d’Hydro-Québec. C’est grâce à
lui que la problématique d’interpolation de données de précipitations dans un cadre de
modélisation hydrologique a pu être traitée. Il a collaboré activement à mon mémoire
en fournissant temps, conseils, données et en acceptant d’en être un examinateur.
D’autres personnes ont été impliquées dans ce projet, notamment madame Josée
Fitzback, étudiante au doctorat en sciences géomatiques de l’Université Laval. Je la
remercie d’avoir établi les contacts qui ont rendu le projet possible. Travailler avec elle
fut toujours un plaisir. Merci également à messieurs Jean-François Mahfouf et Pierre
Pellerin, météorologues à Environnement Canada, pour l’intérêt qu’ils ont porté envers
le projet. J’exprime de plus ma gratitude envers le CRSNG, qui a financé mes travaux
de recherche en m’octroyant une bourse d’étude supérieure.
En finissant, je désire bien sûr remercier les membres de ma famille, mes parents et
mes soeurs, qui ont toujours encouragé mes études et grâce à qui j’ai grandi dans un
milieu épanouissant. Mes derniers remerciements s’adressent à mon amoureux, Mathieu,
qui est toujours un stimulant à aller de l’avant. Il m’a écouté avec intérêt parler des
hauts et des bas de ma rédaction, m’a aidé à garder le sourire en temps de stress et a
répondu à mes nombreuses questions LATEX. Je le remercie simplement d’avoir été là.
Table des matières
Résumé ii
Avant-Propos iii
1 Introduction 1
1.1 Objectifs du mémoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Structure du mémoire . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Interpolation spatiale 3
2.1 Définition et notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Revue des méthodes d’interpolation spatiale . . . . . . . . . . . . . . . 5
2.2.1 Méthodes barycentriques . . . . . . . . . . . . . . . . . . . . . . 5
2.2.2 Méthodes d’interpolation par partitionnement de l’espace . . . . 6
2.2.3 Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.4 Régression classique . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.5 Régression locale . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.6 Krigeage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.7 Autres méthodes . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Conclusion du chapitre . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3 Analyse variographique 16
3.1 Hypothèse de stationnarité . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Décomposition de la variation spatiale . . . . . . . . . . . . . . . . . . 18
3.3 Propriétés du semi-variogramme . . . . . . . . . . . . . . . . . . . . . . 19
3.4 Estimation du semi-variogramme . . . . . . . . . . . . . . . . . . . . . 21
3.5 Modélisation du semi-variogramme . . . . . . . . . . . . . . . . . . . . 23
3.6 Conclusion du chapitre . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
v
4 Théorie du krigeage 27
4.1 Démarche générale de résolution des équations du krigeage . . . . . . . 27
4.2 Krigeage simple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.3 Krigeage ordinaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.4 Krigeage avec modèle de tendance . . . . . . . . . . . . . . . . . . . . . 36
4.4.1 Lien entre le krigeage avec modèle de tendance et le krigeage sur
les résidus d’une régression . . . . . . . . . . . . . . . . . . . . . 39
4.4.2 Problème de l’analyse variographique en krigeage avec modèle de
tendance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.5 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.5.1 Normalité des données . . . . . . . . . . . . . . . . . . . . . . . 41
4.5.2 Transformation de données . . . . . . . . . . . . . . . . . . . . . 43
4.5.3 Géostatistique multivariable . . . . . . . . . . . . . . . . . . . . 45
4.5.4 Autres types de krigeage . . . . . . . . . . . . . . . . . . . . . . 48
4.6 Conclusion du chapitre . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
7 Conclusion 89
7.1 Littérature en interpolation de données de précipitations . . . . . . . . 89
7.2 Synthèse du mémoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
7.3 Travaux futurs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
Bibliographie 94
Introduction
Des données à référence spatiale sont de plus en plus exploitées, et ce, dans di-
vers domaines de recherche. Qu’il s’agisse de précipitations mesurées à des stations
météorologiques, de la densité d’un minerai dans des échantillons de sol ou de la concen-
tration de gaz carbonique dans l’air en certains sites, ces données possèdent toutes un
point en commun : elles sont localisées dans l’espace géographique. Le traitement statis-
tique de ce type de données demande une attention particulière car l’hypothèse classique
selon laquelle les observations sont indépendantes et identiquement distribuées est ra-
rement vérifiée. Des méthodes statistiques adaptées à l’analyse de données à référence
spatiale ont été développées (Ripley, 1981; Cressie, 1993). Ce mémoire porte sur l’une
de ces méthodes, le krigeage, développée par Matheron (1962, 1963a,b). Le krigeage sert
à effectuer de l’interpolation spatiale, c’est-à dire qu’il permet de prévoir la valeur prise
par un phénomène naturel un site à partir d’observations ponctuelles de ce phénomène
en des sites voisins.
mise en oeuvre du krigeage. Les fondements du krigeage se basent sur une connaissance
de la structure de dépendance spatiale des données. Cependant, en pratique, cette
structure n’est pas connue. Elle est estimée lors de l’analyse variographique. Cette
étape préalable au krigeage se doit donc aussi d’être exposée. Pour compléter cette
étude du krigeage, nous nous sommes demandé si le krigeage permettait de répondre
efficacement à une problématique d’interpolation spatiale. La problématique traitée ici
concerne l’interpolation de données de précipitations dans le but de fournir ces données
en entrée à un modèle de simulation hydrologique.
Interpolation spatiale
Ce chapitre présente une revue des principales méthodes d’interpolation spatiale afin
de comprendre leurs principes de base et de comparer le krigeage aux autres méthodes.
Il ressort alors que malgré les nombreux points communs entre les méthodes, le krigeage
se distingue par sa prise en compte de la structure de dépendance spatiale des données.
En pratique, on n’a pas une valeur régionalisée en tout point s du champ D. Par
exemple, pour un gisement d’or, la densité du minerai dans le sol n’est mesurée qu’en
quelques sites où une carotte de sol est prélevée. Peu importe l’application en question,
notons ces sites d’observation si avec i = 1, ..., n. L’interpolation spatiale répond au be-
Chapitre 2. Interpolation spatiale 4
La figure 2.1 résume la notation présentée dans cette section. Elle représente du
même coup la modélisation d’un phénomène naturel selon les deux niveaux d’abstraction
décrits précédemment.
Dans cette section, les principales méthodes d’interpolation spatiale sont recensées.
Le livre de Arnaud et Emery (2000) a été le point de départ de cet inventaire. La ter-
minologie proposée dans ce livre est employée ici. Trois grandes classes de méthodes
déterministes ont été dénombrées : les méthodes barycentriques, les méthodes d’inter-
polation par partitionnement de l’espace et les splines. Du côté des méthodes stochas-
tiques, les techniques de régression classique, de régression locale et de krigeage ont
été listées. Contrairement aux méthodes déterministes, les méthodes stochastiques in-
corporent le concept de hasard. Elles proposent toutes un modèle probabiliste incluant
un ou des termes d’erreurs aléatoires pour formaliser le comportement du phénomène
naturel à l’étude. Grâce à cette modélisation, des erreurs de prévision peuvent être cal-
culées. Les paragraphes suivants proposent une brève introduction à chacune de ces six
classes de méthodes.
Les poids λi sont contraints de sommer à la valeur 1 afin que la prévision ne présente
pas de distorsion par rapport à la valeur réelle. Ces poids sont fonction de la distance
euclidienne |si − s0 | entre le site d’observation si et le site de prévision s0 de façon à
ce que les sites les plus proches aient plus d’influence dans l’interpolation. Souvent, un
poids nul est accordé aux observations les plus éloignées. Ainsi, seules les observations
localisées dans un certain voisinage de s0 , noté V (s0 ), sont prises en compte. La façon la
plus simple de déterminer un voisinage est de prendre les n0 sites d’observation les plus
proches ou les sites tombant à l’intérieur d’un cercle centré en s0 de rayon prédéterminé.
Le champ peut aussi être divisé en quadrants ou en octants dont l’origine est s0 . Le
voisinage peut alors comprendre les n∗0 sites d’observation le plus proche de s0 pour
chaque secteur.
valeurs régionalisées observées. La figure 2.3 illustre comment le poids de chacune des
observations est déterminé. Préalablement, les polygones de Thiessen associés aux sites
d’observation sont tracés (e.g. mosaı̈que de droite de la figure 2.3). Ensuite, les polygones
de Thiessen sont reformés en ajoutant au champ le site s0 pour lequel une prévision est
voulue (e.g. mosaı̈que centrale de la figure 2.3). Cette nouvelle mosaı̈que est finalement
superposée aux polygones de Thiessen initiaux sans le site s0 (e.g. mosaı̈que de gauche
de la figure 2.3). Le poids de l’observation en si est alors l’aire Ai de l’intersection entre
le polygone de s0 et le polygone initial de si divisée par l’aire totale du polygone de
P P
s0 . Dans l’exemple de la figure 2.3, la prévision est donc ẑ(s0 ) = 7i=1 Ai z(si )/ 7i=1 Ai
avec A3 , A4 = 0.
critère de voisinage et que les prévisions qu’elles fournissent aux points d’observation
sont égales aux observations. Ces méthodes sont déterministes car elles ne modélisent
pas la variable régionalisée par une fonction aléatoire. Notons que les méthodes d’in-
terpolation par partitionnement fournies en exemple ici sont les plus simplistes. Elles
produisent des surfaces d’interpolation présentant, sauf exception, des changements
abruptes. En outre, l’interpolation linéaire ne permet pas d’extrapoler en dehors de
l’enveloppe convexe des sites. Plusieurs autres méthodes d’interpolation par partition-
nement ne présentent pas ces inconvénients.
Chapitre 2. Interpolation spatiale 9
2.2.3 Splines
où ρ est un paramètre de lissage fixé à priori. L’intégrale incluse dans ce critère à mi-
nimiser est appelée « énergie de flexion » . Elle réfère au concept physique d’énergie de
flexion d’une mince plaque de métal. Les surfaces générées par ces splines peuvent donc
être comparées à des plaques de métal flexible ajustées de façon à minimiser leurs degrés
de courbure tout en passant à proximité des points d’observation. D’autres types de
splines permettent d’effectuer de l’interpolation spatiale, notamment les fonctions mul-
tiquadratiques (Hardy, 1971), aussi appelées fonctions radiales de base (Myers, 1994),
et les splines régularisées avec tension (Mitášová et Mitáš, 1993).
Tout comme les splines, la régression classique permet d’effectuer une interpolation
en ajustant une surface aux valeurs régionalisées observées. Cependant, la régression
est une méthode stochastique. Elle suppose que la variable régionalisée à l’étude est
une fonction aléatoire qui se décompose comme suit :
où µ(·) est la structure déterministe pour l’espérance fonction de la localisation des
observations et ²(·) est une fonction aléatoire normale d’espérance nulle, de variance
homogène et ne présentant pas de structure de dépendance spatiale. Ainsi, ²(·) est un
processus gaussien de bruit blanc représentant des erreurs de mesure indépendantes. La
tendance µ(s) peut prendre plusieurs formes. La plus utilisée est un polynôme de degré
d des coordonnées : X
µ(s) = µ(x, y) = βlk xl y k . (2.2)
l+k≤d
Notons que l’interpolation par régression classique est considérée déterministe par
certains auteurs (Arnaud et Emery, 2000, p.73 ; Burrough et McDonnell, 1998, p.158).
D’ailleurs, l’emploi de la méthode des moindres carrés rend superflue la modélisation
stochastique de la variable aléatoire. Par contre, cette modélisation, accompagnée des
hypothèses de normalité, d’indépendance et d’homoscédasticité des erreurs du modèle,
permet de calculer des tests de significativité de la surface et des erreurs de prévision.
Chapitre 2. Interpolation spatiale 11
Cependant, ces statistiques ne sont presque jamais fiables en interpolation spatiale car
le postulat d’indépendance est rarement vérifié sur des données spatiales. Malgré tout,
l’interpolation par régression classique est ici considérée stochastique, comme le font
Klinkenberg et Waters (1990), car cette technique fut à l’origine développée dans un
cadre statistique.
De plus, l’interpolation par régression classique possède les propriétés d’être ap-
proximative et globale. La surface d’interpolation qu’elle génère ne passe donc pas
nécessairement par les points d’observation et toutes les observations, mêmes celles
éloignées, vont influencer avec le même poids la prévision en n’importe quel point du
champ. En raison de ces caractéristiques, la méthode produit des surfaces plutôt lisses.
où ẑ(si )|s0 est la valeur prédite de z(·) au point d’observation si à partir de la surface
ajustée en s0 et z(si ) est la valeur régionalisée observée en si . Le poids de l’observation
au point si prend typiquement la forme :
µ ¶
|s0 − si |
λi (s0 ) = K
h(s0 )
où K(·) est une fonction de poids, |s0 − si | la distance euclidienne entre le point de
prévision s0 et le point d’observation si et h(·) est la taille du voisinage. Plusieurs
fonctions de poids sont proposées dans la littérature (Loader, 1999, p.48 ; Wand et
Jones, 1995, p.31). Par exemple, une de ces fonctions est la fonction Epanechnikov
définie par K(u) = 1 − u2 pour 0 ≤ u < 1, et 0 sinon. En outre, le paramètre h(s0 ) est
habituellement déterminé en choisissant d’abord la fraction α des points d’observation à
inclure dans l’ajustement du modèle. Ensuite, on attribue à h(s0 ) la valeur de la distance
Chapitre 2. Interpolation spatiale 12
entre s0 et le [nα]ième plus proche point d’observation de s0 (Loader, 1999, p.20), où
[nα] représente la valeur entière du produit nα. Ainsi, plus une valeur régionalisée à été
mesurée loin du point de prévision s0 , moins elle a d’impact sur l’interpolation en s0 .
Au-delà d’une certaine distance elle n’a plus d’impact du tout.
Afin de retomber sur une régression classique en partant d’une régression locale, il
suffit de prendre une fonction de poids uniforme et une fraction de voisinage égale 1.
De plus, une regression locale avec une tendance polynomiale de degré 0 ajustée par la
méthode des moindres carrés pondérés revient à une méthode barycentrique (Loader,
2004). Dans ce cas, la seule distinction entre les deux méthodes est que la modélisation
stochastique de la régression locale permet le calcul d’erreurs de prévision. Cependant,
pour la même raison qu’en régression classique (voir section 2.2.4), il faut douter de la
fiabilité de ces erreurs.
2.2.6 Krigeage
Les deux méthodes stochastiques présentées dans les sections précédentes ne pren-
nent pas en considération la structure de dépendance spatiale des données. Le krigeage
est la première méthode d’interpolation spatiale à en avoir tenu compte. Les travaux de
l’ingénieur minier sud-africain Krige (1951) sont précurseurs de la méthode. Cependant,
le terme krigeage et le formalisme de cette méthode sont dus au français Matheron
(1962, 1963a,b), qui en a aussi assuré le développement au Centre de Géostatistique de
l’École des Mines de Paris. En fait, les fondements de la méthode ont été développés
parallèlement par d’autres chercheurs, notamment le météorologue Gandin (1963) de
l’ex-URSS, mais c’est aujourd’hui sous la terminologie proposée par Matheron qu’elle
Chapitre 2. Interpolation spatiale 13
Les poids λi associés à chacune des valeurs régionalisées observées sont choisis de façon
à obtenir une prévision non biaisée et de variance minimale. Ces poids dépendent de
la localisation des observations et de leur structure de dépendance spatiale. En fait, le
krigeage est le nom donné à la meilleure prévision linéaire sans biais, en anglais « best
linear unbiased predictor » ou « BLUP », dans un cadre spatial.
Le krigeage est une méthode d’interpolation très souple. Il peut être global ou local
dépendamment du voisinage choisi. De plus, selon le développement classique du kri-
geage qui sera suivi dans ce mémoire, il s’agit d’une méthode d’interpolation exacte.
Chapitre 2. Interpolation spatiale 14
Il restitue donc les valeurs régionalisées mesurées aux sites d’observation. Cependant,
il est aussi possible d’effectuer un krigeage dit avec erreurs de mesure (Cressie, 1993,
p.128) qui est approximatif. Finalement, comme en régression classique ou locale, des
variables régionalisées auxiliaires peuvent être intégrées au krigeage en les ajoutant à
la tendance générale µ(·) du modèle. Le krigeage incorporant une telle tendance est
nommé « krigeage avec dérive externe » (Goovaerts, 1997). Le krigeage possède aussi
d’autres extensions multivariables, notamment le « cokrigeage » (Wackernagel, 2003).
Par cette méthode, ẑ(s0 ) prend la forme d’une combinaison linéaire pondérée des obser-
vations de la variable régionalisée à interpoler et des variables régionalisées auxiliaires
notées {wj (s), s ∈ D} avec j = 1, ..., q :
q
X X X
ẑ(s0 ) = a + λi,0 z(si ) + λi,j wj (si,j ).
i∈V (s0 ) j=i i∈V (s0 )
Les sections précédentes ont introduit les méthodes d’interpolation les plus connues
et utilisées. Mentionnons tout de même ici quelques-unes des autres méthodes exis-
tantes. Certaines faiblesses du krigeage ont mené au développement de méthodes bayé-
siennes (Gaudard et al., 1999; Banerjee et al., 2004). Ces faiblesses sont principalement
liées à la présence de biais dans l’estimation du variogramme en krigeage universel
(Cressie, 1993, p.165) et au fait que l’incertitude d’une prévision calculée en krigeage
n’incorpore pas l’incertitude associée à l’estimation du variogramme. Le modèle bayésien
typique en interpolation spatiale effectue un krigeage avec une fonction de covariance
de paramètres inconnus. Ces modèles peuvent aussi permettre d’intégrer à l’interpo-
lation de l’information à priori concernant la tendance µ(·). Un des grands avantages
de l’approche bayésienne est la perspective conditionnelle qui permet un raisonnement
probabiliste, après observation des données, sans recourir à l’idée d’une répétition de
l’expérience, condition souvent difficile à concevoir en géostatistique. Certains utilisent
aussi les réseaux de neurones artificiels pour interpoler des données (Bryan et Adams,
2002; Rigol et al., 2001). Les réseaux de neurones agissent alors comme une méthode
de régression non paramétrique. Dans ce mémoire, aucune de ces méthodes ne sera
approfondie.
Chapitre 2. Interpolation spatiale 15
Analyse variographique
Le processus naturel étudié étant unique, une seule réalisation de la fonction aléatoire
Z(·) est observable. Afin de rendre possible l’inférence statistique malgré l’unicité de la
réalisation disponible, une hypothèse de stationnarité est émise concernant la fonction
aléatoire δ(·). Au sens strict, la stationnarité signifie que la loi de probabilité de la
fonction aléatoire est invariante par translation, c’est-à-dire qu’elle ne dépend pas de
l’origine du champ. Par contre, en krigeage, la stationnarité postulée est faible : elle ne
concernent que les moments d’ordre 1 et 2 de la fonction aléatoire ou de ses accroisse-
ments plutôt que sa distribution entière. La théorie du krigeage peut être développée en
postulant la stationnarité de second ordre ou la stationnarité de second ordre intrinsèque
(plus simplement nommée stationnarité intrinsèque) de δ(·) (Cressie, 1993, p.40 et 53 ;
Arnaud et Emery, 2000, p.106). Ces deux types de stationnarité se définissent ainsi :
Chapitre 3. Analyse variographique 17
Cette généralité est utile en krigeage car certains phénomènes régionalisés présentent
une dispersion infinie. C’est le cas des gisements d’or selon Krige (1951). Postuler
l’existence de la variance de la fonction aléatoire δ(·) est donc parfois inadéquat. Pour
cette raison, l’hypothèse de stationnarité intrinsèque sera privilégiée dans ce mémoire,
tel que le font la majorité des auteurs sur le sujet. Ainsi, une dépendance spatiale
sera représentée par un variogramme plutôt qu’un covariogramme. D’ailleurs, le vario-
gramme est plus facile à estimer car il ne requiert pas d’estimation de l’espérance de
δ(·). Pour ces raisons, le variogramme est l’outil privilégié en analyse variographique.
En fait, la plupart des auteurs travaillent avec γ(·), la demie du variogramme, appelée
semi-variogramme. C’est aussi ce qui sera fait dans le reste de ce mémoire.
Chapitre 3. Analyse variographique 18
où
– µ(·) = E(Z(·)) : variation à grande échelle, structure déterministe pour l’espérance
de Z(·).
– ω(·) : variation lisse à petite échelle (échelle plus grande que la distance minimale
entre deux sites d’échantillonnage), structure stochastique de fluctuations autour
de µ(·) dépendantes spatialement.
– η(·) : variation micro-échelle (échelle plus petite que la distance minimale entre
deux sites d’échantillonnage), structure stochastique présentant une dépendance
spatiale.
– ²(·) : erreur de mesure, structure stochastique sans dépendance spatiale (bruit
blanc).
Ainsi, la variable régionalisée z(·) est supposée être composée de variations à différentes
échelles emboı̂tées les une dans les autres (Oliver, 2001). La fonction aléatoire δ(·) du
modèle 2.4 est formée par le regroupement des termes ω(·), η(·) et ²(·) du modèle
précédent.
Tel que défini à la section 3.1, le semi-variogramme est une fonction paire, i.e.
γ(h) = γ(−h), nulle en h = 0 et positive partout ailleurs. Toutes ces caractéristiques
découlent en fait d’une propriété générale des semi-variogrammes : ce sont des fonctions
de type négatif conditionnel (Christakos, 1984), c’est-à-dire que :
l X
X l
ai aj γ(si − sj ) ≤ 0
i=1 j=1
pour n’importe quel ensemble fini de points {si : i = 1, ..., l} et n’importe quels nombres
P
réels {ai : i = 1, ..., l} tel que li=1 ai = 0. Ainsi, une fonction continue peut représenter
un semi-variogramme si et seulement si elle respecte cette propriété qui assure la posi-
tivité de la variance de toute combinaison linéaire de variables aléatoires issues de δ(·).
Cette propriété est en fait l’extension aux semi-variogrammes du caractère semi-défini
positif bien connu des fonctions de covariance.
Chapitre 3. Analyse variographique 20
où |h| est la norme du vecteur de h. Ces fonctions sont tracées dans les graphiques
de la figure 3.2. Il s’agit de semi-variogrammes isotropes, d’effet de pépite valant 0.2,
de palier exact ou asymptotique valant 1 et de portée exacte ou pratique valant 6. Les
paragraphes suivants définissent ces caractéristiques des semi-variogrammes.
Isotropie
1
γ(r) = V ar(δ(s) − δ(s + h)) ∀ h de norme r et ∀ s, s + h ∈ D (3.2)
2
Chapitre 3. Analyse variographique 21
Effet de pépite
Portée et palier
Par contre, lorsque µ(·) n’est pas une fonction constante tel qu’en krigeage universel
ou en krigeage avec dérive externe, l’estimation doit se baser sur la fonction δ(·). Des
observations de δ(·) sont disponibles si la fonction µ(·) est connue. Le semi-variogramme
peut alors être estimé à partir des valeurs z(si )−µ(si ). Cependant, en pratique, µ(·) n’est
pas connue. Cette fonction est donc estimée, puis le semi-variogramme expérimental est
calculé à partir des résidus e(si ) = z(si ) − µ̂(si ). Toutefois, Matheron (1970, p.155) a
lui même soulevé que cette estimation est biaisée. Une discussion plus approfondie de
cette problématique se trouve à la section 4.4.
où N (r) = {(i, j) tel que |si − sj | = r} et |N (r)| est le nombre de paires distinctes de
l’ensemble N (r).
Si les données disponibles sont réparties sur une grille régulière couvrant le champ
D, γ(·) est estimable pour un petit nombre de distances, chacune associée à plusieurs
couples de données. Toutefois, les données sont plus souvent réparties irrégulièrement
sur D. Dans ce cas, plus de distances sont représentées, mais avec un très petit nombre de
paires de données. Afin de rendre γ̂(·) plus robuste, des tolérances sont introduites sur les
distances. Le semi-variogramme expérimental est donc calculé seulement pour certaines
distances rk , avec k = 1, ..., K, en considérant les couples de données éloignées d’une
Chapitre 3. Analyse variographique 23
distance à peu près égale à rk . L’ensemble N (rk ) est ainsi redéfini par {(i, j) tel que |si −
sj | = rk ± bk } où bk est une tolérance déterminée par l’utilisateur.
Tester le caractère négatif conditionnel d’une fonction est une tâche assez difficile.
Christakos (1984) explique une procédure pour y arriver. Cependant, en pratique, il est
plus simple d’utiliser un modèle variographique classique, qui a été démontré valide.
Voici donc certains modèles variographiques couramment utilisés (Arnaud et Emery,
2000, p.133), qui seront d’ailleurs employés dans les chapitres 5 et 6 :
Portée exacte :
Portée asymptotique :
γ(r) = c0 + mr pour r ≥ 0
où γ̂(·) est le semi-variogramme empirique, γ(·, θ) est le modèle variographique de pa-
ramètres θ, wk est le poids associé à la donnée γ̂(rk ), et r1 à rK sont les distances pour
lesquelles une estimation du semi-variogramme a été calculée. Plusieurs poids wk sont
proposés dans la littérature. Ils peuvent évidemment tous être fixés à un, ce qui revient
aux moindres carrés ordinaires. Il est aussi intuitif d’envisager les poids wk = 1/|N (rk )|,
soit l’inverse du nombre de paires de données ayant entré dans le calcul de γ̂(rk ). Fuentes
(2000) discute d’autres poids pertinents et présente différentes méthodes d’estimation,
telles la méthode du maximum de vraisemblance ainsi que des approches bayésiennes.
Chapitre 3. Analyse variographique 26
Les rudiments de l’analyse variographique ont été présentés dans ce chapitre. Nous
avons défini un outil, le semi-variogramme, pour représenter une dépendance spatiale et
nous avons montré comment estimer et modéliser cette fonction. Le chapitre 5 reviendra
sur l’analyse variographique pour expliquer comment elle s’effectue en pratique. Notam-
ment, le chapitre 5 traitera de comment choisir l’estimateur du semi-variogramme et
le modèle variographique. Mais tout d’abord, passons à l’élaboration des équations du
krigeage.
Chapitre 4
Théorie du krigeage
Par la suite, la formule de prévision par krigeage peut être trouvée. Pour effec-
tuer cette étape, la démarche proposée par Chauvet (1999, p.117) sera suivie dans ce
mémoire. Chauvet conseille de procéder par l’écriture des contraintes du krigeage que
voici :
1. Contrainte de linéarité
La contrainte de base du krigeage est que la prévision prenne la forme d’une
combinaison linéaire des données. Elle doit donc s’écrire ainsi :
X
Ẑ(s0 ) = a + λi Z(si )
i∈V (s0 )
– Γ est la matrice n0 xn0 dont l’élément (i, j) est γ(s[i] − s[j] ), soit le semi-
variogramme entre δ(s[i] ) et δ(s[j] ), les éléments i et j de δ
– γ 0 est le vecteur n0 x1 dont l’élément i est γ(s[i] − s0 ), le semi-variogramme
entre δ(s[i] ) et δ(s0 )
Notons que la prévision Ẑ(s0 ) obtenue est une variable aléatoire. Sa valeur numérique
est calculée en remplaçant les variables aléatoires Z(si ) par les valeurs régionalisées
observées z(si ).
avec m constante connue et δ(·) fonction aléatoire stationnaire de second ordre d’espé-
rance nulle et de structure de dépendance connue. La stationnarité de second ordre
implique que le semi-variogramme de δ(·) atteigne un palier.
Ce modèle peut être réécrit sous forme vectorielle en fonction de la variable aléatoire
à prévoir Z(s0 ) et de celles qui serviront à faire la prévision Z = (Z(s[1] ), Z(s[2] ), ...,
Z(s[n0 ] ))t , soient les variables aléatoires associées aux données. Notons Z ∗ = (Z(s0 ), Z t )
et δ ∗ = (δ(s0 ), δ t ), alors le modèle s’écrit :
m constante connue
E[δ ∗ ] = 0
Z ∗ = m1(n0 +1) + δ ∗ avec à ! (4.1)
2 t
∗ σ c0
V ar[δ ] = connus
c0 Σ
Toutes les informations nécessaires à l’écriture des contraintes sont comprises dans
ce modèle. Suivons donc la démarche proposée afin de trouver la formule de prévision
en krigeage simple.
Chapitre 4. Théorie du krigeage 30
1. Contrainte de linéarité
2. Contrainte d’autorisation
3. Contrainte de non-biais
La propriété d’absence de biais de la prévision doit être vérifiée. Il faut donc s’assurer
du respect de l’égalité suivante :
E[Ẑ(s0 ) − Z(s0 )] = E[a + λt Z − Z(s0 )] = a + λt m1n0 − m = 0
Cette égalité implique que a = m(1 − λt 1n0 ). La prévision peut donc être réécrite sous
la forme :
Ẑ(s0 ) = m(1 − λt 1n0 ) + λt Z = m + λt (Z − m1n0 )
4. Contrainte d’optimalité
Cette expression doit être vue comme une fonction des λi , qui peut être notée
f (λ1 , ..., λn0 ) = f (λ). Le gradient de cette fonction, c’est-à-dire le vecteur de ses dérivées
Chapitre 4. Théorie du krigeage 31
partielles, s’écrit :
∂ ∂ ¡ t ¢
f (λ) = λ Σλ + σ 2 − 2λt c0 = 2Σλ − 2c0
∂λ ∂λ
Toutes les composantes de ce vecteur sont nulles au point λ̂ = Σ−1 c0 . De plus, la
2
matrice n0 × n0 dont l’élément (i, j) est ∂λ∂j ∂λi f (λ), appelée hessien, est 2Σ. Cette
matrice est semi-définie positive car il s’agit d’une matrice de variances-covariances
multipliée par une constante positive. Ainsi, la fonction f (λ) est convexe et le point
critique λ̂ = Σ−1 c0 est un minimum global (Khuri, 1993, p.282).
Cette expression est équivalente à E[Ẑ(s0 )|Z] si Z(·) est supposé de loi normale multi-
dimentionnelle. Cette équivalence est discutée à la section 4.5.1. Une autre statistique
d’intérêt est la valeur minimale de la variance de l’erreur de prévision, appelée « variance
de krigeage ». En krigeage simple elle vaut :
signifie que l’espérance n’est pas contrainte à demeurer la même partout dans le champ
D. Elle doit par contre rester constante à l’intérieur de chaque voisinage de krigeage
(Arnaud et Emery, 2000, p.185). Ainsi, une prévision au point s0 se base sur le modèle
vectoriel suivant :
µ constante inconnue
∗ ∗
Z = µ1(n0 +1) + δ avec E[δ ∗ ] = 0
Γ, γ connus
0
Voici, étape par étape, la résolution des équations du krigeage ordinaire se basant sur
Cressie (1993, p.119).
1. Contrainte de linéarité
La prévision Ẑ(s0 ) doit être une combinaison linéaire des variables aléatoires Z(s[1] )
à Z(s[n0 ] ). Elle s’exprime donc sous la forme :
X
Ẑ(s0 ) = a + λi Z(si ) = a + λt Z
i∈V (s0 )
2. Contrainte d’autorisation
il ressort que cette erreur est une combinaison linéaire d’accroissements de δ(·) si et
P
seulement si i∈V (s0 ) λi − 1 = 0. Il faudra donc travailler par la suite avec la contrainte
que la somme des poids λi vaut un pour s’assurer de l’existence des deux premiers
moments de l’erreur de prévision.
3. Contrainte de non-biais
Comme en krigeage simple, la prévision doit être non biaisée, il faut donc que :
X X
E[Ẑ(s0 ) − Z(s0 )] = E[a + λi Z(si ) − Z(s0 )] = a + µ( λi − 1) = 0
i∈V (s0 ) i∈V (s0 )
P
De l’étape précédente, la contrainte i∈V (s0 ) λi = 1 doit être respectée. Ainsi, a doit
être fixé à zéro. La forme de la prévision se simplifie donc à :
X X
Ẑ(s0 ) = λi Z(si ) avec λi = 1
i∈V (s0 ) i∈V (s0 )
4. Contrainte d’optimalité
V ar[Ẑ(s0 ) − Z(s0 )]
= E[(
"µẐ(s0 ) − Z(s0 ))2 ]
X ¶2 #
=E λi Z(si ) − Z(s0 )
"µi∈V (s0 ) ¶2 #
X X
=E λi µ + λi δ(si ) − µ − δ(s0 )
"µi∈V (s0 ) i∈V (s0 )
¶2 #
X
=E λi δ(si ) − δ(s0 )
"µi∈V (s0 ) ¶2 #
X X
=E λi δ(si ) − 2δ(s0 ) λi δ(si ) + δ(s0 )2
" i∈V (s0 ) i∈V (s0 ) #
X X X
=E λi λj δ(si )δ(sj ) − 2δ(s0 ) λi δ(si ) + δ(s0 )2
"i∈V (s0 ) j∈V (s0 ) i∈V (s0 )
X X X
=E λi λj δ(si )δ(sj ) − λi δ(si )2 +
i∈V (s0 ) j∈V (s0 ) i∈V (s0 )
| {z }
1 #
X X
λi δ(si )2 − 2δ(s0 ) λi δ(si ) + δ(s0 )2
i∈V (s0 ) i∈V (s0 )
| {z }
2
Donc
Il faut minimiser cette expression sous la contrainte λt 1n0 = 1. Cette étape sera
effectuée à l’aide d’un Lagrangien noté `. La fonction à minimiser est : f (λ, `) =
−λt Γλ + 2λt γ 0 + 2`(λt 1n0 − 1). Le vecteur de ses dérivées partielles par rapport aux
λi est :
∂
f (λ, `) = −2Γλ + 2γ 0 + 2`1n0
∂λ
∂
La fonction f (λ, `) admet un point critique lorsque ∂λ f (λ, `) = 0, c’est-à-dire au point
−1
λ̂ = Γ (γ 0 + `1n0 ).
Le Lagrangien ` est estimé en utilisant la contrainte que la somme des poids vaille
1, donc que 1tn0 λ̂ = 1tn0 Γ−1 (γ 0 + `1n0 ) = 1. On obtient :
ˆ 1 − 1tn0 Γ−1 γ 0
`=
1tn0 Γ−1 1n0
Chapitre 4. Théorie du krigeage 36
1 − 1tn0 Γ−1 γ 0
λ̂ = Γ−1 (γ 0 + 1n0 )
1tn0 Γ−1 1n0
Matheron (1969, p.35) a prouvé qu’il s’agit bien de l’unique minimum de f (λ, `).
1 − 1tn0 Γ−1 γ 0
Ẑ(s0 ) = (γ 0 + t −1 1n0 )t Γ−1 Z
1n0 Γ 1n0
(1 − 1tn0 Γ−1 γ 0 )2
σ 2 (s0 ) = V ar[Ẑ(s0 ) − Z(s0 )] = γ t0 Γ−1 γ 0 − .
1tn0 Γ−1 1n0
L’hypothèse de stationnarité sur laquelle repose les deux types de krigeage présentés
précédemment peut souvent être mise en doute. En particulier, il semble souvent er-
roné de postuler que l’espérance de la fonction aléatoire étudiée reste constante ou
quasi-constante sur le champ D. En krigeage avec un modèle de tendance, aussi appelé
krigeage en présence d’une dérive, l’espérance est une fonction des coordonnées spa-
tiales ou de variables régionalisées auxiliaires connues exhaustivement. La majorité des
auteurs en géostatistique parlent alors de krigeage universel (Matheron, 1969 ; Cressie,
1993, p.151) et de krigeage avec dérive externe (Goovaerts, 1997, p.194 ; Wackernagel,
2003, p.283) respectivement.
avec fj (s) fonctions de la position s = (x, y), βj paramètres inconnus et δ(·) fonc-
tion aléatoire stationnaire intrinsèque d’espérance nulle et de structure de dépendance
connue. Les fj (s) sont déterminés par l’utilisateur. Les graphiques des z(si ) en fonction
des coordonnées xi et yi sont utilisés pour guider le choix de ces fonctions. Souvent, une
tendance linéaire ou quadratique est choisie. Par exemple, dans le cas d’une tendance
linéaire, les fj (s) sont : f0 (s) = 1, f1 (s) = x et f2 (s) = y.
Chapitre 4. Théorie du krigeage 37
Krigeage avec dérive externe : En krigeage avec dérive externe, le modèle s’écrit :
p
X
Z(s) = fj (w)βj + δ(s), s∈D
j=0
où w = (w1 (s), ..., wq (s)) est le vecteur des valeurs prises par les q variables régionalisées
auxiliaires au point s. Si une tendance linéaire est choisie, les fj (w) seront f0 (w) = 1,
f1 (w) = w1 , ... , fq (w) = wq .
La résolution des équations du krigeage est la même pour ces deux types de krigeage.
Dans les paragraphes qui suivent, nous utilisons la notation du krigeage universel, mais
la démarche est exactement la même en krigeage avec dérive externe en changeant les
fj (s) pour des fj (w). En fait, le modèle de tendance peut aussi être composé de fonctions
des coordonnées spatiales et de variables régionalisées auxiliaires simultanément.
Sous forme matricielle, le modèle utilisé pour prévoir Z(s0 ) s’énonce comme suit :
à ! (
x 0 β E[δ ∗ ] = 0
Z∗ = + δ∗ avec (4.4)
Xβ Γ, γ 0 connus
où Z ∗ = (Z(s0 ), Z), δ ∗ = (δ(s0 ), δ), β = (β0 , β1 , ..., βp ), x0 = (f0 (s0 ), f1 (s0 ), ..., fp (s0 ))
et X est une matrice n0 × (p + 1) dont l’élément (i, j) est fj (si ). Suivons la démarche
des contraintes afin de prévoir Z(s0 ) par krigeage universel.
1. Contrainte de linéarité
Encore une fois, Ẑ(s0 ) doit être une combinaison linéaire des Z(si ). Il s’écrit donc
sous la forme : X
Ẑ(s0 ) = a + λi Z(si ) = a + λt Z
i∈V (s0 )
2. Contrainte d’autorisation
3. Contrainte de non-biais
Afin que cette espérance vaille zéro pour tout βj , j = 0, ..., p, il faut que a = 0 et
P t t
i∈V (s0 ) λi fj (si ) − fj (s0 ) = 0 pour j = 0, ..., p. Ces contraintes s’écrivent λ X = x0
sous forme matricielle.
P
Ainsi, sans oublier la contrainte d’autorisation i∈V (s0 ) λi = 1, il y a au total p +
2 contraintes sur les poids λ en krigeage universel. Toutefois, afin de simplifier les
P
calculs, on suppose que f0 (·) = 1. Ainsi, i∈V (s0 ) λi fj (si ) = fj (s0 ) pour j = 0 revient
P
à i∈V (s0 ) λi = 1. Ce postulat est usuel en krigeage (Cressie, 1993, p.152 ; Goovaerts,
1997, p.140). Il permet d’éliminer une contrainte. L’estimateur devient donc :
X X
Ẑ(s0 ) = λi Z(si ) avec λi fj (si ) = fj (s0 ) pour j = 0, ..., p
i∈V (s0 ) i∈V (s0 )
Chapitre 4. Théorie du krigeage 39
4. Contrainte d’optimalité
Ensuite, ` est remplacé par ˆ` dans λ̂. L’unique point minimum de la fonction f (λ, `)
(Matheron, 1969, p.35) est :
Notons que le krigeage ordinaire peut être vu comme un cas particulier du krigeage
universel avec p = 0 et f0 (·)=1. Ainsi, la prévision en krigeage ordinaire, ainsi que sa
variance, peuvent être retrouvées en remplaçant X par 1n0 et x0 par 1 dans les deux
expressions précédentes.
Dans un cas stationnaire d’ordre deux, une prévision par krigeage avec modèle de
tendance prend la forme :
La même prévision peut être obtenue en effectuant un krigeage simple avec espérance
nulle (m = 0) sur les résidus d’une régression linéaire qui tient compte de la dépendance
spatiale des erreurs (Hengl et al., 2003). Les paramètres du modèle sont donc estimés
par une méthode des moindres carrés généralisés ou par une méthode du maximum de
vraisemblance avec une matrice de variances-covariances non diagonale. Dans ce cas,
une prévision de Z(s0 ) est formée en additionnant la prévision de la tendance générale
par régression à la prévision par krigeage simple des erreurs e = Z − X β̂ gls :
Cependant, certains utilisent cette approche en effectuant une régression qui ne tient pas
compte de la dépendance spatiale des données, par exemple une régression des moindres
carrés ordinaires (β̂ ols = (X t X)−1 X t Z). Dans ce cas, les prévisions ne sont pas tout
à fait les mêmes que celles obtenues par krigeage avec modèle de tendance. Ce type
de krigeage est parfois appelé en anglais « detrended kriging » (Nalder et Wein, 1998;
Phillips et al., 1992). Enfin, une autre approche parfois utilisée consiste à simplement
effectuer un krigeage ordinaire sur des erreurs z(si ) − w(si ). Ce type de krigeage porte
parfois le nom de « krigeage résiduel » (Hessami, 2002, p.29). Ce krigeage est en fait
un krigeage avec dérive externe avec une seule variable régionalisée auxiliaire, donc
q = 1 et w = w(s), en fixant f0 (w) = 1, f1 (w) = w(s) et β1 = 1. Ainsi, β0 est le seul
paramètre inconnu dans la dérive. Le modèle s’écrit donc Z(s) = β0 + w(s) + δ(s) , ou
Z(s) − w(s) = β0 + δ(s) pour s ∈ D. Ce modèle est particulièrement intéressant lorsque
la variable régionalisée auxiliaire mesure le même phénomène naturel que la variable
régionalisée à interpoler, mais avec moins de précision.
Afin de sortir de ce cercle vicieux, une solution consiste à d’abord obtenir un esti-
mateur des moindres carrés ordinaires de β, estimer ensuite un semi-variogramme sur
les résidus, puis calculer un estimateur des moindres carrés généralisés de β, et ainsi de
suite (Cressie, 1993, p.166). Au fil des itérations, les résidus des moindres carrés ordi-
naires s’approchent de plus en plus des résidus des moindres carrés généralisés. Hengl
et al. (2003) affirment qu’en pratique une seule itération suffit pour obtenir des résultats
de krigeage satisfaisants.
Cependant, Cressie (1993, p.167) soulève que même en appliquant cette méthode
itérative, l’estimation du semi-variogramme reste biaisée (Matheron, 1969, p.37, 1970,
p.155). Les conséquences de ce biais sont cependant difficiles à évaluer. Ainsi, le krigeage
avec un modèle de tendance reste une méthode qui peut s’avérer bonne, mais la prudence
est de mise lors de son utilisation. Le biais dans l’estimation du semi-variogramme
en krigeage avec modèle de tendance est une des motivations au développement de
méthodes bayésiennes.
4.5 Discussions
Les trois sections précédentes ont présenté l’essentiel de la théorie du krigeage for-
malisée par Matheron (1962, 1963a,b, 1965, 1969, 1970). Ces sections devraient donc
suffirent à une compréhension de la méthode dans le but d’une utilisation de base. Ce-
pendant, avant de passer à la mise en oeuvre du krigeage, il est intéressant de s’attarder
sur certains aspects théoriques plus poussés ainsi que sur d’autres types de krigeage
développés plus récemment. Les remarques théoriques traitées ici touchent la norma-
lité des données, le travail sur des données transformées, notamment par la fonction
logarithmique, et le cokrigeage, une extension du krigeage dans le cas multivariable.
sur les moments de la distribution de Z(·) par l’intermédiaire de δ(·) et non sur la loi
de la distribution. L’utilisateur du krigeage peut donc croire que la performance de la
méthode est indépendante de la loi de distribution des données. Malheureusement, ce
n’est pas le cas. Le krigeage a tendance à fournir de meilleures prévisions lorsque les
données suivent une loi normale.
Cette remarque peut être vérifiée en pratique, mais elle s’explique aussi théori-
quement. Elle est attribuable à la contrainte en apparence anodine de linéarité de la
prévision. Cette hypothèse est typiquement inadéquate dans un contexte non normal
(Cressie, 1993, p.110). Cependant, si Z(·) est une fonction aléatoire gaussienne, cette
linéarité est appropriée.
Pour justifier cette affirmation, rappelons d’abord que la prévision par krigeage est
linéaire, sans biais et qu’elle minimise la variance de l’erreur de prévision. Elle minimise
donc l’erreur quadratique moyenne :
Une prévision non contrainte à être linéaire et minimisant l’erreur quadratique moyenne
conditionnelle E[(Ẑ(s0 ) − Z(s0 ))2 |Z] serait E[Z(s0 )|Z] (Cressie, 1993, p.108). Thé-
oriquement, cette prévision spatiale est équivalente ou meilleure à celle obtenue par
krigeage. Cependant, elle dépend de la loi de distribution de Z(·). Supposons donc que
Z(·) est une fonction aléatoire gaussienne telle que :
à ! Ãà ! à !!
Z(s0 ) x0 t β σ 2 ct0
∼ N(n+1) ,
Z Xβ c0 Σ
Cette expression est bien une combinaison linéaire des Z(si ). L’hypothèse de linéarité
est par conséquent appropriée dans le cas normal. De plus, cette prévision est en
fait la même que celle du krigeage universel dans un cadre stationnaire de second
ordre si β est remplacée par son estimation du maximum de vraisemblance β̂ MV =
(X t Σ−1 X)−1 X t Σ−1 Z. En effet, nous obtenons alors :
Il faut donc garder en tête que la prévision spatiale par krigeage peut mal performer
lorsque Z(·) est loin d’être gaussienne. Dans un tel cas, une transformation des données
ou encore l’utilisation de méthodes de la géostatistique non-linéaire peuvent être envi-
sagées. Ces deux solutions sont détaillées dans les prochains points de discussion.
Il peut parfois être avantageux de travailler sur une transformation des données
plutôt que sur les données brutes. Par transformation de données, nous entendons ici
une fonction mathématique (par exemple le logarithme, la racine, le carré, etc.) ap-
pliquée aux données. Une transformation peut aider à rendre la distribution des données
plus près d’une distribution gaussienne. Elle peut aussi être employée afin de s’assurer
que les prévisions par krigeage prennent des valeurs à l’intérieur d’un certain intervalle.
Par exemple, si une variable régionalisée telle la concentration d’un polluant dans l’air
est étudiée, ses valeurs régionalisées seront obligatoirement comprises entre 0 et 100.
Cependant, rien n’empêche une prévision par krigeage de cette variable de prendre une
valeur à l’extérieur de cet intervalle. Le même problème se présente pour des variables
régionalisées positives telle une quantité de précipitations tombées. Afin de régler ce
problème, une solution possible est d’ajouter des contraintes aux poids de krigeage.
Par exemple, en forçant les poids du krigeage à être positifs, les prévisions seront as-
surément positives si les valeurs régionalisées sont toutes positives. Une autre solution
est d’appliquer une transformation aux données.
Krigeage lognormal
où Ŷ (s0 ) est la prévision de Y (s0 ) par krigeage ordinaire, σY2 (s0 ) est la variance de
cette prévision et `Y est le lagrangien intervenant dans la résolution des équations du
krigeage ordinaire sur les Y (si ), i ∈ V (s0 ). Le facteur exp (σY2 (s0 )/2 − `Y ) sert donc
à corriger pour le biais de la prévision. La variance de cet estimateur est donnée par
Cressie (1993, p.136).
Il faut cependant mentionner que Roth (1998) s’interroge sur la pertinence du kri-
geage lognormal pour la prévision spatiale. Il soutient que la propriété d’absence de biais
de sa prévision provient d’une surestimation de la majorité des valeurs régionalisées
compensée par une sous-estimation de quelques valeurs plus grandes. Il conseille donc
d’utiliser avec précaution cette méthode.
Chapitre 4. Théorie du krigeage 45
Krigeage trans-gaussien
Plus généralement, Cressie (1993, p.137) propose une formule de prévision qui sup-
pose que la fonction aléatoire transformée est Y (s) = φ−1 (Z(s)) pour tout s ∈ D
où φ est une fonction deux fois différentiable choisie telle que Y (·) soit gausienne. La
prévision s’écrit :
où Ŷ (s0 ) est la prévision de Y (s0 ) par krigeage ordinaire, σY2 (s0 ) est la variance de cette
prévision, `Y est le lagrangien intervenant dans la résolution des équations du krigeage
ordinaire sur les Y (si ) = φ−1 (Z(si )), i ∈ V (s0 ) et µ̂Y = 1tn0 Σ−1 t −1
Y Y /(1n0 ΣY 1n0 ). Il est
bien mis en évidence par cette formule que l’analyse variographique doit être effectuée
sur les données transformées. Cette prévision se base sur un développement en série
Taylor de φ(Y (s)) autour de µY , l’espérance de Y (s). Elle approximativement sans
biais pour σY2 (s0 ) petit. Les détails de la démarche menant à cette prévision peuvent
être trouvés dans Kozintseva (1999, p.17). Notons que pour φ(·) = exp (·), les prévisions
des krigeages trans-gaussien et lognormal sont approximativement les mêmes.
Dans la même ligne d’idées, De Olivieira et al. (1997) ont proposé un modèle bayésien
trans-gaussien. Ce modèle ne requiert pas la détermination préalable de la transforma-
tion, mais plutôt d’une famille paramétrique de transformations.
Jusqu’à maintenant, dans ce chapitre, l’accent a été mis sur l’aspect stationnaire ou
non des différents types de krigeage. Deux branches de la géostatistique se sont des-
sinées : la géostatistique stationnaire et la géostatistique non-stationnaire. Le krigeage
simple et le krigeage ordinaire (sections 4.2 et 4.3) se classent dans la division station-
naire et le krigeage universel et le krigeage avec dérive externe (section 4.4) dans la
division non-stationnaire. Un autre classement possible des techniques géostatistiques
concerne le nombre de variables régionalisées incluses dans l’analyse. Cet aspect a
d’ailleurs été abordé à la section 2.2.6. En géostatistique univariable, une seule va-
riable régionalisée intervient dans les calculs tandis qu’en géostatistique multivariable,
plus d’une variable régionalisées est exploitée. Ainsi, le krigeage simple, le krigeage or-
dinaire et le krigeage universel se classent dans la catégorie univariable et le krigeage
avec dérive externe dans la branche multivariable.
cokrigeage, brièvement présenté à la section 2.2.6. Dans cette section, nous revenons sur
cette méthode pour l’approfondir davantage. Mais avant, mentionnons qu’une interpo-
lation multivariable peut parfois être effectuée en utilisant une technique univariable.
Par exemple, des variables auxiliaires catégoriques, telles qu’un type de roche, peuvent
permettre de partitionner le champ D en strates distinctes. Dans ce cas, s’il y a suf-
fisamment de points d’observation, le krigeage peut s’effectuer à l’intérieur des strates
(Goovaerts, 1997, p.187) plutôt que sur tout le champ. Ainsi, un semi-variogramme dis-
tinct est calculé par strate et une prévision en un point est effectuée uniquement à partir
des points d’observation appartenant à la même strate. En ce qui concerne les variables
auxiliaires continues, une façon simple de les exploiter est d’effectuer un krigeage simple
avec une espérance variant localement (Goovaerts, 1997, p.190). Cette espérance, sup-
posée connue, est alors estimée en chaque point de prévision par un modèle de régression
incluant les variables régionalisées auxiliaires.
Cokrigeage
Les poids de cette combinaison linéaire sont choisis de façon à minimiser la variance de
l’erreur de prévision sous une contrainte de non-biais, comme en krigeage. Pour ce faire,
toutes les variables régionalisées sont modélisées par des fonctions aléatoires, même les
variables régionalisées auxiliaires. C’est d’ailleurs ce qui distingue le cokrigeage des
autres techniques de la géostatistique multivariable. Cette modélisation stochastique
permet d’exploiter, en plus des relations entre les variables régionalisées, la comporte-
ment spatial de chacune de ces variables.
Chapitre 4. Théorie du krigeage 47
Z(s) = µ0 (s) + δ0 (s) et Wj (s) = µj (s) + δj (s) pour tout j ∈ {1, ..., q}, s∈D
où les µj (·) sont des structures déterministes pour les espérances et les δj (·) sont
des fonctions aléatoires, chacune d’espérance nulle et conjointement stationnaires (de
deuxième ordre ou intrinsèquement) de structure de dépendance connue. La stationna-
rité conjointe (Arnaud et Emery, 2000, p.157) implique que :
Dans le cadre stationnaire d’ordre deux : Cov(δj (s), δk (s + h)) soit uniquement
fonction de h, le vecteur de translation entre les points s et s+h, pour tout couple
j, k ∈ {0, ..., q} et toute paire de points s, s + h ∈ D. Ces fonctions de covariance,
notées Cjk (h), sont nommées covariogrammes simples (j = k) et croisés (j 6= k) ;
Dans le cadre stationnaire intrinsèque : Cov[δj (s+h)−δj (s), δk (s+h)−δk (s)] soit
uniquement fonction de h, le vecteur de translation entre les points s et s+h, pour
tout couple j, k ∈ {0, ..., q} et toute paire de points s, s + h ∈ D. Ces fonctions de
covariance, notées 2γjk (h), sont nommés variogrammes simples (j = k) et croisés
(j 6= k).
La structure de dépendance spatiale de l’ensemble des fonctions aléatoires résiduelles est
donc représentée par les Cjk (h) ou par les γjk (h). Cressie (1993, p.66) note qu’il existe en
fait deux généralisations du semi-variogramme au cas multivariable. Au lieu de la fonc-
tion γjk (·) introduite précédemment, certains utilisent la fonction πjk (h) = 21 V ar(δj (s +
h) − δk (s)), que Wackernagel (2003, p.149) nomme pseudo semi-variogramme croisé.
L’analyse variographique précédant un cokrigeage doit inclure l’estimation et la mo-
délisation des covariogrammes, ou des (pseudo) semi-variogrammes, pour tout j, k ∈
{0, ..., q}. Cette tâche s’avère souvent complexe. Comme dans le cas univariable, seuls les
modèles assurant la positivité de la variance de toute combinaison linéaire de variables
aléatoires sont admissibles. Les « modèles linéaires de corégionalisation » garantissent
cette positivité. Ils sont donc fréquemment utilisés dans la pratique du cokrigeage. Le
lecteur est référé à Goovaerts (1997, p.109), Wackernagel (2003, p.145) ou Arnaud et
Emery (2000, p.157) pour plus d’informations sur ces modèles ou sur tout autre aspect
de l’analyse variographique menant au cokrigeage.
Comme pour le krigeage, les versions simple, ordinaire et avec modèle de tendance du
cokrigeage ont été développées. Le type de cokrigeage dépend de la forme des fonctions
µj (·) pour j ∈ {0, ..., q} de la même façon qu’en krigeage. La résolution des équations
de cokrigeage ne sera pas présentée dans ce mémoire. Hoef et Cressie (1993) ainsi que
Matheron (1970, p.187) le font pour le cokrigeage avec modèle de tendance. Goovaerts
(1997, p.204) soulève cependant que ce type de cokrigeage est très difficile à implan-
ter, particulièrement en raison du biais dans l’estimation des covariogrammes ou des
semi-variogrammes sur les résidus. Ce problème, traité dans la section 4.4 pour le kri-
geage, devient encore plus sévère dans le cas multivariable. L’utilisation du cokrigeage
Chapitre 4. Théorie du krigeage 48
simple ou ordinaire est donc plus recommandé. Goovaerts (1997, p.203), Wackernagel
(2003, p.159) ainsi que Arnaud et Emery (2000, p.182) présentent ces types de cokri-
geage. Notons que plusieurs variantes de formules de prévision peuvent être obtenues
en cokrigeage ordinaire dépendamment des contraintes de non-biais choisies (Isaaks et
Srivastava, 1989, p.403).
Le krigeage classique étant assez sensible aux données extrêmes, des versions ro-
bustes du krigeage ordinaire et du krigeage universel ont été développées.
Krigeage robuste : Proposé par Hawkins et Cressie (1984), le krigeage robuste repré-
sente une solution de rechange au krigeage ordinaire. L’idée derrière cette méthode
Chapitre 4. Théorie du krigeage 49
est de diminuer le poids des données extrêmes attribuables à une distribution non
gaussienne des données. Cressie (1993, p.144) présente cette technique.
Krigeage « median polish » : Le krigeage universel possède aussi sa version ro-
buste : le krigeage « median polish » Cressie (1984, 1986). Cette méthode est
spécialement adaptée aux données se présentant sur une grille dans le plan, qui
sont alors vues comme un tableau de fréquences.
Le krigeage avec modèle de tendance n’est pas le seul à entrer dans la catégorie non
stationnaire. Des généralisations de ce krigeage ont fait place à d’autres techniques.
Krigeage bayésien : Une modélisation autre que celle du krigeage avec modèle de
tendance est de considérer la fonction de variation à grande échelle µ(·) aléatoire
plutôt que déterministe. Les paramètres de cette fonction deviennent alors des
variables aléatoires qui sont associées à des loi à priori. C’est l’idée de base du
krigeage bayésien, proposé par Omre (1987) afin de permettre l’introduction de
connaissances à priori dans la prévision par krigeage.
Krigeage intrinsèque généralisé : Le biais dans l’estimation du semi-variogramme
en krigeage avec modèle de tendance a mené Matheron (1973) à proposer le kri-
geage intrinsèque généralisé. Chauvet (1999, p.177) présente cette méthode, qui
fait intervenir des fonctions aléatoires plus complexes que celles stationnaires de
second ordre ou intrinsèque ainsi qu’un outil structural autre que le semi-vario-
gramme appelé covariance généralisée.
Tel qu’expliqué à la section 4.5.1, la linéarité des prévisions par krigeage implique
que celui-ci accomplisse parfois de mauvaises performances lorsque les données sont
loin d’être normales. La géostatistique non linéaire (Cressie, 1993, p.278) propose des
solutions à ce problème.
Ainsi, le krigeage consiste simplement en une prévision linéaire sans biais à variance
minimale. La résolution d’un système d’équations conforme aux hypothèses du modèle
permet d’obtenir une prévision par krigeage. Cependant, cette prévision a tendance a
moins bien performer dans certaines circonstances telles la non normalité des données
ou lorsque les données sont peu nombreuses.
Seule l’interpolation spatiale ponctuelle a été traitée dans ce chapitre. Notons que
la théorie pour effectuer par krigeage une prévision sur une sous-région du champ D est
développée dans la littérature. Il s’agit du krigeage par bloc (Goovaerts, 1997, p.152).
Tel que l’indique Cressie (1993, p.289), la mise en oeuvre du krigeage s’effectue
en suivant certaines étapes. La première est l’analyse exploratoire qui permet de se
familiariser avec les données et d’éclairer le choix du modèle. Ensuite, le modèle peut
être énoncé. Cette deuxième étape inclue le choix de la forme de la tendance déterministe
pour l’espérance de Z(·), l’analyse variographique, puis une validation croisée si désirée.
Cette dernière sous-étape permet de comparer la performance de différents modèles
afin de sélectionner celui susceptible de mener aux meilleures prévisions. Finalement,
l’interpolation est effectuée par krigeage. La figure 5.1 schématise cette démarche. La
flèche courbe reliant la validation croisée à la détermination de µ(·) illustre le fait que
les étapes de la formulation du modèle peuvent être reprises afin de comparer plusieurs
modèles. Chacune des étapes sont vues plus en détails dans les paragraphes suivants.
Chapitre 5. Mise en oeuvre du krigeage 52
1- ANALYSE EXPLORATOIRE
2- FORMULATION DU MODÈLE
a) Détermination de µ(.)
b) Analyse variographique
c) Validation croisée
3- KRIGEAGE
Comme dans toute analyse statistique, il est bon de débuter nos travaux par une
analyse exploratoire des données. Isaaks et Srivastava soulignent l’importance de cette
étape en lui consacrant pratiquement le tiers de leur livre Applied Geostatistics (1989,
p.10 à p.183). Tel que le font Isaaks et Srivastava, cette étape préliminaire sera présentée
ici en la divisant en trois volets : un volet univarié, un volet bivarié et un volet spatial.
Description univariée :
L’analyse exploratoire vise à donner une idée de la distribution des données. Rap-
pelons qu’en géostatistique, les données sont considérées comme une réalisation par-
tielle d’une fonction aléatoire. Ainsi, chaque donnée serait une observation d’une va-
riable aléatoire différente et chaque variable aléatoire aurait sa propre distribution.
L’échantillon de données serait donc multivarié. Toutefois, l’analyse exploratoire est ha-
bituellement débutée en oubliant cette hypothèse de base et en considérant les données
comme plusieurs réalisations d’une même variable aléatoire. Des statistiques descriptives
univariées telles la moyenne et la variance expérimentale ainsi que des outils graphiques
Chapitre 5. Mise en oeuvre du krigeage 53
tels l’histogramme et le diagramme en boı̂te peuvent alors être utilisés afin de se fa-
miliariser avec la distribution de la variable aléatoire. Si cette distribution ne semble
vraiment pas être gausienne, une transformation de variable peut être envisagée. La
transformation la plus adéquate sera donc recherchée. Cette tâche peut être accomplie
par une méthode d’essais/erreurs ou encore en utilisant des techniques proposées par
Box et Cox (1964). Si une transformation est choisie, le reste de l’analyse exploratoire
devrait être effectuée sur les données transformées.
Description bivariée :
Toutes ces techniques requièrent des couples de données aux points d’observation
s1 à sn . Si certaines variables régionalisées auxiliaires n’ont pas été observées en ces
points, une interpolation spatiale est préalablement effectuée afin de prévoir ces valeurs
régionalisées. Une simple technique déterministe, telle une méthode barycentrique, est
habituellement employée pour effectuer cette tâche. Les données interpolées obtenues
seront ensuite nécessaires à la résolution des équations du krigeage si le krigeage avec
dérive externe est choisi.
Description spatiale :
En se rappelant maintenant que les observations sont en fait multivariées, les va-
riables régionalisées peuvent être explorées spatialement. Cette analyse descriptive est
particulièrement importante pour la variable régionalisée à interpoler. Les données
peuvent d’abord être visualisées à l’aide de graphiques 3D ou de cartes de courbes
de niveaux. Ces outils graphiques aident notamment à juger si l’espérance de la fonc-
tion aléatoire modélisant la variable régionalisée peut être vue comme une fonction
Chapitre 5. Mise en oeuvre du krigeage 54
constante. La présence de valeurs fortement différentes les unes des autres indiquerait
que ce n’est pas le cas. Ainsi, la stationnarité du premier moment est étudiée par ces
graphiques.
Cette stationnarité peut aussi être examinée à l’aide d’une fenêtre mobile s’il y a
assez de données (Isaaks et Srivastava, 1989, p.46). Cette fenêtre mobile est simplement
un outil qui permet de délimiter des sous-régions du champ qui peuvent se chevaucher.
Pour chacune de ces sous-régions, des statistiques descriptives de base sont calculées,
telles la moyenne et la variance expérimentale. Ainsi, la stationnarité du deuxième
moment peut aussi être étudiée avec une fenêtre mobile. Par exemple, une variance
expérimentale qui augmenterait toujours en agrandissant la fenêtre suggérerait une
variance infinie, donc une hypothèse de stationnarité intrinsèque serait plus appropriée
qu’une hypothèse de stationnarité de deuxième d’ordre.
La structure de dépendance spatiale que présente les données sera étudiée en profon-
deur à l’étape de l’analyse variographique. Elle peut aussi être explorée préliminairement
à l’aide de h-scatterplots (Isaaks et Srivastava, 1989, p.52 ; Arnaud et Emery, 2000,
p.117), nommés nuages de corrélation différée en français. Il s’agit de nuages de points
formés par les couples (z(si ), z(si + h)). Ils indiquent, pour un vecteur de transla-
tion donné h, si les valeurs prises aux sites séparés par h sont semblables. Lorsque les
sites d’observation ne sont pas répartis de façon régulière sur le territoire et que la
dépendance spatiale est supposée isotrope, ces graphiques de dispersion peuvent être
faits pour les couples (z(si ), z(sj )) avec |si − sj | approximativement égale à certaines
distances prédéterminées.
Une bonne exploration des données est essentielle au choix du modèle. En se ba-
sant sur les résultats de cette analyse descriptive, le type de krigeage à employer peut
être choisi judicieusement et l’analyse variographique menée de façon plus éclairée.
Si plusieurs modèles sont envisagés, la validation croisée peut être utilisée afin d’en
sélectionner un. Mais avant tout, il faut déterminer si l’interpolation sera faite directe-
Chapitre 5. Mise en oeuvre du krigeage 55
ment sur la fonction aléatoire d’intérêt ou sur une transformation de celle-ci. Le plus
simple est bien sûr de travailler sur les données non transformées. Cette option est
toujours privilégiée. Cependant, un utilisateur peut choisir de travailler sur une trans-
formation des données pour améliorer la qualité de l’interpolation dans le cas d’un
important écart à la normalité. Dans ce cas, il doit toutefois s’assurer d’être en mesure
d’effectuer la transformation inverse après le krigeage.
a) Détermination de µ(·) :
Dans le cas où aucune variable régionalisée auxiliaire n’est disponible, ou que celles-
ci ne sont pas assez reliées à la variable régionalisée d’intérêt, un krigeage univariable est
sélectionné. Si l’analyse descriptive a clairement mis en évidence une non stationnarité
de l’espérance de la fonction aléatoire à interpoler, le krigeage universel est choisi. Les
fonctions fi (·) composant le modèle de tendance sont alors déterminées à l’aide des gra-
phiques de comportement directionnel. Par contre, lorsque la stationnarité du premier
moment peut être acceptée, un krigeage stationnaire est sélectionné. Si l’espérance de
la fonction aléatoire à interpoler est connue, le krigeage simple est employé, sinon le
krigeage ordinaire est utilisé.
b) Analyse variographique :
c) Validation croisée :
modèles (Isaaks et Srivastava, 1989, p.351 ; Arnaud et Emery, 2000, p.152 ; Cressie, 1993,
p.101 ; Wackernagel, 2003, p.87). Cette technique peut être utilisée non seulement pour
comparer les modèles, mais pour en choisir un, c’est-à-dire choisir le type de krigeage
et le modèle variographique employé. La validation croisée de type « leave-one-out »
consiste à retirer une à une les observations pour ensuite les prévoir par krigeage à
partir des autres données. La prévision par validation croisée d’une observation Z(si )
sera notée Ẑ−i (si ). Des erreurs de validation croisée sont calculées en soustrayant la
valeur observée à la la valeur estimée : e−i (si ) = Ẑ−i (si ) − Z(si ). À partir de ces
erreurs, un ou des indices sont calculés (voir Myers, 1993, pour une liste d’indices
P
possibles). Souvent, l’indice de la moyenne des erreurs au carré n1 ni=1 e−i (si )2 , nommé
EQM (Erreur Quadratique Moyenne), est employé. Les avantages de cet indice sont
qu’il incorpore le biais et la variance de la distribution des erreurs et que l’absence de
division par l’écart type de krigeage évite la dépendance des résultats par rapport au
type de krigeage effectué. Un indice équivalent mais plus robuste est la moyenne ou la
médiane des erreurs en valeur absolue. Cet indice peut donc être préféré à l’EQM lors
de la présence de valeurs extrêmes. Un indice tel que l’EQM est vu comme une erreur
de prévision. On cherche donc à le minimiser. Ainsi, les étapes a) à c) de la formulation
du modèle sont effectuées pour plusieurs modèles et ceux obtenant les plus petits EQM
sont jugés meilleurs. Une façon d’automatiser le choix du modèle est donc de retenir
celui qui minimise l’EQM.
d’estimer les paramètres du semi-variogramme par une méthode des moindres carrés
(voir section 3.5).
5.1.3 Krigeage
Après avoir sélectionné un modèle, l’interpolation peut être effectuée. Cette étape
en est une de calculs seulement. Les prévisions et les variances de krigeage sont évaluées
pour tous les points désirés. Pour ce faire, il suffit de d’abord calculer Γ et γ 0 , ou Σ
et c0 , à partir du modèle variographique, puis de remplacer ces valeurs ainsi que les
observations de la variable régionalisée d’intérêt et des variables régionalisées auxiliaires,
s’il y a lieu, dans les expressions trouvées au chapitre 4 pour le type de krigeage choisi.
Un des logiciels les plus connus en géostatistique est GSLIB de Deutsch et Jour-
nel (1998). Le code source de ce logiciel peut être téléchargé gratuitement sur le site
web http://www.gslib.com/. Cependant, le livre GSLIB : Geostatistical Software Li-
Chapitre 5. Mise en oeuvre du krigeage 59
brary and User’s Guide est indispensable pour son utilisation. Ce logiciel permet d’ef-
fectuer plusieurs types de krigeage. Il est fréquemment utilisé par les chercheurs en
géostatistique. Par contre, dans ce mémoire, les calculs géostatistiques ont plutôt été
effectués avec le logiciel S-Plus (Insightful, 2004) bien connu des statisticiens. Ce logiciel
propose un module géostatistique, S+SpatialStats, qui se base sur le livre de Cressie
(1993). Ce module s’utilise facilement, mais il est assez limité, ajustant uniquement
des modèles variographiques bornés et effectuant seulement des prévisions par krigeage
ordinaire ou universel. Le logiciel Gstat (Pebesma et Wesseling, 1998), offert en librairie
S-Plus (Pebesma, 2004a), est beaucoup plus complet. Il propose plusieurs modèles va-
riographiques bornés ou non et effectue du krigeage ou du cokrigeage simple, ordinaire
ou avec modèle de tendance. Il permet aussi de faire des simulations géostatistiques et
de la validation croisée. La librairie Gstat se télécharge gratuitement à partir du site
web http://www.gstat.org. Elle est aussi offerte en extension pour le logiciel R (R
project, 2004), qui peut être vu comme une version gratuite de S-Plus. De plus, la do-
cumentation très complète de Gstat (Pebesma, 2004b, 1999) rend facile son utilisation
même si celle-ci s’effectue par la saisie de commandes. Gstat est donc un très bon logi-
ciel géostatistique accessible à tous. C’est ce logiciel, sous sa forme de librairie S-Plus,
qui a été utilisé pour effectuer les analyses de la section suivante et du chapitre 6. Ce
mémoire contient d’ailleurs à l’annexe B le programme S-Plus qui a permis d’effectuer
les interpolations au chapitre 6.
À ce point-ci du mémoire, assez d’informations ont été fournies pour savoir comment
effectuer un krigeage. Le chapitre suivant traitera de l’utilisation du krigeage afin de
résoudre une problématique d’interpolation spatiale. Celle-ci concerne la préparation
d’intrants pour des modèles de simulation hydrologique. Ces modèles, qui produisent
des prévisions de débit des cours d’eau, basent leurs simulations sur les quantités de
pluie tombées. Ces données observées doivent être fournies aux modèles sur une cer-
taine grille. Par contre, la principale source de données de précipitations, soit les sta-
tions météorologiques, fournie des données ponctuelles irrégulièrement réparties sur
le territoire. Une interpolation spatiale doit donc être effectuée afin de d’estimer les
précipitations en chaque noeud de la grille.
Les données traitées au chapitre 6 ont été fournies par Hydro-Québec et Environne-
ment Canada. Il s’agit de mesures de quantités de pluie tombées provenant de 16 sta-
Chapitre 5. Mise en oeuvre du krigeage 60
Notons que dans ce mémoire, la localisation d’un site est toujours exprimée en
coordonnées UTM (Tranverse de Mercator Universelle) plutôt qu’en latitude-longitude.
Les coordonnées latitude-longitude sont des mesures d’angle se basant sur la forme
sphérique de la Terre, tandis que les coordonnées UTM sont obtenues en projetant
la surface sphérique du globe sur un plan. En conséquence, l’emploi d’une distance
euclidienne est correct pour des coordonnées UTM, mais n’est pas du tout recommandé
pour des coordonnées latitude-longitude.
Les précipitations observées proviennent des 16 stations localisées sur la figure 5.3.
Bien que certaines de ces stations se retrouvent légèrement à l’extérieur du champ
d’étude du bassin de la rivière Gatineau, elles sont incluses dans l’interpolation pour
augmenter le nombre de stations et parce que leur proximité par rapport au champ
d’étude signifie qu’elles permettront très probablement d’améliorer la qualité de l’inter-
polation.
fourni les observations analysées dans ce mémoire ont généralement une précision de
lecture de cumul de précipitations égale à 0.2 mm. Seuls 4 précipitomètre font excep-
tion, ceux des stations Mercier Amont, Barrière Amont, Cabonga Amont et Dozois
Est, avec une précision de lecture de cumul de 3.75 mm. Ces précipitomètres sont donc
beaucoup moins précis que les autres. De plus, Yang et al. (1999) ont mis en évidence
que le vent cause une sous-estimation des précipitations par les précipitomètres. Cette
sous-estimation, qui augmente en fonction de la force des vents, est évidemment plus
importante pour des précipitomètres très exposés aux vents. Ainsi, la qualité des ob-
servations dépend notamment du type de précipitomètre, de son exposition, de son
entretien et des conditions de vent et de température au moment de la mesure. Malheu-
reusement, cette variabilité dans la qualité des observations n’est pas prise en compte
lors des interpolations spatiales effectuées dans ce mémoire. L’emploi de techniques
de krigeage avec erreurs de mesures permettrait par contre de le faire (Cressie, 1993,
p.128).
Notons que pour certaines périodes de 6 heures, au moins une des 16 stations n’a
pas rapporté d’observations pour une raison quelconque. En fait, jusqu’à 4 données sur
Chapitre 5. Mise en oeuvre du krigeage 62
Fig. 5.3 – Localisation des 16 stations météorologiques par rapport aux limites du
bassin versant de la rivière Gatineau
16 sont manquantes pour quelques périodes. Les interpolations sont donc effectuées à
partir de 12 à 16 mesures de la variable régionalisée à l’étude. Webster et Oliver (1993)
affirment que 150 données sont nécessaires à une estimation vraiment fiable du semi-
variogramme et plus encore pour décrire l’anisotropie. Pour les données utilisées ici, il
est donc impensable d’étudier l’anisotropie de la variable régionalisée et les estimations
de semi-variogrammes sont peu précises. En conséquence, il est difficile de modéliser vi-
suellement les semi-variogrammes expérimentaux et l’emploi de la procédure de sélection
automatique du modèle variographique par validation croisée est justifiable.
24 km par 24 km. Pour obtenir les données sur la grille de 10 km par 10 km mentionnée
à la section 5.4, Environnement Canada a interpolé les sorties du modèle.
Le modèle GEM utilise comme système de référence temporel le temps UTC (Temps
Universel Coordonné). Le temps UTC est la base légale de l’heure dans le monde. Il
s’agit en fait de l’heure du méridien d’origine de Greenwich. Ce système de référence
facilite le travail sur de grands territoires traversant plusieurs fuseaux horaire. Dans la
vie de tous les jours au Québec, c’est l’heure normale de l’est (HNE) en hiver ou l’heure
avancée de l’est (HAE) en été qui est utilisé. Pour se ramener à l’heure HNE à partir
de l’heure UTC il suffit de soustraire 5 heures et pour se rammener à l’heure HAE il
faut soustraire 4 heures.
Le modèle GEM fournit des prévisions deux fois par jour : à 00 UTC et à 12 UTC.
Les prévisions du modèle GEM régional concernent les 48 prochaines heures. Ces 48
heures sont divisées en 8 périodes de 6 heures. Les variables météorologiques prévues
par GEM sont nombreuses. Une seule de ces variable est utilisée ici : les précipitations
cumulées sur des périodes de 6 heures. Étant donné que le modèle GEM fournit à toutes
les 12 heures des prévisions pour un intervalle de temps de 48 heures, des prévisions
GEM se chevauchent dans le temps. Pour chaque période de 6 heures, 4 prévisions GEM
sont disponibles. Pour illustrer cet énoncé, considérons la période entre 18 UTC le 4
août 2003 et 00 UTC le 5 août 2003. Cette période est localisée par le rectangle hachuré
sur la figure 5.4. Quatre prévisions GEM concernent cette période : la prévision émise à
00 UTC le 3 août 2003 pour la huitième période de 6 heures plus tard (P8), la prévision
émise à 12 UTC le 3 août 2003 pour la sixième période de 6 heures plus tard (P6), la
prévision émise à 00 UTC le 4 août 2003 pour la quatrième période de 6 heures plus tard
(P4), la prévision émise à 00 UTC le 4 août 2003 pour la deuxième période de 6 heures
plus tard (P2). Ainsi, un choix doit être fait parmi ces prévisions. Des spécialistes en
prévision numérique du temps avec le modèle GEM ont jugé que les prévisions les plus
exactes étaient celles effectuées pour la deuxième et la troisième période de 6 heures sui-
vant le moment de l’émission des prévisions. Les prévisions exploitées dans ce mémoire
sont donc celles associées aux intervalles P2 et P3 de chaque émission de prévisions,
soit les segments bleus de la figure 5.4. Un météorologue dirait donc que les données
auxiliaires employées dans le chapitre 6 sont les prévisions de précipitations du modèle
GEM régional 24 km 6h-18h ramenées à 10 km. Des informations sur le fonctionnement
du modèle atmosphérique GEM peuvent être trouvées sur le site internet du Service
météorologique du Canada (SMC, 2002), une division d’Environnement Canada.
Chapitre 5. Mise en oeuvre du krigeage 64
....
.
Prévisions émises à P1 P2 P3 P4 P5 P6 P7 P8
.
00 UTC le 3 août 2003
Prévisions émises à . P1 P2 P3 P4 P5 P6 P7 P8
.
12 UTC le 3 août 2003
Prévisions émises à
. P1 P2 P3 P4 P5 P6 P7 P8
.
00 UTC le 4 août 2003
Prévisions émises à . P1 P2 P3 P4 P5 P6 P7 P8
.
12 UTC le 4 août 2003
-0.03. De plus, les prévisions GEM sont ici nettement supérieures aux observations des
stations avec des précipitations moyennes de 5.40 mm contrairement à 1.53 pour les
stations. Malgré cette surestimation de GEM par rapport aux stations, les données
des deux sources sont relativement reliées, avec une corrélation de 0.35. Notons cepen-
dant que cette corrélation est plus ou moins précise car elle est calculée à partir de 16
couples de données seulement. Pour obtenir ces couples de données, les valeurs générées
par le modèle GEM ont dû être interpolées aux sites où sont localisées les stations
météorologiques. Cette interpolation a été effectuée par la méthode de l’inverse de la
distance au carré, en considérant pour chaque site d’observation un voisinage d’inter-
polation constitué des quatre points de la grille du modèle GEM formant le carré dans
lequel tombe le site. Les couples de données stations-GEM seront aussi nécessaires pour
déterminer les poids de la prévision par krigeage avec dérive externe.
La figure 5.5 présente les histogrammes des données. Les données des deux sources
présentent des distributions asymétriques avec des fréquences plus élevées pour les va-
leurs faibles. Cette observation est particulièrement critique pour les observations de
stations. La normalité de ces données peut être sérieusement mise en doute. Une trans-
formation logarithmique aiderait certainement à rendre plus symétrique la distribution
des données, mais cette transformation est déconseillée par plusieurs (Roth, 1998). De
plus, la librairie gstat utilisée ici pour mettre en oeuvre l’interpolation par krigeage
ne permet pas d’effectuer la transformation inverse des données après l’interpolation.
Dans l’étude présentée ici, les calculs seront donc effectués sur les données non trans-
formées. Cependant, étant donné que le krigeage est spécialement adéquat pour des
données normales (voir section 4.5.1), il est possible que ses performances laisse ici à
désirer. Notons que la distribution de données de précipitations cumulées sur un petit
intervalle de temps est souvent asymétrique en raison de la masse de données à zéro.
Pour des pas de temps annuel ou mensuel un tel problème se pose rarement, mais
plus le pas de temps est petit, plus ce problème est prononcé. Pour l’interpolation de
données de précipitations à une fine résolution temporelle comme celle de 6 heures, il
serait intéressant d’étudier la modélisation proposée par Velarde et al. (2004) qui est
spécialement adaptée aux distributions de variables non négatives avec une masse de
probabilité à zéro.
Chapitre 5. Mise en oeuvre du krigeage 66
Voyons maintenant comment les données se répartissent dans l’espace. Les figures
5.6 et 5.7 présentent cette répartition spatiale en 3 dimensions et en 2 dimensions
respectivement. Afin de tracer les courbes de niveaux de la figure 5.7, les observations
des stations ont d’abord été interpolées sur les points de la grille GEM par la méthode
de l’inverse de la distance. Le graphique de gauche de la figure 5.6 contient les valeurs
des observations de stations. Parmi les 16 stations, six ont mesuré plus de 0.2 mm de
pluie. Les précipitations sont donc très localisées. Les deux plus fortes observations de
stations sont celles des stations Maniwaki et La Vérendrye. De ces figures, il ressort
que les observations de pluie les plus élevées se situent au nord-est du bassin versant.
Les deux sources de données sont en accord sur ce point. Par contre, l’autre pic de
précipitations n’est pas tout à fait situé au même endroit par les deux sources. Les
stations le situe plutôt au sud tandis que le modèle GEM le situe plutôt à l’ouest
du bassin. À partir des graphiques de gauche des figures 5.6 et 5.7, il est difficile de
tirer des conclusions sur la stationnarité de la variable régionalisée à interpoler. De
plus, le petit nombre de stations rend inutile l’étude de la stationnarité au moyen
d’une fenêtre mobile. Des graphiques de comportement directionnel pourraient peut-
Chapitre 5. Mise en oeuvre du krigeage 67
être aider. La figure 5.8 contient ces graphiques sur lesquels une droite de régression a
été tracée. Comme il avait été observé sur les figures 5.6 et 5.7, les précipitations ont
tendance à augmenter en se déplaçant vers le nord et l’est. Cependant, en moyenne
cette augmentation n’est pas si forte. À la limite, on pourrait donc juger plausible le
postulat de stationnarité de ces données.
Quatre types de krigeages sont effectués ici : le krigeage ordinaire, le krigeage uni-
versel, le krigeage avec dérive externe et le cokrigeage ordinaire. Tous ces krigeages
ont été effectués avec un voisinage composé de toutes les données en raison de leur
faible nombre. Pour le krigeage et le cokrigeage ordinaire, aucun choix ne doit être fait
concernant les tendances des modèles. Il s’agit toujours de constantes inconnues. Pour
les krigeages avec modèle de tendance, des tendances linéaires ont été choisies pour mi-
nimiser le nombre de paramètres du modèle. Ainsi, la tendance est µ(s) = β0 +β1 x+β2 y
en krigeage universel et µ(s) = β0 + β1 w en krigeage avec dérive externe.
Chapitre 5. Mise en oeuvre du krigeage 69
Ensuite, l’analyse variographique est effectuée. La figure 5.9 présente les semi-
variogrammes expérimentaux pour le krigeage ordinaire, le krigeage universel et le
krigeage avec dérive externe. Ces variogrammes sont très semblables. Ainsi, il semble
qu’estimer le semi-variogramme sur les résidus ou sur les données brutes n’influence pas
beaucoup l’allure du semi-variogramme. En cokrigeage, en plus des semi-variogrammes
simples pour les données de stations et les données du modèle GEM, le semi-variogramme
croisé entre ces données doit être calculé. La figure 5.10 présente ces semi-variogrammes.
Celui des données GEM est beaucoup plus structuré que celui des stations. Cette obser-
vation est peu étonnante car le modèle GEM est déterministe. La figure 5.10 contient
aussi les modèles variographiques ajustés aux semi-variogrammes expérimentaux. Ces
modèles sont de la forme linéaire sans palier et sans effet de pépite car des modèles
plus complexes sont difficiles à ajuster avec la librairie gstat de S-Plus. Sur le semi-
variogramme croisé de la figure 5.10, il ressort cependant qu’un modèle avec effet de
pépite aurait été plus juste. Les trois modèles forment un modèle linéaire de corégiona-
lisation, ce qui assure la positivité des variances tel que mentionné à la section 4.5.3.
Pour les autres types de krigeage, le modèle variographique est choisi par valida-
tion croisée. Ainsi, plusieurs modèles sont ajustés et ensuite comparés. Pour illustrer
cette étape, la figure 5.11 présente les ajustements des sept modèles testés pour le kri-
geage ordinaire. Ces sept modèles sont ceux présentés à la section 3.5. Il est difficile de
sélectionner à l’oeil le meilleur modèle. La validation croisée facilite le travail. Elle a ici
sélectionné le modèle sphérique.
La figure 5.12 représente les surfaces d’interpolation qui sont obtenues avec les
différents modèles. Le modèle pépitique donne une surface plane en krigeage ordinaire
car il revient à une régression classique et le modèle de tendance est ici une constante.
Les modèles sphérique, exponentiel et pépitique ont une allure relativement semblable
dans la figure 5.11. En conséquence, ils produisent des surfaces d’interpolations qui se
Chapitre 5. Mise en oeuvre du krigeage 71
Fig. 5.12 – Comparaison des surfaces d’interpolation obtenues avec les différents
modèles variographiques
Chapitre 5. Mise en oeuvre du krigeage 72
ressemblent. La surface produite avec le modèle gaussien ressemble aussi à ces surfaces,
mais elle comporte des pics plus élevés. Finalement, l’utilisation d’un modèle linéaire
sans palier plutôt qu’avec palier a pour effet de lisser davantage la surface d’interpola-
tion.
5.4.3 Krigeage
Une fois le modèle spécifié, l’interpolation peut être effectuée. Les surfaces d’interpo-
lation obtenues des différents types de krigeage effectués sont présentées dans la figure
5.13 et le tableau 5.2 contient des statistiques descriptives sur ces valeurs interpolées.
On remarque tout de suite que les méthodes multivariables produisent de plus grandes
valeurs interpolées, ce à quoi on s’attendait étant donné que les données GEM ont
tendance à être plus élevées que les données de stations. Le cokrigeage accorde plus
d’importance aux données GEM que le krigeage avec dérive externe. En effet, la sur-
face d’interpolation qu’il produit ressemble beaucoup aux données GEM tandis que la
surface du krigeage avec dérive externe se rapproche plus des données de stations. En
outre, le krigeage ordinaire et le krigeage universel produisent des valeurs interpolées
très semblables.
Fig. 5.13 – Comparaison des surfaces d’interpolation obtenues des différents types de
krigeage
Interpolation statistique
multivariable de données de
précipitations dans un cadre de
modélisation hydrologique
Résumé
Divers organismes québécois responsables de la gestion de l’eau se servent
de modèles hydrologiques pour simuler et prévoir le débit des rivières. Les
champs de précipitations observés donnés en entrée à ces modèles sont ob-
tenus par interpolation spatiale des données de stations météorologiques.
Cependant, la faible densité du réseau de stations rend les valeurs inter-
polées incertaines. Pour améliorer la qualité de cet intrant, l’intégration
de nouvelles sources de données est envisagée. Dans le cadre d’un projet
de maı̂trise, des méthodes statistiques d’interpolation spatiale multivariable
ont été étudiées. Une première comparaison théorique a permis d’établir que
les techniques de régression locale et de krigeage étaient les plus promet-
teuses pour produire l’intrant des précipitations observées. Des variantes
de ces techniques ont été testées par validation croisée et à l’aide de si-
mulations hydrologiques. Les données de test provenaient de 16 stations
météorologiques localisées sur le bassin versant de la rivière Gatineau pour
Chapitre 6. Interpolation statistique multivariable de données de précipitations 76
Abstract
Various organizations in charge of water management in Quebec use hydro-
logical models to simulate and forecast rivers discharge. The observed preci-
pitation fields given as input for those models are obtained by spatial inter-
polation of weather stations data. However, the low density of the station’s
network causes noisy interpolated values. To improve the quality of this
input, the integration of new data sources is considered. For the purposes
of a master’s degree, multivariable statistical methods for spatial interpola-
tion were studied. A preliminary theoretical comparison led to the conclu-
sion that the local regression techniques and the kriging techniques were
the most promising to produce the input of observed precipitation fields.
Consequently, variants of these techniques were tested by cross-validation
and hydrological simulations. The test data came from 16 weather stations
located in the Gatineau watershed for August 2003. Precipitation forecasts
generated with the numerical weather prediction model GEM (Environment
Canada) were used as an auxiliary variable. Finally, the proposed approach
consists in interpolating data, for a given period, by the technique among
the examined ones which produces the smallest cross-validation errors. The
tests suggest that this approach performs as well as the inverse distance
method.
Dans le cadre d’un projet de maı̂trise réalisé à l’Université Laval, des méthodes
statistiques d’interpolation spatiale multivariable pouvant résoudre cette problématique
ont été étudiées. Ce projet visait d’abord à recenser et comparer théoriquement les
méthodes. Les techniques de régression locale et de krigeage sont ressorties comme
étant prometteuses et assez simples à implanter. La section 2 de cet article comporte
une brève introduction à ces méthodes. Ensuite, le projet avait pour but de comparer
la performance des méthodes choisies pour interpoler des données de précipitations.
Les travaux visant l’atteinte de cet objectif ont finalement mené au développement
d’une démarche qui semble fournir des prévisions de précipitations un peu plus justes
que celles obtenues de la méthode témoin de l’inverse de la distance. De plus, cette
approche permet d’intégrer à l’interpolation des données provenant d’autres sources
que les stations météorologiques. Les données sur lesquelles cette étude a été menée
sont décrites à la section 3. Ensuite, la méthodologie suivie est détaillée à la section 4,
Chapitre 6. Interpolation statistique multivariable de données de précipitations 78
La régression locale est une méthode de lissage qui permet d’ajuster une surface de
régression (Cleveland et Devlin, 1988). Pour répondre à une problématique d’interpo-
lation spatiale, la variable d’intérêt est modélisée par une fonction linéaire ou quadra-
tique des coordonnées spatiales x et y. Cette fonction est ajustée par la méthode des
moindres carrés pondérés. Ce qui différencie la régression locale des surfaces de tendance
est cette pondération des données, qui est fonction de la distance géographique entre
les sites d’observation et le site pour lequel une prévision est voulue. Certains nomment
« régression pondérée géographiquement » ce type de régression locale permettant de
faire de l’interpolation spatiale (Fotheringham et al., 2002).
Pour utiliser cette technique, il faut préalablement choisir une fonction de poids et
Chapitre 6. Interpolation statistique multivariable de données de précipitations 79
spécifier la taille du voisinage. La fonction de poids détermine dans quelle mesure les
observations les plus proches ont plus d’importance dans l’interpolation. Dans ce projet,
la fonction Epanechnikov est toujours employée. Ainsi, pour effectuer une prévision au
2
point s0 = (x0 , y0 ), le poids de l’observation prise au point si = (xi , yi ) est 1 − |sh(s
0 −si |
0)
pour 0 ≤ |s0 − si |/h(s0 ) < 1, et 0 sinon. L’expression |s0 − si | représente la distance
euclidienne entre les points s0 et si , et h(s0 ) est la distance au-delà de laquelle les
observations se voient accorder un poids nul. Ainsi, h(s0 ) constitue la limite du voisinage
de s0 . Une pratique courante en régression locale est de spécifier ce voisinage par la
fraction des points d’observation que l’utilisateur désire inclure dans la prévision. Plus
cette fraction est grande, plus la surface ajustée est lisse.
6.2.2 Krigeage
En krigeage, la valeur de la variable d’intérêt est prévue en un point par une somme
pondérée des observations ponctuelles disponibles. Les poids des données sont choisis
de façon à ce que l’interpolation soit sans biais et à variance minimale. Cette technique
a été introduite par le Français G. Matheron en 1962 (Matheron, 1962, 1963b). Il s’agit
de la première méthode d’interpolation à tenir compte de la structure de dépendance
spatiale des données. Notons que le krigeage repose sur les mêmes bases théoriques que
l’« interpolation optimale » employée en météorologie (Gandin, 1963).
Il existe plusieurs types de krigeage, qui diffèrent selon la forme postulée pour
l’espérance de la variable d’intérêt. Par exemple, lorsqu’il est supposé que l’espérance
soit constante et connue, on parle de krigeage simple. S’il est postulé que l’espérance
soit constante mais inconnue, il s’agit de krigeage ordinaire. Enfin, le krigeage univer-
sel repose sur l’hypothèse que cette espérance soit une fonction des coordonnées spa-
tiales. Ainsi, ce dernier type de krigeage n’est pas stationnaire par rapport à l’espérance
contrairement aux deux autres.
En plus des stations météorologiques, une autre source de données a été considérée :
le modèle numérique de prévision météorologique, ou modèle atmosphérique, GEM
(Global Environnemental Multi-échelles) d’Environnement Canada. Les données géné-
rées par ce modèle qui ont été employées dans l’interpolation sont des prévisions de
précipitations cumulées sur 6 heures. Ces données se présentent sur une grille de 10km
par 10km dont les points sont représentés sur le graphique de droite de la figure 6.1.
Ce graphique présente aussi le champ moyen des précipitations 6h prévues par GEM
en août 2003. Il indique que le modèle GEM a tendance à surestimer les précipitations
comparativement aux stations.
Les données du modèle atmosphérique GEM ont été choisies parce qu’elles étaient
facilement accessibles et qu’elles présentaient un bon potentiel. Notons cependant que la
Chapitre 6. Interpolation statistique multivariable de données de précipitations 81
Fig. 6.1 – Emplacement des stations météorologiques ainsi que des points de la grille
du modèle GEM sur le bassin versant de la rivière Gatineau et champs moyens de
précipitations observées et prévues pour une période de 6 heures en août 2003 (coor-
données spatiales UTM sur la zone 18T en km)
variable auxiliaire aurait pu provenir d’autres sources, telles que des radars météorolo-
giques au sol ou des capteurs satellitaires. Bien sûr, des traitements préalables adéquats
auraient dû être appliqués à ces données avant de les combiner aux observations des
stations.
6.4 Méthodologie
krigeage ont été appliquées aux données. La moitié de ces techniques sont univariables,
elles n’utilisent donc que les données des stations. Les autres sont multivariables ; elles
intègrent les données du modèle atmosphérique GEM à l’interpolation. L’utilisation de
ces deux groupes de techniques a permis d’évaluer l’apport de l’intégration des données
GEM à l’interpolation. Le tableau suivant décrit les méthodes employées :
Paramètre
Identifiant Méthode Type (krigeage : modèle variographique,
régression locale : fraction voisinage)
Régression locale par
RLxy rapport aux coordonnées univariable Sélection par validation croisée
spatiales x et y
Régression locale par
RLxyw rapport à x et y, ainsi multivariable Sélection par validation croisée
qu’à la variable auxiliaire
KO Krigeage ordinaire univariable Sélection par validation croisée
KU Krigeage universel univariable Sélection par validation croisée
Krigeage
KDE multivariable Sélection par validation croisée
avec dérive externe
CKO Cokrigeage ordinaire multivariable Linéaire sans palier
Les six méthodes du tableau 6.1 ont été employées pour interpoler les données à
chaque période de temps pour laquelle il a plu. En août 2003, cela représente 92 périodes
de 6 heures sur 124. En fait, le krigeage avec dérive externe (KDE) et la régression locale
multivariable (RLxyw) n’ont pu être appliqués lorsque les données GEM étaient trop
peu variables. Ainsi, ces méthodes n’ont été utilisées que pour 63 des 92 périodes de
pluie du mois. En outre, les techniques de régression locale et de krigeage examinées
dans ce projet peuvent produire des prévisions inférieures à zéro, ce qui est inapproprié
pour une variable non négative telle qu’une quantité de précipitations. Les prévisions
négatives ont donc toujours été ramenées à zéro.
Les calculs ont tous été effectués avec le logiciel S-Plus. Les fonctions de la librairie
Chapitre 6. Interpolation statistique multivariable de données de précipitations 84
Gstat (Pebesma, 2004a) ont permis d’exécuter tous les types de krigeage.
6.5 Résultats
La figure 6.2 permet de comparer les erreurs de validation croisée des différentes
procédures d’interpolation. L’identifiant « InvDist » désigne la méthode témoin de
l’inverse de la distance. Les autres méthodes sont identifiées comme dans la section
précédente. Pour les diagrammes en boı̂te, les moustaches représentent le premier et le
neuvième décile, les contours de la boı̂te le premier et le troisième quartile, la ligne cen-
trale la médiane et le cercle la moyenne. L’indice EQM (erreur quadratique moyenne)
a aussi été ajouté à la figure 6.2. Toutes ces statistiques sont calculées à partir des
erreurs de validation croisée de toutes les stations et de toutes les périodes de temps
d’utilisation des méthodes, ce qui représente entre 900 et 1400 données environ.
Fig. 6.2 – Diagrammes en boı̂te des erreurs de validation croisée pour les techniques
d’interpolation spatiale examinées
Les lignes pointillées de la figure 6.2 mettent en évidence que les résidus de l’approche
Sélection ont tendance à être moins extrêmes que ceux de toute autre méthode. En effet,
Chapitre 6. Interpolation statistique multivariable de données de précipitations 85
le diagramme en boı̂te de Sélection possède les plus courtes moustaches. De plus, son
erreur quadratique moyenne est inférieure à toutes les autres. Par contre, la boı̂te du
diagramme en boı̂te de la méthode de l’inverse de la distance est plus petite que celle de
l’approche Sélection. Les erreurs de validation croisée de la méthode témoin sont donc
plus concentrées autour de zéro. En fait, ce grand nombre d’erreurs presque nulles est
associé aux précipitations observées de zéro millimètre de pluie. La méthode de l’inverse
de la distance reproduit bien ces observations, possiblement à cause du petit voisinage
de prévision qu’elle emploie. Dans ce projet, en krigeage, toutes les données étaient
incluses dans l’interpolation. Il serait intéressant de tester l’utilisation des techniques
de krigeage avec un voisinage restreint.
Fig. 6.3 – Débit de la rivière Gatineau aux Rapides Ceizur en août 2003
Cette figure indique que les méthodes Sélection et Inverse de la distance accom-
plissent des performances semblables. Les simulations qu’elles permettent d’obtenir
d’HYDROTEL sont toujours très rapprochées. Ainsi, pour les données de test em-
ployées, la méthode Sélection proposée ne semble pas donner de résultats significative-
ment meilleurs que la méthode de l’inverse de la distance programmée dans HYDRO-
TEL. Par contre, étant donné que les tests ont été effectués sur une courte période et un
seul bassin versant, les résultats ne sont pas généralisables à l’ensemble du territoire du
Québec. Il serait donc intéressant de refaire ces tests pour d’autres périodes présentant
Chapitre 6. Interpolation statistique multivariable de données de précipitations 86
Fig. 6.4 – Fréquences de sélection des méthodes par l’approche Sélection pour les 92
périodes de pluie d’août 2003 catégorisées selon l’intensité des précipitations
Pour la majorité des périodes (au total 57 sur 92), l’approche Sélection a choisi
d’appliquer le krigeage ordinaire (KO) aux données. Il n’est donc pas étonnant que les
erreurs de validation croisée de la procédure Sélection se distinguent peu de celles du
krigeage ordinaire dans la figure 6.2. Ces deux méthodes présentent des diagrammes
en boı̂te très semblables. Elles sont notamment les deux méthodes qui surestiment le
moins les précipitations, présentant les résidus négatifs les moins gros. En outre, il est
pertinent de se demander si l’erreur quadratique moyenne de 7.59 pour la méthode
Sélection est significativement plus petite que celle de 8.05 du krigeage ordinaire. En
fait, il n’est pas étonnant que l’approche Sélection ait un indice EQM plus petit que
toutes les méthodes du tableau 1 considérant le fait qu’à chaque période de temps la
procédure Sélection applique justement la méthode qui minimise le EQM. Ainsi, il ne se-
rait pas vraiment exact de conclure ici que l’approche Sélection accomplit de meilleures
performance que le krigeage ordinaire employé systématique. Cependant, la procédure
Sélection est définitivement plus souple, permettant d’inclure des variables auxiliaires
dans l’interpolation. Elle offre donc un bon potentiel de performance, notamment pour
des bassins versants contenant moins de stations que celui de la rivière Gatineau.
Chapitre 6. Interpolation statistique multivariable de données de précipitations 87
Comme pour les méthodes RLxy et KU, l’emploi du krigeage avec dérive externe
(KDE) ou du cokrigeage ordinaire (CKO) pour interpoler les données n’est avantageux
que s’il est effectué au bon moment. Le cokrigeage ordinaire se distingue surtout pour les
périodes de pluie très faible (précipitations moyennes inférieures à 0.1 mm). L’approche
Sélection l’utilise 22% des périodes de cette classe (figure 6.4). À l’opposé, le krigeage
avec dérive externe semble mieux performer pour des périodes de pluies fortes. La figure
4 indique que cette méthode produit les plus petites erreurs de validation croisée selon
l’indice EQM pour 18% des périodes de précipitations moyennes supérieures à 1 mm.
Cependant, il est étonnant de remarquer que les périodes de temps pour lesquelles une
de ces méthodes multivariables est choisie ne présentent pas nécessairement une forte
corrélation entre les données de stations et les données GEM. En fait, ces corrélations
sont inférieures à 0.5 pour 86% des 92 périodes de pluie. Les grandes étendues des
diagrammes en boı̂te de CKO et de KDE sur la figure 6.2 s’expliquent sûrement par ces
faibles corrélations. L’intégration à l’interpolation de variables auxiliaires plus corrélées
aux observations de stations donnerait possiblement de meilleurs résultats.
spatiale par validation croisée à chaque période de temps. Cette approche est plus
intéressante que de toujours employer la même technique, car elle combine les forces
des différentes méthodes. En particulier, les données de variables auxiliaires peuvent
être intégrées à la prévision par le biais d’une méthode multivariable d’interpolation
si la validation croisée indique que ce serait avantageux de le faire. Pour les données
de test du bassin versant de la Gatineau en août 2003, la méthode proposée accomplit
d’aussi bons résultats que la méthode de l’inverse de la distance. Ces méthodes ont été
comparées par validation croisée et à l’aide du modèle hydrologique de prévision et de
simulation HYDROTEL.
Rappelons qu’une seule variable auxiliaire a été employée ici : les prévisions générées
par le modèle atmosphérique GEM. Il serait intéressant d’étudier l’intégration d’autres
variables auxiliaires, qui fourniraient possiblement des résultats plus notables. L’utilisa-
tion d’une technique permettant d’exploiter la dépendance temporelle entre les périodes
de temps serait aussi une avenue intéressante, qui n’a pas été explorée dans ce projet.
Finalement, il faut noter que les méthodes examinées ici ne sont pas vraiment adaptées
à l’interpolation de variables non négatives prenant souvent des valeurs nulles, comme
les données de précipitations aux 6 heures. Une version lognormale du krigeage assure
des prévisions positives, mais cette technique est déconseillée en raison des difficultés
que la transformation inverse des prévisions présente (Roth, 1998). L’emploi d’autres
approches, notamment des méthodes bayésiennes (Banerjee et al., 2004; Velarde et al.,
2004), serait sûrement indiqué pour des données de précipitations.
Chapitre 7
Conclusion
aussi pour interpoler des données de précipitations des techniques simples de parti-
tionnement de l’espace (e.g. Dirks et al., 1998; Haberlandt et Kite, 1998; Goovaerts,
2000) et de régression classique (e.g. Kruizinga et Yperlaan, 1978; Tabios et Salas,
1985; Abtew et al., 1985), ainsi que des méthodes bayésiennes plus complexes (e.g.
Gaudard et al., 1999; Johns et al., 2001; Velarde et al., 2004). Bref, toutes les classes de
méthodes présentées au chapitre 2 se retrouvent dans la littérature sur l’interpolation de
précipitations. Les techniques qui semblent les plus fréquemment utilisées en pratique
sont les méthodes barycentriques d’inverse de la distance et les techniques classiques de
krigeage tel que le krigeage ordinaire.
Dans la littérature, les résultats des comparaisons diffèrent d’une étude à l’autre. Les
performances des méthodes dépendent de plusieurs facteurs, notamment des résolutions
Chapitre 7. Conclusion 91
temporelle et spatiale des données, ainsi que des paramètres des modèles comme le
semi-variogramme en krigeage. Les articles répertoriés ici traitent de l’analyse de pré-
cipitations annuelles, mensuelles, journalières, horaires ou totales pour des événements
de précipitations de durées quelconques. Les densités des réseaux de stations utilisés
varient d’environ une station pour 3 km2 à une station pour 50 000 km2 . Il est donc
impossible de tirer une conclusion générale et aucune méthode ne ressort comme étant
universellement la meilleure. D’ailleurs, même à l’intérieur d’une étude, les performances
d’une méthode se démarquent rarement nettement des autres. Par exemple, Kruizinga
et Yperlaan (1978) ne préconisent l’emploi d’aucune méthode spécifique car ils n’ob-
servent pas de différences marquées entre les méthodes. D’autres auteurs recommandent
la méthode parmi les meilleures qu’ils jugent la plus pratique. Par exemple, Tabios et
Salas (1985), Abtew et al. (1985) et Syed et al. (2003) notent tous des performances
relativement équivalentes entre le krigeage ordinaire et les fonction multiquadratiques
(type de splines). Cependant, Tabios et Salas (1985) ainsi que Abtew et al. (1985) re-
commandent l’utilisation du krigeage ordinaire car il permet de calculer des erreurs de
prévision, tandis que Syed et al. (2003) choisissent d’employer les fonctions multiquadra-
tiques car ils les considèrent plus faciles d’utilisation. De même, sur un réseau dense de
stations, Dirks et al. (1998) n’obtiennent pas d’améliorations significatives des résultats
en employant un krigeage ordinaire plutôt qu’une méthode de l’inverse de la distance.
Ils recommandent donc la méthode plus simple de l’inverse de la distance. Ainsi, des
résultats comme ceux obtenus au chapitre 6 qui ne mettent pas en évidence la supériorité
de méthodes élaborées comme le krigeage sur des méthodes simples tel l’inverse de la
distance sont usuels avec des données de précipitations. En fait, cette remarque semble
se généraliser aussi à d’autres types de données. Dans une étude expérimentale sur des
données simulées, Zimmerman et al. (1999) obtiennent une meilleure interpolation avec
le krigeage ordinaire ou universel qu’avec la méthode de l’inverse de la distance seule-
ment lorsque le plan d’échantillonnage des données est régulier, le bruit est faible et la
corrélation spatiale est forte. Dans l’étude du chapitre 6, aucune de ces caractéristiques
n’est rencontrée.
défini, approfondi, puis illustré. C’est ainsi que les objectifs du mémoire de comprendre
et de savoir utiliser le krigeage ont été rencontrés.
Il ressort de ce mémoire que le krigeage repose sur des bases mathématiques so-
lides. Ces bases sont plus rigoureuses que celles de la plupart des autres méthodes
d’interpolation. La mise en oeuvre du krigeage est cependant complexifiée par l’étape
de l’analyse variographique. Cette étape, lorsqu’effectuée tel que recommandé, requiert
l’intervention humaine. Dans une problématique telle que celle traitée au chapitre 6 où
l’interpolation doit être complètement automatisée, le krigeage pose donc un problème
méthodologique. La solution proposée ici pour résoudre ce problème est d’utiliser la
validation croisée pour sélectionner le modèle variographique. Par contre, les résultats
du chapitre 6 montrent que la complexité du krigeage n’est pas nécessairement garante
de meilleures performances en pratique par rapport à des méthodes d’interpolation plus
simples.
Afin de poursuivre les travaux de recherche présentés dans ce mémoire, les méthodes
bayésiennes d’interpolation spatiale ainsi que les différentes méthodes d’interpolation
spatio-temporelle pourraient être étudiées. Toutes ces méthodes seraient particulière-
ment intéressante à utiliser pour résoudre la problématique d’interpolation de données
de précipitations traitée au chapitre 6. L’interpolation spatio-temporelle surpasserait
possiblement l’interpolation spatiale pour de fines résolutions temporelles. En effet,
dans l’analyse d’observations d’un phénomène naturel évoluant dans le temps, plus
les intervalles de temps entre les collectes de données sont courts, plus la dépendance
temporelle entre les données a tendance à être forte. Étant donné que les outils de
collecte de données travaillent à des résolutions temporelles de plus en plus fines, la
dépendance temporelle entre les observations est de plus en plus présente dans les
données. Kyriakidis et Journel (1999) proposent une revue très complète des modèles
géostatistiques permettant d’exploiter la dimension temporelle d’un jeu de données.
Bibliographie
Abtew, W., Obeysekera, J. et Shih, G. (1985). Spatial analysis for monthly rainfall
in south florida. Water Resources Bulletin, 29(2):179–188.
Cressie, N. (1984). Towards resistant geostatistics. Dans Verly, G., David, M.,
Journel, A. et Marechal, A., éditeurs : Geostatistics for Natural Resources Cha-
racterization, Part 1, pages 21–44. Reidel, Dordrecht.
Cressie, N. A. C. (1993). Statistics for spatial data. Wiley Series in Probability and
Mathematical Statistics : Applied Probability and Statistics. John Wiley & Sons Inc.,
New York. Revised reprint of the 1991 edition, A Wiley-Interscience Publication.
Côté, J., Gravel, S., Méthot, A., Patoine, A., Roch, M. et Staniforth, A.
(1998). The operational cmc-mrb global environmental multiscale (gem) model : Part
i - design considerations and formulation. Monthly Weather Review, 126:1373–1395.
Fuentes, M. (2000). Fitting algorithms for the semivariogram, and examples. Notes
du cours ST 810M : Spatial Statistics, Lectures 4 à 8, North Carolina State Univer-
sity, Raleigh. Disponibles en ligne : http://www4.stat.ncsu.edu/∼fuentes/st810/
lectures/lectures.html (Page consultée le 28 avril 2005).
BIBLIOGRAPHIE 97
Gaudard, M., Karson, M., Linder, E. et Sinha, D. (1999). Bayesian spatial pre-
diction. Environmental and Ecological Statistics, 6:147–171.
Insightful (2004). Statistics software and data mining software. [En ligne], http:
//www.insightful.com/ (Page consultée le 28 avril 2005).
Johns, C., Nychka, D., Kittel, T. et Daly, C. (2001). Infilling sparse records of
precipitation fields. Journal of the American Statistical Association, 98(464):796–806.
Khuri, A. (1993). Advanced Calculus with Applications in Statistics. John Wiley &
Sons, Inc., New York.
Krige, D. (1951). A statistical approach to some basic mine valuation problems on the
witwatersrand. Journal of the Chemical, Metallurgical and Mining Society, 52:119–
139.
Lee, S., Cho, S. et Wong, P. M. (1998). Rainfall prediction using artificial neural
networks. Journal of Geographic Information and Decision Analysis, 2(2):233–242.
Disponible en ligne : http://www.geodec.org/gida 4.htm (Page consultée le 28
avril 2005).
Loader, C. (1999). Local regression and likelihood. Statistics and Computing. Springer-
Verlag, New York.
Marcotte, D. (1988). Trend surface analysis as a special case of irf-k kriging. Ma-
thematical Geology, 20(7):821–824.
Matheron, G. (1973). The intrinsic random functions and their applications. Advances
in Applied Pobability, 5:439–468.
R project (2004). The r project for statistical computing. [En ligne], http://www.
r-project.org/ (Page consultée le 28 avril 2005).
Ripley, B. D. (1981). Spatial statistics. John Wiley & Sons Inc., New York. Wiley
Series in Probability and Mathematical Statistics.
SMC (2002). Service Météorologique du Canada. [En ligne], version actuelle de GEM
http://www.msc-smc.ec.gc.ca/cmc/op systems/regional forecast f.html, an-
cienne version de GEM http://www.msc-smc.ec.gc.ca/cmc/GEWEX/regional
forecast f.html , présentation de GEM http://www.cmc.ec.gc.ca/rpn/gef
html public/ (Pages consultées le 28 avril 2005).
Wahba, G. (1990). Spline models for observational data, volume 59 de CBMS-NSF Re-
gional Conference Series in Applied Mathematics. Society for Industrial and Applied
Mathematics (SIAM), Philadelphia, PA.
Yang, D., Goodison, B., Metcalfe, J., Louie, P., Leavesley, G., Emerson, D.,
Golubev, V., Elomaa, E., Gunther, T., Hanson, C., Pangburn, T., Kang, E.
et Milkovic, J. (1999). Quantification of precipitation measurement discontinuity
induced by wind shields on national gauge. Water Resources Research, 35(2):491–508.
Disponible en ligne : http://www.uaf.edu/water/faculty/yang/1998WR900042.
pdf (Page consultée le 28 avril 2005).
Dans cette annexe, quelques résultats supplémentaires intéressants relatifs aux tra-
vaux effectués dans le chapitre 6 sont présentés. Tout d’abord, le tableau A.1 indique
le nombre de fois que chaque fraction de voisinage a été choisie par la validation croisée
pour les deux types de régression locale effectués. Pour la régression locale incluant la
variable auxiliaire, seules les fractions 0.8, 0.9 et 1 ont été considérées car le modèle
contient un paramètre de plus qu’en régression locale univariée. Un plus grand nombre
de données sont donc requises dans l’ajustement du modèle. Probablement en raison
du petit nombre d’observations analysées, la fraction choisie le plus fréquemment est 1
pour les deux types de régression locale.
Les fréquences de sélection des modèles variographiques par la validation croisée sont
présentées dans le tableau A.2 en fonction du type de krigeage. Le modèle pépitique est
celui le plus souvent choisi pour les trois types de krigeage. Afin de mieux comprendre
ce résultat, regardons les fréquences de sélection des modèles en fonction de l’intensité
Annexe A. Complément d’analyse des résultats du chapitre 6 104
Modèle variographique
Exponentiel
avec palier
sans palier
Sphérique
Puissance
Pépitique
Gaussien
Linéaire
Linéaire
Type de krigeage Total
36 20 8 4 15 6 3
ordinaire 92
0.39 0.22 0.09 0.04 0.16 0.07 0.03
36 23 6 8 8 9 2
universel 92
0.39 0.25 0.07 0.09 0.09 0.10 0.02
22 16 2 3 15 4 1
avec dérive externe 63
0.35 0.25 0.03 0.05 0.25 0.06 0.02
des précipitations. Le tableau A.3 présente ces fréquences pour le krigeage ordinaire. On
remarque que le modèle pépitique est surtout retenu lorsque les précipitations sont de
faible ou moyenne intensité. Lorsque les précipitations moyennes sur le bassin versant
sont supérieures à 1 mm, les modèles les plus souvent sélectionnés sont plutôt le modèle
linéaire avec palier et le modèle gaussien.
Modèle variographique
Exponentiel
avec palier
sans palier
Sphérique
Puissance
Pépitique
Gaussien
Linéaire
Linéaire
Précipitations moyennes 24 5 2 0 4 2 0
37
inférieures à 0.1 mm 0.65 0.14 0.05 0.00 0.11 0.05 0
Précipitations moyennes 11 7 5 2 5 1 2
33
entre 0.1 mm et 1 mm 0.33 0.21 0.15 0.06 0.15 0.03 0.06
Précipitations moyennes 1 8 1 2 6 3 1
22
supérieures à 1 mm 0.05 0.36 0.05 0.09 0.27 0.14 0.05
Total 36 20 8 4 15 6 3 92
des erreurs de validation croisée par intensité de précipitations pour la méthode témoin
de l’inverse de la distance et pour la méthode Sélection qui a été développée. Bien que la
méthode Sélection ne semble pas donner des prévisions vraiment plus justes que celles
générées par la méthode de l’inverse de la distance lorsque les précipitations sont de
faible ou moyenne intensité, on note une amélioration des prévisions par la méthode
Sélection dans le cas de précipitations plus fortes. En effet, pour des précipitations
moyennes sur le bassin versant supérieures à 1 mm, le diagramme en boı̂te des er-
reurs de validation croisée de la méthode Sélection a des moustaches plus courtes et
une boı̂te moins longue que le diagramme en boı̂te pour la méthode de l’inverse de
la distance. Dans une perspective de prévision hydrologique, ce sont ces périodes de
fortes précipitations qui ont le plus d’importance. Ainsi, la méthode Sélection proposée
améliore l’interpolation spatiale comparativement à la méthode de l’inverse de la dis-
tance lorsque ça importe le plus.
Fig. A.1 – Diagrammes en boı̂te des erreurs de validation croisée catégorisées selon
l’intensité des précipitations pour la méthode de l’inverse de la distance et la méthode
Sélection
Annexe A. Complément d’analyse des résultats du chapitre 6 106
Tab. A.4 – Tests des rangs signés de Wilcoxon pour comparer les EQM de validation
croisée par station en fonction de l’intensité des précipitations
signifie que les EQM pour la méthode Sélection sont plus petits que pour la méthode de
l’inverse de la distance. Peu de tests sont significatifs au seuil de 5%. Les tests significa-
tifs pour des précipitations de faible intensité indiquent tous que la méthode de l’inverse
de la distance performe significativement mieux que la méthode Sélection, tandis que
le test significatif pour des précipitations de forte intensité indique le contraire.
Annexe B
Programme S-Plus
Voici le programme S-Plus qui a permis d’interpoler les données pour chaque période
de 6 heures d’août 2003 par la méthode de l’inverse de la distance, la régression locale par
rapport à x et y, la régression locale par rapport à x, y et les données GEM, le krigeage
ordinaire, le krigeage universel, le krigeage avec dérive externe et le cokrigeage ordinaire.
Ce programme appelle les fichiers Donnees et Grille. Le premier fichier contient les
valeurs régionalisées mesurées aux stations et générées par le modèle GEM pour toutes
les périodes du mois d’août 2003 et tous les sites d’observation. Les coordonnées spatiales
UTM sur la zone 18T en km des sites d’observation sont aussi incluses dans le fichier. Le
deuxième fichier comprend plutôt les coordonnées spatiales UTM des sites de prévision
ainsi que les données fournies par GEM en ces sites. Afin de décrire la forme de ces
fichiers, les tableaux B.1 et B.2 en contiennent des extraits.
Notons que ce programme SPLUS est fonctionnel, mais non optimisé. Il est donc
assez long et ce principalement en raison de la technique de sélection des modèles
variographiques par validation croisée utilisée.
Periode x y GEM
1 430.000 5050.000 0.00
1 440.000 5050.000 0.00
1 410.000 5060.000 0.00
1 420.000 5060.000 0.00
.. .. .. ..
. . . .
Programme :
# Prévisions :
Prev=data.frame(matrix(rep(NA,dim(Grille)[1]*32),nrow=dim(Grille)[1]))
dimnames(Prev)<-list(1:dim(Prev)[1],c("Periode","x","y","InvDist","KO1","KO2","KO3","KO4","KO5",
"KO6","KO7","KO","KU1","KU2","KU3","KU4","KU5","KU6","KU7","KU","KDE1","KDE2",
"KDE3","KDE4","KDE5","KDE6","KDE7","KDE","CKO","RLxy","RLxyw","Selection"))
Prev[,1:3]=Grille[,1:3]
# Variances de prévision :
Var=data.frame(matrix(rep(NA,dim(Grille)[1]*32),nrow=dim(Grille)[1]))
dimnames(Var)<-list(1:dim(Var)[1],c("Periode","x","y","InvDist","KO1","KO2","KO3","KO4","KO5",
"KO6","KO7","KO","KU1","KU2","KU3","KU4","KU5","KU6","KU7","KU","KDE1","KDE2",
"KDE3","KDE4","KDE5","KDE6","KDE7","KDE","CKO","RLxy","RLxyw","Selection"))
Var[,1:3]=Grille[,1:3]
# Informations Diverses :
Info=data.frame(matrix(rep(NA,123*137),nrow=123))
dimnames(Info)<-list(1:dim(Info)[1],c("NbStations","Nbplus0","Classe","MoyStations",
"EtypeStations","MoyInvDist","EtypeInvDist","MoyGEM","EtypeGEM","InvDistEQM","InvDistEAM",
"KO1pepite","KO1EQM","KO1EAM","KO2pepite","KO2ppalier","KO2portee","KO2EQM","KO2EAM",
"KO3pepite","KO3ppalier","KO3portee","KO3EQM","KO3EAM","KO4pepite","KO4ppalier","KO4porteep",
"KO4EQM","KO4EAM","KO5pepite","KO5ppalier","KO5porteep","KO5EQM","KO5EAM","KO6pepite",
"KO6pente","KO6EQM","KO6EAM","KO7pepite","KO7echelle","KO7puissance","KO7EQM","KO7EAM",
"KOmodele","KOEQM","KOEAM","KUcorrX","KUpvalueX","KUcorrY","KUpvalueY","KU1pepite","KU1EQM",
"KU1EAM","KU2pepite","KU2ppalier","KU2portee","KU2EQM","KU2EAM","KU3pepite","KU3ppalier",
"KU3portee","KU3EQM","KU3EAM","KU4pepite","KU4ppalier","KU4porteep","KU4EQM","KU4EAM",
"KU5pepite","KU5ppalier","KU5porteep","KU5EQM","KU5EAM","KU6pepite","KU6pente","KU6EQM",
"KU6EAM","KU7pepite","KU7echelle","KU7puissance","KU7EQM","KU7EAM","KUmodele","KUEQM","KUEAM",
"KDEpossible","KDEcorrGEM","KDEpvalueGEM","KDE1pepite","KDE1EQM","KDE1EAM","KDE2pepite",
"KDE2ppalier","KDE2portee","KDE2EQM","KDE2EAM","KDE3pepite","KDE3ppalier","KDE3portee",
"KDE3EQM","KDE3EAM","KDE4pepite","KDE4ppalier","KDE4porteep","KDE4EQM","KDE4EAM","KDE5pepite",
"KDE5ppalier","KDE5porteep","KDE5EQM","KDE5EAM","KDE6pepite","KDE6pente","KDE6EQM","KDE6EAM",
Annexe B. Programme S-Plus 109
"KDE7pepite","KDE7echelle","KDE7puissance","KDE7EQM","KDE7EAM","KDEmodele","KDEEQM","KDEEAM",
"CKOpenteStations","CKOpenteGEM","CKOpenteStaGEM","CKOEQM","CKOEAM","RLxyfen","RLxyEQM",
"RLxyEAM","RLxywfen","RLxywEQM","RLxywEAM","Selection","SelectionEQM","SelectionEAM"))
library(gstat)
# FILTRE : Mise de c^
oté des périodes de temps pour lesquelles aucune station n’a mesuré de pluie
Pluie=na.omit(ifelse(Info$Classe=="NA",row(Info),NA))
# INVERSE DE LA DISTANCE
for (t in Pluie) {
Periodet=na.omit(ifelse(Donnees$Periode==t,row(Donnees),NA))
# Validation croisée :
for (i in 1:length(Periodet)){
distances=sqrt((Donnees$x[Periodet[i]]-Donnees$x[Periodet[-i]])**2
+(Donnees$y[Periodet[i]]-Donnees$y[Periodet[-i]])**2)
voisinage=na.omit(ifelse(rank(distances)<=3,Periodet[-i],NA))
distances=na.omit(ifelse(rank(distances)<=3,distances,NA))
poids=as.vector((1/distances)/sum(1/distances))
valid=sum(poids*Donnees$Stations[voisinage],na.rm=T)
Err$InvDist[Periodet[i]]=Donnees$Stations[Periodet[i]]-valid
}
ErrValid=data.frame(matrix(rep(NA,dim(Donnees)[1]*7),nrow=dim(Donnees)[1]))
for (t in Pluie) {
ErrValid=data.frame(matrix(rep(NA,dim(Donnees)[1]*3),nrow=dim(Donnees)[1]))
for (t in KDEpossible) {
Periodet=na.omit(ifelse(Donnees$Periode==t,row(Donnees),NA))
j=1
for(a in c(1,0.9,0.8)){
k=trunc(a*Info$NbStations[t])
for (i in 1:length(Periodet)){
distances=sqrt((Donnees$x[Periodet[i]]-Donnees$x[Periodet[-i]])**2+(Donnees$y[Periodet[i]]
-Donnees$y[Periodet[-i]])**2)
voisinage=na.omit(ifelse(rank(distances)<=k,Periodet[-i],NA))
distances=na.omit(ifelse(rank(distances)<=k,distances,NA))
poids=as.vector(1-(distances/ceiling(max(distances)))**2)
reg=lm(Stations~1+x+y+GEM,data=Donnees[voisinage,],weights=poids)
valid=predict(reg,Donnees[Periodet[i],])
ErrValid[Periodet[i],j]=Donnees$Stations[Periodet[i]]-ifelse(valid<0,0,valid)
}
j=j+1
}
indice=data.frame(apply(ErrValid[Periodet,]**2,2,mean))
minimums=na.omit(ifelse(indice==min(indice,na.rm=T),row(indice),NA))
Info$RLxywfen[t]=if(minimums[1,1]==1)1 else if(minimums[1,1]==2)0.9 else if(minimums[1,1]==3)0.8
Err$RLxyw[Periodet]=ErrValid[Periodet,minimums[1,1]]
Info$RLxywEQM[t]=mean((Err$RLxyw[Err$Periode==t])**2,na.rm=T)
Info$RLxywEAM[t]=median(abs(Err$RLxyw[Err$Periode==t]),na.rm=T)
# KRIGEAGE ORDINAIRE
for(t in Pluie){
# 1- Modèle pépitique :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele1=fit.variogram(vario,vgm(max(vario[,3]),"Nug"),fit.method=1)
Info$KO1pepite[t]= modele1$psill[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele1[dim(modele1)[1],2]>0) {
valid=krige.cv(Stations~1,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,])
Err$KO1[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,],
Annexe B. Programme S-Plus 112
newdata=Grille[Grille$Periode==t,])
Prev$KO1[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KO1[Var$Periode==t]=Interpo[,4]
} else {
Err$KO1[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KO1[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KO1[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KO1EQM[t]=mean((Err$KO1[Err$Periode==t])**2,na.rm=T)
Info$KO1EAM[t]=median(abs(Err$KO1[Err$Periode==t]),na.rm=T)
# 3- Modèle sphérique :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele3=fit.variogram(vario,vgm(max(vario[,3]),"Sph",140,0.1*max(vario[,3])),fit.method=1)
modele3=if (modele3[1,2]<0 || modele3[2,2]<0 || modele3[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Sph",140),fit.method=1) else modele3
Info$KO3pepite[t]=if(dim(modele3)[1]==2) modele3$psill[1] else 0
Info$KO3ppalier[t]=if(dim(modele3)[1]==2) modele3$psill[2] else modele3$psill[1]
Info$KO3portee[t]=if(dim(modele3)[1]==2) modele3$range[2] else modele3$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele3[dim(modele3)[1],2]>0 && modele3[dim(modele3)[1],3]>0) {
valid=krige.cv(Stations~1,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,])
Err$KO3[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KO3[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KO3[Var$Periode==t]=Interpo[,4]
} else {
Err$KO3[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KO3[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KO3[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Annexe B. Programme S-Plus 113
}
Info$KO3EQM[t]=mean((Err$KO3[Err$Periode==t])**2,na.rm=T)
Info$KO3EAM[t]=median(abs(Err$KO3[Err$Periode==t]),na.rm=T)
# 4- Modèle exponentiel :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele4=fit.variogram(vario,vgm(max(vario[,3]),"Exp",140,0.1*max(vario[,3])),fit.method=1)
modele4=if (modele4[1,2]<0 || modele4[2,2]<0 || modele4[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Exp",140),fit.method=1) else modele4
Info$KO4pepite[t]=if(dim(modele4)[1]==2) modele4$psill[1] else 0
Info$KO4ppalier[t]=if(dim(modele4)[1]==2) modele4$psill[2] else modele4$psill[1]
Info$KO4porteep[t]=if(dim(modele4)[1]==2) 3*modele4$range[2] else 3*modele4$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele4[dim(modele4)[1],2]>0 && modele4[dim(modele4)[1],3]>0) {
valid=krige.cv(Stations~1,~x+y,model=modele4,data=Donnees[Donnees$Periode==t,])
Err$KO4[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1,~x+y,model=modele4,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KO4[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KO4[Var$Periode==t]=Interpo[,4]
} else {
Err$KO4[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KO4[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KO4[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KO4EQM[t]=mean((Err$KO4[Err$Periode==t])**2,na.rm=T)
Info$KO4EAM[t]=median(abs(Err$KO4[Err$Periode==t]),na.rm=T)
# 5- Modèle gaussien :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele5=fit.variogram(vario,vgm(max(vario[,3]),"Gau",140,0.1*max(vario[,3])),fit.method=1)
modele5=if (modele5[1,2]<0 || modele5[2,2]<0 || modele5[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Gau",140),fit.method=1) else modele5
Info$KO5pepite[t]=if(dim(modele5)[1]==2) modele5$psill[1] else 0
Info$KO5ppalier[t]=if(dim(modele5)[1]==2) modele5$psill[2] else modele5$psill[1]
Info$KO5porteep[t]=if(dim(modele5)[1]==2) sqrt(3)*modele5$range[2] else sqrt(3)*modele5$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele5[dim(modele5)[1],2]>0 && modele5[dim(modele5)[1],3]>0) {
valid=krige.cv(Stations~1,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,])
Err$KO5[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KO5[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KO5[Var$Periode==t]=Interpo[,4]
} else {
Err$KO5[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KO5[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KO5[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KO5EQM[t]=mean((Err$KO5[Err$Periode==t])**2,na.rm=T)
Info$KO5EAM[t]=median(abs(Err$KO5[Err$Periode==t]),na.rm=T)
# 7- Modèle puissance :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele7=fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1,0.1*max(vario[,3])),fit.method=1)
modele7=if (modele7[1,2]<0 || modele7[2,2]<0 || modele7[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1),fit.method=1) else modele7
Info$KO7pepite[t]=if(dim(modele7)[1]==2) modele7$psill[1] else 0
Info$KO7echelle[t]=if(dim(modele7)[1]==2) modele7$psill[2] else modele7$psill[1]
Info$KO7puissance[t]=if(dim(modele7)[1]==2) modele7$range[2] else modele7$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if(modele7[dim(modele7)[1],2]>0&&modele7[dim(modele7)[1],3]>0&&modele7[dim(modele7)[1],3]<2){
valid=krige.cv(Stations~1,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,])
Err$KO7[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KO7[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KO7[Var$Periode==t]=Interpo[,4]
} else {
Err$KO7[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KO7[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KO7[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KO7EQM[t]=mean((Err$KO7[Err$Periode==t])**2,na.rm=T)
Info$KO7EAM[t]=median(abs(Err$KO7[Err$Periode==t]),na.rm=T)
# Choix du modèle :
indice=data.frame(c(Info$KO1EQM[t],Info$KO2EQM[t],Info$KO3EQM[t],Info$KO4EQM[t],
Info$KO5EQM[t],Info$KO6EQM[t],Info$KO7EQM[t]))
minimums=na.omit(ifelse(indice==min(indice,na.rm=T),row(indice),NA))
Info$KOmodele[t]=minimums[1,1]
if (Info$KOmodele[t]==1) {
Err$KO[Err$Periode==t]=Err$KO1[Err$Periode==t]
Prev$KO[Prev$Periode==t]=Prev$KO1[Prev$Periode==t]
Var$KO[Var$Periode==t]=Var$KO1[Var$Periode==t]
} else
if (Info$KOmodele[t]==2) {
Annexe B. Programme S-Plus 115
Err$KO[Err$Periode==t]=Err$KO2[Err$Periode==t]
Prev$KO[Prev$Periode==t]=Prev$KO2[Prev$Periode==t]
Var$KO[Var$Periode==t]=Var$KO2[Var$Periode==t]
} else
if (Info$KOmodele[t]==3) {
Err$KO[Err$Periode==t]=Err$KO3[Err$Periode==t]
Prev$KO[Prev$Periode==t]=Prev$KO3[Prev$Periode==t]
Var$KO[Var$Periode==t]=Var$KO3[Var$Periode==t]
} else
if (Info$KOmodele[t]==4) {
Err$KO[Err$Periode==t]=Err$KO4[Err$Periode==t]
Prev$KO[Prev$Periode==t]=Prev$KO4[Prev$Periode==t]
Var$KO[Var$Periode==t]=Var$KO4[Var$Periode==t]
} else
if (Info$KOmodele[t]==5) {
Err$KO[Err$Periode==t]=Err$KO5[Err$Periode==t]
Prev$KO[Prev$Periode==t]=Prev$KO5[Prev$Periode==t]
Var$KO[Var$Periode==t]=Var$KO5[Var$Periode==t]
} else
if (Info$KOmodele[t]==6) {
Err$KO[Err$Periode==t]=Err$KO6[Err$Periode==t]
Prev$KO[Prev$Periode==t]=Prev$KO6[Prev$Periode==t]
Var$KO[Var$Periode==t]=Var$KO6[Var$Periode==t]
} else {
Err$KO[Err$Periode==t]=Err$KO7[Err$Periode==t]
Prev$KO[Prev$Periode==t]=Prev$KO7[Prev$Periode==t]
Var$KO[Var$Periode==t]=Var$KO7[Var$Periode==t]
}
Info$KOEQM[t]=mean((Err$KO[Err$Periode==t])**2,na.rm=T)
Info$KOEAM[t]=median(abs(Err$KO[Err$Periode==t]),na.rm=T)
}
# KRIGEAGE UNIVERSEL :
for(t in Pluie){
corr=cor.test(Donnees$Stations[Donnees$Periode==t],Donnees$x[Donnees$Periode==t],
alternative="two.sided",method="pearson")
Info$KUcorrX[t]=corr$estimate
Info$KUpvalueX[t]=corr$p.value
corr=cor.test(Donnees$Stations[Donnees$Periode==t],Donnees$y[Donnees$Periode==t],
alternative="two.sided",method="pearson")
Info$KUcorrY[t]=corr$estimate
Info$KUpvalueY[t]=corr$p.value
# Première itération :
# 1- Modèle pépitique :
modele1=fit.variogram(vario,vgm(max(vario[,3]),"Nug"),fit.method=1)
if (modele1[dim(modele1)[1],2]>0) {
valid=krige.cv(Stations~1+x+y,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,])
Annexe B. Programme S-Plus 116
Err$KU1[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KU1[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KU1EQM[t]=mean((Err$KU1[Err$Periode==t])**2,na.rm=T)
Info$KU1EAM[t]=median(abs(Err$KU1[Err$Periode==t]),na.rm=T)
# 3- Modèle sphérique :
modele3=fit.variogram(vario,vgm(max(vario[,3]),"Sph",140,0.1*max(vario[,3])),fit.method=1)
modele3=if (modele3[1,2]<0 || modele3[2,2]<0 || modele3[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Sph",140),fit.method=1) else modele3
if (modele3[dim(modele3)[1],2]>0 && modele3[dim(modele3)[1],3]>0) {
valid=krige.cv(Stations~1+x+y,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,])
Err$KU3[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KU3[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KU3EQM[t]=mean((Err$KU3[Err$Periode==t])**2,na.rm=T)
Info$KU3EAM[t]=median(abs(Err$KU3[Err$Periode==t]),na.rm=T)
# 4- Modèle exponentiel :
modele4=fit.variogram(vario,vgm(max(vario[,3]),"Exp",140,0.1*max(vario[,3])),fit.method=1)
modele4=if (modele4[1,2]<0 || modele4[2,2]<0 || modele4[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Exp",140),fit.method=1) else modele4
if (modele4[dim(modele4)[1],2]>0 && modele4[dim(modele4)[1],3]>0) {
valid=krige.cv(Stations~1+x+y,~x+y,model=modele4,data=Donnees[Donnees$Periode==t,])
Err$KU4[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KU4[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KU4EQM[t]=mean((Err$KU4[Err$Periode==t])**2,na.rm=T)
Info$KU4EAM[t]=median(abs(Err$KU4[Err$Periode==t]),na.rm=T)
# 5- Modèle gaussien :
modele5=fit.variogram(vario,vgm(max(vario[,3]),"Gau",140,0.1*max(vario[,3])),fit.method=1)
modele5=if (modele5[1,2]<0 || modele5[2,2]<0 || modele5[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Gau",140),fit.method=1) else modele5
if (modele5[dim(modele5)[1],2]>0 && modele5[dim(modele5)[1],3]>0) {
valid=krige.cv(Stations~1+x+y,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,])
Err$KU5[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KU5[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KU5EQM[t]=mean((Err$KU5[Err$Periode==t])**2,na.rm=T)
Info$KU5EAM[t]=median(abs(Err$KU5[Err$Periode==t]),na.rm=T)
Err$KU6[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KU6[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KU6EQM[t]=mean((Err$KU6[Err$Periode==t])**2,na.rm=T)
Info$KU6EAM[t]=median(abs(Err$KU6[Err$Periode==t]),na.rm=T)
# 7- Modèle puissance :
modele7=fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1,0.1*max(vario[,3])),fit.method=1)
modele7=if (modele7[1,2]<0 || modele7[2,2]<0 || modele7[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1),fit.method=1) else modele7
if(modele7[dim(modele7)[1],2]>0&&modele7[dim(modele7)[1],3]>0&&modele7[dim(modele7)[1],3]<2){
valid=krige.cv(Stations~1+x+y,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,])
Err$KU7[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KU7[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KU7EQM[t]=mean((Err$KU7[Err$Periode==t])**2,na.rm=T)
Info$KU7EAM[t]=median(abs(Err$KU7[Err$Periode==t]),na.rm=T)
# Choix du modèle :
indice=data.frame(c(Info$KU1EQM[t],Info$KU2EQM[t],Info$KU3EQM[t],Info$KU4EQM[t],
Info$KU5EQM[t],Info$KU6EQM[t],Info$KU7EQM[t]))
minimums=na.omit(ifelse(indice==min(indice,na.rm=T),row(indice),NA))
Info$KUmodele[t]=minimums[1,1]
modele = if (Info$KUmodele[t]==1) modele1 else if (Info$KUmodele[t]==2) modele2 else
if (Info$KUmodele[t]==3) modele3 else if (Info$KUmodele[t]==4) modele4 else
if (Info$KUmodele[t]==5) modele5 else if (Info$KUmodele[t]==6) modele6 else modele7
# Deuxième itération :
# 1- Modèle pépitique :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele1=fit.variogram(vario,vgm(max(vario[,3]),"Nug"),fit.method=1)
Info$KU1pepite[t]= modele1$psill[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele1[dim(modele1)[1],2]>0) {
valid=krige.cv(Stations~1+x+y,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,])
Err$KU1[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+x+y,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KU1[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KU1[Var$Periode==t]=Interpo[,4]
} else {
Err$KU1[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KU1[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KU1[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KU1EQM[t]=mean((Err$KU1[Err$Periode==t])**2,na.rm=T)
Info$KU1EAM[t]=median(abs(Err$KU1[Err$Periode==t]),na.rm=T)
Annexe B. Programme S-Plus 118
# 3- Modèle sphérique :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele3=fit.variogram(vario,vgm(max(vario[,3]),"Sph",140,0.1*max(vario[,3])),fit.method=1)
modele3=if (modele3[1,2]<0 || modele3[2,2]<0 || modele3[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Sph",140),fit.method=1) else modele3
Info$KU3pepite[t]=if(dim(modele3)[1]==2) modele3$psill[1] else 0
Info$KU3ppalier[t]=if(dim(modele3)[1]==2) modele3$psill[2] else modele3$psill[1]
Info$KU3portee[t]=if(dim(modele3)[1]==2) modele3$range[2] else modele3$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele3[dim(modele3)[1],2]>0 && modele3[dim(modele3)[1],3]>0) {
valid=krige.cv(Stations~1+x+y,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,])
Err$KU3[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+x+y,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KU3[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KU3[Var$Periode==t]=Interpo[,4]
} else {
Err$KU3[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KU3[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KU3[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KU3EQM[t]=mean((Err$KU3[Err$Periode==t])**2,na.rm=T)
Info$KU3EAM[t]=median(abs(Err$KU3[Err$Periode==t]),na.rm=T)
# 4- Modèle exponentiel :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele4=fit.variogram(vario,vgm(max(vario[,3]),"Exp",140,0.1*max(vario[,3])),fit.method=1)
modele4=if (modele4[1,2]<0 || modele4[2,2]<0 || modele4[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Exp",140),fit.method=1) else modele4
Info$KU4pepite[t]=if(dim(modele4)[1]==2) modele4$psill[1] else 0
Annexe B. Programme S-Plus 119
# 5- Modèle gaussien :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele5=fit.variogram(vario,vgm(max(vario[,3]),"Gau",140,0.1*max(vario[,3])),fit.method=1)
modele5=if (modele5[1,2]<0 || modele5[2,2]<0 || modele5[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Gau",140),fit.method=1) else modele5
Info$KU5pepite[t]=if(dim(modele5)[1]==2) modele5$psill[1] else 0
Info$KU5ppalier[t]=if(dim(modele5)[1]==2) modele5$psill[2] else modele5$psill[1]
Info$KU5porteep[t]=if(dim(modele5)[1]==2) sqrt(3)*modele5$range[2] else sqrt(3)*modele5$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele5[dim(modele5)[1],2]>0 && modele5[dim(modele5)[1],3]>0) {
valid=krige.cv(Stations~1+x+y,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,])
Err$KU5[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+x+y,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KU5[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KU5[Var$Periode==t]=Interpo[,4]
} else {
Err$KU5[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KU5[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KU5[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KU5EQM[t]=mean((Err$KU5[Err$Periode==t])**2,na.rm=T)
Info$KU5EAM[t]=median(abs(Err$KU5[Err$Periode==t]),na.rm=T)
Prev$KU6[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KU6[Var$Periode==t]=Interpo[,4]
} else {
Err$KU6[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KU6[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KU6[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KU6EQM[t]=mean((Err$KU6[Err$Periode==t])**2,na.rm=T)
Info$KU6EAM[t]=median(abs(Err$KU6[Err$Periode==t]),na.rm=T)
# 7- Modèle puissance :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele7=fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1,0.1*max(vario[,3])),fit.method=1)
modele7=if (modele7[1,2]<0 || modele7[2,2]<0 || modele7[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1),fit.method=1) else modele7
Info$KU7pepite[t]=if(dim(modele7)[1]==2) modele7$psill[1] else 0
Info$KU7echelle[t]=if(dim(modele7)[1]==2) modele7$psill[2] else modele7$psill[1]
Info$KU7puissance[t]=if(dim(modele7)[1]==2) modele7$range[2] else modele7$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if(modele7[dim(modele7)[1],2]>0&&modele7[dim(modele7)[1],3]>0&&modele7[dim(modele7)[1],3]<2){
valid=krige.cv(Stations~1+x+y,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,])
Err$KU7[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+x+y,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KU7[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KU7[Var$Periode==t]=Interpo[,4]
} else {
Err$KU7[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KU7[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KU7[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KU7EQM[t]=mean((Err$KU7[Err$Periode==t])**2,na.rm=T)
Info$KU7EAM[t]=median(abs(Err$KU7[Err$Periode==t]),na.rm=T)
# Choix du modèle :
indice=data.frame(c(Info$KU1EQM[t],Info$KU2EQM[t],Info$KU3EQM[t],Info$KU4EQM[t],
Info$KU5EQM[t],Info$KU6EQM[t],Info$KU7EQM[t]))
minimums=na.omit(ifelse(indice==min(indice,na.rm=T),row(indice),NA))
Info$KUmodele[t]=minimums[1,1]
if (Info$KUmodele[t]==1) {
Err$KU[Err$Periode==t]=Err$KU1[Err$Periode==t]
Prev$KU[Prev$Periode==t]=Prev$KU1[Prev$Periode==t]
Var$KU[Var$Periode==t]=Var$KU1[Var$Periode==t]
} else
if (Info$KUmodele[t]==2) {
Err$KU[Err$Periode==t]=Err$KU2[Err$Periode==t]
Prev$KU[Prev$Periode==t]=Prev$KU2[Prev$Periode==t]
Var$KU[Var$Periode==t]=Var$KU2[Var$Periode==t]
} else
if (Info$KUmodele[t]==3) {
Err$KU[Err$Periode==t]=Err$KU3[Err$Periode==t]
Prev$KU[Prev$Periode==t]=Prev$KU3[Prev$Periode==t]
Var$KU[Var$Periode==t]=Var$KU3[Var$Periode==t]
} else
if (Info$KUmodele[t]==4) {
Annexe B. Programme S-Plus 121
Err$KU[Err$Periode==t]=Err$KU4[Err$Periode==t]
Prev$KU[Prev$Periode==t]=Prev$KU4[Prev$Periode==t]
Var$KU[Var$Periode==t]=Var$KU4[Var$Periode==t]
} else
if (Info$KUmodele[t]==5) {
Err$KU[Err$Periode==t]=Err$KU5[Err$Periode==t]
Prev$KU[Prev$Periode==t]=Prev$KU5[Prev$Periode==t]
Var$KU[Var$Periode==t]=Var$KU5[Var$Periode==t]
} else
if (Info$KUmodele[t]==6) {
Err$KU[Err$Periode==t]=Err$KU6[Err$Periode==t]
Prev$KU[Prev$Periode==t]=Prev$KU6[Prev$Periode==t]
Var$KU[Var$Periode==t]=Var$KU6[Var$Periode==t]
} else {
Err$KU[Err$Periode==t]=Err$KU7[Err$Periode==t]
Prev$KU[Prev$Periode==t]=Prev$KU7[Prev$Periode==t]
Var$KU[Var$Periode==t]=Var$KU7[Var$Periode==t]
}
Info$KUEQM[t]=mean((Err$KU[Err$Periode==t])**2,na.rm=T)
Info$KUEAM[t]=median(abs(Err$KU[Err$Periode==t]),na.rm=T)
}
# FILTRE : Mise de c^oté des périodes pour lesquelles les variances leave-one-out
# des prévisions GEM ne sont pas assez grandes pour KDE
for(t in Pluie){
IDlignes=na.omit(ifelse(Donnees$Periode==t,row(Donnees),NA))
PetiteVar=rep(NA,length(IDlignes))
for (i in 1:length(IDlignes)){
PetiteVar[i]=ifelse(var(Donnees$GEM[IDlignes[-i]])<0.0001,1,0)
}
Info$KDEpossible[t]= if (sum(PetiteVar)>0) "non" else "oui"
}
KDEpossible=na.omit(ifelse(Info$KDEpossible=="oui",row(Info),NA))
for(t in KDEpossible){
corr=cor.test(Donnees$Stations[Donnees$Periode==t],Donnees$GEM[Donnees$Periode==t],
alternative="two.sided", method="pearson")
Info$KDEcorrGEM[t]=corr$estimate
Info$KDEpvalueGEM[t]=corr$p.value
# Première itération :
# 1- Modèle pépitique :
modele1=fit.variogram(vario,vgm(max(vario[,3]),"Nug"),fit.method=1)
if (modele1[dim(modele1)[1],2]>0) {
valid=krige.cv(Stations~1+GEM,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,])
Err$KDE1[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Annexe B. Programme S-Plus 122
} else Err$KDE1[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KDE1EQM[t]=mean((Err$KDE1[Err$Periode==t])**2,na.rm=T)
Info$KDE1EAM[t]=median(abs(Err$KDE1[Err$Periode==t]),na.rm=T)
# 3- Modèle sphérique :
modele3=fit.variogram(vario,vgm(max(vario[,3]),"Sph",140,0.1*max(vario[,3])),fit.method=1)
modele3=if (modele3[1,2]<0 || modele3[2,2]<0 || modele3[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Sph",140),fit.method=1) else modele3
if (modele3[dim(modele3)[1],2]>0 && modele3[dim(modele3)[1],3]>0) {
valid=krige.cv(Stations~1+GEM,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,])
Err$KDE3[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KDE3[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KDE3EQM[t]=mean((Err$KDE3[Err$Periode==t])**2,na.rm=T)
Info$KDE3EAM[t]=median(abs(Err$KDE3[Err$Periode==t]),na.rm=T)
# 4- Modèle exponentiel :
modele4=fit.variogram(vario,vgm(max(vario[,3]),"Exp",140,0.1*max(vario[,3])),fit.method=1)
modele4=if (modele4[1,2]<0 || modele4[2,2]<0 || modele4[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Exp",140),fit.method=1) else modele4
if (modele4[dim(modele4)[1],2]>0 && modele4[dim(modele4)[1],3]>0) {
valid=krige.cv(Stations~1+GEM,~x+y,model=modele4,data=Donnees[Donnees$Periode==t,])
Err$KDE4[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KDE4[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KDE4EQM[t]=mean((Err$KDE4[Err$Periode==t])**2,na.rm=T)
Info$KDE4EAM[t]=median(abs(Err$KDE4[Err$Periode==t]),na.rm=T)
# 5- Modèle gaussien :
modele5=fit.variogram(vario,vgm(max(vario[,3]),"Gau",140,0.1*max(vario[,3])),fit.method=1)
modele5=if (modele5[1,2]<0 || modele5[2,2]<0 || modele5[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Gau",140),fit.method=1) else modele5
if (modele5[dim(modele5)[1],2]>0 && modele5[dim(modele5)[1],3]>0) {
valid=krige.cv(Stations~1+GEM,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,])
Err$KDE5[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KDE5[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KDE5EQM[t]=mean((Err$KDE5[Err$Periode==t])**2,na.rm=T)
Info$KDE5EAM[t]=median(abs(Err$KDE5[Err$Periode==t]),na.rm=T)
} else Err$KDE6[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KDE6EQM[t]=mean((Err$KDE6[Err$Periode==t])**2,na.rm=T)
Info$KDE6EAM[t]=median(abs(Err$KDE6[Err$Periode==t]),na.rm=T)
# 7- Modèle puissance :
modele7=fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1,0.1*max(vario[,3])),fit.method=1)
modele7=if (modele7[1,2]<0 || modele7[2,2]<0 || modele7[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1),fit.method=1) else modele7
if(modele7[dim(modele7)[1],2]>0&&modele7[dim(modele7)[1],3]>0&&modele7[dim(modele7)[1],3]<2){
valid=krige.cv(Stations~1+GEM,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,])
Err$KDE7[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
} else Err$KDE7[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Info$KDE7EQM[t]=mean((Err$KDE7[Err$Periode==t])**2,na.rm=T)
Info$KDE7EAM[t]=median(abs(Err$KDE7[Err$Periode==t]),na.rm=T)
# Choix du modèle :
indice=data.frame(c(Info$KDE1EQM[t],Info$KDE2EQM[t],Info$KDE3EQM[t],Info$KDE4EQM[t],
Info$KDE5EQM[t],Info$KDE6EQM[t],Info$KDE7EQM[t]))
minimums=na.omit(ifelse(indice==min(indice,na.rm=T),row(indice),NA))
Info$KDEmodele[t]=minimums[1,1]
modele = if (Info$KDEmodele[t]==1) modele1 else if (Info$KDEmodele[t]==2) modele2 else
if (Info$KDEmodele[t]==3) modele3 else if (Info$KDEmodele[t]==4) modele4 else
if (Info$KDEmodele[t]==5) modele5 else if (Info$KDEmodele[t]==6) modele6 else modele7
# Deuxième itération :
# 1- Modèle pépitique :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele1=fit.variogram(vario,vgm(max(vario[,3]),"Nug"),fit.method=1)
Info$KDE1pepite[t]= modele1$psill[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele1[dim(modele1)[1],2]>0) {
valid=krige.cv(Stations~1+GEM,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,])
Err$KDE1[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+GEM,~x+y,model=modele1,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KDE1[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KDE1[Var$Periode==t]=Interpo[,4]
} else {
Err$KDE1[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KDE1[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KDE1[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KDE1EQM[t]=mean((Err$KDE1[Err$Periode==t])**2,na.rm=T)
Info$KDE1EAM[t]=median(abs(Err$KDE1[Err$Periode==t]),na.rm=T)
Annexe B. Programme S-Plus 124
# 3- Modèle sphérique :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele3=fit.variogram(vario,vgm(max(vario[,3]),"Sph",140,0.1*max(vario[,3])),fit.method=1)
modele3=if (modele3[1,2]<0 || modele3[2,2]<0 || modele3[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Sph",140),fit.method=1) else modele3
Info$KDE3pepite[t]=if(dim(modele3)[1]==2) modele3$psill[1] else 0
Info$KDE3ppalier[t]=if(dim(modele3)[1]==2) modele3$psill[2] else modele3$psill[1]
Info$KDE3portee[t]=if(dim(modele3)[1]==2) modele3$range[2] else modele3$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele3[dim(modele3)[1],2]>0 && modele3[dim(modele3)[1],3]>0) {
valid=krige.cv(Stations~1+GEM,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,])
Err$KDE3[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+GEM,~x+y,model=modele3,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KDE3[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KDE3[Var$Periode==t]=Interpo[,4]
} else {
Err$KDE3[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KDE3[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KDE3[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KDE3EQM[t]=mean((Err$KDE3[Err$Periode==t])**2,na.rm=T)
Info$KDE3EAM[t]=median(abs(Err$KDE3[Err$Periode==t]),na.rm=T)
# 4- Modèle exponentiel :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele4=fit.variogram(vario,vgm(max(vario[,3]),"Exp",140,0.1*max(vario[,3])),fit.method=1)
modele4=if (modele4[1,2]<0 || modele4[2,2]<0 || modele4[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Exp",140),fit.method=1) else modele4
Info$KDE4pepite[t]=if(dim(modele4)[1]==2) modele4$psill[1] else 0
Info$KDE4ppalier[t]=if(dim(modele4)[1]==2) modele4$psill[2] else modele4$psill[1]
Annexe B. Programme S-Plus 125
# 5- Modèle gaussien :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele5=fit.variogram(vario,vgm(max(vario[,3]),"Gau",140,0.1*max(vario[,3])),fit.method=1)
modele5=if (modele5[1,2]<0 || modele5[2,2]<0 || modele5[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3]),"Gau",140),fit.method=1) else modele5
Info$KDE5pepite[t]=if(dim(modele5)[1]==2) modele5$psill[1] else 0
Info$KDE5ppalier[t]=if(dim(modele5)[1]==2) modele5$psill[2] else modele5$psill[1]
Info$KDE5porteep[t]=if(dim(modele5)[1]==2) sqrt(3)*modele5$range[2] else sqrt(3)*modele5$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if (modele5[dim(modele5)[1],2]>0 && modele5[dim(modele5)[1],3]>0) {
valid=krige.cv(Stations~1+GEM,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,])
Err$KDE5[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+GEM,~x+y,model=modele5,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KDE5[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KDE5[Var$Periode==t]=Interpo[,4]
} else {
Err$KDE5[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KDE5[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KDE5[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KDE5EQM[t]=mean((Err$KDE5[Err$Periode==t])**2,na.rm=T)
Info$KDE5EAM[t]=median(abs(Err$KDE5[Err$Periode==t]),na.rm=T)
Var$KDE6[Var$Periode==t]=Interpo[,4]
} else {
Err$KDE6[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KDE6[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KDE6[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KDE6EQM[t]=mean((Err$KDE6[Err$Periode==t])**2,na.rm=T)
Info$KDE6EAM[t]=median(abs(Err$KDE6[Err$Periode==t]),na.rm=T)
# 7- Modèle puissance :
# Pour ajuster le modèle et noter les info par rapport à ce modèle :
modele7=fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1,0.1*max(vario[,3])),fit.method=1)
modele7=if (modele7[1,2]<0 || modele7[2,2]<0 || modele7[2,3]<0)
fit.variogram(vario,vgm(max(vario[,3])/150,"Pow",1),fit.method=1) else modele7
Info$KDE7pepite[t]=if(dim(modele7)[1]==2) modele7$psill[1] else 0
Info$KDE7echelle[t]=if(dim(modele7)[1]==2) modele7$psill[2] else modele7$psill[1]
Info$KDE7puissance[t]=if(dim(modele7)[1]==2) modele7$range[2] else modele7$range[1]
# Pour faire validation croisée, interpolation, sauver résultats et noter informations :
if(modele7[dim(modele7)[1],2]>0&&modele7[dim(modele7)[1],3]>0&&modele7[dim(modele7)[1],3]<2){
valid=krige.cv(Stations~1+GEM,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,])
Err$KDE7[Err$Periode==t]=valid[,5]-ifelse(valid[,3]<0,0,valid[,3])
Interpo=krige(Stations~1+GEM,~x+y,model=modele7,data=Donnees[Donnees$Periode==t,],
newdata=Grille[Grille$Periode==t,])
Prev$KDE7[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$KDE7[Var$Periode==t]=Interpo[,4]
} else {
Err$KDE7[Err$Periode==t]=rep(NA,dim(Donnees[Donnees$Periode==t,])[1])
Prev$KDE7[Prev$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
Var$KDE7[Var$Periode==t]=rep(NA,dim(Grille[Grille$Periode==t,])[1])
}
Info$KDE7EQM[t]=mean((Err$KDE7[Err$Periode==t])**2,na.rm=T)
Info$KDE7EAM[t]=median(abs(Err$KDE7[Err$Periode==t]),na.rm=T)
# Choix du modèle :
indice=data.frame(c(Info$KDE1EQM[t],Info$KDE2EQM[t],Info$KDE3EQM[t],Info$KDE4EQM[t],
Info$KDE5EQM[t],Info$KDE6EQM[t],Info$KDE7EQM[t]))
minimums=na.omit(ifelse(indice==min(indice,na.rm=T),row(indice),NA))
Info$KDEmodele[t]=minimums[1,1]
if (Info$KDEmodele[t]==1) {
Err$KDE[Err$Periode==t]=Err$KDE1[Err$Periode==t]
Prev$KDE[Prev$Periode==t]=Prev$KDE1[Prev$Periode==t]
Var$KDE[Var$Periode==t]=Var$KDE1[Var$Periode==t]
} else
if (Info$KDEmodele[t]==2) {
Err$KDE[Err$Periode==t]=Err$KDE2[Err$Periode==t]
Prev$KDE[Prev$Periode==t]=Prev$KDE2[Prev$Periode==t]
Var$KDE[Var$Periode==t]=Var$KDE2[Var$Periode==t]
} else
if (Info$KDEmodele[t]==3) {
Err$KDE[Err$Periode==t]=Err$KDE3[Err$Periode==t]
Prev$KDE[Prev$Periode==t]=Prev$KDE3[Prev$Periode==t]
Var$KDE[Var$Periode==t]=Var$KDE3[Var$Periode==t]
} else
if (Info$KDEmodele[t]==4) {
Err$KDE[Err$Periode==t]=Err$KDE4[Err$Periode==t]
Annexe B. Programme S-Plus 127
Prev$KDE[Prev$Periode==t]=Prev$KDE4[Prev$Periode==t]
Var$KDE[Var$Periode==t]=Var$KDE4[Var$Periode==t]
} else
if (Info$KDEmodele[t]==5) {
Err$KDE[Err$Periode==t]=Err$KDE5[Err$Periode==t]
Prev$KDE[Prev$Periode==t]=Prev$KDE5[Prev$Periode==t]
Var$KDE[Var$Periode==t]=Var$KDE5[Var$Periode==t]
} else
if (Info$KDEmodele[t]==6) {
Err$KDE[Err$Periode==t]=Err$KDE6[Err$Periode==t]
Prev$KDE[Prev$Periode==t]=Prev$KDE6[Prev$Periode==t]
Var$KDE[Var$Periode==t]=Var$KDE6[Var$Periode==t]
} else {
Err$KDE[Err$Periode==t]=Err$KDE7[Err$Periode==t]
Prev$KDE[Prev$Periode==t]=Prev$KDE7[Prev$Periode==t]
Var$KDE[Var$Periode==t]=Var$KDE7[Var$Periode==t]
}
Info$KDEEQM[t]=mean((Err$KDE[Err$Periode==t])**2,na.rm=T)
Info$KDEEAM[t]=median(abs(Err$KDE[Err$Periode==t]),na.rm=T)
}
# COKRIGEAGE ORDINAIRE
for (t in Pluie) {
# Validation croisée
Periodet=na.omit(ifelse(Donnees$Periode==t,row(Donnees),NA))
for (i in 1:length(Periodet)){
g.valid=gstat(id="Stations",formula=Stations~1,loc=~x+y,data=Donnees[Periodet[-i],],
model=g.interpo$model$Stations)
g.valid=gstat(g.valid,id="GEM",formula=GEM~1,locations=~x+y,data=Grille[Grille$Periode==t,],
model=g.interpo$model$GEM)
g.valid=gstat(g.valid,id=c("Stations","GEM"),model=g.interpo$model$Stations.GEM)
valid=predict(g.valid,Donnees[Periodet[i],])
Err$CKO[Periodet[i]]=Donnees$Stations[Periodet[i]]
-ifelse(valid$Stations.pred<0,0,valid$Stations.pred)
}
Info$CKOEQM[t]=mean((Err$CKO[Err$Periode==t])**2,na.rm=T)
Info$CKOEAM[t]=median(abs(Err$CKO[Err$Periode==t]),na.rm=T)
Prev$CKO[Prev$Periode==t]=ifelse(Interpo[,3]<0,0,Interpo[,3])
Var$CKO[Var$Periode==t]=Interpo[,4]
}
for( t in Pluie) {
indice=data.frame(c(Info$KOEQM[t],Info$KUEQM[t],Info$KDEEQM[t],
Info$CKOEQM[t],Info$RLxyEQM[t],Info$RLxywEQM[t]))
minimums=na.omit(ifelse(indice==min(indice,na.rm=T),row(indice),NA))
if (minimums[1,1]==1) {
Info$Selection[t]="KO"
Err$Selection[Err$Periode==t]=Err$KO[Err$Periode==t]
Prev$Selection[Prev$Periode==t]=Prev$KO[Prev$Periode==t]
Var$Selection[Var$Periode==t]=Var$KO[Var$Periode==t]
} else
if (minimums[1,1]==2) {
Info$Selection[t]="KU"
Err$Selection[Err$Periode==t]=Err$KU[Err$Periode==t]
Prev$Selection[Prev$Periode==t]=Prev$KU[Prev$Periode==t]
Var$Selection[Var$Periode==t]=Var$KU[Var$Periode==t]
} else
if (minimums[1,1]==3) {
Info$Selection[t]="KDE"
Err$Selection[Err$Periode==t]=Err$KDE[Err$Periode==t]
Prev$Selection[Prev$Periode==t]=Prev$KDE[Prev$Periode==t]
Var$Selection[Var$Periode==t]=Var$KDE[Var$Periode==t]
} else
if (minimums[1,1]==4) {
Info$Selection[t]="CKO"
Err$Selection[Err$Periode==t]=Err$CKO[Err$Periode==t]
Prev$Selection[Prev$Periode==t]=Prev$CKO[Prev$Periode==t]
Var$Selection[Var$Periode==t]=Var$CKO[Var$Periode==t]
} else
if (minimums[1,1]==5) {
Info$Selection[t]="RLxy"
Err$Selection[Err$Periode==t]=Err$RLxy[Err$Periode==t]
Prev$Selection[Prev$Periode==t]=Prev$RLxy[Prev$Periode==t]
Var$Selection[Var$Periode==t]=Var$RLxy[Var$Periode==t]
} else {
Info$Selection[t]="RLxyw"
Err$Selection[Err$Periode==t]=Err$RLxyw[Err$Periode==t]
Prev$Selection[Prev$Periode==t]=Prev$RLxyw[Prev$Periode==t]
Var$Selection[Var$Periode==t]=Var$RLxyw[Var$Periode==t]
}
Info$SelectionEQM[t]=mean((Err$Selection[Err$Periode==t])**2,na.rm=T)
Info$SelectionEAM[t]=median(abs(Err$Selection[Err$Periode==t]),na.rm=T)
}