Skip to navigation – Site map

HomeSectionsSystems, Modelling, Geostatistics2004Modélisations prospectives de l’o...

2004
295

Modélisations prospectives de l’occupation du sol. Le cas d’une montagne méditerranéenne

Prospective modelling of land cover in Mediterranean mountain areas
Martin Paegelow, Nathalie Villa, Laurence Cornez, Frédéric Ferraty, Louis Ferré and Pascal Sarda

Abstracts

The authors apply three methods of prospective modelling to high resolution georeferenced land cover data in a Mediterranean mountain area: GIS approach, generalized linear model and neural networks. Land cover prediction to the latest known date is used to validate the models. In the frame of spatial-temporal dynamics in open systems, results are encouraging and comparable. Correct prediction scores average around 73 %. The results analysis focuses on geographic location, land cover categories and parametric distance to reality of the residuals. Comparing the three models shows the high degree of convergence and a relative similitude of the results obtained by the two statistical approaches compared to the GIS supervised model. Further steps are the application of the models to other test areas and the identification of respective advantages to develop an integrated model.

Top of page

Full text

Problématique et objectifs

1L’objet de notre recherche est la modélisation de dynamiques environnementales dans le cadre de systèmes complexes et ouverts. Dans ce cadre, la variable étudiée – ainsi que les variables d’environnement, susceptibles d’expliquer son évolution dans l’espace et dans le temps – contient une part d’incertitude ou d’aléa ce qui exclut, de fait, une approche déterministe. Ainsi, nous utilisons une approche stochastique (probabiliste) tenant compte de la dépendance dans le temps (effet mémoire) et dans l’espace : les outils probabilistes utilisés sont notamment la distribution multinomiale et l’analyse de Markov. En outre, notre approche fait également appel à la logique floue et à un automate cellulaire. L’ensemble de ces méthodes est mis en œuvre dans trois modèles à but prévisionnel différents afin de comparer leurs performances respectives. Plus précisément, il s’agit de comparer un modèle géomatique combiné de simulation prospective dont la mise en œuvre est possible en utilisant les fonctions de logiciels SIG disponibles sur le marché à deux modèles statistiques dont la mise en œuvre, plus longue, est extérieure au SIG. En contrepartie, l’intérêt des deux approches statistiques réside dans leur caractère automatique tandis que le modèle SIG nécessite une analyse thématique experte.

2Notre choix en matière de modélisation statistique a porté sur deux approches classiques, l’une basée sur le maximum de vraisemblance (modèle linéaire généralisé) et l’autre utilisant un réseau de neurones. Ces deux méthodes sont proches du point de vue de la modélisation et diffèrent essentiellement en ce qui concerne les algorithmes de mise en œuvre.

3Un des défis actuels de la recherche en géomatique est celui de la modélisation prospective de données géographiques à haute résolution. Des méthodes géostatistiques éprouvées pour l’interpolation et l’extrapolation spatiales existent depuis plusieurs décennies et sont implémentées dans nombre de logiciels géomatiques commercialisés. Par contre, des outils de modélisation temporelle et d’aide à la décision ne furent implémentés dans les SIG que récemment et doivent être considérés plutôt comme des algorithmes expérimentaux intéressants que des techniques opérationnelles.

4Depuis les années 1990, la demande sociale en outils d’aide à la décision et de modélisation capables d’assister différentes tâches de gestion environnementale (notamment la prévention de risques) et d’aménagement des territoires s’est fortement accrue.

5Cet article illustre les premiers résultats d’une étude comparative de trois méthodes de modélisation prospective appliquée à l’occupation du sol dans des anthroposystèmes montagnards de l’Europe du sud. Nous considérons l’occupation du sol comme un indicateur pertinent, disponible à haute résolution, d’une combinaison d’activités humaines que les sociétés déploient dans l’espace – et auxquelles l’occupation du sol réagit avec une certaine inertie – et de facteurs naturels.

6Les montagnes méditerranéennes font l’objet d’une profonde restructuration socio-économique qui se manifeste, entre autres, dans de spectaculaires changements paysagers. Cette réorganisation commença dans les Garrotxes (Pyrénées Orientales, zone d’études) à la fin de la première moitié du 19ème siècle par le déclin du système agropastoral traditionnel provoquant l’exode rural.

7Une base de données géoréférencées matérialise les connaissances des dynamiques passées et actuelles de l’occupation du sol ainsi que des facteurs d’environnement potentiellement explicatifs. Elle alimente trois méthodes de modélisation prospective : l’une, supervisée, met en œuvre des algorithmes implémentés dans des logiciels SIG; les deux autres approches –modèle linéaire généralisé et réseau neuronal – peuvent être qualifiées de non supervisées dans la mesure où les règles du comportement spatio-temporel sont automatiquement détectées par l’outil. L’objectif principal étant la mise en œuvre et l’optimisation de chacune des trois approches sur le même jeu de données. L’interprétation critique des résultats, notamment des résidus de la prédiction, permet de cerner avantages et inconvénients respectifs. A partir de cette analyse comparative, peuvent être envisagées la possibilité de construction d’un modèle optimisé intégrant les points forts de chacune des méthodes ainsi que les modalités de transposition et les limites de généralisation.

  • 1 GEODE UMR 5602 CNRS, Groupe SMASH –EA 3686 GRIMM, UTM et Instituto de Désarrollo Régional – Univers (...)

8Afin de valider et d’optimiser les modèles ceux-ci sont appliqués, dans un premier temps, à prédire l’occupation du sol à la dernière date connue avant de proposer des scénarii prospectifs. Cette étude est menée dans le cadre d’une coopération entre trois équipes de recherche1 travaillant sur deux sites : les Garrotxes dont nous présentons ici les premiers résultats et la Alta Alpujarra Granadina (Andalousie, Espagne – travaux en cours).

Zone d’étude et base de données

Les Garrotxes

9Les Garrotxes, situées dans le département des Pyrénées Orientales, à l’extrémité NO du Conflent forment un ensemble géographique constitué de cinq communes et d’une superficie de 8 570 ha. Ce bassin versant présente une rive droite granitique, à modelé géomorphologique relativement lourd, où sont localisées la quasi-totalité des anciennes terrasses de cultures et des forêts de pins à crochet (Pinus uncinata) et de pins sylvestre (Pinus sylvestris) ; un espace à dynamique végétale très rapide. La rive gauche du Cabrils, cours d’eau collecteur se jetant dans la Têt à Olette, est une large soulane orientée SO sur substrat schisteux avec un métamorphisme de contact dans les zones les plus basses et occupée par des landes majoritairement ligneuses (à base de Genista purgans et, dans une moindre mesure, de Calluna vulgaris), fortement embroussaillée aux altitudes les plus basses par des chênes verts (Quercus ilex). La particularité des Garrotxes est leur enclavement : le bassin versant, à l’écart des grandes routes, est délimité au nord par le massif du Madrès (2469 m), à l’ouest (Cami Ramader) et au sud (Puig de la Tossa, Serrat del Cortal) par des chaînes culminant entre 1600 et 2000 m d’altitude et à l’est par la crête (Lloumet) de la soulane rejoignant le Madrès. La vallée du Cabrils présente une dégradation progressive du climat méditerranéen ; la remontée de l’influence méditerranéenne au cœur des Pyrénées Orientales étant assurée par la vallée de la Têt modifiant ainsi la rudesse du climat montagnard.

Image 10000000000001FE00000131CD054FCA01BF6F56.jpg

Fig. 1 : Localisation des Garrotxes à l’intérieur du département des Pyrénées Orientales

10Autrefois modèle d’organisation agropastorale traditionnelle, de nos jours l’agriculture a quasiment disparu tandis que l’activité pastorale, longtemps en déclin, donne des signes d’un renouveau suite à une profonde réorganisation entamée durant les années 1980 ([Met04], [Pae04b]).

11Le maximum démographique au début du 19ème siècle se traduisait par une mise en valeur de toutes les ressources montagnardes traditionnelles mobilisables (agriculture, élevage, sylviculture). Ainsi en 1826 (cadastre napoléonien) un quart de la surface totale était cultivé. Le support quasi exclusif de l’agriculture a été les terrasses de cultures (feixtes). Le déclin démographique (de 1 832 habitants en 1830 à 90 en 1999) et la reconversion des terrasses de culture en pâturages, broussailles et forêts allaient de pair.

12Parmi les agents externes considérés responsables du déclin de cette société locale à faible degré d’insertion dans l’économie nationale on peut citer, outre les processus d’industrialisation et de mise en valeur agricole de plaines au cours du 19ème siècle, une variabilité interannuelle accrue des précipitations, observée au milieu du 19ème siècle ([Tab03]), qui eût pu contribuer à la rupture d’un système poussé à bout par la pression anthropique sur le milieu. Deux évènements ponctuels – l’arrivée du chemin de fer à Olette (1911) et la Première Guerre Mondiale – ont accéléré l’exode rural.

13Ainsi il est probable que l’avenir proche se jouera en termes de gestion - ou de non gestion - pastorale se matérialisant par divers moyens de blocage, voire d’inversion, de l’embroussaillement et du reboisement spontané des espaces pastoraux (berger guidant le bétail, clôtures, écobuage). L’instauration de groupements pastoraux (GP) et d’associations foncières et pastorales (AFP) à partir des années 1980 a effectivement conduit à une reprise de l’activité pastorale avec un remplacement partiel du cheptel ovin par des bovins et des équins. Les signes de reconversion économique sont récents (années 1990) mais d’une portée limitée (ouverture d’un gîte d’étape à Sansa, tentatives de valorisation en tourisme vert) malgré la concrétisation prévue prochainement du Projet de Parc Naturel Régional des Pyrénées Catalanes.

La base de données et l’évolution de l’occupation du sol

14La base de données géoréférencées consiste en une série de cartes d’occupation du sol assorties de plans d’information représentant des facteurs environnementaux et sociaux. Les cartes existent – selon les traitements envisagés – soit en mode image (résolution du pixel d’environ 18 m), soit en mode objet. Ainsi les principaux traitements pour la modélisation font appel à une logique image (analyse spatiale) tandis que le mode objet offre une plus grande souplesse pour réaliser des requêtes attributaires.

15Les cartes d’occupation du sol ont la même résolution spatiale mais ont des origines et des légendes variables. La première carte est basée sur le cadastre napoléonien (1826) – un support permettant de distinguer entre forêts, espaces pastoraux, prairies, champs agricoles et le bâti (villages). La première mission aérienne disponible (1942) rend possible le renseignement de la catégorie broussailles (landes très embroussaillées contenant de groupes d'arbres ou un nombre important d'arbres isolés) – maillon manquant en 1826 entre les formations arborescentes denses et arbustives (landes ligneuses). La carte de 1962 conserve la même légende tandis que l’échelle et la qualité des missions aériennes plus récentes (1980 et 1989), également panchromatiques, facilitent la distinction entre forêts de conifères et formations boisées de type feuillus. Il en va de même pour une meilleure discrimination des espaces pastoraux : landes ligneuses (notamment à base de Genista purgans) et landes à graminées. La carte d’occupation du sol la plus récente (2000) est basée sur des observations de terrain et est, par conséquent, nettement plus détaillée (20 catégories).

16En raison de la nature des sources, la classification de l'occupation du sol, sous forme de trois nomenclatures emboîtées, est surtout d’ordre physionomique.

Image 100000000000025D00000145B76B37D7798C496C.jpg

17Figure 2 : Répartition de l’occupation du sol entre 1826 et 2000 (Garrotxes)

18L’évolution de l’occupation du sol (cf. fig. 2) est classique. Les terres labourées délaissées ont été d’abord utilisées comme zones de pâture avant un embroussaillement menant souvent à la reconquête par la forêt.

19La base de données géoréférencées contient nombre de plans d’information ayant trait à l’occupation du sol :

  • Plans issus du MNT : carte altitudinale, carte des pentes, carte d’occupation du sol ;

  • Plans d’accessibilité calculés à partir du réseau routier et l’emplacement des villages (habitat groupé) : indice d’accessibilité (analyse de distance-coût) selon la date ;

  • Plans relatifs à la gestion pastorale : unités pastorales (UP), associations foncières et pastorales (AFP), pression pastorale ;

  • Plans tenant compte du statut particulier de certains zones : forêts domaniales et zone militaire ;

  • Limites administratives ;

  • Réseau hydrographique.

20Les données purement attributaires (recensement de la population, recensement général agricole, …) sont, selon leur degré de confidentialité, connus, mais apportent peu de connaissances compte tenu de leur unité spatiale de rattachement (la commune) incompatible avec une analyse haute résolution.

Méthodologie et mise en œuvre

21Il s’agit de construire trois modèles prédictifs de l’occupation du sol (variable qualitative) à haute résolution, de comparer leurs efficacités respectives sur la dernière date connue en utilisant la même base de données. La première approche fait appel aux techniques disponibles dans les SIG et peut être qualifiée de supervisée dans la mesure où l’analyse du géographe conduit à l’implémentation des règles nécessaire au calcul des cartes de probabilité. A l’opposé pour les deux autres approches statistiques, le rôle du mathématicien, peu familier avec la thématique, consiste plutôt à optimiser l’algorithme que d’y introduire les conclusions de sa propre analyse du comportement spatio-temporel du milieu considéré.

Approche supervisée par SIG

22Le modèle que nous présentons se veut être simple à deux égards : simple dans le lever des données d’entrée (limitation à quelques données facilement disponibles, cf. paragraphe 2.2.) et simple dans sa mise en œuvre informatique (recours à des algorithmes implémentés dans un logiciel SIG commercialisé). Ce modèle géomatique combiné :

  • Fait appel à la logique floue afin d’ajuster les données environnementales dans l’évaluation multicritère.

  • Est stochastique pour l’aspect prédictif de la simulation à événements discrets et états finis – chaînes de Markov d’ordre 2 (deux dates antérieures).

  • Remédie aux limites de l’analyse markovienne en recourant à une évaluation multicritère optimisant l’affectation spatiale des probabilités markoviennes (prise en compte de la rugosité de l’espace) par constitution d’une base de connaissances et de règles d’inférence (variables d’environnement) relatives à la phase d’apprentissage (calibration) du modèle.

  • Utilise un automate cellulaire simple pour favoriser l’émergence de zones de probabilités d’états à extension surfacique réaliste.

23Mise en œuvre sous le logiciel Idrisi 32, la modélisation se découpe en trois phases :

  • La constitution de la base de connaissances de la dynamique spatio-temporelle de l’occupation du sol par évaluation multicritère (EMC) des variables d’environnement. Les variables d’environnement d’origine sont transformées, pour chacun des types d’occupation du sol, par traitements statistiques et par logique floue en plans de probabilité d’occurrence de chacune des catégories d’occupation du sol. Ces plans de probabilité résultants, se basant sur la période d’apprentissage (1980 à 1989) servent à l’allocation spatiale des probabilités de transition.

  • Le calcul des probabilités de transition par analyse de chaînes de Markov (ACM) entre les dates de la phase d’apprentissage et la date simulée (2000 – dernière date connue).

  • L’allocation spatiale des probabilités de transition markoviennes : cette dernière étape utilise les résultats catégoriels de l’EMC. Ceux-ci sont intégrés, par évaluation multiobjectif (EMO), en une seule carte de l’occupation du sol simulée laquelle est traitée par un automate cellulaire (AC) basé sur un filtre de contiguïté spatiale.

24Le calibrage du modèle est obtenue en modélisant l’état de l’occupation du sol à la dernière date connue (2000) sur la base de l’information provenant d’une période d’apprentissage englobant les deux dates précédentes (1980 et 1989). Bien que nous disposons d’une connaissance historique plus approfondie, il est évident que les conditions socio-économiques actuelles ne s’appliquent pas aux états de l’occupation du sol du 19ème et du 20ème siècle jusque dans les années 1960.

Construction d’une base de connaissances des dynamiques de l’occupation du sol

25La connaissance des dynamiques récentes est essentielle pour appréhender l’évolution future et sa modélisation. Nous entendons par connaissance des mesures statistiquement significatives du comportement spatial et temporel de l’occupation du sol en relation avec des critères environnementaux considérés explicatifs d’une partie de sa variabilité. Dans une évaluation multicritère (EMC), on distingue entre des critères binaires, les contraintes, et les critères ayant une aptitude variable dans l’espace, les facteurs. Les contraintes booléennes masquent certaines zones (occurrence possible ou non) ; ils peuvent s’appliquer à toutes formes d’occupation du sol (exclusion des espaces bâtis) ou être spécifique à certaines formations (limite altitudinale des conifères). Les facteurs traduisent une connaissance graduée pour l’objectif en question (une forme d’occupation du sol) – ils peuvent être pondérés et on peut définir leur degré de compensation.

26Pour chaque catégorie d’occupation du sol la connaissance de son comportement spatio-temporel provient d’une analyse diachronique des dynamiques et de la friction géographique en comparant la répartition théorique (espace homogène) et la répartition réelle (niveaux de probabilité 99% et 99.9%). La rugosité géographique est exprimée par les facteurs d’environnement cartographiés (altitude, pente, exposition, accessibilité, proximité aux entités de même nature, statut de gestion particulière de certaines zones et probabilité de changement) et disponibles à la même résolution spatiale.

27Les facteurs sont standardisés par recodage manuel ou par l’emploi de fonctions fuzzy et pondérés par l’emploi de la matrice de Saaty ([Saa77]) qui renvoie le vecteur propre de chaque facteur.

28L’approche de l’EMC développée ([Eas93]) inclut des poids d’ordre (ordered weighted averaging – OWA) permettant le choix du niveau de risque et de compensation entre facteurs. Ces poids d’ordre classent les aptitudes par rang croissant et aboutissent à un classement spécifique à chaque pixel. Nous avons opté pour une approche conservatrice (peu de risques et niveau de compensation limité) comme l’exprime la fig. 3.

Image 1000000000000148000000CBBFEAA6D3F912093C.jpg

Figure 3 : Espace de décision et approche EMC-OWA choisie.

Calcul des probabilités de transition

29Le calcul prédictif de l’occupation du sol est opéré par une analyse des chaînes de Markov, un processus discret avec des pas temporels discrets et dont les valeurs à la date prédite dépendent des valeurs à des dates antérieures. La prédiction est exprimée par une estimation des probabilités de transition.

30Le test consiste à prédire l’occupation du sol en 2000 (dernière date connue) sur la base de 1980 et de 1989. Le résultat se présente sous forme d’une matrice dans laquelle sont codées les probabilités de changement de chaque catégorie d’occupation du sol ainsi que le nombre de pixels affectées entre la dernière date d’apprentissage et la date projetée. La fonction calcule également une carte de probabilité conditionnelle pour chaque catégorie d’occupation du sol indiquant la probabilité markovienne par pixel de la modalité en question à la date projetée. Une intégration stochastique de l’ensemble des cartes par modalité en une seule est possible mais aboutit à un résultat relativement éloigné de la réalité (image bruitée) car cette procédure purement statistique ne tient pas compte ni des règles de connaissances établies par EMC, ni de la continuité spatiale.

Allocation spatiale des probabilités prédites

31L’allocation spatiale des probabilités markoviennes intègre la connaissance sur la répartition probable de l’occupation du sol (EMC), une évaluation multiobjectif (EMO) tenant compte des objectifs concurrents (chaque modalité d’occupation du sol étant un objectif) et un automate cellulaire, basé sur un filtre de contiguïté spatiale. La fonction utilisée sous Idrisi est CA_Markov dont l’algorithme est itératif afin de tenir compte des distances temporelles entre les deux dates d’apprentissage et la dernière date d’apprentissage et la date de projection. Elle donne en sortie une carte prospective de l’occupation du sol probable.

32La fig. 4 résume les principales étapes de modélisation par SIG ([Pae03], [Pae04a]).

Image 1000000000000283000001AF35DC7A6DBE90C289.png

Figure 4 : Modélisation de l’occupation du sol par SIG. Aperçu des principales fonctions et de leur enchaînement

Approche par réseaux de neurones

33L’idée de l’utilisation de réseaux de neurones a été principalement motivée par leurs remarquables capacités d’adaptation et de souplesse face à un très grand nombre de problèmes, notamment lorsque ceux-ci présentent des aspects non linéaires ou lorsque les variables explicatives sont fortement corrélées ; ces deux aspects se retrouvent dans le travail de prédiction que nous cherchons à effectuer. Aussi, les réseaux de neurones ont récemment connu une grande popularité et ont très favorablement concurrencé les méthodes statistiques classiques. On les retrouve notamment dans la prédiction de séries chronologiques (cf. [Bis95], [Lai01] et [Par00]). Le cadre dans lequel nous sommes amenés à travailler est encore plus étendu puisqu’il s’agit ici d’un processus spatio-temporel auquel s’ajoutent des variables explicatives comme nous le détaillerons plus loin.

Réseaux de neurones multi-couches (perceptrons)

34Nous avons travaillé avec une classe particulière de réseaux de neurones, les réseaux multi-couches ou perceptrons. Ceux-ci ont été les premiers à connaître un essor important ; leur création est issue des premières tentatives de modélisation des principes de base régissant le fonctionnement du cerveau même si leur champ d’application s’est, depuis, considérablement élargi, notamment au traitement de données statistiques (pour plus de détails, consulter [Dav69] ).

35Lorsque l’on parle de réseaux à couches, le nombre de couches est à définir par l’utilisateur mais doit comprendre au minimum une couche d’entrée et une couche de sortie ; les autres couches dont le nombre varie sont appelées couches cachées ; pour leurs remarquables propriétés d’approximation, nous avons choisi d’utiliser un réseau à une couche cachée dont l’architecture générale est celle de la figure 5.

Image 10000000000001D100000116D3888A63AD9C77F0.png

Figure 5 : Architecture d’un réseau à une couche cachée

36Détaillons un peu le fonctionnement de ce réseau : en entrée, la valeur des neurones est celle des variables explicatives du modèle ; chacune de ces valeurs numériques est multipliée par un certain nombre de poids pour être, finalement, additionnée et transformée par une fonction de lien au niveau des neurones de la couche cachée. Enfin, les valeurs numériques des neurones de la couche cachée subissent à leur tour une multiplication par des poids et leur addition donne la valeur des neurones de sortie qui modélisent la variable expliquée. Les poids, généralement notés w, sont choisis lors d’une phase dite d’apprentissage sur un jeu de données test et minimisent l’erreur quadratique de ce jeu de données. Finalement, les réseaux de neurones à une couche cachée sont les fonctions de la forme :

Image 10000000000000C2000000311B8B5266D9E90596.png

x est le vecteur des variables explicatives du modèle, q2 le nombre de neurones sur la couche cachée, g la fonction de lien de la couche cachée (typiquement g est la sigmoïde g : x Image 10000000000000150000001031BBC8CAF4A3DABA.pngImage 1000000000000030000000285DE549678AAD0C0C.png), w (1) sont les poids entre la couche d’entrée et la couche cachée et w (2) les poids entre la couche cachée et la couche de sortie.

37L’intérêt de ce type de réseau est expliqué par le résultat suivant ([Hor93]) : les réseaux de neurones à une couche cachée permettent d’approcher, avec la précision souhaitée, n’importe quelle fonction continue (ou d’autres fonctions qui ne sont pas nécessairement continues) : c’est ce que l’on appelle la capacité d’approximateur universel et c’est aussi ce qui leur permet de s’appliquer avec une grande efficacité à un grand nombre de modèles.

Modèle pour les données de Garrotxes

38Dans l’exemple des données de Garrotxes, plusieurs facteurs ont été retenus comme pouvant influencer l’occupation du sol d’un pixel donné à la date t :

  • pour l’aspect temporel (processus d’ordre 1) : la valeur du pixel considéré à la date précédente (t – 1) que l’on exprime sous forme disjonctive (par exemple, si l’on dispose de 8 catégories d’occupation du sol, la première sera codée sous la forme d’un vecteur à 8 coordonnées : (1 0 0 0 0 0 0 0), la seconde : (0 1 0 0 0 0 0 0), etc) ;

  • pour l’aspect spatial : la fréquence de chaque type d’occupation du sol dans le voisinage du pixel considéré à la date précédente (t – 1). Se pose alors le problème du choix du voisinage (taille et forme) : pour la forme, diverses possibilités s’offrent à nous, de la plus simple (voisinage carré ou en étoile) à des voisinages plus sophistiquées (voisinage suivant la pente pour mieux tenir compte des influences morphologiques du terrain). Quant à la taille du voisinage, il s’agira de déterminer jusqu’à quelle distance un pixel est susceptible d’influencer le pixel considéré. Afin de respecter la spatialisation de la carte, on pondèrera l’influence d’un pixel par une fonction décroissante de la distance au pixel considéré (cf. fig. 6) ;

Image 10000000000001DD0000011F7E6D31238AB794C4.png

Figure 6 : Un exemple de voisinage

  • des variables environnementales (pente, altitude, …).

Mise en œuvre

39A l’issue d’une phase exploratoire qui nous a permis de cerner les variables pertinentes et divers paramètres du modèle comme la forme du voisinage (que, dans un soucis de simplicité, nous avons choisi carré), sa taille (finalement fixée à 3 ; cf. fig. 6) ou le nombre de neurones optimal sur la couche cachée, l’architecture choisie compte, en entrée, 19 neurones :

  • pour l’aspect temporel : 7 neurones pour le codage disjonctif de la valeur du pixel (le bâti, constant, ayant été retiré de l’étude) à la date précédente (t – 1) ;

  • pour l’aspect spatial : 8 neurones pour la fréquence des divers types d’occupation du sol dans le voisinage (fréquence pondérée par une fonction décroissante de la distance) ;

  • 4 neurones pour les variables d’environnement pente, altitude, exposition et distance aux infrastructures (préalablement centrées et réduites).

40Le réseau dispose aussi de 8 neurones sur la couche cachée et de 7 neurones en sortie, chacun estimant la probabilité d’appartenance du pixel à un type d’occupation du sol (hors bâti).

41Pour la phase d’apprentissage, nous avons utilisé comme jeu de données la carte de 1980 avec comme cible à prévoir la carte de 1989 et la carte de 1989 avec comme cible à prévoir la carte de 2000. Après avoir repéré de larges zones dans lesquelles l’occupation du sol était stable, nous avons considéré uniquement les pixels dont un voisin au moins avait une valeur d’occupation du sol différente (exploitant ainsi pleinement la spatialisation des données) : ces pixels seront, dans la suite, nommés pixels frontières. Les autres pixels ont été considérés comme stables dans un intervalle de temps de 10 ans.

42D’un point de vue calculatoire, les programmes ont été réalisés à l’aide du logiciel Matlab (cf [Bea98] ) et sont disponibles sur demande.

Approche par modèle linéaire généralisé

43Le modèle de régression logistique est un modèle linéaire généralisé (et donc paramétrique) dans lequel la variable réponse est qualitative et qui permet d'obtenir une prédiction de celle-ci en tenant compte d'un ensemble d'informations issues de variables explicatives. Lorsque la réponse possède plus de deux modalités, on parle de modèle de régression logistique multiple ou modèle de régression polychotomique ([Hos89]). D'autres développements plus récents concernant ce modèle ont été réalisés par Kooperberg et al. [Koo97]. Ce type de modèle logistique est particulièrement bien adapté à notre problématique puisqu'il s'agit de prédire pour chaque pixel de la carte un type d’occupation du sol (variable réponse ayant 8 modalités codées par un entier v, v = 1, 2, …, 8). La spécificité de notre étude vient du fait que, outre la nature topographique du pixel (pente, altitude, ...), le modèle choisi doit tenir compte d'un effet spatial (état de la végétation dans l'environnement du pixel) et d'un effet temporel (évolution du type d’occupation du sol du pixel et de son environnement). En ce sens il s'agit d'adapter le modèle de régression logistique à notre cadre, un des enjeux les plus importants étant le choix de la forme et de la taille de l'environnement pris en compte par le modèle.

44De façon générale, le modèle de régression logistique permet de modéliser, en fonction d'un certain nombre de paramètres, la probabilité pour que le type d’occupation du sol d'un pixel au temps t (c'est-à-dire la variable réponse) soit égal à un des 8 catégories d’occupation du sol. Il s'agit donc d'estimer les paramètres inconnus du modèle, et ensuite les probabilités a posteriori de type d’occupation du sol sachant les valeurs des différentes variables explicatives. On utilise ensuite une règle de type bayésien consistant à affecter au temps t à un pixel donné l'indice de végétation ayant la plus forte probabilité a posteriori.

Régression logistique multiple spatio-temporelle

45Indexons par i = 1, 2, …, N les pixels de la carte d’occupation du sol et notons Ii l'ensemble des informations dont on dispose concernant le pixel n° i. D'un point de vue formel, le modèle de régression logistique multiple que nous adoptons peut se présenter sous la forme générale suivante :

Image 1000000000000210000000422A98DB5BBE89C70E.png

46αv est un paramètre associé au type d’occupation du sol v que l'on souhaite prédire pour le ième pixel et 𝛾v,Ii un ensemble de paramètres lies à v ainsi qu'aux informations concernant toujours ce pixel n° i. Ainsi, le nombre total de paramètres mis en jeu dans ce modèle dépend uniquement du nombre de types d’occupation du sol et du nombre de variables explicatives. Dans l'expression (1) Prob(pixeli = v‌│Ii) représente la probabilité que l’occupation du sol du pixel i soit du type v lorsque les variables explicatives prennent les valeurs décrites par l'ensemble Ii. Notons que l'expression (1) modélise le rapport (son logarithme) de la probabilité qu'un pixel prenne la modalité v sur la probabilité que ce pixel prenne la modalité codée 8 ce qui permet d'intégrer la contrainte que la somme des huit probabilités est égale à 1.

47Dans l'expression (1), nous devons intégrer l'effet temporel : celui-ci est pris en compte en faisant dépendre le type d’occupation du sol du pixel i du temps t c'est-à-dire en posant pixeli = pixeli(t). Par ailleurs l'information (ou plus exactement une partie de cette information) dépend du temps t – 1 : Ii = Ii(t – 1). L'idée consiste donc à calculer la probabilité qu'un pixel prenne un type d’occupation du sol v à l'instant t en fonction de l'information que l'on possède sur ce même pixel à l'instant précédent t – 1; on répète cette procédure pour tous les pixels de la carte. Connaissant les cartes 1980 (t – 1) et 1989 (t), on peut estimer l'ensemble des paramètres de notre modèle de sorte à ajuster au mieux la carte 1989. Il suffit alors d'incrémenter le temps dans notre modèle pour prédire la carte à l'instant futur t + 1 (2000) à partir de la carte observée à l'instant t (1989).

48Quant à l'effet spatial, il est pris en compte de la même façon que dans l'approche par réseau de neurones. Il est en effet naturel de penser que l’évolution de l’occupation du sol du pixel i dépend de celle des pixels environnants. Pour cela on considère un voisinage carré Vi autour du pixel n° i que l'on souhaite prédire et on extrait comme information du voisinage Vi le nombre de pixels prenant le type n°1 d’occupation du sol, le type n°2 d’occupation du sol, ... Cette façon de procéder revient à supposer une invariance isotrope, c'est-à-dire que le type d’occupation du sol autour du pixel i ne dépend pas de la direction. Dans la mise en oeuvre de la méthode, nous avons privilégié la simplicité de la forme du voisinage (carré). Lors de développements ultérieurs, on pourrait envisager d'autres formes que le carré (étoile, rectangle, ...) et voir s'il en résulte un gain ou pas en terme de prédiction. On peut également envisager une modélisation privilégiant certaines directions c'est-à-dire rompant avec l'hypothèse d'invariance isotrope. Notons cependant qu'il en résulterait un modèle avec un plus grand nombre de paramètres et que de ce point de vue on doit également composer avec la capacité à bien estimer un modèle qui serait trop complexe.

49En combinant effet temporel et effet spatial, on est finalement amené à considérer le modèle suivant

Image 10000000000001B50000003BD8F41CC1EEBA8FB0.png

50où Ii(t – 1) englobe l'information extraite du voisinage Vi(t – 1), c'est-à-dire tient compte de l’occupation du sol observée autour du pixel n° i à l'instant t – 1.

51Enfin, Ii(t – 1) comprend également l'information issue des variables telle que la pente ou l'altitude décrites plus haut.

Mise en oeuvre

52D'un point de vue pratique, la mise en oeuvre se décompose en deux étapes : une étape d'estimation et une étape de validation.

53Etape d'estimation

On estime les paramètres du modèle (αv et ceux contenus dans Image 1000000000000057000000194392DBA89B2FA7EF.jpg). La procédure d'estimation est basée sur la maximisation de la vraisemblance pénalisée, critère bien connu en statistique pour la stabilité des solutions obtenues. L'algorithme d'optimisation utilise est de type Newton-Raphson. Remarquons que la pénalisation introduit un nouveau paramètre, appelé paramètre de pénalisation et noté ε, qu'il faudra choisir. Comme cela a été dit précédemment, on utilise les cartes 1980 et 1989 pour estimer les paramètres, ceci pour différentes tailles de voisinage et valeurs ε.

54Etape de validation

55Il s'agit de déterminer la taille de voisinage et le paramètre de pénalisation optimaux en ce sens que ces choix fourniront une prédiction de la carte 2000 la plus proche possible de celle observée Plus précisément, à l'aide de étape précédente, il est possible de construire plusieurs prédictions de la carte 2000, chacune correspondant à différentes tailles de voisinage et valeurs ε. En comparant les cartes prédîtes pour 2000 avec la carte réelle de 2000, on repère la carte qui possède le plus petit nombre de pixels mal prédits; la taille de voisinage et le paramètre de pénalisation correspondants seront considérés comme étant optimaux.

Résultats et interprétation

56Les trois approches ont été testées par la modélisation de l’occupation du sol à la dernière date connue (2000)

Image 100000000000021A000002724E0A2E29225C20E2.jpg

Figure 7 : résultats des modélisations de l’occupation du sol en 2000 et occupation du sol réelle

57Les résultats globaux (% de surface) de ce test sont proches de la réalité (cf. fig. 7 et 8).

Image 10000000000001CF000000AE20B08FD1B6172CE4.jpg

Figure 8 : Surface en pourcentage de l’occupation du sol en 2000 réelle et modélisée

58Cependant il convient d’analyser la répartition spatiale de ces sommes de surface prédites par catégorie à haute résolution (coté du pixel environ 18 m). Le tableau de la figure 9 compare les résidus par catégorie en % de la surface réelle de chaque catégorie mise à part les cultures dont le nombre de pixels tend vers zéro. Le premier enseignement de ce tableau est la relative concordance des résultats des trois approches. On remarque également que la modélisation de modalités ayant une grande superficie (forêt de conifères, landes à genêt) est plus facile que celle des catégories d’occupation du sol de faible ampleur. Ainsi, quel que soit le modèle, moins d’un pixel sur deux a été correctement prédit pour la catégorie broussailles. Parmi les modalités de faible surface, on note cependant de grandes différences selon leur stabilité spatio-temporelle : les prairies étant plus stables que les landes à graminées, leur taux de prédiction est meilleur.

59Les taux de prédiction globale des trois méthodes sont très proches : 72.8 % (SIG), 74.3 % (réseaux neuronaux) et 72.8 % (modèle linéaire généralisé).

Image 100000000000021D0000008E694546B915CE907B.jpg

Figure 9: Pourcentage de résidus par catégorie d’occupation du sol et approche modélisatrice

60Les modélisations n’ont pas vocation à prédire la réalité mais peuvent nous aider à mieux comprendre des changements spatio-temporels environnementaux et sociaux complexes. Dans ce sens, l’interprétation des résultats des modélisations doit tenir compte des limites des modèles. La modélisation de l’occupation du sol signifie une simulation de ce que la réalité pourrait être, un scénario raisonné et quantifiable dans le contexte d’aide à la décision.

61Cependant une interprétation minutieuse des résultats devrait nous permettre à améliorer le modèle et, par conséquent, le taux de prédiction. Dans ce sens, l’analyse se focalise surtout sur les résidus.

62La catégorie d’occupation du sol la plus représentée (les conifères) obtient un très bon score de prédiction alors que les broussailles, relativement présentes sur le territoire, obtiennent un très mauvais score (plus de la moitié de mal prédits). Diverses remarques permettent d’expliquer ces phénomènes et de penser à des stratégies d’amélioration de la prédiction. Tout d’abord, on peut constater que les broussailles sont la catégorie naturellement la plus dynamique sur un territoire caractérisé par un équilibre entre espace forestier et pastoral régi notamment par la gestion pastorale. Les broussailles sont également soumises le plus à des effets aléatoires : un feu de forêt, une coupe ou bien l’abandon de pâturage transforment en l’espace de 10 ans une parcelle en broussailles ; ce sont des phénomènes complètement incontrôlables.

63Bien que l’occupation du sol soit décrite de manière qualitative, ses différentes catégories s’échelonnent entre des formations fermées (forêt de conifères, forêt de feuillus) et ouvertes (cultures). Ces rangs « paysagers » permettent de quantifier l’erreur de prédiction, exprimée en nombre de catégories (cf. fig. 10). Ainsi, quelle que soit la méthode de modélisation, pour environ la moitié des pixels mal prédits, l’erreur de prédiction n’est que d’une catégorie à une résolution spatiale élevée. Le nombre de résidus décroît fortement avec l’augmentation de l’écart entre la réalité et la projection.

Image 10000000000002220000005BCA391BF2EEBE7975.jpg

Figure 10 : Analyse des résidus de la modélisation par l’écart de catégories entre la réalité et la modélisation (données en % de la surface totale)

64Un autre aspect intéressant est la grande concordance des trois modèles (cf. fig. 11 et 12). Ainsi 66.5 % de la surface totale sont correctement prédites par chacun des modèles, 20.3 % par aucune.

Image 10000000000001FB000001BEF0A6936258F93699.png

Figure 11 : Prédictions correctes croisées par les 3 modèles

Image 1000000000000221000000BA2AC220E304C8816A.jpg

Figure 12 : Mise en relation des scores de prédiction correctes des trois modèles avec l’occupation du sol en 2000. Données en % de la surface totale. RN = modèle par réseaux neuronaux ; MLG = modèle linéaire généralisé ; SIG = modèle SIG

65Un autre enseignement du croisement des pixels correctement prédits par les différents modèles est la similitude des résultats des deux modèles statistiques. Globalement, les surfaces correctement prédites par deux modèles sur trois le sont le plus souvent par les réseaux neuronaux et le modèle linéaire généralisé (3.75 %). Le taux de prédiction combiné de la méthode SIG avec chacun des modèles statistiques est nettement plus faible. La relative distinction du modèle SIG est corroborée par les surfaces correctement prédites par seulement une des approches (3.23 %). Une analyse catégorielle souligne ce constat. Ainsi, pour les surfaces correctement prédites par uniquement deux modèles, cinq catégories d’occupation du sol sur six le sont par les approches statistiques. A l’inverse, les surfaces uniquement prédites correctement par un seul modèle le sont le plus souvent par le modèle SIG (quatre catégories d’occupation du sol sur six). En outre, le modèle SIG se distingue par le fait que son taux de prédiction est meilleur pour des zones affectées par un changement d’état, soit entre la dernière date de la phase d’apprentissage et la date simulée, soit durant la phase d’apprentissage. Les deux méthodes statistiques, au contraire, aboutissent à des résultats légèrement plus proches de l’observation sur les secteurs stables. Ce comportement spécifique sur les marges, s’explique par les procédures d’affectation spatiale différentes des transitions temporelles où le modèle SIG fait appel à une analyse géographique de la rugosité de l’espace par évaluation multicritère. Il s’agit, en l’occurrence, d’indices à approfondir dans la perspective de l’intégration des trois modèles en un seul.

66Les résultats font également ressortir certaines limites des modélisations :

  • La variable modélisée est décrite par un nombre limité de catégories. Ce choix, intentionnel, est motivé par la qualité et l’origine des données source afin de minimiser les erreurs d’interprétation. L’inconvénient est une certaine variabilité à l’intérieur de chaque catégorie qui n’est pas prise en compte par les modèles. Ainsi peut-on voir dans les trois cartes de modélisation (fig. 7) une vaste zone prédite en broussailles dans le secteur sud-est qui est réellement occupée par des bois de feuillus. Pendant la période d’apprentissage cette zone a été photo-interprétée comme broussailles en 1980 et en 1989. Cependant les broussailles se sont densifiées, leur composition floristique a changé au profit de Quercus ilex formant une strate arborescente dominante en 2000 où la zone a été classée forêt de feuillus.

  • La période d’apprentissage est peu fournie. Les modélisations ne se basent que sur deux dates connues. Ceci pose des problèmes par rapport à la source de données à utiliser. En l’occurrence il s’agit de données haute résolution : des missions de photographies aériennes espacées dans le temps au point que l’utilisation de missions plus anciennes nous semble problématique sachant que le contexte socio-économique a beaucoup évolué.

  • Le taux d’explication de la variabilité spatio-temporelle de l’occupation du sol par les variables d’environnement est inégal selon les catégories d’occupation du sol.

  • Finalement la prévision obtenue à l’aide de chaque modèle est affecté par d’éventuels événements exceptionnels tels que des incendies de forêt, du chablis ou des programmes de reforestation. On pourrait envisager une modélisation plus complexe tenant compte d’événements rares. Cependant, dans notre contexte, où on ne dispose que de trois cartes, ceux-ci sont difficilement prévisibles.

Perspectives

67Les résultats exposés sont les premiers au sein d’un projet de recherche portant sur trois sites d’études. Les mêmes modèles seront appliqués à la Montagne de Lure en Haute Provence ainsi qu’à la Alta Alpujarra Granadina formant la partie occidentale du versant sud de la Sierra Nevada (Espagne).

68Cette comparaison croisée entre plusieurs approches appliquées sur plusieurs sites devrait, d’une part, limiter l’effet de singularité propre à chaque terrain d’études et, d’autre part, nous donner à terme des indications afin de proposer un modèle intégré, issu des trois méthodes mises en œuvre actuellement parallèlement.

69En même temps, les modèles restent évolutifs au sens où l’interprétation des premiers résultats réoriente les prochaines étapes. A ce propos il nous semble intéressant de remédier au problème de la variabilité interne à chaque catégorie d’occupation du sol par une approche semi quantitative. Aux modalités qualitatives décrites s’ajoutent des données ordonnées sur le taux de recouvrement de la strate arborescente.

70Enfin, les approches mises en œuvre sans intervention « spécialiste » (du type SIG) donnent des résultats aussi bons que la méthode « supervisée ». La combinaison de méthodes purement mathématiques avec un guidage expert pourrait permettre de gommer les imperfections du modèle. Cette voie est encore à explorer afin de progresser vers un modèle intégré de la simulation prospective de l’occupation du sol.

Top of page

Bibliography

BEAL M. et DEMUTH H. (1998) Neural Network Toolbox User’s Guide, Version 3, The Mathworks Inc.

BISHOP C.M. (1995) Neural Networks for pattern recognition, Oxford University Press, New York.

DAVALO E. et NAIM P. (1969) Des réseaux de neurones, Ed. Eyrolles.

EASTMAN J.R., KYEM P.A.K., JOLEDANO J, JIN W. (1993) Explorations in Geographic Systems Technology. volume 4: GIS and Decision Making, UNITAR, Genève.

HORNIK K. (1993) Some new results on neural network approximation, Neural Networks, 6 (8), p 1069-1072.

HOSMER D., LEMESHOW S. (1989). Applied Logistic Regression, NewYork: Wiley

KOOPERBERG C., BOSE S., STONE, J. (1997). Polychotomous Régression, J. Amer. Statist. Assoc., 92, p 117-127.

LAI T.L. et WONG S. (2001) Stochastic Neural Networks With Applications to Nonlinear Time Series, Journal of the American Statistical Association, vol 96, n° 455, Theory and Methods.

METAILIE J.-P., PAEGELOW M. (2004) La dynamique du pin à crochets (Pinus uncinata Ram.) dans l’est des Pyrénées Françaises: le retour de la forêt en montagne pastorale et métallurgique, Contribution à ouvrage collectif. Ed. Casa de Velasquez; à paraître

PAEGELOW M. (2003) Prospective modelling with GIS of land cover in Mediterranean mountain regions, 6th AGILE Conference on GIScience. 24-26 avril 2003, Lyon

PAEGELOW M., CAMACHO OLMEDO M.T., MENOR TORIBIO J. (2004) Modelizacion prospectiva del paisaje mediante Sistemas de Informacion Geografica., GEOFOCUS n° 3 p 22-44

PAEGELOW M., CAMACHO OLMEDO M.T. (2004) Le processus d’abandon des cultures et la dynamique de reconquête végétale en milieu montagnard méditerranéen : L’exemple des Garrotxes (P.O., France) et de la Alta Alpujarra Granadina (Sierra Nevada, Espagne), Sud Ouest Européen; à paraître

PARLITZ U., MERKWIRTH C. (2000) Nonlinear prediction of spatio-temporal time series, ESANN’2000 proceedings, Bruges, 26-28, p 317-322.

SAATY T.L. (1977)  A Scaling Method for Priorities in Hierarchical Structures,, J. Math. Psychology, 15, p 234-28

TABEAUD M., PECH P., SIMON L.,(2003) Weather hazards, vulnerabilities and risks in Mediterranean hinterlands from the 19th century – the Lure Mountain (France), IAG Working Group on Geochaeology ; Colloque International “dynamiques environnementales et histoire en domaines méditerranéens” ; à paraître

Top of page

Notes

1 GEODE UMR 5602 CNRS, Groupe SMASH –EA 3686 GRIMM, UTM et Instituto de Désarrollo Régional – Universidad de Granada

Top of page

References

Electronic reference

Martin Paegelow, Nathalie Villa, Laurence Cornez, Frédéric Ferraty, Louis Ferré and Pascal Sarda, « Modélisations prospectives de l’occupation du sol. Le cas d’une montagne méditerranéenne », Cybergeo: European Journal of Geography [Online], Systems, Modelling, Geostatistics, document 295, Online since 06 December 2004, connection on 18 April 2024. URL : http://journals.openedition.org/cybergeo/2811 ; DOI : https://doi.org/10.4000/cybergeo.2811

Top of page

About the authors

Martin Paegelow

GEODE UMR 5602 CNRS, Maison de la Recherche, 5, all. A. Machado, F-31058 Toulouse cedex, France
paegelow@univ-tlse2.fr

By this author

Nathalie Villa

GRIMM, Equipe d’accueil 3686, Université de Toulouse – Le Mirail, France

Laurence Cornez

Stagiaire en 2003
GEODE UMR 5602 CNRS, France

Frédéric Ferraty

GRIMM, Equipe d’accueil 3686, Université de Toulouse – Le Mirail, France

Louis Ferré

GRIMM, Equipe d’accueil 3686, Université de Toulouse – Le Mirail, France

Pascal Sarda

GRIMM, Equipe d’accueil 3686, Université de Toulouse – Le Mirail, France

Top of page

Copyright

CC-BY-4.0

The text only may be used under licence CC BY 4.0. All other elements (illustrations, imported files) are “All rights reserved”, unless otherwise stated.

Top of page
Search OpenEdition Search

You will be redirected to OpenEdition Search