This application claims priority to foreign French patent application No. FR 2202800, filed on Mar. 29, 2022, the disclosure of which is incorporated by reference in its entirety.
The invention relates to the field of ground-penetrating radars or georadars which cover all of the techniques making it possible to detect, locate or identify underground targets by means of a radio-frequency system.
Underground targets are, for example, piping of different diameters and types (steel, PVC, cement, concrete, etc.) which may be buried at various depths.
One objective of ground-penetrating radars is to locate such objects with precision in order to be able to correctly map a subsurface, for example for safety needs during works.
The images delivered by a ground-penetrating radar contain unwanted elements that may be of a number of types (effect of antenna coupling, thermal noise, radio interference, etc.). When an object is irradiated by the radar, it reflects energy that may be measured. Since the ground is a non-uniform medium, imaging algorithms generally cause a clutter of points of interest to appear. This clutter consists of reflections of the radar signal from an interface between two segments of the ground of different natures or from a reflective object present in the ground (pebbles, stones).
One problem addressed by the present invention is that of detecting and mapping the presence of piping or more generally of longitudinal structures in the ground while decreasing the clutter of points of interest present in the 3D image reconstructed based on the acquired radar signals.
Clutter-decreasing techniques are widely addressed in the literature. Most existing techniques may be classed into two main categories: methods based on modelling of the clutter, and methods aiming to decrease clutter by filtering.
The main disadvantages of methods for modelling clutter reside in the fact that the performance of these techniques is dependent on the adopted clutter model, on the difficulty of precisely estimating the parameters of the model, and/or on prior knowledge of the response of the terrain without targets. Techniques based on modelling of the (clutter and target) signals are especially described in references [1], [2], [3], [4], and [5].
As regards techniques based on filtering clutter, one drawback is that certain methods of this kind as they filter the clutter degrade the signal corresponding to the targets, whereas others need to make the assumption that the clutter signal is stronger than that of the targets or that the frequency spectrum of the clutter signal is concentrated in a region different from that of the signal of the targets. Examples of such methods are given in references [6], [7] and [8].
Moreover, techniques based on shape detection via Hough transform are described in articles [9] and [10]. These methods are exclusively dedicated to 2D detection of the radargram of a ground-penetrating radar on a uniform grid. These methods are based on recognition of a parabola in a radargram for a radar with a single transmit-receive antenna (SISO device).
The invention provides a method for detecting longitudinal structures in a 3D cloud of points of interest that is determined based on a plurality of radar measurements of a region of ground. The proposed method is based on searching for lines in the 3D space that have a sufficient number of points of the point cloud located in proximity to them. Detection of lines in the space of the 3D point cloud makes it possible to characterize longitudinal shapes corresponding to objects to be detected such as piping, pipelines and pipes.
In one variant of embodiment, the invention also makes provision to pre-process the point cloud so as to allow points belonging to the clutter to be filtered by means of a search for unit vectors that are substantially collinear in the 3D point cloud.
Contrary to methods for modelling clutter, the invention requires no assumption to be made in respect of a statistical model of the clutter.
Contrary to existing solutions based on the Hough transform, the provided solution is based on processing of a non-uniformly sampled 3D point cloud, each point corresponding to the potential detection of a region of interest.
One subject of the invention is a method for detecting buried longitudinal structures using a ground-penetrating radar, the method comprising the steps of: acquiring a plurality of radar signals for a region of ground, determining, based on said radar signals, a 3D point cloud, each point corresponding to one radar detection and being geolocated in space, searching for at least one set of substantially aligned points in the 3D point cloud by:
According to one particular aspect of the invention, the search for at least one set of substantially aligned points is carried out by means of steps of:
According to one particular aspect of the invention, the search for at least one set of substantially aligned points is carried out by means of steps of:
According to one particular aspect of the invention, the straight lines are defined in a plane by a first parameter regarding distance between the origin of a reference coordinate system and the straight line and a second parameter regarding angle between an axis of the reference coordinate system and a line perpendicular to the straight line and passing through the origin of the reference coordinate system.
In one variant of embodiment, the method further comprises a step of selecting lines for which the sum of the amplitudes of the points of the set of points located at a distance from the line smaller than said minimum distance is larger than a predetermined threshold.
According to one particular aspect of the invention, the points of the 3D point cloud are geolocated using a geolocation device of the ground-penetrating radar.
According to one particular aspect of the invention, the radar signals are acquired for a plurality of planes in the region of ground.
Another subject of the invention is a ground-penetrating radar comprising at least one transmit antenna and at least one receive antenna and a device for detecting buried longitudinal structures in a region of ground, which is configured to execute the steps of the detecting method according to the invention.
Another subject of the invention is a computer program comprising code instructions for executing the method according to the invention when the program is executed by a processor.
Another subject of the invention is a processor-readable storage medium on which is stored a program comprising instructions for executing the method according to the invention, when the program is executed by a processor.
Other features and advantages of the present invention will become more clearly apparent on reading the following description with reference to the following appended drawings.
The scene of
One objective of the invention is to decrease the points 101 belonging to the clutter so as to better detect the presence of a longitudinal structure such as a pipe, a pipeline or piping.
On analysis of these three cross-sectional planes, it may be seen that the clutter 210 corresponds to noise that is coherent with the target to be detected, i.e. it may be of equivalent amplitude to a target in a radar image. However, the clutter 210 is not spatially coherent as the radar is moved along the scene. This is especially due to the fact that the clutter results from reflections of radar signals from objects that are very small or from interfaces between two types of ground that may vary spatially.
In contrast, it may be seen that an object such as a pipe is spatially coherent, i.e. the same amplitudes are found in the various images corresponding to various cross-sectional planes. In other words, the radar echoes of such a target are spatially correlated in various images corresponding to various cross-sectional planes.
This property is exploited by the invention to search, in a point cloud of the kind illustrated in
The first step 301 consists in acquiring radar images of a plurality of cross-sectional planes of a region of ground, in order to obtain radar images of the kind illustrated in
The radar images are acquired by means of a ground-penetrating radar comprising at least one transmit antenna and one receive antenna.
The radar is moved over the region to be imaged in order to perform a plurality of successive acquisitions. The raw signals measured are processed in order to generate a radar image each point of which has an intensity characteristic of the reflection of the signals from a buried target.
Step 301 of the method may for example be carried out by means of the acquisition method described in patent application FR 3111994 of the applicant. This application describes use of a so-called MIMO radar (MIMO standing for Multiple Input Multiple Output) having a plurality of transmit and receive antennas coupled to a demodulation algorithm that produces the radar image. Any other radar acquisition method allowing radar images of a plurality of cross-sectional planes of a 3D region may be used.
In step 302, a 3D point cloud is then constructed based on the various radar images of the cross-sectional planes of the scene. This step requires each image to be geolocated, in order to allow a position of the cross-sectional planes and therefore of the points to be deduced. Geolocation may be achieved by means of a device for receiving satellite radio-navigation signals, such as a GPS device, or of an odometer or more generally of any kind of locating device.
In step 303, a method for detecting aligned points in the 3D point cloud 402 is then applied in order to detect lines, corresponding to targets to be detected, and to remove points belonging to the clutter.
Step 303 consists in searching, in the 3D space, for a straight line such that a set of points of the 3D point cloud 402 are located in proximity to this straight line.
To do this, a representation of the straight lines in polar coordinates is used.
A first possible representation of a straight line in 2D space is given by the equation:
y=tan Φ0x+b0
A straight line may therefore be defined in a plane by the pair of parameters (Φ0, b0).
Similarly, the equation of a straight line in space is given by:
P=(x,tan Φ0x+b0,tan Φ1x+b1)
A straight line in space may be defined by the four parameters (Φ0, b0, Φ1, b1).
One drawback of the preceding representations is that the slope tan Φ0 tends to infinity when the straight line tends towards a parallel to the axis (Oy), this preventing representation of this type of straight line.
To solve this drawback, another representation is proposed in
ρ=x·cos(θ)+y·sin(θ)
with [OH] the segment perpendicular to the straight line D of length ρ and making an angle θ to the axis (Ox).
Each straight line D may therefore be represented by its parameters (ρ, θ).
Similarly, the equation of a straight line in space is defined by four parameters (ρ1, θ1, ρ2, θ2) as the intersection of two planes.
Step 303 of the method according to the invention consists, for each straight line defined by a set of four parameters (ρ1, θ1, ρ2, θ2) varying in a predefined set, in computing a distance, for example a Euclidean distance, between each point of the point cloud and the straight line. Next, points of which the distance to the straight line is smaller than a minimum distance are preserved and the number of obtained points is counted.
This process is iterated for all the possible straight lines of the set, according to a predefined sampling scheme.
Next, only straight lines for which the set of points detected to be in their proximity contains a number of points higher than a minimum number allowing a longitudinal structure to be characterized are preserved.
Thus, the method according to the invention allows one or more lines in the 3D point cloud to be detected, these lines approximately joining a set of points of the cloud or being close to this set. Each set of points in proximity to a detected line corresponds to a longitudinal structure that it is desired to detect in the point cloud.
A number of ways of carrying out the step 303 of detecting lines of the method according to the invention in 3D space are conceivable.
According to a first variant of embodiment, lines are detected using an approach in which successive planes are considered, the space defined by a pair of parameters (ρ, θ) being scanned with two orthogonal planes, the planes (XY) and (XZ) for example.
More precisely, a first search is carried out to determine the sets of points of the 3D point cloud that are close to the straight lines of parameters (ρ, θ) in a first plane (for example the plane (XY)). Next, among the points selected in this first search, a second search is carried out for straight lines of the second plane (for example the plane (XZ)).
The parameters of the straight lines in the plane (XY) are denoted (ρi, θj) and the parameters of the straight lines in the plane (XZ) are denoted (ρk, θl). These parameters are varied in the intervals [0; ρmax] and [0; θmax], respectively.
In step 801, a distance between each point of the 3D point cloud and the straight line defined by the parameters (ρi, θj) in the plane (XY) is computed. This distance is, for example, a Euclidean distance between a point and its projection onto the straight line.
In step 802, each distance computed in step 801 is compared with a threshold distance dist_th. If the distance is larger than the threshold distance, the point is not retained, otherwise step 803 is passed to, in which step the points P the distance of which is smaller than the threshold distance dist_th are selected to form a set of points close to the straight line (ρi, θj).
Next, among the set of points P determined in step 803 for the straight line (ρi, θj), the number of points close to each straight line of the plane (XZ) parametrized by the parameters (ρk, θl) are determined and these numbers are saved in a matrix A.
The process is iterated for all the straight lines of the first plane (XY).
At the end of steps 801-804, a matrix A giving, for each pair of straight lines parametrized by (ρi, θj, ρk, θl), the number of points of the point cloud in proximity to these two straight lines is obtained.
Next, in step 805, the maxima of the matrix A are sought with a view to determining the sets of points corresponding to nbmax lines, with nbmax the number of lines to be found in the 3D space.
In step 806, solely the straight lines and associated sets of points for which the number of points of the set is higher than a detection threshold corresponding to a minimum number of points allowing a longitudinal structure to be characterized are preserved.
In step 807, the obtained lines are drawn then, in step 808, the points close to these lines are selected.
To carry out step 808, i.e. to select the points of the point cloud that correspond to the lines obtained in 3D space based on lines determined for two projection planes (XY) and (XZ), for each pair of lines belonging to the two planes, respectively, the number of points in common of these two lines is sought, and if this number is higher than a predetermined threshold Tu, these points are added to a list of retained points.
One drawback of the first variant of embodiment illustrated in
According to a second variant of embodiment, the line-detecting step 303 may be carried out using an approximate algorithm that has a lower cost and complexity than the first variant.
This second variant consists in scanning the three projection planes of the 3D space, namely the planes (XY), (YZ) and (XZ) independently, so as to accelerate the process of line detection.
For each projection plane (indexed by the index n in
The indices of the sets of points close to each straight line are preserved in a matrix V(i,j) and the number of points of each set is preserved in a matrix A(i,j).
As in the first variant, the proximity of the points is measured by a parameter that defines a maximum Euclidean distance between the line defined by (ρi, θj) and the points of the cloud. The indices of points under the maximum distance are stored in the matrix V, which is used later on to allow the lines to be drawn.
The following step 902 consists in constructing, based on the indices saved in the matrix V, a matrix M that contains, for each projection plane, the highest occurrences of the number of indices. An occurrence corresponds to the indices of the points of the cloud that are close to the straight line obtained based on (ρi, θj). The maximum number of occurrences, which corresponds to the maximum number of lines to be detected, is defined by the parameter nbmax. The value of this parameter is such that the existing number of targets is covered. These occurrences are stored by decreasing size of the number of indices close to the lines defined based on (ρi, θj).
This matrix M is constructed using another matrix A which stores, in turn, the number of indices of points close to a line. In order to remove similar occurrences (nearby lines), a tolerance is established on the indices of the matrix A and, said matrix, zeroed for the set of indices in the margin of tolerance.
Thus, at the end of step 902, straight lines (of maximum number equal to nbmax) for which the number of points located in their proximity is highest will have been retained.
The operations 901, 902 are carried out for each of the three planes (XY), (YZ) and (XZ) independently.
Next, a voting system 903 is used to select the lines in the 3D space.
Step 903 first consists in establishing a threshold in order to preserve the highest occurrences among the nbmax occurrences. This prevents, subsequently, short segments of aligned points that do not correspond to a target that it is desired to detect from being drawn. In other words, there is a minimum number of aligned points corresponding to a longitudinal structure.
Next, the matrices M obtained for the 3 planes are combined, for example by retaining the sets of points present in at least two of the three projection planes.
The criterion of selection of the indices of points that define a line finally retained in the 3D space is given by:
S(nb*)=(S(nb*)XY∩S(nb*)XZ)∪(S(nb*)XZ∩S(nb*)YZ)∪(S(nb*)XY∩S(nb*)YZ), with S(nb*)XY, S(nb*)XZ, S(nb*)YZ the sets of indices of points obtained for each of the three projection planes for the occurrence nb*.
Once the set of indices of points has been determined for each occurrence nb* (corresponding to one straight line), the line that best fits the corresponding set of points is drawn, then the points of the 3D point cloud that are closest to this line are drawn by again computing a distance between all the points of the cloud and the line finally retained.
In one variant of embodiment, the method described in
In another variant of embodiment, a step 304 of filtering the points of the 3D point cloud is carried out before the line-detecting step 303.
The filtering step 304 is carried out using a method for detecting aligned points in the 3D point cloud, in order to remove points belonging to the clutter.
Step 304 is a step of searching for similar unit vectors, or in other words unit vectors that are substantially collinear between pairs of points forming thus a line.
The algorithm of
The algorithm of
An index i is initialized to 1 corresponding to the number of lines to be detected in the point cloud. In step 502, the method continues if i is lower than or equal to Nr the maximum number of lines that it is sought to detect.
In step 503, the first point of the list (that of highest amplitude) is selected. This step is illustrated in
In step 504, the set of points that are aligned or substantially aligned with the selected point 601 is sought.
To do this, first all the unit vectors having as start point the selected point 601 and the direction of which is given by one of the other points of the point cloud are determined. This step is illustrated in
The unit vectors are given by the formula:
where pk is the start point 601 and p; is the end point.
Among all the computed unit vectors, those that are substantially collinear with a predetermined margin of tolerance are retained. This step is shown in
There are a number of possible solutions that may be used to determine the set of substantially collinear vectors.
A first solution consists in executing the following steps.
Firstly, a rotation of 180° is applied to unit vectors a coordinate (for example the x-coordinate) of which is negative:
if ujkx<0→ujk−−ujk
This step allows all the unit vectors almost aligned in a common direction to be obtained.
Next, each component of ujk is approximated to n significant digits. For example, n is set equal to 2 or 3. Use of n>3 is recommended only in a very noisy scenario.
Next, for each component, the dominant value among all the approximated unit vectors is computed. In other words, the largest set of vectors having the same approximate component value is sought.
Next, the results obtained for the three components x,y,z are combined in order to determine points that are almost aligned, with a tolerance given by the approximation by the number n.
One possible combination consists in collating all the points present in at least two projection planes. In other words, if the sets of indices of the points of the cloud with the most recurrent value of the components x,y,z are denoted Ix,Iy,Iz, respectively, the set of almost aligned points is given by the following formula:
I=(Ix∩Iy)∪(Ix∩Iz)∪(Iy∩Iz), where ∪ designates the union operator and ∩ designates the intersection operator.
Another possible approach consists in measuring the angle between each pair of unit vectors and in preserving all of the points for which the absolute value of this angle does not exceed a predetermined threshold close to 0 and dimensioned to accept a certain tolerance in the alignment of the points.
A third possible approach consists in converting the unit vectors into spherical coordinates (r, θ, φ) and in introducing a tolerance into the angular variations of the angular components (θ, φ).
To do this, a histogram of the values of each angular component (θ, φ) is computed. The histogram is defined by an increment that gives the desired tolerance in the angular variation.
Next, the most recurrent values in these two histograms are determined by taking into account phase ambiguities (modulo π).
The sought set of almost aligned points corresponds to the intersection of the points having the most recurrent values of the two angular components, respectively.
At the end of step 504, a set of points that are almost aligned with the point selected in step 503 is obtained.
Next, in step 505, the number of points of the obtained set is compared with a threshold Nmin that is a minimum number of points that may be considered to belong to a target. This threshold is set depending on the type of structure that it is desired to detect. In the case of pipes, pipelines or piping, this threshold allows detection of objects of small size that may rather belong to the clutter to be excluded. The threshold value Nmin especially depends on the resolution of the movement of the radar. If the measurements carried out by the radar are very spaced apart spatially, the value of the threshold Nmin may be set very low. If in contrast the measurements are not very spaced apart, this value may be set higher. The value of the threshold Nmin also depends on the size of the 3D region scanned by the radar and on the length of the structures to be detected. The value of the threshold Nmin is for example set in the interval [10; 50].
If the number of points of the set is strictly higher than Nmin then, in step 506, all the aligned points are preserved, these being associated with a detected structure. These points are then withdrawn from the list to be processed and the index i of the list is incremented (step 508) to pass to the following point of highest amplitude among the remaining points.
If the number of points of the set is lower than or equal to Nmin then the point selected in step 503 is removed (step 507) from the list to be processed and is considered to belong to the clutter—it is therefore filtered from the 3D point cloud.
Steps 505 to 507 are illustrated in
The various modules DET,SB,NP,SVU may be produced in software and/or hardware form, notably using one or more processors and one or more memories. The processor may be a generic processor, a specific processor, an application-specific integrated circuit (ASIC) or a field-programmable gate array (FPGA).
Number | Date | Country | Kind |
---|---|---|---|
2202800 | Mar 2022 | FR | national |