Automatic modelling of urban subsurface with ground-penetrating radar using multi-agent classification method

ABSTRACT The subsurface of urban cities is becoming increasingly congested. In-time records of subsurface structures are of vital importance for the maintenance and management of urban infrastructure beneath or above the ground. Ground-penetrating radar (GPR) is a nondestructive testing method that can survey and image the subsurface without excavation. However, the interpretation of GPR relies on the operator’s experience. An automatic workflow was proposed for recognizing and classifying subsurface structures with GPR using computer vision and machine learning techniques. The workflow comprises three stages: first, full-cover GPR measurements are processed to form the C-scans; second, the abnormal areas are extracted from the full-cover C-scans with coefficient of variation-active contour model (CV-ACM); finally, the extracted segments are recognized and classified from the corresponding B-scans with aggregate channel feature (ACF) to produce a semantic map. The selected computer vision methods were validated by a controlled test in the laboratory, and the entire workflow was evaluated with a real, on-site case study. The results of the controlled and on-site case were both promising. This study establishes the necessity of a full-cover 3D GPR survey, illustrating the feasibility of integrating advanced computer vision techniques to analyze a large amount of 3D GPR survey data, and paves the way for automating subsurface modeling with GPR.


Introduction
The rapid expansion of cities has led to an increasing number of congested subsurface infrastructures. In modern cities, there are underground utilities (drainage, sewer, power supply, telecommunication, or gas), building foundations, archeological sites, and other subsurface structures. These structures belong to different organizations due to which the record drawings are prone to being misplaced or becoming outdated. Construction projects without detailed information about subsurface structures may cause urban hazards. Therefore, to maintain functioning cities, comprehensive records of subsurface structures are of vital importance, necessitating point-in-time subsurface surveys and modeling techniques.
Among many nondestructive testing (NDT) methods, ground-penetrating radar (GPR) has the advantages of penetrability, high resolution, high efficiency, and noncontact measurements. It has gradually gained popularity and has become the primary method for shallow subsurface imaging. Electromagnetic waves are produced by GPR instrumentation and then emitted from a transmitting antenna. Subsequently, the signal spreads out into the ground in a downward conical form and penetrates through subsurface layers. When the EM wave encounters any electrical parameter contrast in the ground, it is backscattered and received by a receiving antenna. The main electrical parameters of any medium, that is, permittivity, permeability, and conductivity, are generally a function of frequency. Contrasts in permittivity is the main cause of reflections of the radiated EM waves. The change in conductivity affects the absorption of the radar signal by the medium, while the variation in permittivity determines the variations in the impedance characteristics of the medium (Annan 2004;Benedetto and Pajewski 2015). Full-cover GPR measurements can depict the survey area in a 3D manner and supply a large volume of measurement data. Thus, full-cover GPR surveys substantially improve the interpretation of the subsurface environment.
A-, B-, and C-scans are three forms of GPR data presentation from one to three dimensions, respectively. An A-scan denotes a single waveform recorded by the GPR at a particular position. When moving antennas along a traverse, a set of A-scans form a vertical section through the ground, and this section is called a B-scan, as shown in Figure 1. C-scans map 3D data volumes comprising multiple data "slices", each of which is a horizontal plan for a plane at a certain depth. A C-scan is formed by stacking multiple B-scans collected in this horizontal plane, and then plotting the amplitudes of the recorded data at a given time (Goodman and Piro 2013). C-scans show reflection intensities from an overall perspective, while B-scans offer full waveform information in vertical sections. Integrated analyses of C-and B-scans can facilitate the complete utilization of the information they offer. However, data interpretation is always difficult. Currently, theinterpretation of GPR data is based on the waveform information presented in the radargram. The recorded GPR reflections consistof complex scattering noises, which distort the GPR reflections of the targeted objects, making it difficult to distinguish targetobjects from inhomogeneous subsurface environments. Thus, the GPR data interpretation relies heavily on the operators'experience and their knowledge, such that the degree of automation is less than satisfactory. Furthermore, the large data volumeproduced by full-cover surveys contribute to the demand for interpretation efficiency. Hence, automating the GPR datainterpretation process is an important research concern.
Many studies have been conducted on automating the recognition of GPR reflections from B-scans. These studies can be categorized into two groups: machine learning and deep learning based. Some machine learningbased studies and their research objectives are listed as follows: Ayala-Cabrera and Herrera (2014) combined a multi-subject and clustering approach to extract strong scattering regions to identify soil looseness; Pasolli and Melgani (Pasolli, Melgani, and Donelli 2009) extracted hyperbolic reflections with edge detection and used the Hough transform method to identify small cavities; Xie et al. (2013) extracted typical hyperbolic reflections using a support vector machine (SVM) to identify small cavities in road structures; and Luo and Lai (2020b) extracted typical reverberation patterns with the pyramid method to identify larger-scale cavities under roads. However, in reality, owing to the numerous sources of disturbance in the subsurface, the response of subsurface objects is often distorted and nontypical, such that previous studies mainly validated their methods with synthetic data which are cleaner and simpler. In recent years, many studies have investigated the application of deep learning in B-scan based object recognition, for example, Kilic and Eren (2018), Zhang et al. (2020), andHou et al. (2021). Supervised deep-learning techniques require a large amount of tagging data to train robust models. Owing to the lack of validation techniques, it is difficult to collect sufficient GPR ground truths. Hence, there is substantial scope for improvements in the application of supervised deep-learning techniques for GPR detection.
Few studies have explored feature extraction from C-scans. In a C-scan, the target reflectors were clearly visible only when the reflection intensities of the target reflectors were significantly higher than that of the background. Using the image segmentation technique, the abnormal area can be distinguished from the subsurface environment. For instance, at a certain depth, in a C-scan image, cavities under roads are usually irregular in shape and presented as local but discontinuous objects. In contrast, a continuous reflection in a C-scan is most likely generated by an underground utility (Liu et al. 2021). Therefore, the buried reflector can be roughly distinguished from C-scans based on their geometry. However, unlike traditional remote sensing images composed of multiple bands, GPR Cscans present only single-band reflection intensities. In addition, owing to the diffraction of GPR signals at the edges of the reflector, the object boundary is not necessarily sharp in GPR C-scans. Scattering caused by aggregates in subsurface environment leads to heavily distributed noise. Therefore, the segmentation of abnormal areas from the GPR C-scan is comparatively more difficult than that of traditional images.
Numerous methods have been developed for the segmentation of single-band images. These methods include Laplacian and gradient counts of grayscale values focused on the identification of maximum degrees of difference (Gou et al. 2013), Gaussian determination of edges based on image frequency (Permuter et al. 2006), and K-means and Otsu's clustering-based threshold definition (Lee et al. 1990). Originally proposed by Kass et al. (1988), ACM is a computer vision framework for delineating an object outline from a possibly noisy 2D image, making it suitable for analyzing blurry and noise C-scans of GPR surveys. Among these methods, ACM was adopted in this study for its computational efficiency.
Currently, GPR C-scans only present the reflection intensities at a certain depth, the information provided by them is less than sufficient for determining a target object. As discussed in the previous section, cavities are presented as local high reflections in C-scans, but many subsurface objects (e.g. stones, manmade structures) generate strong GPR reflections, as long as the dielectric contrast is strong. Thus, further illustration with GPR waveform (the B-scan) is critical before a judgment on reflector types can be made. Many studies have indicated the necessity and feasibility of the full-cover GPR survey: it can not only narrow down the interpretation scope, but also contribute to object recognition, as the 3D images present the geometry of the buried reflectors in a straightforward manner. Therefore, in this study, we proposed a "coarse-to-fine" approach for subsurface modeling with GPR. First, the anomalies from the fullcover images were identified using the active contour model (ACM) and an image segmentation technique; then the reflections were recognized and classified in abnormal areas using a pattern recognition technique known as aggregate channel feature (ACF). The proposed workflow was evaluated with controlled and onsite case studies, and the results showed that over 70% of the subsurface reflectors were successfully identified and classified automatically. However, the GPR survey efficiency could still be improved through further research.

Workflow of automatic subsurface modeling with GPR
The proposed workflow consists of three stages: 1) producing 2D and 3D GPR representation (B-scan and Cscan); 2) roughly extracting the abnormal area from Cscans; and 3) recognizing and classifying the reflector type from B-scans. A custom-built MATLAB program was developed to realize the proposed workflow.

Production of 2D and 3D GPR representations
Collected GPR raw data are firstly processed to enhance the signal to noise ratio. According to Jol (2009), basic signal processing includes de-wow, time-varying gain, time zero correction, frequency filtering, and migration. Afterward, processed B-scans are geo-referenced based on the survey grid. The average of the absolute amplitude of the reflected signal over a vertical time window is applied as the time slice parameter to be mapped. Last but not least, the time-averaged amplitudes are gridded to create a spatially average and full cover image -Cscan, using information from surrounding points.

Abnormal area extraction from C-scans with CV-ACM
The basic principle of ACM is to define the initial curve (C), obtain the energy function based on the image data, trigger a curve change by minimizing the energy function to make it approach the object edge gradually, and finally, find the object edge. The edge curve obtained by this dynamic approximation method has the advantages of being closed and smooth (Yang 2014). Therefore, ACM works automatically as it autonomously and adaptively searches for the minimum state where minimum human interaction is required. While many objects are irregular and blurred in GPR C-scans, ACM as a deformable model can find illusory contours in the image by ignoring missing boundary information.
There are many types of ACM, such as the edge-basedACM (Mignotte and Meunier 2001) and the region-based one (Zhang et al. 2010), level-set based (Osher and Sethian 1988), and Chan-Vese (Chan and Vese 2001). In this study, Chan-Vese ACM (CV-ACM) was selected because it is based on the image energy distribution in which the evolution curve is driven toward the object edge by the minimum value of the energy function. Hence, the CV-ACM is less dependent on the image gradient, which improves its performance on objects with blurred edges whose gradients are relatively smooth.
An ACM is an energy-minimizing, deformable spline influenced by constraints and image forces that pull it toward object contours and internal forces that resist deformation. As illustrated in the equation below, the CV-ACM operates based on the image energy distribution. In a GPR C-scan, the grayscale usually describes the reflection energy strength. The reflection energy is normalized to a range of 0-255, and assigned to each pixel of a C-scan. Thus the grayscale of a C-scan is 8bit.
where C is the iterative curve in the GPR C-scan, C 0 is the mean gray value inside curve C, and C b denotes the mean gray value outside curve C. Combining the energy calculation function, the detailed equation of CV-ACM is as follows: where I(x,y) denotes the pixels in a GPR C-scan. In general, CV-ACM may be understood as a special case of the general technique of matching a deformable model to an image through energy minimization (Equation (1)). It intuitively integrates the external and internal image forces enabling it to perform well in terms of noise exclusion. The effectiveness and efficiency of CV-ACM in noise exclusion and blurred edge delineation make it specifically suitable for interpreting noisy and blurry GPR C-scans.

Anomaly classification from B-scans with ACF
According to Topp and Davis (1980), and Snell's Law, the significant permittivity contrast between the soil and the material of subsurface structures lead to a powerful reflection intensity (Jol 2009). The GPR reflection patterns of different objects have been investigated in previous studies to derive automatic recognition methods. Kofman and Ronen (Kofman, Ronen, and Frydman 2006) illustrated that the inner boundary of an air-filled cavity causes reverberation or "ringing" of the electromagnetic waves. Lai and Chang (Lai et al. 2017) found that although many kinds of buried infrastructure may cause reverberation patterns in B-scans, for example, manholes and metal plates, the reverberation signals from these objects start from time zero, whereas those from air-fill cavities continue to attenuate in a time/depth window that is not close to the surface. Luo and Lai (2020) proved that the size of the cavity determines the morphology of the GPR reflection. Hence, each type of reflector may yield specific reflection patterns on the radar.
After reviewing the strengths and weaknesses of multiple methods, the ACF was applied in this study because it is computationally efficient in lowdimension problems (Yang et al. 2014 ;Luo and Lai 2020a), making it suitable for recognizing irregular reflections of subsurface object from a large amount of data collected 3D GPR surveys.
The ACF trains the detector from a set of training set data by computing registered maps with gradients and histograms with oriented gradients, and then extracts features on these extended channels (Yang et al. 2014;Kuzmenko et al. 2019;Yi et al. 2019). Therefore, it offers a rich representation capacity, whereas the simple feature form guarantees fast computation. The general equation of ACF is as follows: where Ω denotes the function to calculate the features, F denotes the features, ∑ is the aggregation operation, and F^1 is the aggregated feature.
For the B-scan image F = I, the channel feature is an intensity map, while the intensity (grayscale) is the simplest feature channel. After inputting the detection image, the gradient amplitude and histogram pyramid of the gradient direction were calculated. These three features were then combined to form an ACF. The decision trees were applied to the weak classifier to train the boosting classifier and finally obtain the initial detection window and results. A detailed clarification of the ACF is presented by Yang et al. (2014).
This study focuses on four typical objects which include: utility, manhole, cavities, and loose soil. The ground truths for training the classifier were accumulated from previous GPR surveys. The training set consists of approximately 100 sample radargrams collected in various situations. The image size of radargram samples were normalized to 450 × 320. The four typical objects were described by approximately 25 validated radargrams. Figure 2 present the reflection pattern of typical subsurface objects. Then, a multiagent ACF detector was trained with the above radargrams. The amount of ground truth was not sufficient for deep learning-based recognition because the excavation verification is not always available. However, the lack of ground truth still acts as the bottleneck of automatic NDT inspection.

Validation with case studies
It is difficult to obtain ground truths because the subsurface is invisible, and excavation is not always possible. Without the ground truth, the validity of the proposed method cannot be proved. Thus, the adopted algorithms (CV-ACM and ACF) were first tested with a smaller dataset acquired in the laboratory. Subsequently, the proposed workflow was validated using an underground dataset acquired in Hong Kong, China.

Controlled test
In the small-scale controlled case, the GPR dataset was collected in the laboratory at the Hong Kong Polytechnic University. A surface void (0.2 m diameter) was dug in the sand tank. GPR measurements were collected using a GSSI system and a 900 MHz antenna in the common-offset survey mode. The survey area was 1.2 m × 1.6 m, while the profile spacing was 0.1 m and the interval between traces was 0.004 m. C-scans were generated according to the standardized workflow proposed by Luo et al. (2019). As it is suggested that the slice thickness should be similar to the radar wavelength when the object size is unknown, thus in this case, the slice thickness was 2ns with a velocity of 0.12 m/ns. And Inverse distance weighted (IDW) interpolation was applied to fill in the blank areas among GPR traverses, because IDW is a precise interpolator that only consider limited adjacent measurements. The intensity of each C-scan is relatively normalized to 8bit grayscale, which means the color of a C-scan only depicts the intensity contrast of a certain depth. The survey grid is shown in Figure 3(a).

Abnormal area extraction from C-scan
As shown in Figure 3.b, four slices at different depths were selected to illustrate the performance of the CV-ACM. The previous three slices depict the top, middle, bottom of the void dug, and the last slice displays the bottom of the survey tank. Areas with stronger reflection intensities were identified and circled with red lines. The small scatterings caused by the sand particles were successfully excluded, even though their intensities were similar. There is no particular area bounded in the 4th C-scan because the contrast between the intensity contrast of abnormal areas and that of the surrounding environment was subtle.
The dimensions of the areas extracted by CV-ACM were compared with the actual ground truth. It was found that the sizes of the surface void (A) and the manhole (B) are 10% larger than the actual size, while the utility (C) is 20% larger than the actual diameter, and the diffractions generated by the edges of the reflectors contribute to the size of the abnormal areas. Although migration was applied in B-scan processing, considering the footprint of a GPR trace and the interval between traces, the C-scan does not provide sharp boundaries. In addition, the diffractions and interference of the reflected waves result in blurry boundaries of the anomalies. The controlled case proved that the CV-ACM could define the edges without solid parameters.
However, the C-scan can only depict areas that yield different GPR intensities, which means that the materials and defect types cannot be distinguished from C-scans without supplementary information. For instance, in the controlled case, both the manhole and the surface void present as local high reflections in the C-scan, they were considered as cavities. Further analysis using B-scans is required.

Reflector classification from B-scans
The B-scans that cut across the center of the segments were indexed and analyzed with the detector trained by the ACF. Considering the polarization of the GPR, for each segment, two B-scans perpendicular to each other were indexed. The analysis results of the ACF are displayed in Figure 4.
The ACF requires fine sensitivity tuning, which is controlled by the number of negative samples to be used at each stage. After multiple attempts, it was observed that when the sensitivity was normalized to 20, adjacent scatterings could be excluded. The resolution of the training dataset also affected the ACF performance because it determined the scale of the GPR images. As the resolution of the GPR images was controlled by the samples per scan and the scan per meter, the resolution of different GPR datasets was relatively stable. The positions of anomalies identified from B-scans coincided with those of reflectors extracted from C-scans, validating the feasibility of reflector extractions by ACF. However, the anomalies in Figure 4 were all classified as cavities, whereas anomaly A was in fact a manhole. As shown in Figure 4, the GPR reflection pattern of the concrete manhole was similar to that of the surface void, as they both were air-filled local reflectors. Therefore, the ACF classified the GPR reflection patterns rather than the type of reflectors; it could not distinguish reflectors that yielded similar reflection patterns. Furthermore, no other types of subsurface object were found in the controlled case, demonstrating the ACF's performance in terms of classification precision.
The controlled case validates the ACF performance with a small training set. The controlled case was conducted in a laboratory environment in which the surrounding material was relatively homogeneous. In reality, the scattering generated by the subsurface particles may distort the GPR reflections of defects and further affect the defect classification results.

On-site case
To evaluate the validity of the proposed workflow, an on-site case that can represent complex and inhomogeneous realities is necessary. The GPR data of a vehicle lane at the Hong Kong Polytechnic University were used for the uncontrolled case. In the survey area of the on-site case dataset, two underground utilities (UU) was buried at a depth of 1 m. The trial GPR dataset was collected by a common-offset GPR system (GSSI 3000 control unit, 400 MHz antenna), based on a predefined orthogonal grid (20 m × 5 m, profile spacing of 0.5 m). The survey area was covered by buildings and tress adjacent, thus GNSS slocalisation was unavailable and gridding survey was chosen ( Figure 5). The collected GPR data were pre-processed according to the typical signal processing flow introduced by Jol (2009). Signal processing included de-wow, time-varying gain, time zero correction, and frequency filtering. Pre-processing enhanced the signal-to-noise ratio of the GPR data. Detail parameters and values are display in Table 1.
C-scans for various depths were generated from the processed B-scan according to the imaging criteria developed by Luo et al. (2019). The vertical slice thickness was defined as 0.1 m (wave velocity = 0.095 m/ ns), considering the diameter of major underground utilities. IDW interpolation was applied to construct full-cover C-scans. As shown in Figure 5, the generated C-scan offers a top view of the survey area. A brighter color denotes a higher reflection intensity, while darker areas process stronger signal attenuation. The home-built program was tested in a workstation with i7 3 GHz CPU and 32GB RAM.

Results of abnormal area extraction from the on-site case C-scans
The grayscale of the overall C-scans was normalized globally. Because the entire survey area was relatively large, and the discrepancies of GPR reflections caused by particle aggregates in the subgrade within the survey area were significant, the C-scans could only Figure 5. Illustration of the on-site survey site and the generated C-scan at 1 m depth.
display the GPR reflections that were strong enough. Therefore, the tiny anomalies may be hidden by the GPR responses of the buried infrastructure in the subsurface.
Adapted CV-ACM was applied to the C-scan set. The areas with relatively high intensity were randomly distributed among the survey area (bounded by red lines). There are multiple reflectors, including the pipe, manholes and unknown clusters, were extracted. As shown in Figure 6(c), a suspected utility is buried a meter deep. The small scattering caused by aggregates in road structures was excluded; they are determined as background noise. However, ACM segmentation could not distinguish between different reflectors. Further constraints based on reflector geometry and reflection intensities are needed to classify the defects from the subsurface infrastructure.
As shown in Figure 6, the proposed CV-ACM successfully excluded the small scattering and only focused on the larger area that yielded more obvious reflection differences. However, owing to the special mechanism of the GPR survey, the reflector edges were not sharp. The CV-ACM could still define the boundaries of anomalies by striving to balance the inner and outside image energy, proving its feasibility in GPR C-scan interpretation.
The geometry of the segmented abnormal areas was estimated using the ratio of its length and width. The long linear shape segment (E) was considered as the UU, while the remaining irregular-shaped local anomalies (A, C, F) were suspected to be defects under surface. A short linear segment (B, D) was underdefined. Therefore, C-scans can provide an initial classification based on the geometry of the reflectors. Considering the orientation of the linear segment, only the vertical GPR traverse was indexed and labeled A-F (Figure 7).

Results of the reflector classification from the on-site case B-scans
The B-scans that cut across the abnormal areas were indexed for reflector classification. The sensitivity was set to be the same as that used in the control test. The E radargram cut across the centroid of the E segment (suspected utility), and the hyperbolic-shape reflection was extracted and classified, further confirming the initial   classification result from the C-scan. In terms of round anomalies (suspected defects), manhole-related patterns were found in radargram A and F. Radargram B and D crossed the centroid of the underdefined B and D segments. Multiple hyperbolic patterns and loose soil related patterns were observed in these two radargrams. By co-registering the position of the C-and B-scans, it was found that the yellow segment was formed by hyperbolic reflection, which was likely a utility ( Figure 8).
Unfortunately, the survey area was a vehicle lane in use; thus, we could not validate the recognition result with excavation testing. The sensitivity and recall analyses of the anomaly classification results were conducted with visual inspection based on our experiments; in total, 34 anomalies were extracted by the ACF, and categorized as true positive (TP), false positive (FP), true negative (TN), and false negative (FN), according to their waveform information.  Figure 7. Remarks: yellow, blue, red, and green denote cavity, utility, manhole, and loose soil, respectively.
The sensitivity and recall analysis results of the test case were as follows. The accuracy of the abnormal identification is estimated as 76.47%, and the precision of the anomalies-classification is calculated as 70.59%. In conclusion, the ACF reclassified the extracted segments: E, B and D linear segments were utility sections, and A, C and F local segments were manholes. In this manner, the subsurface was depicted automatically.
Finally, the modeled result of the proposed method is shown as Figure 9. The segments extracted from Cscans were re-labeled with the classification results from B-scans. A semantic map was generated. In this way, the model can be communicated among GIS or BIM systems. It took the program less than 2 minutes to finish the whole subsurface modeling process, from C-scan extraction to B-scan classification.

Discussion
Through the case studies, the strengths and weaknesses of the proposed workflow were observed, and suggestions for improving vulnerabilities were proposed. In general, the integration of CV-ACM and ACF identified and interpreted GPR data with a limited training dataset, and the performance was promising.

Necessity of full-cover 3D-survey
Both the controlled test and site case illustrated the necessity of a 3D GPR survey. The 3D GPR images (Cscans) provide a straightforward view of the survey area; the scale and shape of the reflectors can be roughly depicted. Thus, the abnormal area can be identified from the C-can by "a glance", and the interpretations on the B-scans line by line is not necessary anymore: the efficiency of 3D GPR survey is significantly improved. The horizontal and vertical resolution of the antenna depends on the wavelength, wavenumber, beamwidth, etc. Furthermore, the traditional 2D GPR can only image the vertical section along the survey traverse; undetected defects may exist when the defect size is smaller than the gap between adjacent GPR traverses. However, the denser the traverse, the greater the time and labor consumption. Efficiency improvements in 3D GPR surveys have now made large-scale subsurface modeling a possibility. When the 3D GPR provides full-cover measurements in the geographic coordinate system, imaging results can be communicated among management organs and support timely actions. Therefore, the 3D GPR survey contributes to the safe operation of urban infrastructure.
However, as C-scans are currently generated from the GPR reflection intensities, the information provided by C-scans is not sufficient to distinguish different reflectors, and further analysis from the GPR waveform is essential. For instance, utilities present as hyperbola in the B-scans, and in on-site case, the ACF extracted the hyperbolic-like reflections, confirming the utility of the linear segment found in the C-scans. In contrast, the local anomalies that were initially defined as suspected defects were reclassified as manholes in B-scans with ACF, illustrating the necessity of B-scan analysis. Therefore, Bscans that present waveform attribute information contribute to the confidence of the classification. In addition, it is suggested that further research can enrich the information of 3D GPR representations. The workload of a B-scan analysis can be reduced when more GPR waveform information can be learned from C-scans.

Performance of the proposed workflow
This study proposed a 3D GPR data interpretation strategy. In the site case, the proposed workflow achieved promising results and proved the suitability of utilizing CV-ACM and ACF in GPR data interpretation. Since these two algorithms were originally developed for analyzing optical images and the wavelength of GPR is much larger than that of visible light, tailored refinements and adaption were applied to these two methods to enhance their capability to analyze blurred GPR images.
Specifically, the constraint of CV-ACM for eliminating scattering in C-scans was determined by empirical experiments in which subjective biases were inevitable. In future studies, iterative convergence could be introduced into the ACM constraints to find the balance between excluding the scatterings while maintaining the integral part of suspected anomalies. In terms of reflector classification from B-scans, even though the patterns of cavities and manholes are relatively similar, the ACF successfully distinguished them through the range of reverberation. Although there were some errors, the entire interpretation was conducted in a fully automatic manner. Hence, with the proposed workflow, the efficiency of the GPR survey can be improved significantly.
The performance of ACF relies on the fine-tuning of parameters such as sliding window size, training dataset size, and the number of bootstrapping stages. In this study, the suitable parameters used in ACF detection were also empirical: they can be automated by deep learning for further improvements. Moreover, the training dataset size was insufficient, limiting the sensitivity of the recognition and resulting in false alarms and missing data. However, this situation is unlikely to be solved in the near future, as it is challenging to validate subsurface inspection. Therefore, more attention should be paid to the characterization of the GPR responses of different objects to improve the accuracy of GPR pattern recognition. The continued accumulation of the training dataset will contribute to automatic GPR data interpretation using artificial intelligence.

Conclusion
This study proposed and validated an automatic workflow for identifying and classifying subsurface objects in GPR data. The workflow is a "coarse to fine" approach that make use of the information of both C-scans and B-scans. First, full cover GPR measurements are collected to form 2D and 3D GPR representations. Afterward, the abnormal areas with relative high reflection intensity are extracted to segments by the CV-ACM. Finally, the B-scan that across the abnormal areas are indexed to classify the segments with ACF based on a pre-trained model. A laboratory and an on-site case were carried out to validate the proposed method, and it was proved that a semantic map was automatically constructed from GPR data based on the proposed method. The study proves the necessity of a full-cover 3D GPR survey, illustrates the feasibility of integrating advanced computer vision techniques to analyze a large amount of 3D GPR data, and paves the way for automating the process of subsurface modeling with GPR.