Memoire Online - Méthodes géostatistique pour l'interpolation et la modélisation en 2d/3d des données spatiales

C'est un plaisir pour moi de remercier l'ensemble des chercheurs du groupe LEMEL et plus particulièrement Evelyne GOUBERT, David MENIER et Valérie MONBET. Ces derniers m'ont accordé beaucoup de leur temps, m'ont fait bénéficier de leur savoir et de leur expérience tout au long de ce travail. Ils ont été d'une aide précieuse pour comprendre la problématique, la cerner et y répondre.

Résumé

Les travaux de recherche présentés dans ce rapport ont été effectués au sein du laboratoire LEMEL (Laboratoire d'Etude et Modélisation des Environnements Littoraux). Il parcourt différentes techniques d'interpolations spatiales. A partir d'un échantillon de données sédimentologiques, géophysiques, et bathymétriques, il s'agit de réaliser des cartes 2D/3D indiquant la profondeur au toit du socle rocheux à l'est de l'île aux Moines dans le Golfe du Morbihan. Les méthodes d'interpolations sont également utilisées pour localiser des zones à forte présence en crépidules au nord ouest de l'île de Bailleron. Ces deux applications concrètes permettent d'opposer les méthodes déterministes aux méthodes géostatistiques. Les fondements mathématiques des techniques d'interpolations spatiales sont analysés. Des préconisations générales sont données, sur l'usage des modèles d'interpolation pour la cartographie.

Mots-clés : coordonnées géographiques, interpolation, extrapolation, lissage, estimation, sectorisation, variable régionalisé, géostatistique, variogramme, krigeage, validation croisée, toit de substratum.

The research tasks presented in this memory were carried out within laboratory LEMEL (Laboratory for Study and Modelling of the Littoral Environments). It goes through various techniques of spatial interpolations. From a sample of `sedimentologics', geophysics, and bathymetric data, it involves of realizing 2D / 3D maps modelling the depth to the roof of the rocky in the Morbihan's bay. The methods of interpolations are also used to localize zones with strong presence of `crépidules' in the western North of Bailleron's island. These two concrete applications will oppose the determinist methods against the géostatistics methods. The mathematical bases as well as their demonstrations are examined. General recommendations are giving, on the usage of the models of interpolation for the mapping.

Keywords: geographic coordinated, interpolation, extrapolation, smoothing, estimation, division into sectors, regionalized variable, geostatistic, variogram, model fitting, kriging, model validation, roof of substratum.

Introduction : présentation de l'étude

Les travaux de recherche présentés dans ce mémoire ont été effectués au sein du laboratoire LEMEL (Laboratoire d'Etude et Modélisation des Environnements Littoraux). Il est rattaché à la « Faculté de Science et Science de l'Ingénieur » de l'Université de Bretagne Sud. C'est un groupe pluridisciplinaire (géologues, biologistes, mathématiciens, statisticiens...) qui s'attache à modéliser les écosystèmes à l'interface terre/mer. Il passe pour cela par l'intermédiaire de SIG (Système d'Information Géographique) afin d'acquérir, de consulter et de gérer l'information géographique.

Ce mémoire de stage présente l'élaboration d'une carte de la profondeur au toit du socle rocheux à l'ouest de l'île aux Moines et une carte localisant les crépidules au nord ouest de l'île de Bailleron (golfe du Morbihan). Une image 3D au du toit du substratum fournit un outil pour comprendre comment se répartissent les corps sédimentaires dans l'espace. Localiser les zones à présence de crédidules permet de lutter contre cet envahisseur. La connaissance, des profondeurs au toit du substratum^* ou la quantité de crépidules^**, en tout point d'une zone, nécessite :

et sédimentologique que l'on rapporte à des coordonnées géographiques précises (latitude / longitude)

données de mesure, les profondeurs en tout point de la zone : il s'agit alors d'une estimation mathématique de cette profondeur.

Elles ont toutes leurs avantages et inconvénients. Mais dans tous les cas, le principe reste le même : une valeur est estimée pour des sites non échantillonnés, sur la base d'observations existantes.

Le premier chapitre du rapport recense différentes approches déterministes pour interpoler un ensemble de données. Une description très générale en est fournit.

Le deuxième chapitre, propose une analyse des techniques géostatistiques et des conditions de leur mise en oeuvre. Les fondements mathématiques ainsi que leurs démonstrations sont examinées. Le krigeage et l'analyse variographique, une étape préalable du krigeage, sont analysés et commentés. Elles s'appuient sur une synthèse bibliographique.

Les chapitres III à VI présentent une application des méthodes géostatistiques pour aboutir un modèle numérique de terrain au toit du substratum. Ils fournissent également des conseils sur la nature de l'information qu'il est nécessaire de recueillir pour utiliser efficacement ces méthodes. Ils préconisent une démarche à suivre pour interpréter les résultats. Ils mettent en évidence les avantages et les limites des méthodes géostatistiques.

Le dernier chapitre dévoilera une application des méthodes d'analyse spatiale concernant
uniquement la position géographique des objets, et des méthodes non stochastiques présentées

* Substratum : base (souvent un socle rocheux ou cristallin) sur laquelle repose les sédiments. ** Crépidule : mollusque gastéropode

en première partie. Elles ont pour objectif d'identifier des zones géographiques occupées par les crépidules.

L'interpolation spatiale est une procédure qui consiste à estimer la valeur d'un attribut pour des sites non échantillonnés. Dans le contexte de ce stage, nous utilisons l'interpolation spatiale pour mettre en oeuvre des algorithmes mathématiques ou probabilistes afin d'estimer la profondeur de la roche entre les points d'échantillonnages.

Il existe de nombreuses méthodes d'interpolation parmi lesquels il faut faire un choix. Nous distinguerons deux familles :

- les méthodes d'interpolation classiques basées sur des algorithmes purement déterministes.

- Les méthodes d'estimation géostatistiques qui s'appuient sur une modélisation probabiliste du phénomène étudié.

Les méthodes s'appliquent à des variables régionalisées, c'est-à-dire des fonctions numériques qui prennent leurs valeurs dans des régions délimitées de l'espace appelées champ. Dans notre cas, la variable régionalisée est la profondeur au toit du substratum. Le champ est la zone d'étude décrite dans le chapitre III.

1.1/ Forme générale de l'interpolation linéaire

Pour estimer la valeur ponctuelle en un site, nous utilisons des combinaisons linéaires pondérées. La distance entre les lieux d'observations a une influence sur ce que l'on observe. Autrement dit, les valeurs dans deux localités voisines sont souvent plus semblables que dans deux localités éloignées. Les poids doivent donc tenir compte non seulement de la disposition des observations les unes par rapport aux autres, mais aussi de la distance entre le site à estimer et les sites observés.

1.1.1/ Méthodes d'interpolation par partitionnement de l'espace 1.1.1.1/ Polygones de Thiessen

La méthode consiste à partitionner l'espace géographique en polygones, puis à attribuer une valeur à chacun des polygones. Elle permet de déterminer un zonage où la valeur de la variable à prédire est à priori la même que celle du site d'observation^*.

Supposons que l'on désire estimer la valeur en un point S0 du champ (cf. figure 1.1). Ce point appartient nécessairement à l'un des polygones d'influence des sites d'observations.

Un polygone d'influence pour l'observation C s'obtient en traçant les médiatrices des segments joignant C aux sites voisins.

On attribue à S0 et à tout point appartenant au même polygone, la valeur du site ayant ce polygone d'influence. Sur la figure 1.1, la valeur estimée au point S0 sera identique à celle du site C.

La méthode de Thiessen peut s'améliorer tout en gardant le même principe. L'idée de Sibson est de construire un polygone de Thiessen autour du site à estimer. Dans un deuxième temps calculer les surfaces d'intersection (P(S0,i)) entre le polygone précèdent et ceux des sites voisins. Enfin, Sibson propose une estimation de S0 à l'aide d'une combinaison linéaire des valeurs des sites voisins (Z(S0)) pondéré par les surfaces P(S0,i).

P(S0,i) : surface de l'intersection du polygone de l'observation i et celle de S0

La méthode de Sibson a les propriétés suivantes : - l'estimation en S0 est unique

La méthode d'interpolation par triangulation consiste à diviser le champ en triangles disjoints, dont les sommets sont les sites échantillonnés, puis à interpoler à l'intérieur de chaque triangle.

La construction des triangles n'est pas unique, différentes approches sont donc proposées.

On peut à partir du diagramme de Thiessen, en construire le dual. C'est à dire construire un nouveau diagramme où cette fois, on relie par un segment toutes les paires de sites dont les régions de Thiessen correspondantes sont adjacentes (les points séparés par une arête de Thiessen).

On observe que la triangulation de Delaunay opère seulement dans l'enveloppe convexe des sites d'observations et ne recouvre donc pas entièrement le champ.

Supposons que le point S à estimer se trouve à l'intérieur du triangle formé par les sites S1,S2, S3 (fi_gure 1.3). L'estimation de la valeur de la variable régionalisée au point S par interpolation linéaire s'écrit :

L'interpolation linéaire a les propriétés suivantes : - l'estimation en S0 est unique

- la triangulation ne permet pas d'extrapoler les valeurs au-delà de l'enveloppe convexe des sites d'observation.

Les méthodes d'interpolations précédentes ne considèrent, pour estimer la valeur d'un site S0, que les sites d'observation immédiatement voisins. Elles ignorent par conséquent une grande partie de l'information disponible. Les méthodes barycentriques permettent de prendre en compte un nombre plus important de données.

1.1.2.1/ Plus proches voisins

Le plan sur lequel se situe le point à interpoler est découpé en octants. La méthode de calcul consiste en une interpolation sur les huit points de référence les plus proches du point à interpoler, répartis dans les octants. L'importance d'un point de référence est d'autant plus grande dans le résultat du calcul que sa distance au point à interpoler est faible.

L'algorithme est efficace lorsque l'on dispose de beaucoup de points bien répartis dans la fenêtre de calcul.

On peut imaginer découper l'espace en plus ou moins huit parties et prendre plus qu'un point observé par partie. C'est en cela que l'on se rapproche de la méthode qui suit.

1.1.2.2/ Méthode de l'inverse des distances

La première étape est d'effectuer une recherche des sites qui vont intervenir dans l'estimation. On peut par exemple se fixer un rayon de recherche dont le centre est la localisation de la valeur à estimer. On ne retiendra que les sites Si appartenant au cercle.

Dans un deuxième temps on attribue à chaque site Si retenu un poids inversement proportionnel à la distance entre ce site et le point à estimer S0. Considérons que l'on a retenu n0 données dans le cercle alors on obtient comme estimateur :

Les deux méthodes précédentes ont l'inconvénient de prendre en compte uniquement la distance qui sépare les sites entre eux. Cela donne un poids important au groupement de données, alors qu'il n'est pas nécessairement justifié.

1.1.3/ Résumé

Les méthodes ci-dessus ont la caractéristique de traiter uniquement les données de la variable étudiée. Toutes définissent la valeur recherchée en un point comme une combinaison linéaire pondérée des mesures disponibles.

Ce sont des méthodes implémentées dans la plus part des logiciels de Système d'Information Géographique (SIG).

Ces techniques présentent néanmoins des défauts. Elles ignorent la structure spatiale de la variable et produisent du coup des surfaces interpolées très lisses. Des situations locales très spécifiques peuvent alors être omises (zones de fortes ou de très faibles valeurs). Nous prenons le risque d'aboutir à des cartes peu réalistes. Enfin, aucun critère statistique pour juger de la précision de ces cartes n'est formulé.

Si l'on veut optimiser la précision des estimations, il faudra utiliser d'autres outils qui feront appel à des modèles probabilistes. La géostatistique et le krigeage en font partie. Nous allons les étudier dans la suite.

1.2/ Modèles d'estimation géostatistique

Le mot géostatistique est un néologisme forgé à l'École des Mines. La géostatistique est née des problèmes rencontrés dans le secteur de la mine : contrôle des teneurs, optimisation de maille, cartographie des ressources, prévision des réserves récupérables, étude de scénarios d'exploitation...

Daniel Kriege, géologue dans les mines d'or, proposa dans les années 60 une méthode statistique pour estimer la teneur d'un bloc de minerai à partir d'échantillons pris autour du bloc à exploiter. Dix ans plus tard, Georges Matheron développa un outil pour analyser la continuité spatiale des teneurs appelé le « variogramme » et une méthode d'estimation basée sur le « variogramme » appelée le « krigeage ». Nous étudierons ces deux méthodes.

Aujourd'hui, la géostatistique s'exprime dans des champs d'applications comme l'océanographie, la météorologie, le génie civil, l'environnement, la géologie, la qualité de l'air et des sols, la santé, etc.

Techniquement, la géostatistique utilise également une combinaison linéaire des données
observées, mais à la différence des méthodes classiques d'interpolation, elle tient compte à la
fois de l'information relative à leur position et du caractère aléatoire du phénomène étudié. De

plus, elle permet d'intégrer des informations auxiliaires dans l'estimation. Ces avantages font considérablement améliorer les estimations dans le contexte spatial.

La mise en oeuvre de la méthode de krigeage passe par une étape d'analyse des données. Elle est destinée à décrire la structure spatiale de la variable régionalisée. La carte de l'interpolation est alors accompagnée d'indicateurs de fiabilité des résultats.

Chapitre II : Méthodes de la géostatistique linéaire

La géostatistique se réfère aux méthodes d'analyse probabiliste pour étudier des phénomènes corrélés dans l'espace appelés phénomènes régionalisés. A ce titre elle fournit différents outils pour répondre au problème posé par la cartographie du socle rocheux dans le golfe du Morbihan.

1/ Notations et définitions

Une variable régionalisée quantifie des grandeurs mesurées sur l'espace géographique. L'espace dans lequel cette variable prend ses valeurs est appelé champ.

Exemple de variable régionalisée : la profondeur du substratum sous le niveau zéro de la mer, mesurée par des campagnes sismiques dans une zone géographique située à l'est de l'île aux Moines, dans le golfe du Morbihan.

Exemple de champ : la zone géographique située entre la côte et l'est de l'île aux Moines. Nous pourrons estimer les valeurs prises par la variable régionalisée dans cette zone.

1.2/ Notations

D : le champ (domaine sur lequel la variable régionalisée est définie) s ? D : une position dans le champ

Chaque profondeur de roche est une variable aléatoire. Ensemble elles forment une fonction aléatoire de s.

La mesure faite au point si est une réalisation particulière de la fonction aléatoire Z(si).

Définissons à présent les hypothèses indispensables pour utiliser les techniques de la géostatistique linéaire.

2/ Hypothèses de base

Une fonction aléatoire {Z(s), s ? D} est caractérisé par sa loi spatiale F. Elle correspond à la loi de probabilité conjointe de (Z(s1), Z(s2), _Z(s3), ..., Z(s_n)).

Or cette fonction est très complexe par l'infinité des combinaisons possibles. Nous n'allons donc pas pouvoir estimer la fonction de distribution conjointe. La géostatistique linéaire se limite à la fonction de distribution d'ordre un F(v) et d'ordre deux FZ(si),Z(sj)(vi,vj) .

La première nous permet de calculer l'espérance de la fonction aléatoire Z en un point s.

La fonction de distribution d'ordre deux, fourni la loi de probabilité entre les valeurs prises en deux sites si et sj. On utilise la covariance pour quantifier le degré de ressemblance entre les valeurs prises en si et sj et le variogramme pour mesurer la dissemblance entre les valeurs prises aux sites si et sj.

cov( s _i ,s ) = E(Z ( s _i ) Z ( s _j )) - E(Z ( s _i )) E ( Z (s j ))

Le problème qui se pose est que la variable régionalisée est observée qu'une seule fois à un endroit précis. En d'autres termes, nous n'avons qu'une seule réalisation de la variable aléatoire. Or pour estimer les moments d'ordre un et deux il nous faudrait un grand nombre (>30) d'observations. Ce problème ne concerne pas la quantité d'information disponible mais le fait que l'on essaie de décrire un phénomène unique (profondeur au toit du socle rocheux), qui ne se répète pas, à l'aide de lois de probabilités.

Pour palier à ce problème nous posons comme hypothèse que la variable régionalisée est stationnaire. Concrètement cela veut dire que deux paires de points espacés d'un même vecteur h ont des caractéristiques (moyenne et covariance) semblables. Ou encore, la variable régionalisée ne dépend pas de sa position dans l'espace, elle garde les mêmes caractéristiques

où que l'on se place. Cela nous permet de nous détacher de la localisation et de nous restreindre uniquement à la distance qui sépare les points d'observations.

Une fonction aléatoire Z(s) est stationnaire du 2^nd ordre quand l'espérance mathématique existe et ne dépend pas du point s ; et que la covariance entre chaque paire ( Z(s+h) , Z(s)) existe et ne dépend que de h (distance).

- L'espérance mathématique ne dépend pas de s : ? s , E(Z ( s )) = m constante indépendante de s

- La variance existe en tout site s et est une constante indépendante du site s :

Remarque : L'hypothèse de stationnarité d'ordre deux ne peut être validé de manière rigoureuse et infaillible à l'aide d'un test statistique sur les données expérimentales (Arnaud et Emery, 2000 p.107).

On dit _qu'une fonction aléatoire Z(s) est intrinsè_que _quand ses accroissements Z(s+h)-Z(s) sont stationnaires d'ordre 2. C'est-à-dire que

Cette hypothèse permet de dire que la variabilité entre les valeurs prises en deux points différents ne dépend que de h (la distance entre ces points).

Toute fonction aléatoire stationnaire d'ordre deux est également intrinsèque (la réciproque est
fausse). Autrement dit, l'hypothèse de stationnarité intrinsèque est moins restrictive que la

stationnarité du second ordre. L'hypothèse intrinsèque ne requiert pas de connaître l'espérance ni sa covariance de la variable aléatoire.

La fonction la plus utilisée en géostatistique pour décrire la continuité spatiale est le variogramme. La continuité spatiale est réalisée lorsque les valeurs prises entre deux sites proches l'un de l'autre sont similaires.

3/ Variogramme

3.1/ Variogramme théorique et variogramme expérimental

D'après Marcotte (cours de «géostatistique minière »), la nature n'est pas entièrement imprévisible. Deux observations situées l'une près de l'autre devraient en moyenne se ressembler davantage que deux observations éloignées. La différence entre les valeurs prises par deux variables aléatoires est Z(s)-Z(s+h). C'est également une variable aléatoire dont on peut calculer la variance. Cette variance devrait être plus petite lorsque les points sont rapprochés (les valeurs se ressemblent plus en moyenne) et plus grande lorsque les points sont éloignés. On appelle variogramme la demi-variance de cette différence.

L'outil mesure la « variabilité spatiale », c'est-à-dire la dissemblance entre les valeurs en fonction de leurs séparations. Il décrit la continuité spatiale de la variable régionalisée.

- La covariance mesure la ressemblance entre les valeurs en fonction de leur éloignement alors que le variogramme mesure la dissemblance entre les valeurs en fonction de leur éloignement.

- Dans l'hypothèse de stationnarité d'ordre 2, covariance et variogramme existent et sont liés par la relation ã ( h ) = C(0) - C(h ) . Dans l'hypothèse intrinsèque, seul le variogramme existe. C'est pourquoi il est généralement préféré à la covariance pour décrire et interprété la structure spatiale du phénomène étudié.

- Le variogramme réel d'une fonction aléatoire est généralement inconnu, mais il peut être évalué à partir des données d'échantillonnages. On obtient ainsi le variogramme expérimental proposé par Matheron (1962).

N(h) est le nombre de paires dans la classe de distance h. Z(si) est la profondeur de la roche au point de mesure si.

3.2/ Propriétés du variogramme

Le variogramme est une fonction de h, croissante et souvent caractérisé par quatre paramètres :

L'effet de pépite : Le comportement à l'origine du variogramme reflète le degré de régularité spatiale de la variable régionalisée. Si le variogramme présente un saut abrupt à l'origine (effet de pépite), cela indique une absence partielle de corrélation entre les valeurs prises en deux sites très proches. C'est-à-dire qu'il y a une faible ressemblance entre les valeurs régionalisées très voisines.

La portée : Distance où deux observations ne se ressemblent plus du tout. Leur covariance est nulle.

- Si le variogramme est borné alors la covariance existe et l'on peut présumer une stationnarité du second ordre (Journel et Huijbergts, 1993 p.37).

- Si la variable régionalisée est stationnaire du second ordre alors le palier est égal à la variance de cette même variable.

- Inversement, si un variogramme est non borné, il ne possède ni portée, ni palier. La variance de la fonction n'est pas définie et elle n'est donc pas stationnaire du second ordre.

- Arnaud et Emery (2000, p.126) affirment que le variogramme expérimental n'est pas fiable pour des distances supérieures à la moitié du diamètre du champ D. Elles sont peu nombreuses et augmentent la variance de la valeur ãà (.) estimée.

- Les variogrammes directionnels permettent de détecter une éventuelle anisotropie. (cf. § 3.4)

3.3/ Modélisation du variogramme

Le variogramme expérimental n'est pas défini partout, notamment aux distances h pour lesquelles il n'existe pas de paire de points de mesures. Ainsi lui est-il ajusté une fonction mathématique appelée modèle de variogramme. Marcotte recommande d'utiliser des modèles éprouvés ou des modèles construit à partir de modèles éprouvés.

Le variogramme ne dépend que de h, c'est-à-dire le vecteur de déplacement entre les points s et s+h. Ce vecteur contient de l'information sur la distance entre ces deux points, par l'intermédiaire de sa norme, ainsi que sur l'orientation de h. Si le variogramme ne dépend en fait que de la norme de h, il est dit isotrope. S'il dépend aussi de la direction (è) du vecteur de translation, il est dit anisotrope.

Bien qu'il existe une très grande variété d'anisotropie, la plupart des ouvrages de géostatistique montrent uniquement comment modéliser les anisotropies géométriques.

- Les variogrammes des différentes directions ont le même palier et même effet pépite mais des portées différentes.

- Les portées maximales et minimales s'observent selon deux directions orthogonales.

Pour revenir à une situation isotrope, le principe consiste à effectuer une transformation linéaire des coordonnées spatiales c'est-à-dire une rotation en suivant les directions de plus petite et plus grande continuité (application chap. IV § 1.3).

3.5/ Cas multivarié : calcul des variogrammes croisés

Nous avons vu qu'un des avantages des méthodes géostatistique (chap. I § 1.2) était de pouvoir inclure dans le modèle d'interpolation des variables auxiliaires qui apportent de l'information sur la variable cible. Pour décrire la structure de dépendance entre la variable cible et les covariables il convient de calculer le variogramme croisé.

L'analyse du variogramme croisé se fait de la même façon que celle du variogramme simple. C'est-à-dire que l'on relève ses propriétés (effet de pépite, palier, porté, amplitude) et on y ajuste une fonction.

4/ Krigeage

Cette section expose l'une des techniques de géostatistique d'estimation locale, connue sous le nom de krigeage et cokrigeage ordinaire.

Nous cherchons à estimer la valeur d'une variable régionalisée z (profondeur au substratum) en un point s0 quelconque du champ à partir des mesures observées z(si), i=1,..,n (n : nombre de points observés). Le krigeage est un interpolateur exact (la valeur estimée sur un point de mesure est égale à la valeur du point de mesure) et optimal (il minimise la variance sur l'erreur d'estimation).

Il existe trois types de krigeage : le krigeage simple, le krigeage ordinaire et le krigeage universel. Le krigeage ordinaire est le plus fréquemment utilisé en pratique car les hypothèses de départ sont moins contraignantes que celle du krigeage simple. Seul le krigeage ordinaire sera développé ici car il répond aux besoins de notre problématique.

4.1/ Hypothèses et contraintes

Le krigeage ordinaire ne requiert pas la connaissance de l'espérance de la variable régionalisée. Autrement dit, on se place dans l'hypothèse d'une stationnarité intrinsèque. Il n'en reste pas moins que la moyenne doit rester constante à l'échelle du voisinage de krigeage.

Définition du voisinage de krigeage : Domaine du champ qui contient le site à estimer et les données utilisées dans l'estimation.

Comme pour les méthodes barycentriques (chap. I § 1.1.3), le voisinage de krigeage se résume à un cercle ou à une ellipse autour du point à estimer.

avec s D (le champ), _ì(s)=E(Z(s)) quasi-constante inconnue et ä(s) fonction aléatoire stationnaire intrinsèque d'espérance nulle (E(ä(s))=0) et de structure de dépendance connue (ã(|si-sj|) connu). Dans son mémoire (2005, p.31), S. Baillargeon explique que le caractère quasi-constant de _ì(s) 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.

En outre, l'estimateur du krigeage ordinaire doit vérifier les propriétés suivantes :

La linéarité : Z ( s₀ ) est une combinaison linéaire des données Z(si) i=1,..,n.

Le non-biais : Z^à ( s₀ ) est sans biais : E[ Z à ( s ₀ ) - Z(s ₀) ] = 0

Or _ì étant inconnue, le seul moyen de garantir le non-biais de l'erreur est de poser :

Optimalité : c'est-à-dire que Var ( Z^à ( s₀ ) - Z(s ₀)) est minimale.

Var Z(s ₀ ) ) = (2( s ₀ ) - Z(s₀ )) 1 = E( 2( s₀ ) ²) - 2E[ 2( s ₀ ) Z ( s ₀ )] + E( Z (s0) ²⁾

E( Z ( s ₀ ) ²⁾ =E0),Z(s_i) ? E w_iw_iE(Z (s_i)Z (s_j)) =? ? ùiùj cov(Z(s _i ), Z ( s _j))+

car dans le cas d'une variable aléatoire stationnaire, E(Z(si))=E(Z(sj))=_ì, d'où E(Z(si))E(Z(sj))=_ì2

2( s ₀ ) Z ( s ₀ ) ] = i Z(s _i ) Z ( s₀ ) -?1 = _iE( Z (s )Z ( s ₀ ) ) = i cov( Z(s _i), Z ( s ₀)) + ì2

Var ( 2( s ₀ ) - Z(s ₀ )) = Var( Z(s₀ ) ) - ù _i cov(Z(s _i ), Z ( s ₀ )) + ? ? ùiù j cov(Z(s ), Z ( s ))

Le but est de trouver le vecteur ù qui minimise cette variance sous la contrainte ?=ù =

Pour y parvenir, on peut utiliser la méthode de Lagrange. La fonction de Lagrange :

Cette fonction est appelée variance de krigeage. C'est la valeur minimale de la variance de l'erreur de prévision.

Le krigeage est un interpolateur linéaire sans biais. Il prend en compte la géométrie des données, les caractéristiques de la régionalisation et de la variance. Son but est de lisser les données.

Ayant calculé les poids ùi par la résolution du système précédent, nous pouvons écrire l'estimateur de la manière suivante :

Il est possible d'améliorer les estimations obtenues par krigeage en ajoutant à la variable à estimer, l'information fournie par d'autres variables (variables secondaires). La démarche à suivre est la suivante.

Pour estimer la variable Z1 au point s0, on utilise une combinaison linéaire pondérée des mesures concernant toutes les autres variables (Z2, ..., Zp).

- La variance d'estimation (ou de krigeage) obtenue par cokrigeage est toujours inférieure à celle obtenue par krigeage.

- En pratique, on observe que le gain est moindre lorsque les variables sont peu corrélées entre elles (Arnaud et Emery, 2000, p.195).

- L'estimation ponctuelle par (co)krigeage en un site d'observation est égale à la

4.4/ Conclusion

La méthode de krigeage est précédée par l'estimation d'une fonction variographique. C'est cette fonction qui va tenir compte à la fois de la géométrie des données, des caractéristiques de la régionalisation et de la précision de l'estimation. Il est donc important de souligner que la qualité de l'estimation et l'appréciation de sa précision reposent uniquement sur le modèle variographique utilisé. Il doit par conséquent, être le plus cohérent possible avec ce qui a été observé.

Après avoir décrit les fondements mathématiques des méthodes d'interpolation spatiale, nous nous efforçons, dans les chapitres suivants, à les appliquer sur des cas concrets.

Chapitre III : Recueil et analyse préalable de l'information

Les données observées constituent la matière première des méthodes d'interpolation. En géostatistique elles ont la particularité d'être géo-référencées. L'objet de ce chapitre et de présenter la zone du golfe du Morbihan retenue pour appliquer les méthodes d'interpolation géostatistique, puis de présenter les données utilisées pour aboutir à une cartographie du socle cristallin.

1/ Localisation de la zone d'étude

La zone d'étude se localise entre l'île aux Moines et Baden dans le golfe du morbihan. Le passage étroit entre la pointe de l'île aux Moines et le continent provoque des courants importants (figure 6.2). Ils sont à l'origine de formations sédimentaires hétérogènes. La charge de sédiments et la force du courant génèrent de l'érosion et modèlent le fond marin. Les courants assurent la mise en mouvement des particules sédimentaires, ce qui entraînent l'érosion et l'entretient des chenaux. La proximité de la zone d'étude à l'atlantique fait que les conditions hydrodynamiques y sont importantes par rapport aux secteurs à l'est de l'île aux Moines. Elles expliquent les profondeurs bathymétriques de l'ordre de 10 à 15 mètres et des fonds marins rocheux ou sableux.

2/Source des données

Le golfe du Morbihan a fait l'objet de campagnes bathymétriques et sismiques. L'évaluation de la profondeur de la roche dans le milieu marin, s'appuie sur ces relevés. Une partie du travail consiste à relever l'information apportée par les données sismiques (support papier), et de l'intégrer dans un Système d'Information Géographique. Les données bathymétriques sont quant à elles ajoutées à la base d'information. Le tout constitue un échantillon de travail.

2.1/ Informations provenant des campagnes sismiques

La sismique est une technique de mesure qui consiste à enregistrer en surface des échos issus de la propagation dans le sous-sol d'une onde sismique (figure 3.2). Ainsi les hétérogénéités du sous-sol sont relevées. Le passage par exemple d'une couche rocheuse à une couche sédimentaire va se traduire par la présence d'un réflecteur sur les enregistrements (figure 3.3).

Le temps d'arrivée de l'écho permet d'enregistrer la profondeur des obstacles rencontrés par l'onde. En résumé, les études sismiques fournissent une image (profil) de la structure du soussol. On peut y reconnaître la superposition de certaines couches de roche et/ou de sédiment (figure 3.3).

La profondeur du fond marin et du substratum est relevée (au décimètre) tous les 70 mètres (figure 3.4). Le profil fourni les profondeurs en « seconde temps double ». C'est-à-dire le temps que l'onde envoyée par l'appareil met pour revenir à ce dernier. Cette unité peut être convertie en mètre. La vitesse de propagation de l'onde est de 1500 m/s dans l'eau et de 1700 m/s dans les sédiments.

Remarque : les données récoltées peuvent manquer de précision. Les problèmes suivants sont rencontrés :

- Lors de la saisie sur papier avec un décimètre, une erreur de un à deux millimètres est courante (1mm saisi correspond à 18cm de profondeur réelle).

- Le niveau zéro autrement dit le niveau de la mer n'est pas connu pour les profils en notre possession. Il faudra admettre qu'il est approximativement égal au niveau de référence utilisé par le SHOM (Service Hydrographique et Océanographique de la Marine).

Une campagne sismique permet d'obtenir des informations sur la nature et la géométrie de la couverture sédimentaire. Pour rendre une collecte de données efficace, le rapport INERIS (2003, p.52) conseille de disposer les échantillons selon un maillage régulier. A nombre d'échantillons égal, la régularité de la maille facilite l'inférence de la structure spatiale et garantie une meilleure précision. Toujours d'après le même rapport, il est prouvé que la variance de l'erreur d'estimation d'un échantillonnage aléatoire stratifié est inférieure à celle d'un échantillonnage aléatoire uniforme. (On peut retrouver ces affirmations dans l'ouvrage de M. Arnaud et X. Emery, p113.)

En pratique, il n'est pas toujours possible de relever l'information suivant un maillage régulier comme le préconise la théorie. En effet, les nombreux bancs de sable du golfe ne permettent pas aux bateaux de naviguer partout.

Dans tous les cas, il est préconisé de resserrer la maille d'échantillonnage en quelques endroits, afin de faciliter la modélisation du variogramme aux petites distances. C'est pourquoi nous ajoutons aux données sismiques :

- des points de la base bathymétrique (chap. III § 2.2) mesurés sur des zones de roche (figure 3.5).

- des points d'experts, estimés visuellement par des personnes familiarisées avec le terrain (figure 3.5).

Le nombre de points a également son importance. Un nombre de points trop faible rend le variogramme « instable » et les résultats de l'estimation peu fiables. Selon Cressie (1993, p.237) 30 paires de points par classe de distance seraient nécessaires à la construction du variogramme. Dans le cas étudié nous disposons en moyenne de plus de 100 paires de points par classes de distances (cf. chap.IV, § 1).

2.2/ Bathymétrie : Information auxiliaire

Un modèle numérique de terrain (MNT) est une représentation de la topographie (altimétrie et/ou bathymétrie) d'une zone terrestre. En cartographie les altitudes sont habituellement représentées par des courbes de niveaux. Le MNT en notre possession (source SHOM) utilise un maillage régulier carré de 50 mètres. Il permet ainsi de reconstituer une vue en image de synthèse du terrain (figure 3.6, image LIDAR : LIght Detection And Ranging).

La morphobathymétrie reflète très souvent la morphologie au toit du socle cristallin (D. Menier, géologue). En effet la corrélation linéaire entre la profondeur des fonds marins et celle au substratum est de 0,84. Le MNT peut donc fournir une information complémentaire lors de la mise en oeuvre du modèle de cokrigeage.

3/ Etude exploratoire de l'échantillon de la variable d'intérêt

La mise en oeuvre des techniques d'estimation géostatistique exige une analyse préalable des données expérimentales. L'étude exploratoire a pour but d'apprécier la distribution des données dans l'espace, d'appréhender leur degré d'homogénéité, de rechercher et de visualiser les observations atypiques ou tout simplement de se familiariser avec la variable.

Il est ici rappelé que la variable d'intérêt étudiée est la profondeur à laquelle nous atteignons la roche en milieu marin (le toit du socle cristallin ou substratum).

3.1/ Histogrammes et statistiques descriptives

La droite d'Henry : si les probabilités sont alignées le long de la droite d'Henry, nous pouvons supposer que la variable suit une distribution gaussienne.

Le box plot est une autre façon d'appréhender la distribution de la variable. Il a l'avantage de mettre en valeur les quartiles. Les quartiles Q1, Q2 et Q3 sont les éléments essentiels de ce graphique. Le quartile Q2 (ou médiane) partage la série de données en deux groupes d'effectifs égaux. Le quartile Q1 partage le groupe dont les valeurs sont inférieures à la médiane, en deux groupes d'effectifs égaux. Q3 fait de même pour le groupe dont les valeurs sont supérieures à la médiane.

La distribution des profondeurs au toit du socle cristallin (rocheux) est représentée par les trois graphiques de gauche de la figure 3.8. Ils montrent que les valeurs sont comprises entre 2 mètres au-dessus et 25 mètres en dessous du niveau zéro (référence SHOM). L'étendue de 27 mètres peut nous paraître importante. Mais en y regardant de plus près, on s'aperçoit que la moyenne et la médiane sont toutes deux égale à -10 mètres. Il en résulte que 1/4 des valeurs sont comprises entre -5 et -10 mètres et qu'un autre quart entre -10 et -15 mètres. Sachant qu'il n'y a qu'un très petit nombre de valeurs supérieures à zéro ou inférieures à -20 mètres, on peut parler d'une distribution équitablement répartie sur son étendue (25% des valeurs sur 5 mètres d'étendue).

La droite d'Henry et un test de normalité (Kolmogorov-Smirnov) montrent que la distribution des profondeurs au toit du substratum ne suit pas une loi Gaussienne. Par contre le coefficient d'asymétrie (cf. tableau du bas) et l'histogramme permettent de penser qu'elle s'en approche.

Si la propriété de normalité a l'avantage de définir complètement la distribution par sa
moyenne et son écart type, elle n'est pas une condition obligatoire en géostatistique linéaire.
Celle-ci ne fait aucune hypothèse particulière sur la distribution des valeurs. Toutefois d'après

C. Cressie (Statistics for Spatial Data, 1993 p.110), le krigeage a tendance à fournir de meilleures prévisions lorsque les données suivent une loi normale.

A titre d'exemple, l'histogramme des valeurs du MNT n'est pas symétrique. D'un côté, 75% des profondeurs se regroupent entre +2 mètres et -8 mètres (10 m d'amplitude) et de l'autre coté, 25% des valeurs s'étendent entre -8 mètres et -22 mètres (14 m d'amplitude). Certaines valeurs extrêmes apparaissent autour des -20 mètres de profondeur. Alors pour estimer le MNT, Cressie recommande soit une transformation des données par une fonction mathématique pour rendre la distribution gaussienne, soit discrétiser en zones géographiques afin que les profondeurs dans chaque zone soient proches les une des autres. Chaque zone géographique est alors traitée individuellement.

Le coefficient d'asymétrie conforte le fait que la variable d'intérêt a une distribution presque symétrique. Enfin, le coefficient de corrélation linéaire entre les deux variables est égal à 0,84 ce qui indique une forte dépendance entre les variables.

3.2/ Comportements directionnels : Etude de stationnarité

Pour étudier le comportement dans l'espace des variables observées, un outil exploratoire synthétique consiste à visualiser les valeurs des variables le long des directions de coordonnées.

Comportement directionnel de la variable d'intérêt Comportement directionnel de la variable secondaire

Le tracé de l'évolution des profondeurs en fonction de la longitude (axe Est-Ouest) et de la latitude (axe Sud-Nord) est un moyen de détecter la présence d'une dérive spatiale qui, si elle existe, devrait se traduire graphiquement par une croissance ou une décroissance du nuage de points.

Le nuage de points de la variable d'intérêt (figure 3.9 de gauche), est relativement homogène autour de la moyenne dans la direction SN, en revanche les valeurs ont une légère tendance à augmenter lorsque l'on se dirige de l'est vers l'ouest. Une hypothèse de stationnarité serait adéquate dans la direction sud-nord (tendance linéaire de pente nulle), mais pas forcément dans la direction est-ouest, pour laquelle la moyenne n'est pas constante. On pourrait avoir recours à l'hypothèse de stationnarité intrinsèque pour décrire son comportement le long de cette direction.

La variable MNT semble quant à elle être stationnaire au regard de son comportement directionnel (nuage de point relativement homogène et dispersion à peu près constante le long des directions principales).

Les valeurs atypiques indiquées sur la figure 3.9 sont pour certaines de fortes profondeurs dues au canal causé par le passage étroit entre l'île aux Moines et le continent au nord de la zone d'étude. Nous décidons de conserver ces données pour la suite de l'étude. Mais d'autres sont des incohérences entre la source des données sismiques et celle du MNT. Elles se localisent autour du profil 34 (figure 3.5) situé dans une zone rocheuse. Les points proches à valeurs discordantes sont supprimés.

3.3/ Nuage de corrélation différée ou « h-scatter plots »

Pour un vecteur de séparation h donné, le nuage de corrélation différée est le nuage des points formés par les couples de valeurs (z(si), z(si+h)) ; z(si) représente la profondeur de la roche au point si.

Pour que cet outil soit pertinent, il faut que les observations soient équitablement réparties sur
l'ensemble du champ. Plus les valeurs prises aux sites si et si+h sont semblables, plus le nuage

est proche de la première bissectrice ; au contraire, plus ces valeurs sont différentes, plus le nuage est diffus.

Le nuage de corrélation différé est donc un outil qui représente la structuration des valeurs pour une certaine séparation (h).

Dans notre exemple la distance moyenne au voisin le plus proche est d'environ 50 mètres. La figure 3.10 nous montre les nuages de corrélation différée pour un vecteur de séparation de longueur 50 mètres (plus ou moins 5 mètres) orienté dans la direction est-ouest et sud-nord. Les valeurs prises aux sites séparés par 50 mètres dans l'axe est-ouest se ressemblent plus que celles prisent dans l'axe sud-nord. La corrélation est de 0,9 contre 0,8 pour l'axe SN.

Certains points se détachent particulièrement du nuage ; un examen détaillé à l'aide de la figure 3.6 (carte LIDAR : LIght Detection And Ranging), révèle que ces valeurs proviennent de zones de profondeur très localisées comme dans le chenal.

Nous pouvons déduire des « h-scatter plots », que l'écart moyen entre la valeur prise en un site si et un autre en si+50 mètres se rapproche de zéro sur l'axe horizontal comme sur l'axe vertical. L'analyse variographique nous confortera dans cette hypothèse.

3.4/ Description bivariée

Il est important d'examiner la relation entre la variable auxiliaire et la variable d'intérêt. Pour ce faire, les techniques disponibles (coefficient de corrélation, covariance, graphiques) requièrent des couples de données aux points d'observations s1 à s_n. Etant donnée que la variable régionalisé auxiliaire n'a pas été observée en ces points, une interpolation spatiale doit être préalablement effectuée afin de prévoir ces valeurs régionalisées.

Le MNT (variable auxiliaire) est connu sur un quadrillage d'un espacement de 50 mètres. Il
en résulte que la distance maximum qui sépare une observation de la variable principale à la

plus proche observation de la variable secondaire est de moins de 50 mètres. Cette distance est très petite à l'échelle du champ. C'est pourquoi nous avons décidé d'utiliser la méthode d'interpolation des plus proches voisins. Les voisins les plus proche seront les plus représentatifs de la donnée à estimer.

Ceci fait la figure 3.10 montre la dépendance linéaire des deux variables. Les points alignés sur la bissectrice sont les points correspondant aux zones de roche (non couvert d'une couche sédimentaire). Nous leur avions effectivement attribué la profondeur donnée par le MNT.

Le coefficient de corrélation entre la variable d'intérêt et la variable auxiliaire est de 0.84. Mais si on ne prend pas en compte les points des zones de roches ce coefficient décent à 0,75.

Y : Profondeur au toit du socle cristallin X : Profondeur à la surface du fond marin

4/ Transformation des coordonnées

Les coordonnées des points sont en longitude, latitude de projection WGS84. C'est un système géodésique dont l'unité est exprimée en degrés.

Projection WGS84 : une projection est une correspondance entre les points de la terre et ceux de la représentation plane. La projection WGS84 est un système universel, mais il en existe un grand nombre suivant l'endroit où l'on se place sur le globe.

Pour revenir à un système de coordonnées planes (une unité en ordonné égale une unité en abscisse), dans le but de calculer une distance euclidienne nous appliquons une transformation sur les coordonnées.

5/ Conclusion

L'étude exploratoire nous a permis de nous familiariser avec les données. Elle montre comment la distribution de la profondeur au toit du socle cristallin se rapproche d'une loi gaussienne. Il existe bien une relation forte entre cette profondeur et celle du fond marin. Enfin, les nuages des valeurs le long des directions de coordonnées et les nuages de corrélations différées permettent de vérifier la ressemblance des profondeurs dans un rayon de 50 mètres autour d'une observation. Il faut à présent étudier plus précisément la structure spatiale de la variable d'intérêt. C'est l'analyse variographique qui va y répondre.

Chapitre IV : Analyse variographique

L'analyse variographique est l'étape préalable au krigeage. Le manque d'observations de la variable à prédire, l'abondance d'information de la variable auxiliaire et l'information de position (longitude/latitude) pour chacune des observations, nous ont emmené à choisir l'interpolation par « krigeage ».

L'analyse variographique est menée afin d'estimer la fonction de la continuité spatiale de la variable d'intérêt. Dans ce chapitre nous présentons les variogrammes expérimentaux selon diverses directions, l'étude du comportement à courte distance et enfin l'ajustement d'un modèle.

1/ Variogrammes expérimentaux

Le variogramme caractérise la continuité spatiale de la variable régionalisée. Sur la figure 4.1, on trouve le variogramme expérimental omnidirectionnel, calculé pour des distances multiples d'un pas de 50 mètres. L'intervalle de confiance à 95% accordé à chaque point est indiqué par une barre verticale. Cet intervalle de confiance est fonction de la variance et du nombre de paires de points intervenant dans le calcul des points expérimentaux.

En supposant que la distribution Z(si)-Z(si+h) pour i=1 à n, suit une loi gaussienne de moyenne nulle et de variance 2_ã(h), on a 2 ( )

libertés. L'intervalle de confiance à á=95% de la variance est alors donné par :

n : nombre de paire de points intervenant dans le calcul des points expérimentaux

ã à ( h ) : valeur du point expérimental à la distance h (cf. formule estimation du variogramme) _÷2_á : fractile d'ordre á de la loi du ÷² à n degrés de liberté

Avant même de regarder plus précisément les intervalles de confiances, Arnaud et Emery (2000 p .126) recommandent de ne tenir compte dans le calcul du variogramme expérimental, que les distances allant jusqu'à la moitié de la distance maximale rencontrée entre deux points du champ (cf. figure 4.1 de droite). Au-delà, le nombre de couples de points intervenant dans le calcul du variogramme décroît et lui fait perdre en robustesse (cf. chap. II § 3.2).

Le variogramme omnidirectionnel atteint un palier à une distance de 1100 mètres. On peut en déduire que la variable régionalisée (profondeur au toit du substratum) est une réalisation d'une fonction aléatoire stationnaire d'ordre 2 (M. Arnaud et X. Emery, 2000 p.122 / Journel et Huijbergts, 1993 p.37).

1.2/ Variogrammes directionnels

La figure 4.3 montre les variogrammes directionnels, le long des quatre directions principales du plan, d'orientation -35°, 55°, 100° 145° par rapport à la direction est-ouest. Ces directions ont été choisies en fonction de l'orientation des profils 13 et 12 (cf. figure 3.5) et de la géographie observée sur l'image LIDAR (figure 4.2).

La bathymétrie de la zone d'étude montre une continuité spatiale dans la direction 55°. On remarque que si l'on se déplace sur un axe du sud-ouest au nord-est, la profondeur du fond marin reste sensiblement la même, alors que si l'on se déplace du nord-ouest au sud-est on doit faire face à de forte variation de profondeurs. Cette morphologie est le résultat de la baie formé au nord-ouest et de l'écoulement de l'eau provenant en grande force le l'embouchure au nord-est.

Les variogrammes de la figure 4.3 sont sensiblement différents. La différence d'accroissement à l'origine est le symptôme d'une anisotropie. C'est un paramètre qu'il faut prendre en compte.

1.3/ Anisotropie

Une anisotropie caractérise un phénomène qui se déploie préférentiellement dans certaines directions de l'espace.

Les notions d'anisotropie, de continuité maximale et minimale peuvent être schématisées de la manière suivante :

Anisotropie : Bien que la distance entre le point central et celui en direction 55° soit plus grande que celle entre le point central et le point en direction -35°, les valeurs prises en ces points sont très proche (la corrélation est la même).

L'anisotropie géométrique : On observe dans diverses directions des paliers et des composantes pépitiques identiques mais des portées différentes. Les portées maximales (a_g) et minimales (a_p) s'observent selon deux directions orthogonales (cf. chap. II § 3.4).

L'anisotropie zonale : On observe des paliers différents selon les directions. Il n'existe pas de modèle d'ajustement pratique pour traiter ce type d'anisotropie.

En observant la figure 4.3, on remarque que les paliers ne sont pas identiques. Mais si l'on considère que la variable régionalisée est stationnaire d'ordre 2, le palier est approximativement égal à la variance. Les directions de continuités maximale et minimale sont alors respectivement de 55° et -35° par rapport à l'axe ouest-est avec une portée de 1100 mètres et 250 mètres.

De cette façon nous nous replaçons dans un système de coordonnés redressé de 35° dans le sens inverse des aiguilles d'une montre. Pour prendre en compte les portées dans les deux directions on résout le système

C'est cette seconde méthode que l'on retient pour l'estimation par krigeage. 1.4/ Comportement à l'origine du variogramme

Le comportement du variogramme à l'origine, reflète le degré de régularité spatiale de la variable régionalisée (cf. chap. II § 3.2).

- Un comportement de ãà ( h ) à l'origine parabolique indique une grande régularité de la variable régionalisée (continue et différentiable).

- Un comportement à l'origine linéaire montre que la variable régionalisée est moins régulière (elle est continue mais pas différentiable)

- Une discontinuité à l'origine ou effet de pépite, signale une plus grande irrégularité de la variable régionalisée. Il y a absence (partielle ou totale) de corrélation entre les valeurs prises en deux sites très proches. (Journel et Huijbregts, 1993 p.38).

D'après de la figure 4.1 et 4.3, on constate que les variogrammes ont un comportement à
l'origine linéaire ( ãà ( h ) 10*h). La profondeur au toit du socle cristallin est donc une variable

Comme dans le cas monovarié, l'inférence statistique nécessite des hypothèses de stationnarité de la corégionalisation (ensemble des variables régionalisées).

Pour analyser les liens spatiaux entre la profondeur du substratum et la profondeur à la surface du fond marin, la figure 4.5 (de droite) montre les variogrammes croisés selon les directions principales d'anisotropie. A sa lecture, nous faisons l'hypothèse d'une stationnarité d'ordre deux car ils atteignent tous deux un palier environ égal à la covariance.

D'autre part, nous observons que les variogrammes des deux variables régionalisées et le variogramme croisé se ressemblent fortement. Ca nous emmène à choisir le même modèle de base pour leurs ajustements.

2/ Ajustement à un modèle

Le variogrammme expérimental calculé et ses propriétés étudiées, il faut ajuster une courbe théorique. Elle doit être définie pour toutes les distances et toutes les directions de l'espace. Cependant toute fonction mathématique ne peut être utilisée comme modèle. Un modèle variographique doit être une fonction de type négatif conditionnel (Cressie, 1993 p.86) c'est-àdire que :

C'est pourquoi le modèle de variogramme est choisi parmi un ensemble de fonction dont on sait qu'elles sont de type négatif conditionnel (cf. partie II, chap.3.3).

Remarque : quel que soit le mode d'ajustement retenu, la modélisation du variogramme aux courtes distances est particulièrement importante.

2.1/ Cas monovariable

Des anisotropies dans les directions -35° et 55° ont été détectées dans le paragraphe 1.3. Il convient d'ajuster le variogramme dans ces directions. D'autre le palier dans les directions principales est égal à la variance de la variable régionalisée, soit 38,6. N'ayant pas observé d'effet pépite, le comportement à l'origine de la variable aléatoire régionalisée est linéaire.

Ces observations amènent à choisir entre un modèle sphérique et un modèle exponentiel. Le modèle sphérique est préférable, d'une part parce qu'il s'ajuste mieux sur les points de petite distance du variogramme expérimental et d'autre part parce qu'après plusieurs essais il donne de meilleurs résultats.

M. Arnaud et X. Emery montrent que la variance d'estimation (variance de krigeage) est plus élevée dans le cas du modèle exponentiel ; l'explication vient du fait que ce dernier croît plus rapidement que le modèle sphérique, ce qui traduit un phénomène qui se déstructure plus vite, d'où une moins grande précision dans l'estimation.

Remarque : L'estimation du paramètre de la portée se fait en résolvant le problème :

ùi : poids accordé à la distance hi ( les poids sont inversement proportionnels à la distance).

D'autres fonctions variographiques ou paramètres d'ajustement auraient pu être choisis. C'est pour cela qu'il est nécessaire de contrôler la qualité du modèle d'ajustement. On effectue à cette fin une validation croisée.

Pour tester les résultats de la validation croisée, nous comparons ceux obtenus pour le modèle précédent et ceux obtenus pour le modèle variographique qui ne tient pas en compte l'anisotropie (cf. figure 4.7).

D'après le rapport INERIS (2003), la validation croisée doit être réalisée avant d'entreprendre le krigeage. Elle fournit des critères statistiques de sélection dans le choix d'un modèle de variogramme.

La validation croisée consiste à éliminer temporairement un point de l'ensemble des données puis à estimer sa valeur par krigeage à l'aide des données restantes et du modèle de variogramme qui a été ajusté. Cette opération est répétée pour tous les points.

Ainsi en tout point d'observation si, est calculé par krigeage, avec le modèle variographique

retenu, une valeur estimée Z ( s_i ) et un écart-type de krigeage óok(si). Le rapport INERIS

- la moyenne et la variance de l'erreur standardisée ([Z(si)- Z ( s_i ) ]/ óok(si)),

- la moyenne des erreurs d'estimation et des erreurs réduites (standardisées) est plus proche de 0, ce critère assure l'absence de biais,

- la variance des erreurs d'estimation est plus faible, ce critère traduit la robustesse de l'estimateur et renseigne sur la précision de l'estimation,

- la variance des erreurs standardisées est plus proche de 1, ce critère indique que

- la moyenne des erreurs relatives est plus proche de 0, ce critère traduit la bonne

- la corrélation entre Z(si) et Z ( s_i ) est plus proche de 1 et le nuage de corrélation

N=441	Modèle avec anisotropie	Modèle isotrope
Erreur Quadratique Moyenne	3,34	3,18
Moyenne des erreurs d'estimation	-0,01	-0,001
Ecart-type des erreurs d'estimation	1,83	1,78
Moyenne des erreurs relatives	-22,16 %	-21,18%
Moyenne des erreurs standardisées	0,003	0,0007
Ecart-Type des erreurs standardisées	0,66	0,63
Coefficient de corrélation	0,95	0,96
% de données dont l'erreur standardisée est inférieure à 2,5	99	96

Les nuages sont concentrés le long de la première bissectrice, ce qui indique une bonne précision des estimations. Les deux modèles rendent une estimation globalement sans biais, les erreurs d'estimation ne s'écartent pas fortement de 0. D'après Arnaud et Emery (200, p.156), pour valider un modèle on espère au moins 95% de données dont l'erreur standardisée est inférieure à 2,5. C'est le cas pour les deux modèles. Les écarts-types sont proches de 1, les erreurs relatives n'excèdent pas 23%. En revanche quelques fortes profondeurs sont sous estimées.

Il est également intéressant de localiser les erreurs d'estimation. La figure 4.8 nous les montre.

Les cartes montrent que l'on estime mal aux endroits à fort dénivelés (à l'embouchure au nord) et certains points au sud-est. Ces erreurs sont dues à des incohérences entre sources de données différentes (MNT et Sismique, cf. chap. III § 2.1.2).

Après application des deux modèles, le géologue a préféré les résultats obtenus par le modèle à anisotropie.

3/ Conclusion

Dans ce chapitre nous avons exposé les bases de l'analyse variographique et les avons illustrées par une étude de cas réel. La recherche d'un modèle convenable est assez délicat et nécessite une certaine expérience ou de nombreux essais. C'est une procédure qui ne peut être automatisée. Elle nécessite également les connaissances du géologue ou de l'expert sur le terrain pour valider la véracité du phénomène observé (régularité spatiale, variabilité à courtes distances et anisotropie).

Chapitre V : Estimation par krigeage

L'étape de l'estimation par krigeage est relativement rapide si l'étude variographique a été correctement menée. Il ne s'agit ici plus que de définir une grille d'estimation et un voisinage de krigeage (cf. chap. II, § 4.1). L'estimation par krigeage peut alors être entièrement automatisée.

1/ Choix de la grille et du voisinage de krigeage

Nous avons choisi la même grille que celle proposée par le SHOM pour représenter le MNT. Elle facilitera la comparaison entre la morphobathymétrie et la morphologie au toit du socle rocheux. C'est une grille régulière de points espacés d'une distance de 50 mètres. Au regard de l'échelle de la zone d'étude ( 750 ha), le degré de détail de cette grille est élevé. Elle permettra de calculer facilement les courbes de niveaux et de représenter la géographie au toit du socle cristallin en trois dimensions.

Le choix du voisinage de krigeage est un peu plus délicat. Le nombre de points inclus dans le voisinage doit être assez grand pour une estimation de précision. Il dépend également de la continuité spatiale de la variable régionalisée. Le modèle variographique doit être acceptable à l'échelle de ce voisinage. C'est-à-dire qu'il existe une dépendance spatiale. On pourra tester plusieurs tailles de voisinage par validation croisée et choisir celle qui donne les meilleurs résultats. Dans le cas présent nous avons retenu un voisinage de krigeage d'un rayon de 500m. C'est entre la portée de continuité maximale et la portée de continuité minimale du variogramme directionnel.

Enfin, le choix du système de krigeage qui a été fait est le krigeage ordinaire. Comme nous l'avons relaté dans le chapitre II, § 4, il ne nécessite pas la connaissance de l'espérance de la variable régionalisée. Au lieu de l'estimer, il est souvent préférable de la considérer comme inconnue car elle peut varier d'une zone à l'autre du champ. Toutefois, nous rappelons qu'il faut qu'elle reste constante dans le voisinage de krigeage.

2/ Résultats obtenus par krigeage ordinaire

La figure 5.1 montre le résultat de l'estimation par krigeage et l'écart-type de krigeage. Le krigeage utilise le modèle variographique suivant :

avec une portée (a) qui diffère suivant la direction (cf. chap.4 prg. 2.1) L'anisotropie est corrigée par un changement de repère (cf. chap.4 prg. 1.3).

Les profondeurs varient entre 3 et -25 mètres. Une mesure au-dessus de zéro n'est pas aberrante. Cela vient du fait que le niveau zéro de référence est le niveau de l'eau au coefficient le plus fort à marée basse.

L'écart-type de krigeage est un indicateur de la précision de l'estimation. Il est fiable car l'écart-type des erreurs standardisé est proche de 1 (cf. validation croisée Chap. IV § 2.2). Il quantifie la dispersion possible de la valeur vraie, mais inconnue, autour de la valeur estimée. Par contre, il dépend uniquement du modèle de variogramme et de la configuration des données dans le voisinage de krigeage.

A la vue de la figure 5.1, les variations spatiales de l'écart-type de krigeage indiquent la perte de précision lorsqu'on s'éloigne des points de mesure. Elles renseignent sur les zones où l'échantillonnage est suffisamment dense (faible écart-type, bonne précision au centre et au sud de la zone d'étude) et sur celles où l'échantillonnage est trop espacé (grand écart-type, mauvaise précision au nord de la zone d'étude).

3/ Résultats obtenus par cokrigeage ordinaire

Le cokrigeage renforce le krigeage standard en y ajoutant l'information apportée par une deuxième variable. C'est un immense avantage quand l'information de la variable principale vient à manquer. On peut aussi s'imaginer optimiser les coûts des campagnes d'échantillonnage en prélevant la profondeur au toit du socle cristallin uniquement aux endroits à fort volume sédimentaire.

Les principes de base du cokrigeage sont les mêmes que ceux du krigeage. On veut former une estimation linéaire de la variable principale à partir d'observations de la variable principale et de la variable secondaire.

Mais, en plus de créer un variogramme qui représente les variations spatiales de la variable principale entre deux points d'observation, le cokrigeage introduit la notion de variogramme croisé (cf. chap. II § 3.5 ). Tout comme dans le cas du krigeage, le variogramme croisé est modélisé. Il permet alors de résoudre les équations qui fournissent les valeurs des coefficients de pondération en fonction des observations des deux variables utilisés. De la même manière, une évaluation de l'erreur d'estimation est disponible en chaque point d'interpolation. D'autre

part cette variance d'estimation est toujours inférieure à celle obtenue par krigeage (Arnaud et Emery, 200 p.195).

Le variogramme croisé quantifie la corrélation spatiale entre la variable principale et la variable secondaire.

Variogramme directionnel croisé
(Profondeur au toit du socle rocheux / Profondeur à la surface du fond

Le palier est égal à la covariance entre la profondeur au toit du socle cristallin et celle à surface du fond marin.

3.2/ Validation Croisée

La validation croisée est utilisée pour comparer deux approches. La première approche est la même que précédemment, c'est-à-dire que nous avons construit un modèle variographique à partir de données mesurées et de données d'expertise. Ces dernières ont été jugées par un géologue connaissant le terrain (cf. chap. III § 2.1.2). Lors du krigeage standard nous avons gardé ces points pour estimer l'ensemble du champ. Or ce ne sont pas des valeurs exactes. La deuxième approche propose de bénéficier de l'information apportée par la deuxième variable sans prendre en compte les points d'experts lors de l'estimation.

Nous n'avons pas montré la deuxième approche pour l'estimation par krigeage car les résultats n'en valaient pas la peine.

N=441 / N=323	Estimation avec points experts	Estimation sans points experts
Erreur Quadratique Moyenne	2,08	1,20
Moyenne des erreurs d'estimation	-0,02	-0,01
Ecart-type des erreurs d'estimation	1,44	1,09
Moyenne des erreurs relatives	-3,52 %	-1,03%
Moyenne des erreurs standardisées	-0,006	-0,003
Ecart-Type des erreurs standardisées	0,83	0,68
Coefficient de corrélation	0,97	0,98
% de données dont l'erreur standardisée est inférieure à 2,5	98	98

Les indicateurs statistiques des deux méthodes sont proches. Néanmoins l'EQM est bien meilleur quand on ne rajoute pas les points d'expertises. Nous supposons que c'est du aux points d'expertises difficiles à estimer car ils ne reflètent pas forcément la réalité.

Nuage de corrélation entre valeurs vraies et estimées (par cokrigeage)
Figure 5.3

La figure 5.3 compare les nuages de corrélations entre vraies valeurs et des valeurs estimées des deux approches.

Les nuages sont bien plus concentrés le long de la diagonale que ceux de l'estimation par krigeage, ce qui indique une meilleure précision des estimations.

Attardons-nous sur le nuage de gauche, un grand nombre de points d'une profondeur de 15 mètres ont été estimé entre 10 et 20 mètres. Après une recherche de localisation il s'avère que ce sont des points d'expertises ; un élément de plus pour les juger non fiables. Il manque une cohérence entre les points d'expertises et les points qui les entourent.

La figure 5.4 localise les erreurs d'estimation. Elles sont moins fortes que celle faites par krigeage. Au lieu de varier entre -8 et 8 mètres elles ne fluctuent plus qu'entre -8 et 4 mètres. Par contre, les erreurs les plus importantes se localisent toujours aux mêmes endroits.

A la vue de la figure 5.5, on s'aperçoit que les résultats obtenus par cokrigeage sont plus fins que ceux obtenus par krigeage. Le chenal présente une plus forte rugosité, la profondeur diminue en s'approchant de la côte et les monticules rocheux séparant le puits au centre de celui au sud sont plus lisibles.

N'ayant pas de points de mesures au nord ni au niveau de l'embouchure, le cokrigeage n'a pas pu estimer cette zone sans les points d'expertises. En revanche nous supposons que l'estimation sans les points d'expertises, du centre de la zone d'étude, se rapproche plus de la réalité car nous n'avons pas introduit de biais.

L'écart-type du cokrigeage est plus faible que celui du krigeage dans une grande partie de la zone d'étude.

En d'autres termes, l'information apportée par la variable auxiliaire, fortement corrélée avec la variable cible, permet de compenser en grande partie le manque d'information dû à un sous échantillonnage de la variable cible. C'est là l'un des intérêts majeurs du cokrigeage, lorsque les mesures de la variable cible viennent à manquer ou, qu'il est moins coûteux de récupérer une ou plusieurs variables auxiliaires que de mesurer plus largement la variable cible.

Les figures suivantes (figures 5.6) montrent différentes représentations des résultats obtenus.

A titre de comparaison, nous avons dessiné la carte des profondeurs au toit du socle rocheux en utilisant les résultats de quatre méthodes d'interpolation.

Les trois méthodes figurants à côté de celle par cokrigeage sont des techniques d'interpolation automatisées, intégrées dans des logiciels de système d'information géographique (SIG). Leurs procédures de calculs ont été définies dans le chapitre I.

Cette comparaison illustre les résultats des méthodes déterministes. Ils sont différents, mais permettent néanmoins une représentation globale du phénomène à étudier.

Chapitre VI : Les épaisseurs de couches sédimentaires 1/ Epaisseur de sédiment

L'étude doit permettre d'aboutir à l'épaisseur des formations sédimentaires. C'est la différence entre la profondeur des fonds marins et celle au toit du socle rocheux. La figure 6.1 a été obtenue à partir des profondeurs au toit du substratum estimées par cokrigeage et des données bathymétriques délivrées par le SHOM.

La baie est recouverte d'une couche sédimentaire de 5 à 10 mètres d'épaisseur. Des dunes apparaissent le long ouest du chenal. Le socle rocheux du sud-est de la zone est dépourvu de formation sédimentaire.

2/ Courants marins

La répartition, la nature et l'épaisseur des sédiments, mais aussi la profondeur du substratum, dépendent des conditions hydrodynamiques.

Huit heures par jour, il se forme un tourbillon au centre ouest de la zone d'étude (figure 6.2). Les particules de sédiments en suspension dans l'eau se posent alors au fond et attaquent la roche. C'est pourquoi d'importantes masses de sédiments sont constatées.

Les courants en vive eau ramènent les sédiments entassés au nord de l'embouchure (figure 6.1 et 6.2), pour les déposer 500 mètres plus bas lorsque l'intensité du courant s'affaibli. A l'inverse, en morte eau, les courants remontent ces mêmes masses de sédiments mais qui sont pour une partie stoppés par la pointe de l'île aux Moines. Ainsi il se forme des monticules de sédiments juste au sud de cette pointe.

La partie est de la zone est dépourvue de sédiment car constamment parcourue par les courants.

Les courants étant très faibles dans la baie, les sédiments stagnent et ont peu d'effet érosif.

(Source: Modélisation hydrodynamique du GdM - Direction des études et recherches EDF)

La concentration ou la dispersion de l'énergie formée par les courants explique les contrastes des fonds marins.

En sédimentologie, un outil de prélèvement très utilisé est la benne (figure 7.1). Elle permet d'explorer la composition des fonds marins.

L'Ifremer définit la crépidule comme un mollusque gastéropode (crepidula fornicata) originaire d'Amérique du nord et introduit involontairement en Europe. Elle s'étend maintenant de la Suède à la Méditerranée. Son abondance fait qu'elle entre en compétition alimentaire avec d'autres coquillages. Elle a par ailleurs un impact sur les fonds qu'elle colonise (dépôts de coquilles et de déchets, envasement). C'est un mollusque qui a une capacité extraordinaire de prolifération (jusqu'à 10 000 individus au mètre carré). Il a déjà colonisé toute une partie de la côte bretonne

L'objectif de cette présentation est de montrer qu'à partir d'un semis de points on peut analyser une structure spatiale^*. Chacun des points (prélèvement par bennes) nous indique la présence ou non de crépidules. La connaissance de la structure spatiale, nous permet d'utiliser une méthode de lissage pour créer une carte des zones occupées par les crépidules. Nous obtenons alors une idée du peuplement de cet envahisseur dans la zone étudiée.

2/ Zone d'étude

Pour les besoins de cette étude, les relevés par bennes sont localisés dans un secteur à un kilomètre au nord ouest de l'île Tascon. La surface représente un rectangle de 900 mètres d'est en ouest et de 950 mètres du nord au sud. 70 bennes sont relevées sur une surface de 85 ha. Pour chaque prélèvement on note une présence ou non de crépidules.

3/ Statistiques spatiales descriptives

Cette partie propose des méthodes d'analyse concernant uniquement la position des objets, indépendamment de la valeur d'un attribut de ces objets. On présente des indicateurs permettant de décrire la concentration ou la dispersion du semis étudié. Le semis est-il concentré ou réparti avec quel degré de certitude ?

3.1 Position et dispersion

Intéressons-nous aux bennes à présences de crépidules. Compte tenu de leur positionnement, le barycentre se trouve à 200 mètres au sud du point central de la zone d'étude. C'est une première indication pour dire que les bennes à présence de crépidules ne sont pas également réparties sur toute la zone d'étude.

Les trois indicateurs de tendance centrale sont convergents, le point médian et le point de distance minimale se situent à moins de 50 mètres du barycentre.

La dispersion standard des crépidules autour du barycentre est d'environ 350 mètres (distance euclidienne) à vol d'oiseau. L'ellipse de déviation standard (représente le coefficient d'asymétrie) est un tout petit peu plus étirée de l'est à ouest que du nord au sud en raison de la forme de la zone d'étude. Il est intéressant de constater que 60% des bennes à présences de crépidules sont comprises dans l'ellipse. Cela indique une population concentrée dans un espace restreint de la zone d'étude.

3.2 Sectorisation

La sectorisation consiste à rechercher des groupes de sites qui ont la même valeur (présence ou non de crépidules) et qui minimise la distance entre les sites. C'est une façon d'examiner si des groupes homogènes se forment.

En regardant la localisation des bennes à présence de crépidules et celle sans, des frontières apparaissent naturellement (figure 7.3.1).

La figure 7.3.2 montre une classification en 7 groupes par la méthode des K points moyens. C'est une méthode fondée sur la distance (euclidienne).

Les secteurs de la figure 3.2 ont un barycentre qui minimise la dispersion des bennes à présence de crépidules et celles sans crépidules.

3.3 Méthode des quadrats

La répartition des bennes à présence de crépidules peut-elle être assimilée à une répartition aléatoire à l'échelle de la zone d'étude ? Est-elle au contraire plutôt répartie ou plutôt concentrée ?

Soient (Q1 ,Q2 ,..., Qn) les localisations des n quadrats, l'hypothèse d'une répartition spatiale complètement aléatoire est celle issue d'un processus ponctuel homogène de Poisson et reposant sur deux principes :

- un principe d'uniformité, selon lequel chaque localisation de la zone d'étude possède la même probabilité d'être habitée par des crépidules ;

- un principe d'indépendance, en vertu duquel l'occurrence d'un événement en un lieu n'influence en aucune manière l'occurrence d'événements en d'autres lieux.

Autrement dit, on suppose que la zone d'étude est parfaitement homogène et qu'il n'existe aucune interaction entre les individus localisés. Sous ces conditions, ce modèle purement aléatoire peut être formalisé au moyen d'une distribution de probabilités de Poisson.

L'analyse par quadrat dénombre les sites dans les cellules du carroyage (quadrats) et compare les effectifs observés aux effectifs théoriques résultant d'une loi binomiale ou de son approximation par la loi de Poisson (ë=np).

Indice	Valeur calculée (1)	Valeur calculée (2)	Interprétation, la distribution est :
D (rapport de la variance sur la moyenne)	1,34	1,67	1 : aléatoire >1 : concentrée <1 : régulièrement espacée
ICS (indice de dimension des agrégats)	0,34	0,67	0 : aléatoire >0 : concentrée <0 : régulièrement espacée
GI (indice de Green)	0,01	0,02	0 : aléatoire 1 : très concentrée
IMC (indice de voisinage moyen)	2,49	1,31	x : aléatoire > x : concentrée < x : régulièrement espacée
IP (indice d'hétérogénéité)	1,15	2,05	1 : aléatoire >1 : concentrée <1 : régulièrement espacée

Conclusion, on peut considérer que la répartition des bennes à présence de crépidules n'est pas conforme à un processus de Poisson, mais qu'elle manifeste une concentration. On remarque effectivement une densité est plus forte au sud qu'au nord.

Les bennes sans présence de crépidules, sont aussi réparties comme un semis de points concentrés.

- Elle est plus adaptée à l'analyse spatiale d'une population que d'un échantillon. Elle considère que le semis et complet sur la zone d'étude. Exemple `localisation des villes dans une région'.

- Les résultats changent suivant le maillage choisi. Un maillage plus fin conclura à une répartition aléatoire des bennes.

- Elle ne prend pas en compte les observations voisines. Comme indiqué, elle sous-entend qu'en théorie la densité de crépidules est la même dans tout sous-ensemble de la région

(quadrat). C'est-à-dire que si la distribution des sites était répartie de façon régulière il y aurait une densité de 2.16 sites dans chaque carreau.

La méthode suivante tient compte des observations voisines et de la distance qui les sépare. 3.4 Analyse spatiale exploratoire basée sur les distances

Les méthodes suivantes font l'hypothèse que le semis de point étudié est homogène. C'est-àdire qu'on suppose que la densité est constante sur toute la zone d'étude.

3.4.1 Méthode du voisin le plus proche

Elle propose de comparer la distance moyenne d'un point à son plus proche voisin avec une distance théorique attendue sous l'hypothèse nulle de distribution aléatoire.

Soit N observations sur un territoire d'une superficie S. La distance moyenne théorique au

voisin le plus proche est d₀ = =57,83 (S/N : densité du semis de points).

La distance moyenne à vol d'oiseau entre deux prélèvements, à présence de crépidules, les plus proches dans la zone d'étude est de 52 mètres (distance euclidienne). La distance moyenne attendu théorique entre deux présences de crépidules dans un semis de points généré par un processus ponctuel aléatoire de même densité serait de 57 mètres. La statistique R est égale à 0,9 (52/57). Elle est à comparer avec :

- le modèle concentré de référence qui a une statistique R égal à 0, - le modèle aléatoire qui a une statistique R égal à 1,

Avec un R=0,9 nous nous rapprochons fortement d'une répartition aléatoire des crépidules.

On peut connaître la région critique d'acceptation de l'hypothèse nulle. L'intervalle de confiance pour une distribution aléatoire au risque á=5% de se tromper est :

Il en résulte que la distribution est significativement aléatoire, car la distance observée moyenne au voisin le plus proche (52,19) est comprise dans l'intervalle.

Elle serait significativement dispersée si la distance observée moyenne au voisin le plus proche était supérieure à d₀ + (1,96
· ó d0) = 65,88 ;

ou significativement concentré si la distance observée moyenne au voisin le plus proche était inférieure à d₀ - (1,96
· ó d 0 ) = 49,77 pour á=5%.

Remarque : les résultats sont différents de ceux apporté par la méthode des quadrats.

L'extension de la méthode aux 20 premiers rangs de voisinage rend les résultats présentés dans le tableau ci-dessous. La statistique R calculé est croissante avec l'augmentation des ordres de voisinages. Ce qui montre que plus l'ordre de voisinage augmente et plus le semis est régulièrement espacé.

Rang de voisinage	Distance au plus proche voisin	Distance théorique	Statistique R
1	52.19	57.83	0.90
2	88.45	86.74	1.01
3	112.44	108.43	1.03
4	136.01	126.50	1.07
5	155.83	142.31	1.09
6	171.94	156.55	1.09
7	190.72	169.59	1.12
8	212.38	181.71	1.16
9	226.85	193.06	1.17
10	239.60	203.79	1.17
11	254.94	213.98	1.19
12	265.91	223.71	1.18
13	276.83	233.03	1.18
14	287.35	241.99	1.18
15	297.66	250.63	1.18
16	309.94	258.99	1.19
17	320.64	267.08	1.20
18	333.50	274.94	1.21
19	344.34	282.57	1.21
20	354.22	290.01	1.22

D'autre part, si la distance au Kième voisin le plus proche est plus élevée que la distance théorique attendue alors le semis de points est régulièrement espacé à cet ordre de voisinage. Au contraire, une valeur observée faible indique un semis concentré.

La statistique du voisin le plus proche souffre de différents défauts. Elle utilise comme aire d'étude le polygone formé par les points limitrophes et non pas une aire plus large qui serait donnée par l'utilisateur. Selon Cressie (1993), la validité du test décline fortement avec l'élévation de l'ordre K. Il recommande de ne pas dépasser un niveau K<2,5%N. Or dans notre cas cette valeur est égale à 1,35. Les conclusions de l'extension aux 20 premiers rangs sont donc à relativiser.

D'un autre côté la méthode des plus proches voisins est simple à utiliser et ne nécessite pas de connaître toute la population.

3.4.2 Fonction de Ripley modifiée (L de Besag)

La fonction K de Ripley (1977, 1981) est une méthode plus récente pour analyser la structure spatiale d'un semis de points. L'hypothèse de base est que la densité locale à l'intérieur d'un cercle est égale à la densité de l'ensemble de la zone d'étude. Le principe est de calculer la probabilité d'observer un nombre donné de points dans un cercle de rayon d autour d'un point quelconque de la zone d'étude. La règle illustrée en figure 7.4 est répétée pour chacun des points, pour en déduire des résultats moyens.

Sous l'hypothèse nulle de la fonction de Besag, le semis de points est aléatoire. La figure 7.5 représente la courbe de la fonction L de Besag et celles des bornes de l'intervalle de confiance de l'hypothèse nulle.





0	50	100	150	200	250	300	35
0	50			200






0	50	100	150	200	250	300	35
0

Lorsque la fonction de Besag est comprise dans l'intervalle de confiance, on admet que la distribution observée est conforme à une distribution aléatoire de Poisson. Si la fonction de Besag est plus élevée que la borne supérieure de l'intervalle de confiance, alors la distribution est concentrée dans l'espace d'étude. Lorsque la statistique de Besag est plus faible que la borne inférieure de l'intervalle de confiance, alors la distribution est régulièrement dispersée dans l'espace d'étude.

Le calcul de la fonction L de Besag confirme et précise les conclusions tirées de l'analyse du voisin le plus proche. En effet, il y une absence d'organisation spatiale significative sur le semis de points avec présence de crépidules. Toutefois il existe une tendance à la concentration dans les rayons de 10 mètres et de 160 mètres de distance.

Les bennes sans présence de crépidules sont significativement concentrées dans un rayon allant de 130 mètres à 250 mètres.

L'avantage de cette méthode est qu'elle donne de l'information à plusieurs échelles de distance.

3.5 Conclusion de l'analyse spatiale descriptive

Les méthodes précédentes servent à tester une hypothèse globale sur la configuration spatiale de la zone d'étude. Elles permettent de savoir s'il existe un ordre spatial, mais ne sont pas capables de situer les agrégats.

Les méthodes reposent sur une hypothèse d'homogénéité de la variable sur la zone d'étude. Elle est contraignante et ne s'applique pas forcément à notre cas. Au contraire si l'on observe des crépidules sur un site alors il y a de fortes chances qu'il y en ait à côté.

Les méthodes n'ont pas été pensées pour étudier un échantillon mais pour une population connue sur l'ensemble de la zone d'étude. C'est pourquoi nous les utilisons seulement dans le but de localiser les zones susceptibles d'intérêt.

Les méthodes d'interpolations suivantes s'appuient sur une décomposition de la zone d'étude (cf. chap. I). La méthode de partitionnement la plus fréquemment utilisé est celle de Delaunay. Elle découle d'un découpage en polygones de Thiessen.

4.1 Polygones de Thiessen

On attribue à tout point appartenant au même polygone, la valeur du site ayant ce polygone d'influence (cf. chap. I § 1.1).

La validation croisée (cf. chap. IV § 2.2) pour cette méthode rend 81,08% de bonne réponse. 4.2 Interpolation par la méthode des plus proches voisins

L'algorithme se base sur une décomposition de l'espace en polygones de Delaunay (cf. chap. I § 1.1.1.2)

Ayant peu d'observations, et sachant que les crépidules vivent en colonie, si elles sont
présentes à un endroit il y a de forte chance que l'on en trouve juste à coté. D'où l'idée de

découper la zone d'études en une grille régulière à maille rectangulaire de côté mesurant 50 mètres. Elle est alors assez fine pour que chaque benne remontée corresponde à un carreau.

Tout point de la grille appartient à un triangle de Delaunay. La valeur du point que l'on cherche à estimer sera celle du sommet le plus proche.

La validation croisée (cf. chap. IV § 2.2) pour cette méthode rend 82,43% de bonne réponse. 5.3 Interpolation suivant le nombre d'individus comptés

Les méthodes suivantes sont appelées abusivement des méthodes d'interpolation. Elles sont plus précisément des méthodes de lissage. Elles permettent de réaliser une transition entre les régions à valeurs observées différentes. Elles sont normalement utilisées lorsque les sites d'observations sont disposés selon une grille régulière.

Les géologues ont aussi relevé pour certaines bennes remontées, la quantité de crépidule. Le tableau ci-dessous donne quelques statistiques sur ces données.

On remarque que la moyenne est de 29 crépidules par bennes et que 50% des bennes à présence de crépidules remontent plus de 25 crépidules. D'un autre coté, l'amplitude de 83 crépidules est élevée. Il en découle un écart type important.

La figure 8 présente les résultats obtenus suivant quatre méthodes de lissage.

- L'interpolation linéaire (cf. chap. I) se base sur la triangulation de Delaunay. Nous rappelons que l'estimateur pour un point S inclus dans un triangle S1, S2, S3 s'écrit :

- L'interpolation cubique consiste à ajuster à l'intérieur de chaque triangle de Delaunay une surface dont l'équation est un polynôme du troisième rang.

- La méthode du plus proche voisin estime le site S en lui affectant la valeur du site observé le plus proche.

- L'interpolation par splines ne se base pas comme les précédentes sur une triangulation de la

Remarque : les méthodes cubiques et splines produisent des surfaces à aspect très lisse alors que les méthodes linéaires et plus proche voisin présentent des discontinuités.

	Linéaire	Plus Proche Voisin	Cubique	Spline
% de bonnes réponses*	15,78	29,82	14,03	12,28
Moyenne de l'erreur d'estimation	-0,58	-0,54	0,07	-1,07
Variance de l'erreur d'estimation	14,28	17,96	14,27	15,70
Erreur quadratique moyenne	200,33	317,49	199,55	243,51
Moyenne de l'erreur relative	-0,01	-0,01	0,001	-0,02
Coefficient de corrélation	0,728	0,661	0,732	0,723

Le tableau ci-dessus permet de comparer les quatre modèles à partir de critères statistiques
(cf. chap. 4 § 2.2). La moyenne de l'erreur d'estimation des quatre modèles est proche de 0.

Cela montre qu'ils sont sans biais ( E ( Z^à ( S ) ) = Z(S ) ). Par contre ceux ne sont pas des

estimateurs robustes (variance de l'erreur d'estimation élevée). L'avantage d'un estimateur
robuste est qu'il est insensible aux valeurs aberrantes. La moyenne de l'erreur relative,

Si l'on compare les méthodes d'interpolation, la méthode cubique rend les meilleurs résultats globaux. Mais si l'on y regarde de plus près (figure 7.9) on s'aperçoit qu'elle n'estime pas correctement les sites sans présence de crépidules. Ils sont mieux estimés par la méthode du plus proche voisin (PPV).

L'idée qui s'en suit est d'estimer les sites sans présence de crépidules par la méthode du PPV. Les sites restant seront estimés par la méthode cubique. Le résultat se trouve en figure 7.10.

Nous venons de voir quatre méthodes d'interpolations et il en existe encore un grand nombre. Il est difficile de préconiser une méthode plutôt qu'une autre. C'est pourquoi il est conseillé d'utiliser systématiquement la validation croisée qui donne les moyens de comparer les diverses méthodes. Les techniques déterministes ont l'avantage de pouvoir être automatisé, mais à l'inverse du krigeage (méthode stochastique), elles ne prennent pas en compte la structure spatiale du phénomène étudié. Leur but est avant tout d'aboutir à une carte esthétique des valeurs estimées en produisant des surfaces à aspect très lisse.

Conclusion

L'objectif de ce travail a été de créer une carte de la bathymétrie au toit du socle rocheux dans la partie orientale du golfe du Morbihan. A partir d'un semis de 461 points observation, il était demandé d'estimer la valeur (profondeur de la roche sous la couche sédimentaire) des points d'une grille de 3611 points. Une recherche bibliographique a permis d'étudier les bases des techniques d'interpolations spatiales et notamment celle du krigeage ordinaire. Cette dernière a été présentée d'un point de vue théorique et appliquée. Les fondements mathématiques du krigeage et de son étape préalable, l'analyse variographique, ont été examinés et discutés. Le choix de cette méthode a été motivé par le manque de données disponibles de la variable à interpoler. Le cokrigeage (généralisation multivariable du krigeage) est une des rares techniques d'interpolation qui permette d'intégrer des variables auxiliaires dans la modélisation. Nous avons pu ainsi profiter de l'information apportée par le modèle numérique de terrain, fait par le SHOM, pour estimer plus précisément la profondeur de la roche sous la couche sédimentaire. Toute fois, les estimations produites auraient pu être encore améliorées par d'autres variables auxiliaires comme la bathymétrie au début du XX^ème siècle ou la force du courant. Mais à l'heure actuelle, ces données n'existent pas dans un format numérique exploitable.

Cette application concrète a révélé l'importance de l'analyse variographique. Elle fait la force et la faiblesse du krigeage. Elle permet en effet à l'instar des méthodes déterministes, de modéliser la structure spatiale de la variable régionalisée (anisotropie, dépendance et corrélation spatiale entre les observations) ; mais elle ne peut pas être automatisée. Un mauvais ajustement du variogramme ne pourra jamais donner de bonnes estimations et entraîne par ailleurs une mauvaise approximation des erreurs. C'est un point délicat car l'avantage du krigeage par rapport à d'autres méthodes, est qu'il permet justement de mesurer la précision de l'estimation fournie.

Les méthodes déterministes, comme celles d'interpolation par partitionnement de l'espace, sont automatisées et peuvent par conséquent être facilement intégrées aux logiciels de système d'information géographique (SIG). Certaines permettent d'interpoler des variables qualitatives (méthode du plus proche voisin, méthode de Thiesen), ce que ne propose pas le krigeage. C'est ainsi que l'on a pu tracer une carte du territoire occupé par des crépidules sur la base d'un échantillon de prélèvements localisés.

Suivant le nombre de sites d'observation, leur dispersion géographique, les modalités de la variable observée, il se peut qu'une méthode d'interpolation soit préférable à une autre. Il n'en reste pas moins que l'utilisateur ne pourra pas faire l'économie de tester plusieurs méthodes d'interpolation et utiliser la validation croisée pour valider ses résultats.

Ainsi au fil des chapitres, des méthodes déterministes et une méthode stochastique, notamment géostatistique ont été approfondi, illustré et comparé. C'est ainsi que les objectifs du stage, à savoir proposer des méthodes pour créer un modèle numérique de terrain à partir d'un échantillon de données, ont été atteints.

Il aurait aussi été intéressant, d'étudier le suivi spatio-temporel de la couverture sédimentaire, ou encore analyser la propagation des crépidules dans le temps. Des cartes à différentes époques permettraient de mieux comprendre le déplacement des sédiments ou crépidules à l'échelle du golfe et de localiser les zones critiques ainsi que leurs origines. Malheureusement les données en notre possession ne sont pas assez détaillées pour profiter de ces techniques.

Bibliographie

Michel ARNAUD & Xavier EMERY (2000) - Estimation et interpolation spatiale - Hermes

Peter A. BURROUGH & Rachael A. McDONNELLl (1998) - Principles of geographical information systems - 2cd edition - Clarendon Press

Noel A.C.CRESSIE (1993) - Statistics for spatial data - Revised Edition - Wiley Hubert JAYET (1993) - Analyse spatiale quantitative - Economica

A.G. JOURNEL & Ch.J. HUIJBREGTS (1993) - Mining Geostatistics - Academic Press George MATHERON (1963) - Traité de géostatistique appliquée - Technip

Laurent BERTINO et Hans WACKERNAGEL (2004) - Prévoir l'état de l'océan : L'assimilation de données intègre observations et modèles numériques - http://www.larecherche.fr/special/web/wxyz371.pdf

Nicolas JEANNEE (2005) - La Géostatistique : concepts et applications pour l'estimation de phénomènes spatialisés -

Denis MARCOTTE - Cours de géostatistique minière & Cokriging with MATLAB - http://geo.polymtl.ca/