Automatic Subsurface Map Generation Based on GPR Data Processing

In the domain of subsurface scanning, an abundance of GPR data are typically acquired by exploiting technologically advanced Ground Penetrating Radar (GPR) devices. However, the processing and representation of these data in a concise manner that can be easily interpreted by non-experts is still an active research topic. To this end, the paper at hand introduces a new method for automatic registration of sequences of GPR data into a common reference frame. It starts from an elementary GPR scanning session, relied on the existence of a GPR antenna array in a structured device, and proposes a sequential registration approach for multiple elementary scans, by exploiting a similarity measure, that utilizes the multilevel discrete wavelet transform decomposition. The data are processed in the 3D-space to achieve convergence during merging the scans, formulating thus a 3D map of the subsurface. In this map, data of selected depth range are top-down projected to formulate C-Scans in a parametric manner that can be easily interpreted by non-experts. The method has been evaluated in three different datasets and exhibited promising results in the task of automated GPR scans registration and corresponding large-scale subsurface map generation.


I. INTRODUCTION
Subsurface scanning is essential for several application domains including construction, inspection, mine detection and others. In this scope, recently, several technological breakthroughs took place in the field of subsurface scanning machinery [1]. A key milestone in this progress is the development of GPR antenna array systems, which encompass multiple antenna dipoles, aiming to reduce scanning time and human effort involved in the underground inspection procedure. In order to account for the impact of directivity of transmitterreceiver antennas on radar sensing, it is common for such multi-arrays to include dipoles on various configurations. An example of such device is the IDS Stream C [2] which consists of 9 channels on a HH polarization (horizontal transmit and horizontal receive) and 23 channels on a VV polarization (vertical transmit and vertical receive). The former can collect information regarding mostly transversal utilities while the latter regarding mostly longitudinal utilities with respect to system's forward direction. Work done in [3] demonstrates the significance of gathering data from different antenna arrangements towards gaining a better understanding of the inspected area.
In the GPR processing literature, the mapping term refers to the transformation of gathered temporal or frequential data 1 All authors are with the Centre for Research and Technology Hellas -Information Technologies Institute (CERTH / ITI), 6 th Km Charilaou-Thermi Road, Thessaloniki, Greece, 57100 {eskartad}@iti.gr to a spatial arrangement that can be easily interpreted from the experts. Due to specific geometry of scanning sessions, consistent utilities are expected to shape hyperbolic signatures on temporal 2D data. Based on this assumption, the generation of geometrical representations of the underground has been approached in two ways from the GPR community: either as a migration or as a computer vision detection problem. In [4] a comparative study between various migrations methods has been provided, and similarities and equivalence of them has been derived. From the computer vision perspective, gaining spatial understanding has been approached as a hyperbola detection problem. Utilizing hand designed features as robust descriptors of the hyperbola orientation and performing sparse decomposition based on a dictionary have been presented in [5] and [6] respectively. Compared to works originating form this field, approach presented herein introduces, on the top of existing methods, the notion of combining multiple scans in order to generate a global spatial representation of the underground. In the robotics literature the same term, mapping, refers to the problem of a mobile robot that is placed at any unknown location in an unexplored area and incrementally builds a global, consistent map of the environment. GPR data have been only recently involved in the process of producing such a map by a robotic device. Work presented in [7] is to our knowledge the first approach towards a completely autonomous surface/subsurface mapping. A rover navigates on the surface of the area under inspection and tows a GPR antenna. The data acquisition frequency is adjusted by a wheel encoder and acquired signals are georeferenced using GPS positioning. In the work presented in [8], GPR measurements aggregated in a global map are utilized in combination with a dedicated Localization GPR antenna (LGPR) in order to complement existing approaches of LIDAR, camera, and GPS/INS sensors localization. An important limitation of these methods is that they require the existence of surface sensing equipment (i.e. vision sensors, GPS, etc.), the subsurface map is amalgamated with the surface one and, thus, subsurface mapping derives indirectly. Compared to works originating from this context, presented approach introduces generation of a map based solely on GPR measurements, and is therefore a first step towards a completely autonomous robotic underground map. The robotic setup that was utilized in current work is illustrated on Fig. 1.
The paper at hand tackles the issue of registering the data acquired during different scanning sessions, without considering any prior transformation of the GPR antenna motion. At the core of the approach presented herein is the Discrete Wavelet Transform (DWT) decomposition for the matching procedure. The DWT is a powerful tool applied in signal processing domains proved suitable for the frequency analysis of nontransient signals, such as GPR signals. Successful application of DWT for filtering and denoising purposes of GPR data has been reported in [9], [10], while the authors of [11] applied wavelet analysis so as to perform GPR data classification. The proposed method builds upon DWT so as to create consistent 3D maps of the subsurface, where acquired scans are registered and are subsequently top-down projected in order to reveal the existence of buried utilities.

A. Elementary GPR scanning session
A single GPR measurement comprises the transmission of an excitatory pulse at a p = (x, y, z = 0 + ) location above the ground and subsequently recording the back scattered electromagnetic reflections. Since those correspond to the changes in the electromagnetic properties of the materials that the transmitted pulse penetrates, it is possible to infer from these reflections useful knowledge regarding the underground layout. A typical GPR underground inspection process involves movement of the device along a specific trajectory and performance of multiple GPR measurements at discrete, evenly placed locations. This approach resembles the Synthetic Aperture Radar technique, that enables imaging of areas with a surface of many square meters, utilizing an antenna with a diameter of a few centimeters.
The term elementary scanning session introduced herein refers to a GPR inspection procedure on which two geometrical restrictions are imposed regarding the scanning path: (i) the entire path lies on a single plane, the z = 0 plane (i.e. the altitude level of the ground surface remains constant) and (ii) the trajectory is a straight line. When a multi-array system ( e.g. Stream C provided by IDS Georadar) survey meets the aforementioned conditions, the recorded data are arranged on a 2D regular grid. The X axis of the grid coincides with the multi-array main axis while the Y axis is defined by the direction of the scanning straight line. Δx step is the spacing of consecutive channels, while Δy step depends on the spatial sampling frequency of the considered survey. An indicative illustration of a single elementary scanning session is depicted on Fig. 2, where the utilized device comprises 7 channels separated by 4.4cm distance by each other (X axis) and during the survey 21 measurements are performed with a spatial frequency of 3.2cm (Y axis).
Data gathered at any point p = (x, y) that lies on the z = 0 plane are not completely defined by the (x, y) coordinates of the specific points. Taking into consideration that (a) a single GPR measurement is significantly affected by the antenna directivity and (b) the directivity of all antennas remains fixed with respect to the multi-array frame for all measurements along any scanning session, it is inferred that collected data at any point is actually a function of the direction of the straight line defining the elementary scanning session.

B. Registration of two elementary scanning sessions
In the next step, we provide a strategy for registering elementary scanning sessions given their corresponding 2D grids, by determining the relative transformation between them. The registration of elementary scanning sessions is examined here under the restriction that those sessions are performed on the same direction, since data gathered at consecutive sampling points that fall on a straight line depend on the direction of that line, defined by a #» v vector. On the central image of Fig. 2 an example of two 2D grids that organize data collected from two separate sessions performed in the same direction Deviations in scanning directions as illustrated in the right image of Fig. 2 cannot be handled from the proposed method. The relative angle α between the directional vectors of the two sessions is equal to 45 • and consequently the data gathered is not directly comparable to each other, despite the fact that some of the measurements of the two sessions are performed at the exact same locations. The set of points defining a 2D grid can be formulated as: Assuming L is a 2D linear affine transformation, it can be defined as a function f : A 2D → A 2D that maps any point p ∈ A 2D of the 2D space to f (p) ∈ A 2D . The generalized expression for L is given by: where α is the rotation performed around the axes origin and is the translation vector with t x and t y standing for translation along the X and Y axis respectively. The alignment of grids representing parallel scanning sessions does not involve any rotations, consequently α = 0 • . The displacement between the two grids can be calculated in a discrete point-wise manner. Consequently, translation vector itself is a discrete variable and transformation L takes the following form: with N ch standing for the number of channels on the multiarray and N s1 and N s2 standing for the number of scans captured along the trajectories of the two sessions. Estimation of k and l is tackled as a global optimization problem that seeks the optimal transformation between the examined grids that minimizes their resemblance distance (defined in Sect. II-C). Each separate 2D translation, defined by such a (k, l) pair, when applied on the points of the two scanning grids results in a set of matches of GPR measurements captured at the same location during the two considered sessions. Let M be the set of measurement matches between two scanning grids and S be a measure of similarity between two 1D signals. Then (k, l) is derived as follows: Hence, registration is performed by computing the (k, l) values that minimize a loss function. The loss function is defined as the similarity of the resulting matches when the transformation described by the (k, l) values is applied on the G set of points of one of the two grids. There is no indication of the loss function being a convex function and consequently no traditional convex optimization strategies can be utilized. A full search on the domain of (k, l) is carried out with k ranging from −(N ch − 1) to (N ch − 1) and l ranging from −(N s1 − 1) to (N s2 − 1). The registration of multiple scanning sessions with the appropriate transformation on a global grid concerns merging those sessions in a single grid. Herein, we tackle the issue dynamically, by iteratively breaking it down to a pair of sessions on which registration is solved recursively. Given that registration for two sessions G 1 and G 2 is a function of the form f : (G 1 , G 2 ) → G 1,2 , then registering N sessions, G 1 , ..., G n , can be broken down as registration between the first N − 1 sessions G 1 , ..., G n−1 and G n . The same analysis is applied for the registration of the first N − 1, N − 2, ..., 2 sessions always solving a two sessions registration task at each step.

C. Similarity measure based on Multilevel Discrete Wavelet Transform decomposition
Here a measure of the similarity of two GPR measurements (A-Scans) is introduced, which utilizes the multilevel Discrete Wavelet Transform (DWT) decomposition of a 1D signal. In one level of the DWT the entry signal y[n] (A-Scan) is convolved with a high pass filter h[n] and a low pass filter g [n] and is in both cases subsequently down-sampled by two, producing the detail and the approximation coefficients respectively. On a multilevel DWT approximation coefficients are provided to the next level and the same procedure is repeated. The maximum number of levels the above procedure can be repeated is determined by the length of the original signal and the number of coefficients of the applied filters. From one level to the next one, due to the filtering and the subsampling, the frequency band is reduced to half and at the same time the frequency resolution is doubled. The high-pass filters h[n] comprise the child wavelets of the mother wavelet function and are computed from: where j is the scale parameter accounting for the level of resolution and k is the shift parameter. Low-pass g[n] filters are produced by a function which is referred to as scaling function. The wavelet family that is chosen here is Daubechies-6 due to the resemblance of the mother wavelet function to the emitted GPR pulse [12]. Once the decomposition reaches the maximum level, the detail coefficients of all levels along with the approximation coefficients of the maximum level constitute a sparse representation of the signal, which is perfectly invertible and captures information both for frequency and time domain. All coefficients formulate a feature vector. In order to capture the similarity of two GPR measurements, multilevel DWT decomposition is applied on them and Euclidean Distance between derived feature vectors is calculated. The expression of the registration between two scanning sessions takes its final form in Eq. 6. A visualization example of the estimated loss function for all valid transformation pairs is exhibited in Fig. 3.

D. Migration and 3D Map Generation
Once all the measurements captured during multiple scanning sessions have been registered in a common grid, it is possible to produce merged GPR images of the underground (either B-Scans along the scanning direction or T-Scans that correspond to transversal direction, see Fig. 2). Although captured measurements are now arranged on a 2D spatial grid (X, Y axes) they still represent time signals and the resulting GPR images are space-time images. The final step towards a complete 3D map is the transformation of received temporal data to spatial depth information (Z axis), a procedure commonly referred in the literature as migration [4].
Underground utilities are known to produce hyperbolic signatures on GPR space-time images, which are due to the travel time of the transmitted pulse from consecutive scanning positions to underground obstacles and back. In order to produce focused versions of the captured images, migration techniques introduce a reversal of the actual sensing procedure, which can be modelled with the Exploding Source Model (ESM) that introduces the following assumptions: (i) underground scatterers are considered sources of the received measurements, (ii) propagation velocity is set to half the actual velocity and actual back and forth travel is replaced by one way travel from scatterers to GPR receiver. According to the ESM assumption, all underground sources explode simultaneously at t = 0. Migration techniques utilize the distances between the assumed source positions and the actual scanning locations as well as the received temporal measurements on those locations. An extrapolation of the wave field at the explosion time (t = 0) is performed and a spatial image of the underground is obtained. Here the Kirchhoff Migration variation is applied in a 2D fashion, where for a B-Scan image, along Y axis, it assumes scatterers on the (Y, Z) plane and for a T-Scan image, along X axis, assumes scatterers on the (X, Z) plane.
As a result of the proposed registration and migration method, captured data can be finally arranged on a 3D spatial regular grid. In order to provide a meaningful visual representation of the obtained map one last step is required Fig. 4. C-Scans of the three test sites which is related to the top down projection of the focused measurements, originating from scatterers at various specified depth ranges. This is referred to as a C-Scan and constitutes a means of inspecting the underground structure in parametric depth range. The implementation of this last step entails two subordinate routines; The first one is the isolation of the maximum value recorded for all the (x, y) 2D measurement positions at the specified depth range and subsequently the interpolation of the resulting surface at the same step for X and Y axis (Δx ≡ Δy) so as to generate a Cartesian grid. In this way, the C-Scan preserves the ratio of the actual dimensions of the area under inspection as well as the sizes of the underground utilities.

III. EXPERIMENTAL EVALUATION
The test sites and the scanning procedure have been carefully chosen based on the following complementary remarks. On the one hand, for any elementary scanning session performed by a multi-array antenna, such as Stream C, the size of the B-Scans is parametric and for a fixed spatial sampling frequency depends solely on the length of the scanning route. Consequently, gathering more data regarding underground utilities lying transversal to the scanning direction can be achieved simply by increasing scanning path's length. On the other hand, the size of the T-Scans is constant and severely limited, since it is defined by the actual length of the multiarray channels of the utilized GPR antenna. In the case of the Stream C device the 23 VV channels extend in a straight line of 1m, with channel 1 located at 0cm and channel 23 located at 100cm. No matter how long a single elementary session will be, information regarding utilities lying longitudinally to the path direction will be restricted by length of the antenna. Consequently, at a first level, registration, as presented in Sect. II (i.e. in the sense of merging multiple parallel elementary scanning sessions performed along the same lengths, with an offset between them perpendicular to their direction) can significantly benefit information gained from T-Scans. Of course, at a second level, when registration is followed by migration and 3D mapping, information regarding transversal utilities is also enhanced.
The functionality of the proposed method has been tested on three separate cases. In order to demonstrate in a direct way the expected benefits, the scans on the selected sites were  performed so that the underground utility (e.g. pipe) would lie along the scanning direction. At each test site, 7 sessions were performed parallel to each other, with the same start and end point. Moreover, laterally they had an offset of 0.5m to each other ensuring enough overlap of the antenna between consecutive scans. This means that the ground truth discrete translation for two adjacent scans should be considered equal to t x , t y = (11, 0), since the scans are deliberately performed so that there is 100% overlap between them along the scanning direction (Y axis) and 50% overlap between them along the device direction (X axis) and channel 11 of the Stream C antenna is located in the middle of the multi-array. The estimated translations with the proposed registration method between the scans for the three test sites are provided on Table I. It should be noted, that for the last site consecutive scans were performed with a translation between each other in the opposite of the device direction, which is successfully captured by our method, as the negative sign in the estimated t x coordinate indicates. In Fig. 5 Fig. 4 where underground structure is revealed. The 2-D interpolation step that is interleaved guarantees that the size of the C-Scans retains the actual ratio of the dimensions of the area under inspection. As indicated by Fig. 4, automatically generated 3-D representation captures sufficiently the utility lying along the scanning direction in all the three cases, when top-down projection on plane is applied.

IV. CONCLUSIONS
In this work a registration framework of several GPR scanning sessions has been proposed. The method tackles the issue of automatically merging GPR data that have been acquired from a common array of GPR antennas that follow a specific motion pattern. The notion of the elementary scanning session has been introduced and on top of it a matching and registration method, based on a similarity measure, has been presented. The data are processed in the 3D space, however visualizations for the GPR operators are provided as topdown projections, revealing the existence of buried structured utilities in the subsurface. The method has been tested on 3 different test sites isolating utilities, e.g. pipes, in different orientations. The current limitation of the propose method is that the GPR routes during data acquisition should be roughly parallel in order for the method to converge at a valid solution. In our future work, we plan to extend registration of T-Scans from any orientation providing more degrees of freedom to the GPR operators while scanning an area of interest.