Brevet FR2983312A1 — publié le 20130531

Procédé de localisation d'un élément physique compris dans un terrain 1. DOMAINE DE L'INVENTION Le domaine de l'invention est celui de la géophysique appliquée. Plus précisément, l'invention concerne une technique de traitement des données mesurées dans le cadre d'une exploration géophysique. Ainsi, l'invention s'inscrit dans le cadre de la géostatistique multivariable. L'invention a de nombreuses applications, notamment lorsqu'on souhaite générer un modèle numérique de terrain qui soit précis, à des fins de localisation de structures (objets ou bâtiments enfouis, microreliefs, etc...) dans les domaines de l'archéologie terrestre et sous-marine. 2. ARRIÈRE-PLAN TECHNOLOGIQUE L'obtention de données d'altitudes mesurées par un système de type LiDAR (acronyme anglais de « Light Detection And Ranging »), collectées en utilisant conjointement un laser émettant une onde lumineuse (de type pulsé), et un dispositif récepteur permettant de mesurer l'onde lumineuse réfléchie soit par le sol, soit par des éléments présents sur celui-ci (maison, arbres , etc...), ou de données correspondant à la valeur d'un potentiel (mesurées par un magnétomètre, un gravimètre, ...) nécessitent des traitements particuliers pour être pleinement exploitée dans des systèmes d'information géographique (« SIG ») conventionnel. Dans le cas de figure où ce sont des données de type LiDAR qui sont obtenues, on peut souhaiter créer un modèle numérique du terrain (encore noté « MNT ») étudié, qui soit le plus fidèle à la réalité. Dans le cas de figure où ce sont des données comprenant des valeurs liées à un champ magnétique (obtenues par un magnétomètre, un variomètre), ou gravitationnel (obtenues par un gravimètre de fond de mer, un gravimètre de puits ou un instrument mesurant les gradients de pesanteurs), on souhaite détecter, la source des anomalies (magnétique ou gravimétrique) détectées. On s'attache plus particulièrement dans la suite de ce document à décrire l'invention appliquée au traitement de données de type LiDAR, cependant, comme mentionné précédemment, elle ne se limite bien sûr pas à ce domaine particulier d'application, mais présente un intérêt pour toute technique de traitement de données physiques géoréférencées devant faire face à une problématique proche ou similaire. L'obtention d'un modèle numérique de surface (encore appelé « MNS ») peut être réalisé via une simple interpolation des données de type LiDAR brutes. En effet, le modèle numérique de surface permet de visualiser la surface présente sous l'appareil (un avion, un ULM, etc...) comprenant un laser permettant de réaliser l'obtention des données de type LiDAR . Cependant un tel modèle comprend à la fois des informations liées au terrain en lui-même, ainsi que des informations relatives aux objets modélisant la surface tels que, par exemple, des bâtiments, mais aussi la végétation (encore appelée la canopée). Ce modèle est fort utile pour des applications géographiques, mais il est inexploitable en tant que tel lorsqu'on cherche à détecter des éléments précis dans le terrain considéré (comme des objets ou des bâtisses enfouis qui déforment le terrain considéré), et ce même en présence d'une forte canopée (i.e lorsque le couvert forestier du terrain considéré est relativement dense). Deux exemples d'application sont décrits dans le document Dl, correspondant à l'article « Airborne LiDAR, archaeology, and the ancient Maya landscape at Caracol, Belize » de A. F. Chase et al., publié dans la revue : « Journal of Archaeological Science 38 (2011) 387e398 », et dans le document D2, correspondant à l'article « Les apports du laser aéroporté à la documentation de parcellaires anciens fossilisées par la forêt : L'exemple des champs bombés de Rastatt en Pays de Bade » de B. Sittler et al.. Ainsi, lorsqu'on cherche à révéler de tels éléments (notamment à des fins archéologiques), il existe, dans l'état de la technique, plusieurs solutions permettant de filtrer les informations utiles de celles non utiles. En effet, les données de type LiDAR relatifs à une surface donnée sont influencées par les trois paramètres suivants : la hauteur du sol, la hauteur des objets (arbres, maisons, ponts, etc...) présents sur la surface donnée, et du bruit existant (provenant des appareils de mesures, d'éléments perturbateurs comme des oiseaux, etc...). L'Homme du métier pourra consulter notamment les techniques de filtrage décrites dans le document D3 correspondant à l'article « Ground Filtering algorithme for Airborne LiDAR Data : A review of critical issues » de Meng et al., publié dans la revue « Remote Sensing 2010, 2 ». Ainsi, via le filtrage des données de type LiDAR, il est possible de distinguer une altitude « terrain » et une altitude correspondant au haut des objets (toits, arbres, etc...). De même, l'Homme du métier pourra consulter les techniques décrites dans le document D4 correspondant à un mémoire de mastère intitulé « Extracting hidden trails and roads under canopy unsing LiDAR » de A. Karatolios et al., publié par la « Naval Postgraduate School », permettant de révéler des sentiers ainsi que des routes traversant un massif forestier. Cependant, les techniques de filtrage de l'état de la technique, mentionnées précédemment, réduisent certaines informations, or cette réduction peut empêcher la détection et/ou la visualisation de structures d'intérêt pour un utilisateur. Ainsi, il existe un risque non négligeable, avec les techniques mentionnées précédemment, de ne pas faire apparaître, dans le modèle numérique de terrain, certaines structures alors qu'elles existent en réalité. Cependant il semble, dans l'état de la technique, nécessaire de procéder à une réduction des informations, car si les données ne sont pas assez réduites alors le modèle numérique de terrain obtenu sera « chargé » de bruits rendant impossible toute détection de structures fines ou présentera des artefacts de classification révélant des reliefs erronés (par exemple un microrelief présenté comme tel alors que cela ne serait, en réalité, qu'un buisson). Ces inconvénients présents dans l'état de la technique empêchent donc la détection de « petits » élément physiques (par exemple des structures présentent sous le sol et ne déformant que légèrement celui-ci). 3. OBJECTIFS DE L'INVENTION L'invention, dans au moins un mode de réalisation, a notamment pour objectif de pallier ces différents inconvénients de l'état de la technique. Plus précisément, dans au moins un mode de réalisation de l'invention, un objectif est de fournir une technique de traitement de données physiques géoréférencées permettant de révéler des informations liées à des structures fines, n'apparaissant pas via un traitement classique. Au moins un mode de réalisation de l'invention a également pour objectif de fournir une telle technique qui permet de révéler des structures possédant une taille inférieure au mètre. 4. EXPOSÉ DE L'INVENTION Dans un mode de réalisation particulier de l'invention, il est proposé un procédé de localisation d'un élément physique compris dans un terrain, ledit procédé comprenant : - une étape d'obtention d'un ensemble de données physiques géoréférencées dudit terrain, chacune desdites données physiques géoréférencées obtenues étant définie par au moins trois coordonnées x, y, z, correspondant respectivement à une longitude, une latitude et un paramètre physique ; - une étape de traitement desdites données physiques géoréférencées dudit ensemble ; - une étape de localisation d'un élément physique d'intérêt en fonction d'au moins une anomalie détectée par ladite étape de traitement. Un tel procédé est remarquable en ce que l'étape de traitement comprend : - une étape de modification des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble pour ledit terrain; - une étape d'interpolation desdites données physiques géoréférencées dudit ensemble; - une étape de détection et d'amplification d'au moins une anomalie présente au sein d'un modèle numérique de terrain réel en fonction d'une pluralité de modèles numériques de terrain fictifs, ledit modèle numérique de terrain réel étant généré à partir desdites données obtenues en sortie de l'étape d'interpolation. Ainsi, un tel procédé permet de localiser des éléments physiques enfouis dans ledit terrain ainsi que des micro-reliefs qui n'étaient pas conservés lors du traitement des données de type LiDAR effectués dans l'état de la technique. En effet, un tel procédé permet de conserver des informations liées aux données physiques mesurées, qui permettent d'identifier et de localiser, suite à une étape de traitement particulier, des éléments physiques non-localisables via les techniques de l'art antérieur. Un tel procédé permet d'identifier des structures possédant une taille inférieure au mètre (notamment des structures de 60 cm à 80 cm). Selon un aspect particulier de l'invention, l'étape de modification des coordonnées z d'un tel procédé comprend: - une étape d'obtention d'un premier point 0/ qui est un centre d'un premier cercle dans un plan (xy), ledit premier cercle étant une courbe directrice d'un premier cylindre de révolution comprenant l'ensemble desdites données physiques géoréférencées; - une étape d'obtention d'un premier angle a égal à 27rIn1 définissant n1 sections dudit premier cercle, où n1 est un entier supérieur ou égal à 3, et définissant aussi n1 parties dudit premier cylindre de révolution; - une première étape de sélection au sein de chaque partie dudit premier cylindre de révolution, de données physiques géoréférencées de sorte que ladite sélection comprend les pi pourcents de données physiques géoréférencées les plus bas au sein de chaque partie dudit premier cylindre de révolution; - une étape d'interpolation polynomiale de l'ensemble des données physiques géoréférencées sélectionnées, selon la première étape de sélection, par un polynôme bivarié Q1 de degré d avec d un entier ; - une première étape de détermination d'un premier groupe comprenant des coordonnées x,y des données physiques géoréférencées dudit ensemble et des coordonnées z déterminées à partir dudit polynôme bivarié Q1. ; - une étape de modification des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble en fonction des coordonnées z dudit ensemble et dudit premier groupe. Selon une caractéristique particulière, l'étape de modification des coordonnées z d'un tel procédé comprend en outre : - une étape d'obtention d'un deuxième point 02 qui est un centre d'un deuxième cercle dans un plan (xy), ledit deuxième cercle étant une courbe directrice d'un deuxième cylindre de révolution comprenant l'ensemble desdites données physiques géoréférencées; - une étape d'obtention d'un deuxième angle 13 égal à 27r/n2 définissant n2 sections dudit cercle, où n2 est un entier supérieur ou égal à 3, et définissant aussi n2 parties dudit deuxième cylindre de révolution; - une deuxième étape de sélection au sein de chaque partie dudit deuxième cylindre de révolution, de données physiques géoréférencées de sorte que ladite sélection comprend les p2 pourcents de données physiques géoréférencées les plus bas au sein de chaque partie dudit deuxième cylindre de révolution; - une étape d'interpolation polynomiale de l'ensemble des données physiques géoréférencées sélectionnées, selon la deuxième étape de sélection, par un polynôme bivarié Q2 de même degré d; - une deuxième étape de détermination (109) d'un deuxième groupe comprenant des coordonnées x,y des données physiques géoréférencées dudit ensemble et des coordonnées z déterminées à partir dudit polynôme bivarié bivarié Q2 ; et en ce que l'étape de modification (111) des coordonnées z est en outre réalisée en fonction des coordonnées z dudit deuxième groupe. Selon une caractéristique particulière, dans un tel procédé, lesdits premier et deuxième points 0/ et 02 sont confondus, les pourcentages pi et p2 sont identiques, et ledit deuxième angle 13 est déterminé en fonction dudit premier angle a. Selon une caractéristique particulière, dans un tel procédé, l'étape de modification des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble comprend en outre une étape de comparaison de la variance des coordonnées z avec un seuil prédéfini, et dans le cas où ladite variance est supérieure audit seuil, l'étape de modification des coordonnées z est réitérée. Selon une caractéristique particulière, un tel procédé est remarquable en ce que ladite étape d'interpolation desdites données physiques géoréférencées dudit ensemble comprend : - une étape de triangulation de Delaunay utilisant un premier maillage ; - une étape d'interpolation polynomiale pondérée utilisant ledit premier maillage . Selon une caractéristique particulière, un tel procédé est remarquable en ce que que la génération dudit modèle numérique de terrain réel comprend : - une étape de détermination d'une différence des coordonnées z provenant de ladite étape de triangulation et de ladite étape d'interpolation polynomiale pondérée pour chaque point dudit premier maillage; et - une étape de filtrage des extrêmes des coordonnées z déterminées dans l'étape de détermination d'une différence. Selon une caractéristique particulière, un tel procédé est remarquable en ce que ladite étape de détection et d'amplification comprend une étape d'obtention d'une tomographie d'indices de reliefs dudit modèle numérique réel et une étape d'obtention d'une densité tomographique de volumes. Selon une caractéristique particulière, un tel procédé est remarquable en ce que la pluralité de modèles numériques de terrain fictifs est obtenue en utilisant une fonction génératrice gaussienne définie de la manière suivante : f(x,y) = exp(-x2-y2 /(2*sigma2) , où sigma est un paramètre prenant une valeur réelle différente pour chaque modèle numérique de terrain fictif à obtenir. Selon une caractéristique particulière, un tel procédé est remarquable en ce que ladite étape de détection et d'amplification comprend en outre une étape de détermination d'horizons en fonction de ladite tomographie d'indices de reliefs obtenue. Selon une caractéristique particulière, un tel procédé est remarquable en ce que ladite étape de détection et d'amplification comprend en outre une étape de décomposition hyperbolique de relief dudit modèle numérique réel à partir d'une pluralité de modèles numériques de terrain fictifs, ladite étape de décomposition comprenant : - une étape de détermination des dérivées selon l'axe des coordonnées x et l'axe des coordonnées y de ladite pluralité de modèles numériques de terrain fictifs et dudit modèle numérique réel ; - une étape de détermination d'un premier vecteur à partir d'un modèle numérique fictif et des dérivées dudit modèle numérique réel ; - une étape de détermination d'une pluralité de vecteurs à partir des dérivées desdits modèles numériques de terrain fictifs; - une étape de décomposition dudit premier vecteur dans la famille formée par ladite pluralité de vecteurs. Selon une caractéristique particulière, un tel procédé est remarquable en ce que ladite étape de détection et d'amplification comprend en outre une étape de décomposition hyperbolique de relief tournante de relief dudit modèle numérique réel à partir d'une pluralité de modèles numériques de terrain fictifs, ladite étape de décomposition hyperbolique de relief tournante comprenant - une étape de détermination des dérivées selon un premier vecteur de direction el = kos(0), sin(0)] et un deuxième vecteur de direction e2 = [-sin(0), sin (0)1, avec 0 e [0, 27r1, de ladite pluralité de modèles numériques de terrain fictifs et dudit modèle numérique réel ; - une étape de détermination d'un deuxième vecteur à partir d'un modèle numérique fictif et des dérivées dudit modèle numérique réel selon les premier et deuxième vecteur de direction el et e2; - une étape de détermination d'une pluralité de vecteurs à partir des dérivées desdits modèles numériques de terrain fictifs selon les premier et deuxième vecteur de direction el et e2; - une étape de décomposition dudit deuxième vecteur dans la famille formée par ladite pluralité de vecteurs. Selon une caractéristique particulière, pour un tel procédé, ledit paramètre physique correspond à un élément du groupe comprenant : - une hauteur par rapport à un référentiel ; - une intensité d'un champ magnétique pour une hauteur prédéterminée; - une densité. Dans un autre mode de réalisation de l'invention, il est proposé un produit programme d'ordinateur qui comprend des instructions de code de programme pour la mise en oeuvre du procédé précité (dans l'un quelconque de ses différents modes de réalisation), lorsque ledit programme est exécuté sur un ordinateur. Dans un autre mode de réalisation de l'invention, il est proposé un dispositif de localisation d'un élément physique compris dans un terrain, ledit dispositif comprenant : - des moyens d'obtention d'un ensemble de données physiques géoréférencées dudit terrain, chacune desdites données physiques géoréférencées obtenues étant définie par au moins trois coordonnées x, y, z, correspondant respectivement à une longitude, une latitude et un paramètre physique ; - des moyens de traitement desdites données physiques géoréférencées dudit ensemble ; - des moyens de localisation d'un élément physique d'intérêt en fonction d'au moins une anomalie détectée par ladite étape de traitement. Un tel dispositif est remarquable en ce que les moyens de traitement comprennent : - des moyens de modification des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble pour ledit terrain; - des moyens d'interpolation desdites données physiques géoréférencées dudit ensemble; - des moyens de détection et d'amplification d'au moins une anomalie présente au sein d'un modèle numérique de terrain réel en fonction d'une pluralité de modèles numériques de terrain fictifs, ledit modèle numérique de terrain réel étant généré à partir desdites données obtenues en sortie des moyens d'interpolation. 5. LISTE DES FIGURES D'autres caractéristiques et avantages de l'invention apparaîtront à la lecture de la description suivante, donnée à titre d'exemple indicatif et non limitatif, et des dessins annexés, dans lesquels : la figure 1 présente, sous la forme d'un organigramme séquentiel, les différentes étapes réalisées lors d'une étape de modification de données physiques géoréférencées, selon un mode de réalisation particulier de l'invention ; la figure 2 présente un terrain pour lequel des données physiques géoréférencées ont été mesurées, et présente, de manière graphique, une sous étape d'une étape de sélection, selon un mode de réalisation particulier de l'invention ; la figure 3 présente, sous la forme d'un organigramme séquentiel, les différentes étapes réalisées suite à l'obtention de données physiques géoréférencées, selon un mode de réalisation particulier de l'invention; la figure 4 présente une grille régulière telle qu'utilisée dans des étapes de la figure 3 ; la figure 5 présente, sous la forme d'un organigramme séquentiel, les différentes étapes réalisées suite à l'obtention d'un modèle numérique de terrain réel, selon un mode de réalisation particulier de l'invention; la figure 6 présente des images de modèles numériques de terrains obtenus selon des modes de réalisation de l'invention ; la figure 7 présente des images de modèles numériques de terrains obtenus selon des modes de réalisation de l'invention ; la figure 8 présente des images résultant de l'exécution d'une étape décrite en lien avec la figure 5; la figure 9 présente la structure d'un dispositif permettant de réaliser la détection d'une élément physique selon un mode de réalisation particulier de l'invention. 6. DESCRIPTION DÉTAILLÉE Sur toutes les figures du présent document, les éléments et étapes identiques sont désignés par une même référence numérique. La figure 1 présente, sous la forme d'un organigramme séquentiel, les différentes étapes réalisées lors d'une étape de modification de données physiques géoréférencées, selon un mode de réalisation particulier de l'invention. Après avoir collecté avec un dispositif particulier des données de type LiDAR (ou tout type de données physiques géoréférencées), on peut alors représenter celles-ci, dans au moins un espace tridimensionnel, sous la forme d'un ensemble A := { (xi, y' zi) / i C [1 ; N] } , avec N un entier représentant le nombre de points collectés. On instancie ensuite un ensemble Ao := A et on note Ao = (X0, Y0, Zo) où X0, Y0 et Zo sont des vecteurs colonnes comprenant respectivement les xi, les y, et les z,. Dans un mode de réalisation particulier de l'invention, il est réalisé une étape de modification des coordonnées « z » 102 des données présentes dans l'ensemble A. Cette étape de modification 102 comprend : - une étape de détermination 103 d'une variance V relative aux coordonnées « z » de l'ensemble A, permettant ainsi de caractériser la dispersion de l'ensemble A selon l'axe des « z »; - une étape de comparaison 104 de la variance obtenue lors de l'étape 103 avec un seuil donné ; - Si la variance est strictement inférieure à ce seuil, les données de l'ensemble A n'ont plus ou pas à être modifiées, et ces données sont ensuite traitées selon une technique particulière décrite en figure 3 ; - Si la variance est supérieure ou égale à ce seuil, une première 105 et une deuxième 106 étape de sélection sont réalisées. Ainsi, l'étape de sélection 105 comprend : o une étape d'obtention d'un premier point 0/ qui est un centre d'un cercle dans un plan xy, ledit cercle étant une courbe directrice d'un cylindre de révolution comprenant l'ensemble desdites données physiques géoréférencées; o une étape d'obtention d'un premier angle a égal à 27rIn1 où n1 est un entier supérieur ou égal à 3 préalablement établi, ainsi le cercle comprend n1 sections, et définit, par construction, n1parties dudit cylindre de révolution; o une étape de sélection 105 au sein de chaque première partie dudit cylindre de révolution, de données physiques géoréférencées de sorte que ladite sélection comprend les pi pourcents de données physiques géoréférencées les plus bas au sein de ladite première partie dudit cylindre de révolution; puis, une étape d'interpolation polynomiale 107 de l'ensemble des données physiques géoréférencées sélectionnées par un polynôme bivarié Q1 de degré d avec d un entier est réalisée; ensuite une première étape de modification 109 des coordonnées z des données physiques géoréférencées obtenues à partir dudit polynôme bivarié Q1. Cette étape comprend la détermination de l'ensemble C := { yo Q1(x0Y,)) / i C L1 ; L'étape de sélection 106 comprend quant à elle : o une étape d'obtention d'un premier angle 13 égal à 27r/n2 où n2 est un entier supérieur ou égal à 3 préalablement établi, ainsi le cercle comprend n2 sections, et définit, par construction, n2 parties dudit cylindre de révolution; o une étape de sélection 106 au sein de chaque partie dudit cylindre de révolution, de données physiques géoréférencées de sorte que ladite sélection comprend les p2 pourcents de données physiques géoréférencées les plus bas au sein de ladite première partie dudit cylindre de revolution. Dans un mode de réalisation particulier de l'invention, on choisit p1 égal à p2; puis, une étape d'interpolation polynomiale 108 de l'ensemble des données physiques géoréférencées sélectionnées par un polynôme bivarié Q2 de degré d avec d un entier, est réalisée; ensuite une deuxième étape de modification 110 des coordonnées z des données physiques géoréférencées obtenues à partir dudit polynôme bivarié Q2. Cette étape comprend la détermination de l'ensemble B := { yo Q2(x0Y,)) / i C L1 Une fois les ensembles B et C obtenus, on détermine un ensemble D:= { (xi, y' F(z' Qi(x'y,), Q2(x'yi))1 i C[1 ; NI} en utilisant une fonction de sélection particulière F. La fonction de sélection particulière F peut être, par exemple, la composée des deux fonctions F1 et F2 (i.e F = F2 o F1) définies de la manière suivante : FI (z1, Qi(xoY,), Q2(xoY,)) = , min (z1, Qi(xoY,), Q2(xoY,))) ; et F2(a, b) = (a- s1, b) avec s1 un seuil choisi arbitrairement (par exemple, s1 peut, dans des modes de réalisation particulier de l'invention, être compris entre 0,05 m et 0,25 m). Ensuite on effectue l'instanciation suivante A := D, i.e désormais z, := Ftzi, Qi(x'y,), Q2(x'y,)) et on réitère l'étape de détermination de la variance pour les coordonnées « z » du nouvel ensemble A. L'objectif de la modification de l'ensemble des données physiques géoréférencées est de conserver les points M, ayant initialement une altitude « peu élevée», et de transformer des points ayant initialement une altitude élevées en des points ayant une altitude « peu élevée». La figure 2 présente un terrain pour lequel des données physiques géoréférencées ont été mesurées, et présente, de manière graphique, une sous étape d'une étape de sélection, selon un mode de réalisation particulier de l'invention. Ainsi, lors de la recherche des éléments physiques 202 et 201 présents sur ou sous ledit terrain 200, l'étape de sélection 105 ou 106 comprend la détermination d'un cercle 203 ayant un centre donné 0 (par exemple cela peut être le premier point 0/ ou être le deuxième point 02 ) et composé de n sections (par exemple on choisit n = n1 sections comme mentionné en relation avec la description de la figure 1 où on détermine alors un angle a égal à 27rIn1 ). Puis, pour chaque section, en fonction d'un même pourcentage p, on définit un plan 204 qui intersecte la partie 205 d'un cylindre de révolution ayant pour courbe directrice ledit cercle, reliée à ladite section, et sous lequel, dans la partie 205 du cylindre, il y a les p pourcent de point M, ayant les coordonnées les plus petites. Ainsi, à chaque traitement d'une partie du cylindre, on définit un nouveau plan. La figure 3 présente sous la forme d'un organigramme séquentiel, les différentes étapes réalisées suite à l'obtention de données physiques géoréférencées, selon un mode de réalisation particulier de l'invention. Dans un mode de réalisation de l'invention, on instancie 301 un ensemble E de la manière suivante : E := A. Ainsi, l'ensemble E comprend l'ensemble des données physiques géoréférencées modifiées suite à la réalisation des étapes décrites en relation avec la figure 1. Une fois l'ensemble E de données physiques géoréférencées définit, on choisit 302 une grille régulière recouvrant l'ensemble suivant { (x1, y) / i C[1 ; N] }. Une telle grille est, par exemple, représentée en figure 4, où les coordonnées d'un noeud 402 de cette grille régulière 401 (de dimension 15 x 15 dans cet exemple) correspondent au couple suivant : (x''n+ (i-1)*pas x, epas_y), avec pas _x et pas_y, les pas entre deux noeuds de la grille respectivement selon l'axe des x et selon l'axe des y, les coordonnées du point (x''' , y'',) correspondent aux coordonnées du point situé dans le coin de la grille localisé au sud ouest dans le système de géoréférencement considéré, et les indices i, j permettent d'identifier le noeud concerné de la grille. Il convient de remarquer que le pas en x et le pas en y sont dans ce mode de réalisation égaux entre eux, avec des valeurs de 0,25m ou 0,50 m. Selon un mode de réalisation particulier de l'invention, on réalise ensuite une étape de triangulation 303 selon la méthode de Delaunay pour définir la valeur en "z" en chaque noeud de la grille régulière, à partir des données présentes dans l'ensemble E. Remarquons que, dans un mode de réalisation de l'invention, lorsqu'on ne peut pas calculer les valeurs de noeuds qui ne sont pas compris dans un triangle de points de l'ensemble E, on modifie très légèrement les coordonnées x, y des points considérés pour permettre la vérification de la condition relative aux triangles utilisée dans cette méthode, et donc l'obtention d'une valeur en ces noeuds. De même, on réalise, pour la même grille régulière, la détermination de la valeur de chaque noeud, en utilisant par exemple une interpolation du type « par le voisin le plus proche » 304, ou via une technique de pondération particulière. Dans le mode de réalisation utilisant une technique de pondération, on utilise généralement une fonction à deux variables choisies parmi les deux types de fonctions suivants: les fonctions à graphe pointu et base étroite pour localiser les anomalies, ou les fonctions à sommet plat type gaussienne pour adoucir les reliefs, donc les flouter un peu. Dans un autre mode de réalisation, l'étape 304 peut consister à utiliser une technique basée sur un krigeage. Puis, on réalise 305 une interpolation polynomiale pour la même grille, à partir des données présentes dans l'ensemble E, et on moyenne les valeurs obtenues lors des étapes 304 et 305 afin d'obtenir des valeurs moyennées ou pondérées en chaque noeud de la grille régulière. Ensuite, les valeurs obtenues en chaque noeud de la grille sont soustraites 307 afin d'obtenir un écart résiduel, qui est comparé 308, à un seuil de résidu S, et à son opposé, -S. Lorsque la différence en un noeud donné est supérieure ou égale à ce seuil, ou inférieure ou égale à l'opposé de ce seuil, on fixe la différence des noeuds à la valeur du seuil ou à l'opposé de ce seuil. Ainsi, on obtient une grille dite de différence dont les valeurs de chaque noeud sont comprises entre un seuil et l'opposé de ce seuil. Cette grille de différence permet donc de mettre en exergue les différences de valeurs des noeuds suivant la technique d'interpolation utilisée. On peut représenter cette grille sous la forme d'une matrice Delta = (do), i e [1 ; T] , j e [1 ; TJ , avec les coefficients di, qui appartiennent à l'intervalle [ -S ; S [. A partir de cette grille dite de différence, on peut déterminer 309 un modèle numérique de terrain dit réel, que l'on peut représenter sous la forme d'une matrice MR, où les lignes et les colonnes correspondent aux coordonnées des noeuds de la grille régulière, et les valeurs des éléments de cette matrice correspondent aux valeurs des noeuds de la grille régulière. Dans la suite, nous notons MR = (rnro),i C [1 ; TJ , j C [1 ; TJ Ainsi, à l'indice (i,j) de la matrice MR correspond le noeud de la grille régulière ayant pour coordonnées le couple suivant (x''n+ (i-1)*pas x,y,,,,,,+(j-1)*pasy), la valeur rnro correspond à la valeur selon l'axe z en ce noeud. La valeur de z vaut do . Dans un mode de réalisation de l'invention, on réitère les étapes précédentes en choisissant 301 un ensemble E définit de la manière suivante : E := Ao. Ainsi, on obtient un modèle numérique de terrain dit réel, dans lequel, entre autres, les informations liées à la canopée sont présentes. En effectuant une différence des valeurs des noeuds entre les deux modèles numériques de terrain, il ne subsiste que les informations relatives à la canopée. La figure 5 présente, sous la forme d'un organigramme séquentiel, les différentes étapes réalisées suite à l'obtention d'un modèle numérique de terrain réel, selon un mode de réalisation particulier de l'invention. A l'étape 501, on génère un nombre W de modèles de terrains « fictifs » à partir d'une pluralité de fonctions génératrices, i.e pour W grilles régulières identiques à la grille régulière de l'étape 302, on détermine les valeurs en chaque noeud en utilisant, par exemple, la fonction définie par l'équation suivante : f(x,y) = exp(-x2-y2/(2*sigma2)) où sigma prend W valeurs différentes. Ainsi, selon une première approche, on obtient W modèles de terrains « fictifs », chacun de ces modèles pouvant être représenté par une collection de points tridimensionnels noté (X, Y, F,), (X, Y, F2), ..., (X, Y, Fw), où la collection de points (X, Y) représente les coordonnées des noeuds de la grille régulière, et les vecteurs F1, Fw représentent les valeurs obtenues par l'utilisation d'une fonction f telle que définie précédemment en les noeuds de la grille régulière. La notation (X, Y, F'), (X, Y, F2), ..., (X, Y, Fw) sera utilisée ultérieurement. Selon une deuxième approche, chaque grille régulière peut être représentée sous la forme d'une matrice carrée A(k) = ( agi (k)), i e [1 ; TJ , j E [1 ; TJ qui est de dimension T x T, où les a, sont générés par la fonction f mentionnée précédemment, et k un entier parcourant les valeurs de 1 à W. À partir de ces modèles de terrains « fictifs », et du modèle numérique réel obtenu suite à la réalisation de l'étape 309, on cherche à mettre en exergue (i.e détecter) et à amplifier des anomalies présentes au sein du modèle numérique réel lors de la réalisation d'une étape 502. L'étape 502 comprend les étapes suivantes qui peuvent être réalisées indépendamment les unes des autres : l'étape 503 consiste à déterminer un horizon total ainsi que des horizons fractionnés à partir d'une tomographie d'indices de reliefs, et à représenter le résultat sous la forme d'un modèle numérique de terrain (cf. figure 7); l'étape 504 consiste à effectuer une décomposition de relief en fonction des vecteurs de direction ex et es,; l'étape 505 consiste à effectuer une décomposition tournante de relief en fonction des vecteurs de direction el = [cos(0), sin(0)] et e2 = [-sin(0), sin (0)1, avec 0 e [0, 27r1. Ainsi, afin de réaliser l'étape 503, il faut obtenir tomographie d'indices de reliefs. Pour ce faire, il est nécessaire de réaliser la succession d'étapes décrites ci-dessous. On construit, pour chaque matrice, un vecteur colonne u(k) regroupant toutes les éléments de la matrice A(k), en concaténant chaque colonne de la matrice A(k) les unes à la suite des autres. On obtient ainsi W vecteurs colonnes. Ainsi, chaque vecteur colonne u(k) possède une dimension de T2 x 1. Dès lors, on peut calculer la matrice de Gram définit par la relation suivante : G = ( ), i E [1 ; W] , j E [1 ; W] où représente le produit scalaire. Cette matrice est une matrice carrée de dimension W x W. Puis on calcule, une matrice H = (Â.G + u. Id)] où et y sont des nombres réels choisis tels que la relation suivante est vérifiée : + µ = /. On fixe ensuite une matrice tridimensionnelle de pondération L, dont le coefficient d'indice i, j, k, noté Lijk se calcule de la manière suivante : On choisit une fonction r(k) de préférence étant affine et croissante en fonction de k, exprimée par exemple en mètre (dans le mode de réalisation particulier où l'on s'intéresse au gradient du champ magnétique, on choisit r(k) = 0,7 k+ 0,75). On cherche à repérer dans le relief des anomalies dont le diamètre au sol va croissant (famille de gaussienne). On a besoin d'un disque étroit mais plus large que l'anomalie étroite pour repérer une telle anomalie. On a besoin d'un disque large et un peu plus large que l'anomalie pour repérer une anomalie large. On note d la distance du point (i, j) de la grille à une origine (par exemple en mètre, d peut donc être considéré comme un diamètre), et on détermine u(i,j) comme étant la partie positive de 1-(d/r(k)). On détermine ensuite une fonction v(u) croissante (par exemple une fonction de puissance positive, la fonction utilisée est dans un mode de réalisation préférentiel la fonction identité (i.e v(u) = u). Lorsque la fonction v est définie par la relation suivante v(u) = ug avec une valeur de q différente de 1, la distorsion des reliefs mettrait l'accent soit sur les reliefs étroits lorsque la valeur de q est inférieure à 1, soit sur les reliefs larges, lorsque la valeur de q est supérieure à 1. Ainsi, le procédé filtre alors les reliefs repérés). On note cl la somme de tous les coefficients v(u(i,j)) pour i et j parcourant la grille de niveau k. Ainsi, le coefficient d'indices i, j, k, de la matrice tridimensionnelle de pondération L, L uk , a pour valeur v(u(i,j)/c/. À partir de la matrice MR qui correspond au modèle numérique de terrain réel interpolé, on définit une matrice tridimensionnelle, notée DR, en répliquant W fois la grille initiale. Ainsi, les coefficients de cette matrice tridimensionnelle, notés drUk, sont définis par la relation suivante : drUk = mri, pour un indice k, entier, variant de 1 à W. Cette étape de construction de la matrice tridimensionnelle consiste à dupliquer W fois le modèle numérique de terrain qui est identique, niveau par niveau (voir par l'exemple de la figure 8(a)). Ainsi, à chaque niveau k, on retrouve le même modèle numérique de terrain obtenu suite à la réalisation de l'étape 309. On réalise ensuite le produit tensoriel des matrices tridimensionnelles DR et L . Ainsi, on obtient une matrice tridimensionnelle DP ayant pour coefficient dpUk = druk * Lijk. Ainsi, le produit tensoriel permet de créer à chaque niveau un modèle numérique de terrain ayant un relief différent (voir par l'exemple de la figure 8(b)). Dès lors que chaque niveau est devenu différent des autres, on peut chercher un relief d'une taille caractéristique du niveau. Plus précisément, la répétition d'un modèle numérique de terrain (i.e la construction de la matrice DR), puis la multiplication tensorielle de la matrice DR par une matrice de pondération aboutissent à un objet composite. Ainsi, niveau par niveau (i.e pour chaque indice k), le produit tensoriel décrit précédemment force le relief initial explicité dans le modèle de terrain numérique réel MR à « ressembler » à l'anomalie cherchée à ce niveau. Cette « ressemblance » est une corrélation généralisée. Ainsi, il est nécessaire de formuler au préalable des hypothèses relatives à l'anomalie recherchée (par exemple sur les caractéristiques géométriques de l'anomalie recherchée) qui induisent des choix sur la fonction utilisée pour générer les modèles de terrains « fictifs ». Ensuite, les vecteurs colonnes u(k) possédant chacun une dimension de T2 x 1 sont regroupés pour former une matrice rectangulaire G,. = (di) ) de dimension de T2 x W. Puis, on multiplie la matrice rectangulaire G1 par la matrice H pour obtenir une matrice G2 de dimension de T2 x W. Ensuite, on multiplie un vecteur de dimension 1 x T2 , correspondant à un niveau (i.e une sous-matrice rectangulaire) de la matrice DP, dont on met bout à bout les coefficients pour la transformer en vecteur, par la matrice G2 pour obtenir un vecteur de dimension 1 x W. On réalise ensuite une extraction diagonale comprenant les étapes suivantes : pour chaque niveau k, on sélectionne la sous matrice, de la matrice tridimensionnelle DP que l'on transforme en vecteur ligne de taille T2, et l'on multiplie ce vecteur ligne par la matrice G2 afin d'obtenir un premier vecteur ligne de dimension 1 x W. Puis, on sélectionne pour le même indice k, la sous matrice de la matrice DR pour le niveau k, que l'on transforme en vecteur ligne de taille / x T2. Ce vecteur ligne est ensuite multiplié par G2 et l'on obtient alors un deuxième vecteur ligne de dimension 1 x W. L'indice de localisation au niveau k est le cosinus de l'angle de ces deux vecteurs (ou la corrélation de ces deux vecteurs ou le produit scalaire de ces vecteurs divisé par le produit de leurs normes). Ainsi, la tomographie d'indices de reliefs dudit modèle numérique réel correspond à l'ensemble formé par ces k indices de localisation obtenus. Par conséquent, chaque point de coordonnées (x,y) de l'espace se voit attribuer un vecteur de W coordonnées (W étant le nombre de niveau). Ainsi, chaque coordonnée est associée au cosinus d'un angle. En utilisant la fonction h(u) = (11 7r)( 7r -acos(u))) qui transforme le cosinus u en un angle (exprimé en tour), niveau par niveau, une zone spatiale est transformée en une tomographie d'angle. L'étape 503 de détermination d'horizon comprenant une étape consistant à appliquer à chaque indice de reliefs, la fonction (11 7r)*( 7r -acos(.)), où la fonction acos(.) correspond à la fonction arcosinus. Les angles ainsi obtenus sont strictement positifs. A chaque point (x,y) , on associe un vecteur localisation de longueur W dit indice de localisation. Ces indices sont des angles comptés en tours entre 0 et 1. On associe aussi au même point, un vecteur trivial (1 2 3 ... W) de dimension W. Un horizon total, calculé au point (x,y), est défini comme étant le barycentre du vecteur trivial pondérer par le vecteur localisation, soit le produit scalaire du vecteur trivial et du vecteur localisation, divisé par la somme des coordonnées du vecteur localisation. Ainsi, l'horizon correspond, en (x,y), à la profondeur moyenne des indices de localisation transformés en angles par la fonction (1/7r)*(7r -acos(.)). Les horizons « profonds » mettent en valeur les microreliefs invisibles sur l'image de départ. Les horizons exploitent donc de manière spatialement étendue le procédé dans son exploration du sous-sol proche. Ils mettent en valeur des anomalies grandes en horizontal et petite en vertical. Remarquons que ces opérations peuvent être menées avec une sous famille de r niveaux et le produit scalaire s'effectue sur le vecteur [1 2 .... r]. Ce nouvel indice est appelé horizon partiel. L'étape 504 de décomposition hyperbolique de relief permet la détection et l'amplification d'anomalies plus petites que celles qui sont mises en exergue par l'étape 503. Soit la collection de W modèles de reliefs « fictifs » mentionnés précédemment, représentés par les vecteurs suivants : (X, Y, F1), (X, Y, F2), ..., (X, Y, Fw), et un autre modèle de relief « fictif » représenté par les vecteurs suivants : (X, Y, J). Ainsi, on dispose de W+/ modèles de reliefs « fictifs ». On réalise ensuite une première étape de dérivation afin d'obtenir les vecteurs (X, Y, Fa), (X, Y, F,2), ..., (X, Y, F,w) et (X, Y, Jr) qui représentent les modèles de reliefs « fictifs » dérivés selon l'axe des « x » à partir des vecteurs F1, F2, ..., Fw et J. On réalise ensuite une deuxième étape de dérivation afin d'obtenir les vecteurs (X, Y, F), (X, Y, F2), ..., (X, Y, Foy) et (X, Y, J) qui représentent les modèles de reliefs « fictifs » dérivés selon l'axe des «y » à partir des vecteurs F1, F2, ..., Fw et J. Puis, on détermine, pour l'indice j, entier, variant de 1 à W les vecteurs suivants : KID = -(X*(Fv*J-Fj*J,)-Y*(Fv*J-Fj*J,)) et K2,= Kii(X.* Y) Où le symbole « * » désigne le produit tensoriel et le symbole « / » désigne la division tensorielle quand elle existe. En reprenant le modèle du relief interpolé obtenue suite à la réalisation de l'étape 309, qui peut-être représenté par le triplet de vecteurs suivant (X, Y, M), et on peut alors calculer les vecteurs : K3 = -(X*(J,*M-PM,)-Y*(J,*M-PM,)) et K4 = K3/(X.*Y) On décompose ensuite par résolution d'un système linéaire les nouveaux reliefs K3 et K4 dans la famille de vecteurs KIDK21K11. Enfin, on réalise une décomposition hyperbolique tournante du relief en réalisant les étapes précédentes de dérivation, non plus selon l'axe des « x » et des « y » mais selon le vecteur el = kos(0), sin(0)] et le vecteur e2 = [-sin(0), sin (0)] avec 0 e [0, 27r1. La figure 6 présente des images de modèles numériques de terrain obtenues selon des modes de réalisation de l'invention. Plus précisément, la figure 6(a) présente un modèle numérique de terrain d'un terrain boisé. Il convient de remarquer qu'aucun traitement sur les données utilisées pour obtenir un tel modèle n'a été utilisé. La figure 6(b) présente le résultat obtenu après la réalisation de l'étape 309. Cette étape permet de mettre en exergue des structures physiques qui étaient masquées par la présence de la canopée. La figure 6(c) correspond à un accroissement de la zone d'intérêt 601 divulguée dans la figure 6(b). Ainsi, la technique divulguée par ce document permet de manière virtuelle de couper la canopée présente dans une zone géographique sur laquelle on effectue des recherches archéologique par exemple. La figure 7 présente des images de modèles numériques de terrain obtenues selon des modes de réalisation de l'invention. Plus précisément, les figures 7(a), (b) et (c) présente des modèles numériques de terrain, centrés sur la zone d'intérêt 601, obtenus suite à la réalisation de l'étape 503 en utilisant W modèles de terrains fictifs différents. Ainsi, cette technique permet de mettre en exergue des détails qui n'apparaissent pas de prime abord au vu de la figure 6 (c). La figure 8(a) présente un modèle de terrain numérique qui a été dupliqué 5 fois, et la figure 8(b) présente l'équivalent d'une matrice tridimensionnelle DP telle que décrite en relation avec la figure 5. La figure 9 présente la structure d'un dispositif 900 permettant de réaliser la détection d'une élément physique selon un mode de réalisation particulier de l'invention. Un tel dispositif comprend une mémoire RAM 901, une unité de traitement 902, équipée par exemple d'un processeur, et pilotée par un programme d'ordinateur stocké dans une mémoire ROM 903. A l'initialisation, les instructions de code du programme d'ordinateur sont par exemple chargées dans la mémoire RAM 901 avant d'être exécutées par le processeur de l'unité de traitement 902. Cette figure 9 illustre seulement une manière particulière, parmi plusieurs possibles, de réaliser les différents étapes détaillés ci-dessus, en relation avec les figures 1 à 5. En effet, la technique de l'invention se réalise indifféremment : - sur une machine de calcul reprogrammable (un ordinateur PC, un processeur DSP ou un microcontrôleur) exécutant un programme comprenant une séquence d'instructions, ou - sur une machine de calcul dédiée (par exemple un ensemble de portes logiques comme un FPGA ou un ASIC, ou tout autre module matériel). Dans le cas où l'invention est implantée sur une machine de calcul reprogrammable, le programme correspondant (c'est-à-dire la séquence d'instructions) pourra être stocké dans un médium de stockage amovible (tel que par exemple une disquette, un CD-ROM ou un DVD-ROM) ou non, ce médium de stockage étant lisible partiellement ou totalement par un ordinateur ou un processeur. REVENDICATIONS1. Procédé de localisation d'un élément physique (201, 202 ; 601) compris dans un terrain (200), ledit procédé comprenant : - une étape d'obtention (101) d'un ensemble de données physiques géoréférencées dudit terrain, chacune desdites données physiques géoréférencées obtenues étant définie par au moins trois coordonnées x, y, z, correspondant respectivement à une longitude, une latitude et un paramètre physique ; - une étape de traitement desdites données physiques géoréférencées dudit ensemble ; - une étape de localisation d'un élément physique d'intérêt en fonction d'au moins une anomalie détectée par ladite étape de traitement ; caractérisé en ce que l'étape de traitement comprend : - une étape de modification (102) des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble pour ledit terrain; - une étape d'interpolation (303, 306) desdites données physiques géoréférencées dudit ensemble; - une étape de détection et d'amplification (502) d'au moins une anomalie présente au sein d'un modèle numérique de terrain réel en fonction d'une pluralité de modèles numériques de terrain fictifs (501), ledit modèle numérique de terrain réel étant généré (309) à partir desdites données obtenues en sortie de l'étape d'interpolation (303, 306). 2. Procédé selon la revendication 1, caractérisé en ce que l'étape de modification (102) des coordonnées z comprend: - une étape d'obtention d'un premier point 0/ qui est un centre d'un premier cercle dans un plan (xy), ledit premier cercle étant une courbe directrice d'un premier cylindre de révolution comprenant l'ensemble desdites données physiques géoréférencées; - une étape d'obtention d'un premier angle a égal à 27rIn1 définissant n1 sections dudit premier cercle, où n1 est un entier supérieur ou égal à 3, et définissant aussi n1 parties dudit premier cylindre de révolution; - une première étape de sélection (106) au sein de chaque partie dudit premier cylindre de révolution, de données physiques géoréférencées de sorte que ladite sélectioncomprend les pi pourcents de données physiques géoréférencées les plus bas au sein de chaque partie dudit premier cylindre de révolution; - une étape d'interpolation polynomiale (108) de l'ensemble des données physiques géoréférencées sélectionnées, selon la première étape de sélection, par un polynôme bivarié Q1 de degré d avec d un entier ; - une première étape de détermination (110) d'un premier groupe comprenant des coordonnées x,y des données physiques géoréférencées dudit ensemble et des coordonnées z déterminées à partir dudit polynôme bivarié Q1. ; - une étape de modification des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble en fonction des coordonnées z dudit ensemble et dudit premier groupe. 3. Procédé selon la revendication 2, caractérisé en ce que l'étape de modification des coordonnées z comprend en outre : - une étape d'obtention d'un deuxième point 02 qui est un centre d'un deuxième cercle dans un plan (xy), ledit deuxième cercle étant une courbe directrice d'un deuxième cylindre de révolution comprenant l'ensemble desdites données physiques géoréférencées; - une étape d'obtention d'un deuxième angle 13 égal à 27r/n2 définissant n2 sections dudit cercle, où n2 est un entier supérieur ou égal à 3, et définissant aussi n2 parties dudit deuxième cylindre de révolution; - une deuxième étape de sélection (105) au sein de chaque partie dudit deuxième cylindre de révolution, de données physiques géoréférencées de sorte que ladite sélection comprend les p2 pourcents de données physiques géoréférencées les plus bas au sein de chaque partie dudit deuxième cylindre de révolution; - une étape d'interpolation polynomiale (107) de l'ensemble des données physiques géoréférencées sélectionnées, selon la deuxième étape de sélection, par un polynôme bivarié Q2 de même degré d; - une deuxième étape de détermination (109) d'un deuxième groupe comprenant des coordonnées x,y des données physiques géoréférencées dudit ensemble et des coordonnées z déterminées à partir dudit polynôme bivarié bivarié Q2 ;et en ce que l'étape de modification (111) des coordonnées z est en outre réalisée en fonction des coordonnées z dudit deuxième groupe. 4. Procédé selon la revendication 3 caractérisé en ce que lesdits premier et deuxième points 0/ et 02 sont confondus et en ce que les pourcentages pi et p2 sont identiques, et en ce que ledit deuxième angle 13 est déterminé en fonction dudit premier angle a. 5. Procédé selon l'une quelconque des revendications 1 à 4 caractérisé en ce que l'étape de modification (102) des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble comprend en outre une étape de comparaison de la variance des coordonnées z avec un seuil prédéfini, et dans le cas où ladite variance est supérieure audit seuil, l'étape de modification des coordonnées z est réitérée. 6. Procédé selon l'une quelconque des revendications 1 à 5 caractérisé en ce que ladite étape d'interpolation (303, 306, 309) desdites données physiques géoréférencées dudit ensemble comprend : - une étape de triangulation de Delaunay (303) utilisant un premier maillage (302) - une étape d'interpolation polynomiale pondérée utilisant ledit premier maillage (302). 7. Procédé selon la revendication 6 caractérisé en ce que la génération (309) dudit modèle numérique de terrain réel comprend : - une étape de détermination d'une différence (307) des coordonnées z provenant de ladite étape de triangulation (303) et de ladite étape d'interpolation polynomiale pondérée pour chaque point dudit premier maillage (302); et - une étape de filtrage (308) des extrêmes des coordonnées z déterminées dans l'étape de détermination d'une différence (307). 8. Procédé selon l'une quelconque des revendications 1 à 7 caractérisé en ce que ladite étape de détection et d'amplification comprend une étape d'obtention d'une tomographie d'indices de reliefs dudit modèle numérique réel et une étape d'obtention d'une densité tomographique de volumes. 9. Procédé selon la revendication 8 caractérisé en ce que la pluralité de modèles numériques de terrain fictifs est obtenue (501) en utilisant une fonction génératricegaussienne définie de la manière suivante : f(x,y) = exp(-x2-y2 1(2*sigma2) , où sigma est un paramètre prenant une valeur réelle différente pour chaque modèle numérique de terrain fictif à obtenir. 10. Procédé selon l'une quelconque des revendications 8 et 9 caractérisé en ce que ladite étape de détection et d'amplification comprend en outre une étape de détermination (503) d'horizons en fonction de ladite tomographie d'indices de reliefs obtenue. 11. Procédé selon l'une quelconque des revendications 1 à 10 caractérisé en ce que ladite étape de détection et d'amplification comprend en outre une étape de décomposition hyperbolique (504) de relief dudit modèle numérique réel à partir d'une pluralité de modèles numériques de terrain fictifs, ladite étape de décomposition comprenant : - une étape de détermination des dérivées selon l'axe des coordonnées x et l'axe des coordonnées y de ladite pluralité de modèles numériques de terrain fictifs et dudit modèle numérique réel ; - une étape de détermination d'un premier vecteur à partir d'un modèle numérique fictif et des dérivées dudit modèle numérique réel ; - une étape de détermination d'une pluralité de vecteurs à partir des dérivées desdits modèles numériques de terrain fictifs; - une étape de décomposition dudit premier vecteur dans la famille formée par ladite pluralité de vecteurs. 12. Procédé selon l'une quelconque des revendications 1 à 11 caractérisé en ce que ladite étape de détection et d'amplification comprend en outre une étape de décomposition hyperbolique (505) de relief tournante de relief dudit modèle numérique réel à partir d'une pluralité de modèles numériques de terrain fictifs, ladite étape de décomposition hyperbolique de relief tournante comprenant - une étape de détermination des dérivées selon un premier vecteur de direction el = kos(0), sin(0)] et un deuxième vecteur de direction e2 = [-sin(0), sin (0)1, avec 0 e [0, 27r1, de ladite pluralité de modèles numériques de terrain fictifs et dudit modèle numérique réel ;- une étape de détermination d'un deuxième vecteur à partir d'un modèle numérique fictif et des dérivées dudit modèle numérique réel selon les premier et deuxième vecteur de direction el et e2; - une étape de détermination d'une pluralité de vecteurs à partir des dérivées desdits modèles numériques de terrain fictifs selon les premier et deuxième vecteur de direction el et e2; - une étape de décomposition dudit deuxième vecteur dans la famille formée par ladite pluralité de vecteurs. 13. Procédé l'une quelconque des revendications 1 à 12, caractérisé en ce que ledit paramètre physique correspond à un élément du groupe comprenant : - une hauteur par rapport à un référentiel ; - une intensité d'un champ magnétique pour une hauteur prédéterminée; - une densité. 14. Produit programme d'ordinateur, comprenant des instructions de code de programme pour la mise en oeuvre du procédé selon au moins une des revendications 1 à 13, lorsque ledit programme est exécuté sur un ordinateur. 15. Dispositif de localisation d'un élément physique compris dans un terrain, ledit dispositif comprenant : - des moyens d'obtention d'un ensemble de données physiques géoréférencées dudit terrain, chacune desdites données physiques géoréférencées obtenues étant définie par au moins trois coordonnées x, y, z, correspondant respectivement à une longitude, une latitude et un paramètre physique ; - des moyens de traitement desdites données physiques géoréférencées dudit ensemble ; - des moyens de localisation d'un élément physique d'intérêt en fonction d'au moins une anomalie détectée par ladite étape de traitement ; caractérisé en ce que les moyens de traitement comprennent : - des moyens de modification des coordonnées z pour au moins une partie des données physiques géoréférencées dudit ensemble pour ledit terrain; - des moyens d'interpolation desdites données physiques géoréférencées dudit ensemble; - des moyens de détection et d'amplification d'au moins une anomalie présente au sein d'un modèle numérique de terrain réel en fonction d'une pluralité de modèlesnumériques de terrain fictifs, ledit modèle numérique de terrain réel étant généré à partir desdites données obtenues en sortie des moyens d'interpolation.