FR2983312A1 - METHOD FOR LOCATING A PHYSICAL ELEMENT INCLUDED IN A FIELD - Google Patents

Abstract

There is provided a method for locating a physical element (201, 202), included in a terrain, comprising: a step of obtaining a set of georeferenced physical data of said terrain, each of said georeferenced physical data obtained being defined by at least three coordinates corresponding to a longitude, a latitude and a physical parameter; a step of processing said georeferenced physical data of said set; a step of locating a physical element of interest as a function of at least one anomaly detected by said processing step. The processing step is remarkable in that it comprises: a step of modifying the z coordinates for at least a portion of the georeferenced physical data of said set for said terrain; a step of interpolation of said georeferenced physical data of said set; a step of detecting and amplifying at least one anomaly.

Description

FIELD OF THE INVENTION The field of the invention is that of applied geophysics. More specifically, the invention relates to a data processing technique measured in the context of a geophysical exploration. Thus, the invention falls within the framework of multivariable geostatistics. The invention has many applications, particularly when it is desired to generate a digital model of terrain that is accurate, for the purposes of locating structures (buried objects or buildings, microreliefs, etc.) in the archeology fields. terrestrial and underwater. 2. BACKGROUND ART Obtaining altitude data measured by a Light Detection And Ranging (LiDAR) -type system, collected using a light-emitting (pulsed type) laser, and a receiver device for measuring the reflected light wave either from the ground or from elements present on it (house, trees, etc.), or data corresponding to the value of a potential (measured by a magnetometer, a gravimeter, ...) require special treatments to be fully exploited in conventional geographic information systems ("GIS"). In the case where LiDAR type data are obtained, it may be desirable to create a digital model of the terrain (still noted "DTM") studied, which is the most faithful to reality. In the case of figure where these are data comprising values related to a magnetic field (obtained by a magnetometer, a variometer), or gravitational (obtained by a seabed gravimeter, a well gravimeter or an instrument measuring the gradients of gravities), it is desired to detect, the source of the anomalies (magnetic or gravimetric) detected. Of particular interest in the remainder of this document to describe the invention applied to the processing of LiDAR type data, however, as mentioned above, it is of course not limited to this particular field of application, but presents a interest in any technique of georeferenced physical data processing facing a similar or similar problem.

Obtaining a digital surface model (also called "MNS") can be achieved via a simple interpolation of raw LiDAR type data. Indeed, the digital surface model makes it possible to visualize the surface present under the device (an airplane, an ULM, etc ...) including a laser making it possible to obtain the data of the LiDAR type. However, such a model includes both information related to the terrain itself, as well as information relating to objects modeling the surface such as, for example, buildings, but also the vegetation (also called the canopy). This model is very useful for geographical applications, but it is inapplicable as such when one seeks to detect specific elements in the terrain (such as objects or buried buildings that deform the terrain), even in presence of a strong canopy (ie when the forest cover of the considered land is relatively dense). Two examples of application are described in the document Dl, corresponding to the article "Airborne LiDAR, archeology, and the ancient Mayan landscape at Caracol, Belize" by AF Chase et al., Published in the journal: "Journal of Archaeological Science 38 (2011) 387e398 ", and in the document D2, corresponding to the article" The contributions of the airborne laser to the documentation of old parcels fossilized by the forest: The example of the bulging fields of Rastatt in Baden "of B Sittler et al .. Thus, when one seeks to reveal such elements (especially for archaeological purposes), there are, in the state of the art, several solutions making it possible to filter the useful information from the non-useful ones. Indeed, the LiDAR type data relating to a given surface are influenced by the following three parameters: the height of the ground, the height of the objects (trees, houses, bridges, etc.) present on the given surface, and the existing noise (from measuring devices, disturbing elements such as birds, etc.). Those skilled in the art will be able to consult, in particular, the filtering techniques described in document D3 corresponding to the article "Ground Filtering Algorithm for Airborne LiDAR Data: A Review of Critical Issues" by Meng et al., Published in the journal "Remote Sensing" 2010, 2 ". Thus, via the filtering of LiDAR type data, it is possible to distinguish a "terrain" altitude and an altitude corresponding to the top of the objects (roofs, trees, etc.).

Similarly, the skilled person can consult the techniques described in document D4 corresponding to a master thesis entitled "Extracting hidden trails and roads under canopy unsing LiDAR" by A. Karatolios et al., Published by the Naval Postgraduate School ", Revealing trails and roads through a forest. However, the filtering techniques of the state of the art, mentioned above, reduce certain information, but this reduction can prevent the detection and / or visualization of structures of interest for a user. Thus, there is a significant risk, with the aforementioned techniques, not to show, in the digital terrain model, some structures when they actually exist. However it seems, in the state of the art, necessary to reduce the information, because if the data are not reduced enough then the numerical model of terrain obtained will be "loaded" noises making impossible any detection of fine structures or will present classification artifacts revealing erroneous reliefs (for example a microrelief presented as such whereas it would be, in reality, only a bush). These disadvantages present in the state of the art thus prevent the detection of "small" physical elements (for example structures present under the ground and only slightly deforming it). OBJECTIVES OF THE INVENTION The invention, in at least one embodiment, has the particular objective of overcoming these various disadvantages of the state of the art. More specifically, in at least one embodiment of the invention, an objective is to provide a georeferenced physical data processing technique for revealing information related to fine structures, not appearing via conventional processing. At least one embodiment of the invention also aims to provide such a technique that reveals structures having a size less than one meter. 4. DISCLOSURE OF THE INVENTION In a particular embodiment of the invention, there is provided a method for locating a physical element comprised in a field, said method comprising: a step of obtaining a set of georeferenced physical data of said terrain, each of said georeferenced physical data obtained being defined by at least three x, y, z coordinates corresponding respectively to a longitude, a latitude and a physical parameter; a step of processing said georeferenced physical data of said set; a step of locating a physical element of interest as a function of at least one anomaly detected by said processing step.

Such a method is remarkable in that the processing step comprises: a step of modifying the z coordinates for at least part of the georeferenced physical data of said set for said terrain; a step of interpolation of said georeferenced physical data of said set; a step of detecting and amplifying at least one anomaly present within a real-field digital model as a function of a plurality of fictitious digital terrain models, said real-field digital model being generated from said data obtained at the output of the interpolation step. Thus, such a method makes it possible to locate physical elements buried in said terrain as well as micro-reliefs that were not preserved during the processing of LiDAR type data made in the state of the art. Indeed, such a method makes it possible to retain information related to the measured physical data, which makes it possible to identify and locate, following a particular processing step, non-locatable physical elements via the techniques of the prior art. Such a method makes it possible to identify structures having a size less than one meter (in particular structures from 60 cm to 80 cm). According to one particular aspect of the invention, the step of modifying the z coordinates of such a method comprises: a step of obtaining a first point 0 / which is a center of a first circle in a plane ( xy), said first circle being a guide curve of a first revolution cylinder comprising all of said georeferenced physical data; a step of obtaining a first angle equal to 27rIn1 defining n1 sections of said first circle, where n1 is an integer greater than or equal to 3, and also defining n1 parts of said first revolution cylinder; a first selection step within each portion of said first revolution cylinder of georeferenced physical data such that said selection comprises the lowest percentage of georeferenced physical data within each portion of said first revolution cylinder; a step of polynomial interpolation of the set of selected georeferenced physical data, according to the first selection step, with a bivariate polynomial Q1 of degree d with d an integer; a first step of determining a first group comprising x, y coordinates of the georeferenced physical data of said set and z coordinates determined from said bivariate polynomial Q1. ; a step of modifying the z coordinates for at least a portion of the georeferenced physical data of said set as a function of the z coordinates of said set and said first group. According to a particular characteristic, the step of modifying the z-coordinates of such a method further comprises: a step of obtaining a second point 02 which is a center of a second circle in a plane (xy), said second circle being a steering curve of a second revolution cylinder comprising all of said georeferenced physical data; a step of obtaining a second angle 13 equal to 27r / n2 defining n2 sections of said circle, where n2 is an integer greater than or equal to 3, and also defining n2 parts of said second cylinder of revolution; a second selection step within each portion of said second revolution cylinder of georeferenced physical data so that said selection comprises the lowest p2 percent of georeferenced physical data within each portion of said second revolution cylinder; a step of polynomial interpolation of the set of selected georeferenced physical data, according to the second selection step, with a bivariate polynomial Q2 of the same degree d; a second step of determining (109) a second group comprising x, y coordinates of the georeferenced physical data of said set and z coordinates determined from said bivariate bivariate polynomial Q2; and in that the step of modifying (111) the coordinates z is further carried out according to the coordinates z of said second group. According to a particular characteristic, in such a method, said first and second points O / and O2 are the same, the percentages pi and p2 are identical, and said second angle 13 is determined as a function of said first angle a. According to a particular characteristic, in such a method, the step of modifying the z coordinates for at least a portion of the georeferenced physical data of said set further comprises a step of comparing the variance of the z coordinates with a predefined threshold, and in the where said variance is greater than said threshold, the step of modifying the z coordinates is repeated. According to a particular characteristic, such a method is remarkable in that said step of interpolating said georeferenced physical data of said set comprises: a Delaunay triangulation step using a first mesh; a weighted polynomial interpolation step using said first mesh. According to one particular characteristic, such a method is remarkable in that the generation of said real-field digital model comprises: a step of determining a difference of the z-coordinates coming from said triangulation step and said polynomial interpolation step; weighted for each point of said first mesh; and a step of filtering the extremes of the coordinates z determined in the step of determining a difference. According to a particular characteristic, such a method is remarkable in that said detection and amplification step comprises a step of obtaining a tomography of relief indices of said real digital model and a step of obtaining a density tomographic volumes. According to a particular characteristic, such a method is remarkable in that the plurality of fictitious digital terrain models is obtained by using a Gaussian generating function defined as follows: f (x, y) = exp (-x2-y2 / ( 2 * sigma2), where sigma is a parameter taking a different real value for each dummy digital model to be obtained, according to one particular characteristic, such a method is remarkable in that said detection and amplification step further comprises a step of determining horizons according to said tomography of relief indices obtained According to one particular characteristic, such a method is remarkable in that said step of detection and amplification further comprises a step of hyperbolic decomposition of relief of said real digital model from a plurality of fictitious digital terrain models, said decomposition step included nt: - a step of determining the derivatives along the x-coordinate axis and the y-coordinate axis of said plurality of fictitious digital terrain models and said real digital model; a step of determining a first vector from a fictitious digital model and derivatives of said real numerical model; a step of determining a plurality of vectors from the derivatives of said fictitious digital terrain models; a step of decomposing said first vector in the family formed by said plurality of vectors.

According to a particular characteristic, such a method is remarkable in that said step of detecting and amplifying further comprises a step of hyperbolic decomposition of revolving relief relief of said real numerical model from a plurality of fictitious digital terrain models said step of hyperbolic rotational terrain decomposition comprising - a step of determining derivatives according to a first direction vector el = kos (0), sin (0)] and a second direction vector e2 = [-sin (0), sin (0) 1, with 0 e [0, 27r1, from said plurality of dummy digital terrain models and said actual digital model; a step of determining a second vector from a fictitious digital model and derivatives of said real numerical model according to the first and second vectors of direction el and e2; a step of determining a plurality of vectors from the derivatives of said fictitious digital terrain models according to the first and second direction vectors el and e2; a step of decomposing said second vector in the family formed by said plurality of vectors. According to a particular characteristic, for such a method, said physical parameter corresponds to an element of the group comprising: a height relative to a reference frame; an intensity of a magnetic field for a predetermined height; - a density. In another embodiment of the invention there is provided a computer program product which comprises program code instructions for implementing the aforesaid method (in any one of its various embodiments), when said program is run on a computer.

In another embodiment of the invention, there is provided a device for locating a physical element included in a terrain, said device comprising: means for obtaining a set of georeferenced physical data from said terrain, each said georeferenced physical data obtained being defined by at least three x, y, z coordinates corresponding respectively to a longitude, a latitude and a physical parameter; means for processing said georeferenced physical data of said set; means for locating a physical element of interest as a function of at least one anomaly detected by said processing step. Such a device is remarkable in that the processing means comprise: means for modifying the z-coordinates for at least part of the georeferenced physical data of said set for said terrain; means for interpolation of said georeferenced physical data of said set; means for detecting and amplifying at least one anomaly present within a real-field digital model as a function of a plurality of fictitious digital terrain models, said real-field digital model being generated from said data obtained at the output of the interpolation means. 5. LIST OF FIGURES Other features and advantages of the invention will appear on reading the following description, given by way of indicative and nonlimiting example, and the appended drawings, in which: FIG. form of a sequential flow chart, the various steps performed during a step of modifying georeferenced physical data, according to a particular embodiment of the invention; FIG. 2 presents a field for which georeferenced physical data has been measured, and presents, graphically, a substep of a selection step, according to a particular embodiment of the invention; FIG. 3 presents, in the form of a sequential flowchart, the various steps performed following the obtaining of georeferenced physical data, according to one particular embodiment of the invention; Figure 4 shows a regular grid as used in steps of Figure 3; FIG. 5 presents, in the form of a sequential flow chart, the various steps performed following the obtaining of a real digital terrain model, according to a particular embodiment of the invention; FIG. 6 presents digital terrain model images obtained according to embodiments of the invention; FIG. 7 presents digital terrain model images obtained according to embodiments of the invention; FIG. 8 shows images resulting from the execution of a step described with reference to FIG. 5; Figure 9 shows the structure of a device for performing the detection of a physical element according to a particular embodiment of the invention. 6. DETAILED DESCRIPTION In all the figures of this document, the elements and identical steps are designated by the same numerical reference. FIG. 1 presents, in the form of a sequential flowchart, the various steps performed during a step of modifying georeferenced physical data, according to a particular embodiment of the invention. After having collected with a particular device LiDAR type data (or any type of georeferenced physical data), we can then represent them, in at least a three-dimensional space, in the form of a set A: = {(xi , y 'zi) / C [1; N]}, with N an integer representing the number of points collected. We then instantiate a set Ao: = A and we write Ao = (X0, Y0, Zo) where X0, Y0 and Zo are column vectors respectively comprising xi, y, and z ,. In a particular embodiment of the invention, a step of modifying the "z" coordinates 102 of the data present in the set A.

This modifying step 102 comprises: a determination step 103 of a variance V relative to the coordinates "z" of the set A, thus making it possible to characterize the dispersion of the set A along the "z" axis; a comparison step 104 of the variance obtained during step 103 with a given threshold; If the variance is strictly less than this threshold, the data of the set A no longer have to be modified, and these data are then processed according to a particular technique described in FIG. 3; If the variance is greater than or equal to this threshold, a first 105 and a second selection step 106 are performed. Thus, the selection step 105 comprises: a step of obtaining a first point 0 / which is a center of a circle in a plane xy, said circle being a steering curve of a revolution cylinder comprising set of said georeferenced physical data; o a step of obtaining a first angle equal to 27rIn1 where n1 is an integer greater than or equal to 3 previously established, and the circle comprises n1 sections, and defines, by construction, n1parties of said cylinder of revolution; a selection step 105 within each first portion of said revolution cylinder of georeferenced physical data such that said selection comprises the lowest percentage of georeferenced physical data within said first portion of said revolution cylinder; then, a polynomial interpolation step 107 of the set of georeferenced physical data selected by a bivariate polynomial Q1 of degree d with an integer is carried out; then a first step 109 of modifying the z coordinates of the georeferenced physical data obtained from said bivariate polynomial Q1. This step comprises determining the set C: = {yo Q1 (x0Y,)) / i C L1; The selection step 106 comprises: a step of obtaining a first angle 13 equal to 27r / n2 where n2 is an integer greater than or equal to 3 previously established, and the circle comprises n2 sections, and defines by construction, n2 parts of said cylinder of revolution; a selection step 106 within each portion of said revolution cylinder of georeferenced physical data such that said selection comprises the lowest p2 percent of georeferenced physical data within said first portion of said revolution cylinder. In a particular embodiment of the invention, p1 is chosen equal to p2; then, a step of polynomial interpolation 108 of the set of georeferenced physical data selected by a bivariate polynomial Q2 of degree d with an integer, is carried out; then a second step of modifying the coordinates z of the georeferenced physical data obtained from said bivariate polynomial Q2. This step includes the determination of the set B: = {yo Q2 (x0Y,)) / i C L1 Once the sets B and C are obtained, we determine a set D: = {(xi, y 'F (z' Qi (x'y,), Q2 (x'yi)) 1 i C [1; NI] using a particular selection function F.

The particular selection function F can be, for example, the compound of the two functions F1 and F2 (ie F = F2 or F1) defined in the following way: FI (z1, Qi (xoY,), Q2 (xoY,)) =, min (z1, Qi (xoY,), Q2 (xoY,))); and F2 (a, b) = (a-s1, b) with s1 an arbitrarily chosen threshold (for example, s1 may, in particular embodiments of the invention, be between 0.05 m and 0.25 m). Then we perform the following instantiation A: = D, ie now z,: = Ftzi, Qi (x'y,), Q2 (x'y,)) and repeat the variance determination step for the coordinates "Z" of the new set A.

The objective of the modification of all the georeferenced physical data is to keep the M points, initially having a "low" altitude, and to transform points initially having a high altitude into points having a "low" altitude. . FIG. 2 presents a field for which georeferenced physical data have been measured, and presents, graphically, a substep of a selection step, according to a particular embodiment of the invention. Thus, when searching for the physical elements 202 and 201 present on or under said ground 200, the selection step 105 or 106 comprises the determination of a circle 203 having a given center 0 (for example this may be the first point 0 / or be the second point 02) and composed of n sections (for example we choose n = n1 sections as mentioned in connection with the description of Figure 1 where we then determine an angle equal to 27rIn1). Then, for each section, according to the same percentage p, a plane 204 is defined which intersects the portion 205 of a cylinder of revolution whose steering curve is said circle, connected to said section, and under which, in the part 205 of the cylinder, there are the p percent of point M, having the smallest coordinates. Thus, at each treatment of a portion of the cylinder, a new plane is defined. FIG. 3 presents in the form of a sequential flowchart, the various steps performed following the obtaining of georeferenced physical data, according to a particular embodiment of the invention.

In one embodiment of the invention, a set E is instantiated 301 in the following manner: E: = A. Thus, the set E comprises all the georeferenced physical data modified following the realization of the steps described in relation to with Figure 1.

Once the set E of georeferenced physical data defines, one chooses 302 a regular grid covering the following set {(x1, y) / i C [1; NOT] }. Such a grid is, for example, represented in FIG. 4, where the coordinates of a node 402 of this regular grid 401 (of dimension 15 × 15 in this example) correspond to the following pair: (x''n + (i-1 ) * not x, epas_y), with steps _x and pas_y, the steps between two nodes of the grid respectively along the axis of the x and along the axis of the y, the coordinates of the point (x '' ', y' ' ,) correspond to the coordinates of the point situated in the corner of the grid located to the southwest in the georeferencing system considered, and the indices i, j make it possible to identify the node concerned of the grid. It should be noted that the step x and pitch are there in this embodiment equal to each other, with values of 0.25m or 0.50 m. According to a particular embodiment of the invention, a triangulation step 303 is then carried out according to the Delaunay method to define the "z" value at each node of the regular grid, from the data present in the set E Note that, in one embodiment of the invention, when it is not possible to calculate the values of nodes that are not included in a point triangle of the set E, the x, y coordinates are very slightly modified. points considered to allow verification of the condition relating to triangles used in this method, and thus obtaining a value in these nodes. Likewise, the determination of the value of each node is carried out for the same regular grid, for example using an interpolation of the "nearest neighbor" type 304, or via a particular weighting technique. In the embodiment using a weighting technique, a two-variable function is generally chosen from the two following types of functions: the functions with a sharp graph and a narrow base for locating anomalies, or Gaussian flat-top functions for soften the reliefs, so blur them a bit.

In another embodiment, step 304 may involve using a kriging-based technique. Then, a polynomial interpolation is performed for the same grid, from the data present in the set E, and the values obtained during the steps 304 and 305 are averaged in order to obtain averaged or weighted values at each node of the regular grid. Then, the values obtained at each node of the grid are subtracted 307 in order to obtain a residual difference, which is compared 308, with a residual threshold S, and its opposite, -S. When the difference in a given node is greater than or equal to this threshold, or less than or equal to the opposite of this threshold, the difference of the nodes is fixed at the threshold value or at the opposite of this threshold. Thus, we obtain a so-called difference grid whose values of each node are between a threshold and the opposite of this threshold. This difference grid thus makes it possible to highlight the differences in values of the nodes according to the interpolation technique used. We can represent this grid in the form of a matrix Delta = (do), i e [1; T], j e [1; TJ, with the coefficients di, which belong to the interval [-S; S [. From this so-called difference grid, it is possible to determine a real land digital model, which can be represented in the form of a matrix MR, where the rows and columns correspond to the coordinates of the nodes of the grid. regular, and the values of the elements of this matrix correspond to the values of the nodes of the regular grid. In the following, we write MR = (rnro), i C [1; TJ, C [1; TJ Thus, at the index (i, j) of the matrix MR corresponds the node of the regular grid having for coordinates the following pair (x''n + (i-1) * not x, y ,,,,,,, + (j-1) * pasy), the value rnro corresponds to the value along the z axis at this node. The value of z is do. In one embodiment of the invention, the preceding steps are repeated by choosing 301 a set E defined as follows: E: = Ao. Thus, we obtain a digital model of real land, in which, among other things, information related to the canopy is present. By making a difference in the values of the nodes between the two digital terrain models, only the information relating to the canopy remains. FIG. 5 presents, in the form of a sequential flowchart, the various steps performed following the obtaining of a real digital terrain model, according to a particular embodiment of the invention. In step 501, a number W of "fictitious" terrain models is generated from a plurality of generating functions, ie for W regular grids identical to the regular grid of step 302, the values in each node using, for example, the function defined by the following equation: f (x, y) = exp (-x2-y2 / (2 * sigma2)) where sigma takes W different values. Thus, according to a first approach, we obtain W "fictitious" terrain models, each of these models being able to be represented by a collection of three-dimensional points noted (X, Y, F,), (X, Y, F2), .. ., (X, Y, Fw), where the collection of points (X, Y) represents the coordinates of the nodes of the regular grid, and the vectors F1, Fw represent the values obtained by the use of a function f such than previously defined in the nodes of the regular grid. The notation (X, Y, F '), (X, Y, F2), ..., (X, Y, Fw) will be used later. According to a second approach, each regular grid can be represented in the form of a square matrix A (k) = (agi (k)), i e [1; TJ, j E [1; TJ which is of dimension T x T, where the a, are generated by the function f mentioned previously, and k an integer traversing the values of 1 to W. From these models of lands "fictitious", and of the real numerical model obtained following the completion of step 309, it is sought to highlight (ie detect) and amplify anomalies present in the real numerical model when performing a step 502.

Step 502 comprises the following steps that can be performed independently of one another: step 503 consists in determining a total horizon as well as split horizons from a tomography of relief indices, and to represent the result in the form of a digital terrain model (see Figure 7); step 504 is to perform a relief decomposition according to the direction vectors ex and es ,; step 505 consists in performing a rotational decomposition of relief as a function of the direction vectors el = [cos (0), sin (0)] and e2 = [-sin (0), sin (0) 1, with 0 e [0, 27r1. Thus, in order to perform step 503, it is necessary to obtain tomography of relief indices. To do this, it is necessary to carry out the succession of steps described below. For each matrix, we construct a column vector u (k) grouping all the elements of the matrix A (k), concatenating each column of the matrix A (k) one after the other. We thus obtain W column vectors. Thus, each column vector u (k) has a dimension of T2 x 1. From then on, we can calculate the Gram matrix defined by the following relation: G = (<u (i), u (')>), i E [1; W], j E [1; W] where <. > represents the scalar product. This matrix is a square matrix of dimension W x W. Then one calculates, a matrix H = (Â.G + u. Id)] where and y are chosen real numbers such that the following relation is verified: + μ = / .

Then a three-dimensional weighting matrix L, whose index coefficient i, j, k, denoted Lijk, is calculated as follows: A function r (k) is preferably selected being affine and increasing as a function of k, expressed for example in meters (in the particular embodiment where we are interested in the gradient of the magnetic field, we choose r (k) = 0.7 k + 0.75). We seek to locate in the relief anomalies whose ground diameter is increasing (family Gaussian). A narrow but wider disk is needed than the narrow anomaly to identify such an anomaly. We need a wide disc and a little wider than the anomaly to identify a wide anomaly. We denote the distance from the point (i, j) of the grid to an origin (for example in meters, d can therefore be considered as a diameter), and we determine u (i, j) as being the positive part of 1 - (d / r (k)). An increasing function v (u) is then determined (for example a function of positive power, the function used is in a preferential embodiment the identity function (ie v (u) = u) When the function v is defined by the following relation v (u) = ug with a value of q different from 1, the distortion of the reliefs would focus either on the narrow reliefs when the value of q is less than 1, or on the large reliefs, when the value of q is greater than 1. Thus, the method then filters the reliefs identified). We denote by cl the sum of all the coefficients v (u (i, j)) for i and j running through the grid of level k. Thus, the coefficient of indices i, j, k, of the three-dimensional matrix of weighting L, L uk, has for value v (u (i, j) / c /. From the matrix MR which corresponds to the numerical model interpolated real field, we define a three-dimensional matrix, denoted DR, by replicating W times the initial grid.So, the coefficients of this three-dimensional matrix, denoted drUk, are defined by the following relation: drUk = mri, for an index k, integer, ranging from 1 to W.

This step of constructing the three-dimensional matrix consists in duplicating W times the digital terrain model which is identical, level by level (see the example of FIG. 8 (a)). Thus, at each level k, there is the same numerical model of terrain obtained following the completion of step 309. The tensor product of the three-dimensional matrices DR and L is then produced.

Thus, a three-dimensional matrix DP having the coefficient dpUk = druk * Lijk is obtained. Thus, the tensor product allows to create at each level a digital model of terrain having a different relief (see the example of Figure 8 (b)). Since each level has become different from the others, one can look for a relief of a size characteristic of the level. More precisely, the repetition of a digital terrain model (i.e the construction of the matrix DR), then the tensor multiplication of the matrix DR by a weighting matrix result in a composite object. Thus, level by level (i.e for each index k), the tensor product described above forces the initial relief explained in the real digital terrain model MR to "resemble" the anomaly sought at this level. This "resemblance" is a generalized correlation. Thus, it is necessary to first formulate hypotheses relating to the anomaly sought (for example, on the geometric characteristics of the anomaly sought) that induce choices on the function used to generate the "fictitious" terrain models. Then, the column vectors u (k) each having a dimension of T2 x 1 are grouped together to form a rectangular matrix G 1. = (di)) of T2 x W dimension.

Then, one multiplies the rectangular matrix G1 by the matrix H to obtain a matrix G2 of dimension of T2 x W. Then, one multiplies a vector of dimension 1 x T2, corresponding to a level (ie a rectangular sub-matrix) of the DP matrix, which is put end to end coefficients to transform it into a vector, by the matrix G2 to obtain a vector dimension 1 x W. Diagonal extraction is then carried out comprising the following steps: for each level k, we select the under the matrix, of the three-dimensional matrix DP which is transformed into a line vector of size T2, and this line vector is multiplied by the matrix G2 in order to obtain a first line vector of dimension 1 × W. Then, one selects for the same index k, the sub-matrix of the matrix DR for the level k, which one transforms into vector line of size / x T2. This line vector is then multiplied by G2 and a second line vector of dimension 1 × W is then obtained.

The location index at level k is the cosine of the angle of these two vectors (or the correlation of these two vectors or the scalar product of these vectors divided by the product of their norms). Thus, the tomography of relief indices of said real numerical model corresponds to the set formed by these k location indices obtained. Therefore, each coordinate point (x, y) of the space is assigned a vector of W coordinates (where W is the level number). Thus, each coordinate is associated with the cosine of an angle. Using the function h (u) = (11 7r) (7r -acos (u))) which converts the cosine u to an angle (expressed in turns), level by level, a spatial area is transformed into a tomography of angle.

The horizon determination step 503 comprises a step of applying to each index of reliefs the function (11 7r) * (7r -acos (.)), Where the acos (.) Function corresponds to the arcosinus function. The angles thus obtained are strictly positive. At each point (x, y), we associate a location vector of length W said location index. These indices are angles counted in revolutions between 0 and 1. We also associate at the same point, a trivial vector (1 2 3 ... W) of dimension W. A total horizon, calculated at the point (x, y), is defined as the centroid of the trivial vector weighted by the localization vector, ie the scalar product of the trivial vector and the localization vector, divided by the sum of the coordinates of the localization vector. Thus, the horizon corresponds, in (x, y), to the average depth of the localization indices transformed into angles by the function (1 / 7r) * (7r -acos (.)). The "deep" horizons highlight the invisible micro-reliefs on the initial image. The horizons thus exploit the process in a spatially extended way in its exploration of the nearby subsoil. They highlight large anomalies in horizontal and small in vertical.

Note that these operations can be carried out with a sub-family of r levels and the scalar product is carried out on the vector [1 2 .... r]. This new index is called a partial horizon. Step 504 of hyperbolic relief decomposition allows the detection and amplification of anomalies smaller than those highlighted by step 503. Either the collection of W "fictitious" relief models mentioned above, represented by the following vectors: (X, Y, F1), (X, Y, F2), ..., (X, Y, Fw), and another "dummy" relief model represented by the following vectors: (X, Y, J). Thus, we have W + / "fictitious" relief models.

A first derivation step is then performed in order to obtain the vectors (X, Y, Fa), (X, Y, F, 2), ..., (X, Y, F, w) and (X, Y , Jr) which represent the "fictitious" relief models derived along the "x" axis from the vectors F1, F2, ..., Fw and J. A second derivation step is then carried out in order to obtain the vectors (X, Y, F), (X, Y, F2), ..., (X, Y, Foy) and (X, Y, J) which represent the "fictitious" relief models derived along the axis "y" from vectors F1, F2, ..., Fw and J. Then, for the index j, integer, ranging from 1 to W, the following vectors are determined: KID = - (X * (Fv * J-Fj * J,) - Y * (Fv * J-Fj * J,)) and K2, = Kii (X. * Y) Where the symbol "*" designates the tensor product and the symbol "/" designates Tensory division when it exists.

By taking again the model of the interpolated relief obtained following the realization of the step 309, which can be represented by the following triplet of vectors (X, Y, M), and one can then calculate the vectors: K3 = - (X * (J, * M-PM,) - Y * (J, * M-PM,)) and K4 = K3 / (X. * Y) The new K3 and K4 reliefs are then decomposed by solving a linear system. in the family of vectors KIDK21K11. Finally, a rotating hyperbolic decomposition of the relief is carried out by carrying out the preceding steps of derivation, no longer along the axis of "x" and "y" but according to the vector el = kos (0), sin (0)] and the vector e2 = [-sin (0), sin (0)] with 0 e [0, 27r1.

Figure 6 shows digital terrain model images obtained according to embodiments of the invention. Specifically, Figure 6 (a) presents a digital terrain model of a woodlot. It should be noted that no processing on the data used to obtain such a model has been used.

FIG. 6 (b) shows the result obtained after the completion of step 309. This step makes it possible to highlight physical structures that were masked by the presence of the canopy. Figure 6 (c) corresponds to an increase in the area of interest 601 disclosed in Figure 6 (b). Thus, the technique disclosed by this document makes it possible virtually to cut the canopy present in a geographical area on which archaeological research is carried out for example. Figure 7 presents digital terrain model images obtained according to embodiments of the invention. More specifically, FIGS. 7 (a), (b) and (c) present digital terrain models, centered on the area of interest 601, obtained following the realization of step 503 using W fictitious terrain models. different. Thus, this technique makes it possible to highlight details that do not appear at first sight in view of Figure 6 (c). Fig. 8 (a) shows a digital terrain model which has been duplicated 5 times, and Fig. 8 (b) shows the equivalent of a 3-dimensional DP matrix as described in connection with Fig. 5.

Figure 9 shows the structure of a device 900 for performing the detection of a physical element according to a particular embodiment of the invention. Such a device comprises a RAM 901, a processing unit 902, equipped for example with a processor, and controlled by a computer program stored in a ROM 903. At initialization, the code instructions of the program of computer are for example loaded into the RAM 901 before being executed by the processor of the processing unit 902. This FIG. 9 illustrates only one particular way, among several possible, of carrying out the various steps detailed above, in relation to FIGS. 1 to 5. Indeed, the technique of the invention is carried out indifferently: on a reprogrammable calculation machine (a PC computer, a DSP processor or a microcontroller) executing a program comprising a sequence of instructions , or - on a dedicated computing machine (for example a set of logical gates such as an FPGA or an ASIC, or any other hardware module). In the case where the invention is implemented on a reprogrammable calculation machine, the corresponding program (that is to say the sequence of instructions) can be stored in a removable storage medium (such as for example a diskette, a CD-ROM or a DVD-ROM) or not, this storage medium being readable partially or totally by a computer or a processor.

Claims (15)

REVENDICATIONS1. A method of locating a physical element (201, 202, 601) included in a field (200), said method comprising: - a step of obtaining (101) a set of georeferenced physical data from said field, each of said data georeferenced physical features obtained being defined by at least three x, y, z coordinates corresponding respectively to a longitude, a latitude and a physical parameter; a step of processing said georeferenced physical data of said set; a step of locating a physical element of interest as a function of at least one anomaly detected by said processing step; characterized in that the processing step comprises: - a step of modifying (102) the z-coordinates for at least a portion of the georeferenced physical data of said set for said terrain; an interpolation step (303, 306) of said georeferenced physical data of said set; a step of detecting and amplifying (502) at least one anomaly present within a real-field digital model as a function of a plurality of fictitious digital terrain models (501), said digital terrain model real being generated (309) from said data obtained at the output of the interpolation step (303, 306). 2. Method according to claim 1, characterized in that the step of modifying (102) the coordinates z comprises: a step of obtaining a first point 0 / which is a center of a first circle in a plane (xy), said first circle being a guide curve of a first revolution cylinder comprising all of said georeferenced physical data; a step of obtaining a first angle equal to 27rIn1 defining n1 sections of said first circle, where n1 is an integer greater than or equal to 3, and also defining n1 parts of said first revolution cylinder; a first selection step (106) within each portion of said first revolution cylinder of georeferenced physical data such that said selection comprises the lowest percent georeferenced physical data within each portion of said first revolution cylinder; a step of polynomial interpolation (108) of all the selected georeferenced physical data, according to the first selection step, with a bivariate polynomial Q1 of degree d with d an integer; a first step of determining (110) a first group comprising x, y coordinates of the georeferenced physical data of said set and z coordinates determined from said bivariate polynomial Q1. ; a step of modifying the z coordinates for at least a portion of the georeferenced physical data of said set as a function of the z coordinates of said set and said first group. 3. Method according to claim 2, characterized in that the step of modifying the z coordinates further comprises: a step of obtaining a second point 02 which is a center of a second circle in a plane (xy ), said second circle being a steering curve of a second revolution cylinder comprising all of said georeferenced physical data; a step of obtaining a second angle 13 equal to 27r / n2 defining n2 sections of said circle, where n2 is an integer greater than or equal to 3, and also defining n2 parts of said second cylinder of revolution; a second selection step (105) within each portion of said second revolution cylinder, of georeferenced physical data so that said selection comprises the lowest p2 percent of georeferenced physical data within each portion of said second revolution cylinder; ; a step of polynomial interpolation (107) of all the selected georeferenced physical data, according to the second selection step, with a bivariate polynomial Q2 of the same degree d; a second step of determining (109) a second group comprising x, y coordinates of the georeferenced physical data of said set and z coordinates determined from said bivariate bivariate polynomial Q2; and that the modifying step (111) ) z coordinates is further carried out according to the z coordinates of said second group. 4. Method according to claim 3 characterized in that said first and second points 0 / and 02 are combined and in that the percentages pi and p2 are identical, and in that said second angle 13 is determined as a function of said first angle a . 5. Method according to any one of claims 1 to 4 characterized in that the step of modifying (102) z coordinates for at least a portion of the georeferenced physical data of said set further includes a step of comparing the variance of the coordinates z with a predefined threshold, and in the case where said variance is greater than said threshold, the step of modifying the coordinates z is repeated. 6. Method according to any one of claims 1 to 5 characterized in that said step of interpolating (303, 306, 309) said georeferenced physical data of said set comprises: - a Delaunay triangulation step (303) using a first mesh (302) - a weighted polynomial interpolation step using said first mesh (302). 7. Method according to claim 6, characterized in that the generation (309) of said real land digital model comprises: a step of determining a difference (307) of the z-coordinates coming from said triangulation step (303) and said weighted polynomial interpolation step for each point of said first mesh (302); and a step of filtering (308) the extremes of the coordinates z determined in the step of determining a difference (307). 8. Method according to any one of claims 1 to 7 characterized in that said detection and amplification step comprises a step of obtaining a tomography of relief indices of said real digital model and a step of obtaining a tomographic density of volumes. 9. Method according to claim 8, characterized in that the plurality of fictitious digital terrain models is obtained (501) using a Gaussian generating function defined as follows: f (x, y) = exp (-x2-y2 1 ( 2 * sigma2), where sigma is a parameter that takes a different real value for each dummy digital model to obtain. 10. Method according to any one of claims 8 and 9 characterized in that said detection and amplification step further comprises a step of determining (503) horizon according to said tomography of relief indices obtained. 11. Method according to any one of claims 1 to 10 characterized in that said detecting and amplifying step further comprises a step of hyperbolic decomposition (504) of relief of said real numerical model from a plurality of models. fictitious digital terrain maps, said decomposition step comprising: - a step of determining the derivatives along the x-coordinate axis and the y-coordinate axis of said plurality of fictitious digital terrain models and said real digital model; a step of determining a first vector from a fictitious digital model and derivatives of said real numerical model; a step of determining a plurality of vectors from the derivatives of said fictitious digital terrain models; a step of decomposing said first vector in the family formed by said plurality of vectors. 12. Method according to any one of claims 1 to 11 characterized in that said detecting and amplifying step further comprises a step of hyperbolic decomposition (505) of relief of relief of said real digital model from a plurality of fictitious digital terrain models, said hyperbolic rotational terrain decomposition step comprising - a step of determining derivatives according to a first direction vector el = kos (0), sin (0)] and a second direction vector e2 = [-sin (0), sin (0) 1, with 0 e [0, 27r1, from said plurality of fictitious digital terrain models and said actual digital model; -a step of determining a second vector from a fictitious digital model and derivatives of said real numerical model according to the first and second direction vector el and e2; a step of determining a plurality of vectors from the derivatives of said fictitious digital terrain models according to the first and second direction vectors el and e2; a step of decomposing said second vector in the family formed by said plurality of vectors. 13. Method according to any one of claims 1 to 12, characterized in that said physical parameter corresponds to an element of the group comprising: a height relative to a reference frame; an intensity of a magnetic field for a predetermined height; - a density. A computer program product, comprising program code instructions for implementing the method according to at least one of claims 1 to 13, when said program is executed on a computer. 15. Device for locating a physical element included in a terrain, said device comprising: means for obtaining a set of georeferenced physical data of said terrain, each of said georeferenced physical data obtained being defined by at least three coordinates x , y, z, corresponding respectively to a longitude, a latitude and a physical parameter; means for processing said georeferenced physical data of said set; means for locating a physical element of interest as a function of at least one anomaly detected by said processing step; characterized in that the processing means comprises: z-coordinate modification means for at least a portion of the georeferenced physical data of said set for said terrain; means for interpolation of said georeferenced physical data of said set; means for detecting and amplifying at least one anomaly present within a real-field digital model as a function of a plurality of fictitious digital field models, said real-field digital model being generated from said data obtained at the output of the interpolation means.