3D MAPPING OF BURIED PIPES IN MULTI-CHANNEL GPR DATA

: Ground Penetrating Radar (GPR) allows a non-destructive analysis of the subsurface by using electromagnetic waves. GPR is used in application fields such as archeology and civil engineering, where the detection of buried objects is in demand. Such objects, e.g. pipes, lead to a disturbance of the propagation of the radar signal in the underground. A manual detection of disturbing objects can be both time-consuming and tedious, depending on the number of GPR images to be investigated. Hence, automatic methods should be considered. In this study, an object-oriented image analysis approach for the automatic detection and mapping of buried pipes is presented and evaluated. As dataset, measurements of the multi-channel system Stream C are considered.


INTRODUCTION
Techniques, using Ground Penetrating Radar (GPR), allow nondestructive investigations of the subsurface using electromagnetic waves. GPR is used in many different applications such as archeology (Burds et al., 2018), civil engineering (Baek et al., 2018), geophysical surveys and safety (Hu et al., 2018), pipe detection (Zhou et al., 2019;Yamaguchi et al., 2020), forensics and landmine detection (Daniels and Allan, 2009). In detail, GPR systems enable the detection of buried objects, which disturb the propagation of the radar signal in the underground. In essence, all dielectric discontinuities are detected by GPR. These discontinuities can be of natural origin like stones or of human origin such as ruins of historical buildings (Daniels, 2004;Jol, 2008;Goodman and Piro, 2013). In general, GPR images can be generated by a spatial arrangement of measurements of a single-channel system. Consequently, the investigation of the area of interest has to be conducted in a very accurate way. Moreover, in the case of large-scaled areas, this technique can be quite uneconomically. In contrast, multi-channel systems deliver two-dimensional image matrices with each measurement. Focusing on measurements consisting of a large number of GPR images, manual detection of disturbing objects can be both timeconsuming and tedious. Therefore, it seems obvious to consider automatic analysis methods. In this study, an approach for the fully automatic detection of disturbing objects in multi-channel GPR images is presented. Furthermore, as final processing result, the approach delivers visual information about the surveyed underground, given by 3D maps containing locations of buried pipes. For the detection of buried pipes an approach, which bases on the recognition of hyperbolas in GPR images, is used (Dou et al., 2016). The longranged aim of our approach is to develop an algorithm to support the operator in field by calculating a 3D pipe detection map. This algorithm needs to be suitable for commercial multichannel GPR systems, such as Stream C or Mala MIRA. The paper is organized as follows. In Section 2, the test site and the dataset are described. Section 3 contains the pre-processing of the GPR images. The methodology for the detection of buried pipes is given in Section 4. The evaluation of the resulting 3D map takes place in Section 5. The conclusions and an outlook are given in Section 6.

DATA ACQUISITION
In this Section, the Stream C sensor system is outlined. Moreover, the test site and the dataset available are described.

Sensor
The GPR data were recorded in 2020 with the multi-channel GPR system Stream C (IDS GeoRadar), which consists of 34 antennas, 24 of them vertically polarized with a spacing of 4.4 cm and ten horizontally polarized with a spacing of 10 cm. In this study, only the vertically polarized dipoles are used. This means that 23 uniformly spaced channels are available with one scan. Hence, the scan width of the system is 1.0 m. The central antenna frequency is 600 MHz, leading to a suitable depth range for underground utility network detection (Gabrys and Ortyl, 2020). Since Stream C is a multi-channel system, a 2D image can be acquired directly with each measurement.

Location and Data
The test site was built as part of the Detectino project (Detectino GmbH, Hildesheim, Germany) on the grounds of the University of Frankfurt, Germany (Figure 1, left). It enables the non-destructive detection of pipes and cables in the underground. More information about the Detectino project and detailed descriptions of the test area can be found in Naser and Junge (2010). As optical reference data, a digital orthophoto containing the test site is shown. This DOP20 was acquired in 2019 and downloaded from the hessian state office for soil management and geoinformation (https://gds.hessen.de). The test site is divided into two major areas consisting of two troughs, each measuring approximately 13 x 50 m. Trough 1 has a depth of 1.5 m. Trough 2 is 3 m deep and further hydraulically isolated from the surrounding area, allowing analyses with different water saturations. The two areas are filled with seven different subsoil materials such as sand, gravel or construction rubble. In addition, three different surface materials are used (Figure 1, right). In summary, this test site provides a wide variation of geotechnical properties (Naser and Junge, 2010). Concerning the objects, in total 40 different pipes and cables are buried in the subsoil. The pipes are made of different materials such as PE (polyethylene), steel, vitrified clay and concrete. Moreover, the pipelines have different diameters and inclinations. In this paper, a part of the whole test area is investigated (Figure 1, left, zoomed part: section D). The soil of this area consists of sand and the surface material is interlocking plaster (see Figure 1, right). Additionally, some manhole covers are visible on the surface. The examined area has a size of about 15 m x 13 m. This dataset contains 13 adjacent tracks that have been recorded with the GPR sensor ( Figure 1, left, zoomed part). The data acquisition in the sensor moving direction was carried out at walking speed, corresponding to approximately 4 km/h. Each depth measurement contains 512 samples. Depending on the ground conditions, an investigation depth of 3.2 m can be achieved by assuming an average velocity of 0.1 m/ns. For each lane with a length of approximately 15 m, there are about 340 images available, which corresponds to a frame rate of about 4 cm in the sensor moving direction.

Reference Data
As reference data, a site plan of the buried pipelines at the test site is available, which is not in scale ( Figure 2, left). Based on this sketch, a 3D map containing the locations of eight prominent pipes was created ( Figure 2, right). This was manually accomplished by visually analyzing the GPR data.
The hyperbolas were matched with the existing pipelines in the sketch and were manually picked at the apex. Figure 2 shows for example a set of three short pipes consisting of PE (red) perpendicular to the many longitudinal pipelines. All of them are located side by side in a depth of about 0.7 m with diameters between 110 mm and 225 mm. Two longitudinal pipes colored in black and blue are also made of PE with comparable diameters to the red pipes. The blue colored pipe is located in a depth of 1.2 m, whereas the black colored pipe has a depth of 0.7 m. One of the marked longitudinal pipes is made of steel (DN168) and is inclined (green). The depth of this pipe varies from 0.6 m to 2 m (Naser and Junge, 2010). The pipe colored in cyan is located on the edge of the test site and is also inclined, with depths varying from 0.3 m to 1.4 m. The pipeline colored in orange is the only one that runs diagonally at a depth of about 0.9 m. The 3D reference map including eight pipes ( Figure 2, right) is used in Section 5 to evaluate the automatic detection results. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France

PRE-PROCESSING OF GPR IMAGES
As initial processing step, the GPR data has to be prepared both in a geometric and radiometric way.

Geometric Processing
In total, the data cube consists of 512 x 299 x 340 pixels, where the first two dimensions represent the depth and the stack of channels (13 tracks x 23 channels) and the third dimension stands for the measurements per track. A 2D surface image of that cube is shown in Figure 3. The geometric displacements of the individual tracks can easily be observed. Therefore, at first, the even measuring tracks are mirrored. In general, various sources of error, which impose linear and non-linear geometric distortions on the data, have to be corrected. These are, for example:  Varying speed of the sensor during recording;  Deviations from the exact recording position of the sensor, concerning the planned driving track of the sensor array;  Variations in sensor orientation. An elaborate correction of these geometric errors for the system Stream C was investigated by Gabrys and Ortyl (2020). Since, in our study, we could not make use of the recorded coordinates, a simple georeferencing method was used. The shifts between the individual tracks are manually identified and corrected. For this, corresponding points between neighbored tracks are identified and the geometric distance between them is corrected by applying a shift operation. The result is given in Figure 4. In the future, if GPS coordinates are available, it is also possible to utilize open-source software for multi-channel GPR processing, such as presented in Wunderlich (2021). Regarding our long-ranged aim of study, it is significant to apply an automatic registration of the GPR data. For this purpose, the MATLAB based GPR processing software (Wunderlich, 2021) seems very promising.

Image Processing
Geometric correction is followed by a time-zero correction of the data, which is accomplished by using the index of the maximum peak plus half the pulse width. Afterwards, a sigmoid function is applied, which dampens the ground bounce in order to emphasize the signals from the deeper soil layers. Additionally, the data is scaled to amplify the signals from the deeper soil layers (cubic, maximum amplification = 8). The different brightnesses of the individual channels in the data ( Figure 5, left) were eliminated using dewow filtering (Jol, 2008). Further, the ground bounce and the noise dominated deeper parts were cut off. In Figure 6, the result of the pre-    The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France processing is exemplarily shown for one B-Scan. The individual pipes at different depths are clearly visible by several hyperbolas. In Section 4, the automatic detection of such hyperbolas is addressed, whereby the resulting pre-processed GPR data are used as input data.

METHODOLOGY
We favor an object-oriented image analysis approach, in which there is a transition from the iconic level in both 2D and 3D space to more abstract objects.
In GPR measurements, tube-like objects appear as hyperbolas. The first step is to automatically determine the maximum of the hyperbola in the GPR data. This maximum can then be used as a hypothesis for the location of a pipe. The set of 3D hypotheses forms clusters of points, that are approximated to form lines, which in turn are considered as pipe hypotheses. For evaluation, the line approximations are represented as cylinders in 3D space.

Hyperbola Approach
The approach bases on a 2D recognition of hyperbolas. Building on this, a 3D approach, which attempts to discern the linear structure of the buried pipes by clustering clues in 3D space, is applied. The proposed method is presented in Dou et al. (2016) and was implemented fragmentary. The GPR data is arranged in a 3D data cube, which is processed in two directions, as visualized in Figure 7. The reason is that pipes that run parallel to the direction of data acquisition at the same depth are interpreted as background in the GPR data.
Procedures that react to a change (e.g. background removal) fail in this case. If the pipe runs orthogonally to the recording direction, a signal is generated over the entire width of the measuring aperture, as is caused by soil stratification. If the data is processed in both directions, a hyperbola-like signal of a pipe will appear in the GPR data in at least one processing direction. Depending on the processing direction, a B-Scan is defined further on by the number of samples of an A-Scan and the measurements along a sensor track (see Figure 7, processing direction 1), or in the other processing direction by the number of samples of an A-Scan and the number of channels (processing direction 2). An example of B-Scans in different processing directions is shown in Figure 8. The image on the left represents the 40 th plane in driving direction. This B-Scan shows signatures of two pipes orthogonal to the driving direction. The image in the middle shows the plane 143, the right image the result of plane 207 perpendicular to the driving direction. Those corresponding B-Scans of the other processing direction clearly show the corresponding hyperbola signatures. Each B-Scan in the data cube is analyzed to detect hyperbolas as hypotheses of pipes. The necessary threshold is determined by an adaptive operation in the image part below the ground bounce (Bradley and Roth, 2007). In doing so, large-scale fluctuations in the gray values of the background are taken into account. A result of the binarized image of Figure 6 is shown in Figure 9, where different hyperbolas can be identified. To discriminate these hyperbolas from other objects in the image, a filter operation is used. First, a connected component algorithm is applied to generate region objects (Shapiro, 1996) and to fulfil the transition from the iconic to the symbolic level. Afterwards, a set of features is calculated for each region. To reduce the number of candidates, we apply a filter aiming at the number of region pixels. Noise and regions generated by soil layers are eliminated. As a result, we get regions of interest for the position of hyperbolas in the GPR data produced by a pipe. The filter operator is a simple brute forward approach. In the future, it is planned to improve this filter by considering a Support Vector Machine (SVM) approach. The remaining region objects are skeletonized (Zhang and Suen, 1984). As mentioned, pipes located parallel to the driving direction cause hyperbolas in the GPR data. Finally, the coordinates of the maximum of the detected hyperbolas can be used as a hypothesis for the position of the pipe. The minimum value in y-direction of a skeleton line, representing a hyperbola, and the matching x-value build the hypothesis. Migration    The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France effects are not taken into account. A result of two consecutive B-Scans is given in Figure 10. This approach suffers from false suggestions in the case of overlapping hyperbolas. Concerning overlapping hyperbolas, an improvement was investigated that is visualized on synthetic hyperbolas ( Figure 11). First, so-called branch points and endpoints are searched for in the skeletonized data (blue line) of a region object using a morphological operator. In Figure 11, right, the branch points are colored cyan and the endpoints yellow. The skeleton line is split at the branch points and the remaining partial lines (green) are traced from the endpoints. The entire skeleton line is thus split into several nonoverlapping parts. For these, the maxima are now sought as mentioned before. These maxima are highlighted in red in the examples (Figure 11, right).

Generation of 3D Hypothesis
The results of processing each B-Scan as shown in Figure 10 deliver hints to 2D pipe location coordinates. The third dimension is defined by the channel number, respectively the measurement position. The accumulation of all these point coordinates results in a 3D data cube.
To produce a more homogeneous result, a morphological filter called majority is utilized in a next step. Majority is defined as follows: Set a pixel to 1 if five or more pixels in its 3 x 3 neighborhood are 1. Otherwise, set the pixel to 0. The results of the morphological filtering of the data of both processing directions are considered as input for a fusion process basing on a logical OR function. The fusion leads to a 3D data cube containing hypotheses for the position of pipes ( Figure 12, coloring indicates measurement position).

Line Approximation
In further processing steps, this 3D data cube is interpreted as a point cloud. We apply a segmentation of the point cloud into clusters, with a minimum euclidean distance of a defined minimum distance between points from different clusters. Furthermore, a suitable minimum of cluster instances has been found by test runs. A result is shown in Figure 13. Each cluster is marked by a unique color. Subsequently, each cluster is approximated by a line. These lines are determined by multiple linear regression ( Figure 14, red colored lines). The generated line objects are a base for further analyses, for example to aggregate pipe networks.

Visualization of Pipes
In the final step, for a more realistic visualization, the generated lines are transformed into a cylinder representation similar to real pipe geometry. The extracted line represents the center line of the visualized cylinder in the 3D space. Here, the radius of the cylinders is equally chosen for all detected pipes ( Figure  15), since no further information on the diameter of the pipes is available.

EVALUATION
The resulting 3D map of the detected pipes ( Figure 15) is evaluated with the manually created reference map (Figure 2, right). For a better comparability, the cylinder representation is also applied to the reference map ( Figure 16). The result of the evaluation is shown in Figure 17. Detected pipes that could be matched to pipes in the reference map are marked in the same color as the reference (red, green, black, blue and cyan). The two other colors yellow and magenta highlight pipes, which are   The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2022 XXIV ISPRS Congress (2022 edition), 6-11 June 2022, Nice, France very prominent in the detection results, but cannot be found in the reference 3D model. In total, seven of the eight pipes in the reference map can be associated with pipes in the detection result. These include a set of three short pipes consisting of PE (red) perpendicular to the many longitudinal pipelines. In addition, the longitudinal steel pipe colored in green, which is inclined, can be clearly identified. Another longitudinal and inclined pipe at the end of the measurement cube (cyan) is also unmistakably identified. Furthermore, the two longitudinal pipes colored in black and blue can be matched with the reference data. However, in the case of the blue pipe, two detected pipes are in question (both colored blue in Figure 17). The only pipe from the reference map that cannot be assigned is the diagonal pipe colored in orange. One possible reason might be the diagonal orientation, which makes it more difficult for the detection algorithm.
Considering the depth and the material of the pipes, no dependency of the detection results can be recognized on the basis of the examples showed. However, it is clear that deeper located pipes show weaker signals and are therefore more difficult to detect. To further investigate the influences of material and depth of the pipes, additional reference data have to be created. Regarding the different diameters of the pipes, it can be seen that pipes with small diameters are not detected. An example is the additional fourth pipe of the set of three pipes, marked in red (Figure 2, left), which is the thinnest of the whole group with a diameter of 63 mm. This pipe shows only a very weak signal in the data, which made it impossible to pick it for the reference map. The detection result given in Figure 17 shows many additional pipelines, which are not included in the 3D reference map. Four of these pipes, which are all running longitudinal, are highlighted with colors. The two pipes colored in yellow are located at the end of the measurement cube, whereas the two pipes colored in magenta are located at the beginning of the cube. Possibly, these are the pipes in the yellow marked area in the site plan (Figure 2, left). The detection results indicate even many more pipes than are noted in the site plan. Visual inspection of the processed GPR data supports this observation. In summary, it can be stated that the reference pipes are detected in a robust way. The results look quite promising as a basis for ongoing research and optimizations.

CONCLUSION AND OUTLOOK
In this paper, an approach for the automatic 3D mapping of buried pipes in GPR data is presented. It bases on a hyperbola detection method and is applied to multi-channel GPR data acquired with the commercial system Stream C. As a result, a 3D map showing the visualization of the detected pipes is automatically created. This map is evaluated with reference data. Seven out of eight pipes in the reference map are identified with the hyperbola detection method. These detected pipes are orientated parallel and perpendicular to the sensor moving direction. The automatically generated detection result shows many additional pipes, which are not included in the reference. However, this is consistent with the visual impression when looking at the processed GPR data. The results are promising and motivate further investigations.
Concerning an object-oriented image analysis approach, we Figure 14. Approximation of point cloud clusters by lines. Figure 15. Visualization of detected pipes.  plan to optimize the detection of the pipes. The treatment of non-linear pipe structures (e.g. curves, branches) will also be considered. In addition, the possibility of aggregating pipe networks will be examined. Moreover, we will address the following topics: improvement of geometric processing, influence of migration and SVM for filtering hyperbola candidates.