Frequency–Wavenumber Analysis of Deep Learning-based Super Resolution 3D GPR Images

: This paper proposes a frequency–wavenumber ( f–k ) analysis technique through deep learning-based super resolution (SR) ground penetrating radar (GPR) image enhancement. GPR is one of the most popular underground investigation tools owing to its nondestructive and high-speed survey capabilities. However, arbitrary underground medium inhomogeneity and undesired measurement noises often disturb GPR data interpretation. Although the f–k analysis can be a promising technique for GPR data interpretation, the lack of GPR image resolution caused by the fast or coarse spatial scanning mechanism in reality often leads to analysis distortion. To address the technical issue, we propose the f–k analysis technique by a deep learning network in this study. The proposed f–k analysis technique incorporated with the SR GPR images generated by a deep learning network makes it possible to signiﬁcantly reduce the arbitrary underground medium inhomogeneity and undesired measurement noises. Moreover, the GPR-induced electromagnetic waveﬁelds can be decomposed for directivity analysis of wave propagation that is reﬂected from a certain underground object. The e ﬀ ectiveness of the proposed technique is numerically validated through 3D GPR simulation and experimentally demonstrated using in-situ 3D GPR data collected from urban roads in Seoul, Korea.


Introduction
In the past few decades, sinkhole accidents of urban roads have posed a serious hazard to buildings, infrastructures and especially inhabitants of the area [1,2]. Although vision-based road surface inspection techniques have been widely proposed for road degradation evaluation [3,4], early detection of sinkholes, which are typically invisible from the road surface, is still challenging. To effectively detect underground cavities which are most likely extended to sinkholes, various nondestructive testing (NDT) techniques have taken the limelight. Ground penetrating radar (GPR) is one of the widely accepted NDT tools thanks to its high sensitivity to underground media change and rapid inspection capability for broad target areas [5][6][7]. However, the physical interpretation of field GPR data for underground object detection and classification is still challenging in some cases, because electromagnetic waves, which are reflected from a target underground object, are often weaker than underground media's inhomogeneity and undesired measurement noises [8,9]. In general, most of the dominant signals reflected from the road surface often hinder the precise data interpretation of relatively weak signals coming from underground media under air-coupled GPR data acquisition conditions [8,10,11].
To enhance the GPR data interpretability, a number of signal and image processing techniques, such as time-varying gain [10,11], subtraction [8], migration [12], deconvolution [13], basis-pursuit [9], compressive sensing [14], velocity analysis [15], radon transform [16], discrete wavelet transform [17] one of the main distinctive features for underground object identification and classification. However, lack of image resolution often hinders the feature recognition. To obtain high resolution B-and C-scan images, slow scanning speed and dense GPR antenna arrangement are necessary. Unfortunately, the resolution issues are, however, a trade-off with respect to time and cost in reality.
To effectively tackle the lack of resolution issue in existing 3D GPR data without data acquisition condition change, the deep learning-based SR GPR image enhancement network based on a deep residual channel attention network [39] is developed as shown in Figure 1. The deep residual channel attention network is one of the SR image enhancement networks based on CNN, which consists of 500 layers and 1.6 M parameters. This network increases the LR image resolution four times, comprises the four main steps: (1) shallow feature extraction, (2) deep feature extraction, (3) upscaling and (4) reconstruction. First, the shallow feature extraction step, which consists of a single convolution layer with 64 kernels of 3 × 3 size and stride of 1, extracts shallow features from the input low resolution (LR) image.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 18 widely used as one of the main distinctive features for underground object identification and classification. However, lack of image resolution often hinders the feature recognition. To obtain high resolution B-and C-scan images, slow scanning speed and dense GPR antenna arrangement are necessary. Unfortunately, the resolution issues are, however, a trade-off with respect to time and cost in reality.
To effectively tackle the lack of resolution issue in existing 3D GPR data without data acquisition condition change, the deep learning-based SR GPR image enhancement network based on a deep residual channel attention network [39] is developed as shown in Figure 1. The deep residual channel attention network is one of the SR image enhancement networks based on CNN, which consists of 500 layers and 1.6 M parameters. This network increases the LR image resolution four times, comprises the four main steps: (1) shallow feature extraction, (2) deep feature extraction, (3) upscaling and (4) reconstruction. First, the shallow feature extraction step, which consists of a single convolution layer with 64 kernels of 3 × 3 size and stride of 1, extracts shallow features from the input low resolution (LR) image. Subsequently, deep features are extracted through the residual-in-residual structure in the second step of Figure 1. The main purpose of the residual-in-residual structure learning is to allow very deep networks to achieve easy and powerful performance in the training process. Here, the deep feature means high frequency information composed of lines or edges, which is the biggest difference between the SR and LR images. The residual-in-residual structure is the very deep network which consists of 10 residual groups, and each residual group includes 20 residual blocks and 1 convolution layer as shown in Figure 1. To effectively train the deep network, the skip connections are embedded in the residual-in-residual structure in terms of the short and long skip connections. Each residual group is connected by the long skip connection, and the residual blocks inside the residual group are connected by the short skip connection as shown in Figure 1. The multiple skip connections allow enough shallow features to be bypassed, enabling the main network to focus on learning deep features. The skip connections make it possible to stabilize the training performance of the very deep network with residual learning, as the main superiority of the residual-in-residual structure. Each convolution layer inside the residual-in-residual structure similarly works as the convolution layer of the first step. To more efficiently train the high frequency regions of the image, the channel attention mechanism is performed inside each residual block through global average pooling, convolution, rectified linear unit (ReLU) layer and sigmoid function. The channel attention mechanism improves the discriminative learning capabilities by focusing on more useful channels based on the average value extracted from each channel. Each average value is obtained via the global Subsequently, deep features are extracted through the residual-in-residual structure in the second step of Figure 1. The main purpose of the residual-in-residual structure learning is to allow very deep networks to achieve easy and powerful performance in the training process. Here, the deep feature means high frequency information composed of lines or edges, which is the biggest difference between the SR and LR images. The residual-in-residual structure is the very deep network which consists of 10 residual groups, and each residual group includes 20 residual blocks and 1 convolution layer as shown in Figure 1. To effectively train the deep network, the skip connections are embedded in the residual-in-residual structure in terms of the short and long skip connections. Each residual group is connected by the long skip connection, and the residual blocks inside the residual group are connected by the short skip connection as shown in Figure 1. The multiple skip connections allow enough shallow features to be bypassed, enabling the main network to focus on learning deep features. The skip connections make it possible to stabilize the training performance of the very deep network with residual learning, as the main superiority of the residual-in-residual structure. Each convolution layer inside the residual-in-residual structure similarly works as the convolution layer of the first step. To more efficiently train the high frequency regions of the image, the channel attention mechanism is performed inside each residual block through global average pooling, convolution, rectified linear unit (ReLU) layer and sigmoid function. The channel attention mechanism improves the discriminative learning capabilities by focusing on more useful channels based on the average Remote Sens. 2020, 12, 3056 4 of 18 value extracted from each channel. Each average value is obtained via the global average pooling layer. Then, the convolution, ReLU layer and sigmoid function provide non-linearity between the channels and make multiple channel-wise features to be emphasized for a non-mutually exclusive relationship. Next, the upscaling step, as the third step, extends the feature data of the image to the SR resolution size. It consists of the deconvolution layer with 256 kernels of 3 × 3 size and the stride of 1, which increases the size of each pixel by four times in this network. Finally, the SR image is generated through a single convolution layer comprised of three kernels of 3 × 3 size and stride of 1 in the reconstruction step, as the fourth step, of Figure 1.

SR GPR Image-Based f-k Analysis
Once the SR GPR images are generated, the SR C-scan images can be obtained in the t-s domain for the subsequent f-k analysis, as shown in Figure 2. Note that the proposed f-k analysis can be easily extended to B-or D-scan images, although the sequential C-scan images along the depth direction are used in this study. Even though the SR C-scan images contain meaningful wave signals mixed with noise components, it is difficult to remove the undesired noises caused by arbitrary underground medium inhomogeneity and data measurement procedure in the t-s domain. On the other hand, the noise components can be effectively eliminated in the f-k domain. Moreover, wavefield decomposition in the f-k domain is able to precisely analyze wave propagation directivity by reconverting the decomposed wavefield data to the t-s domain data as shown in Figure 2. First, the SR C-scan data in the t-s domain are transformed into the f-k domain through 3D Fourier transform which is given by: where E and U denote the electromagnetic wavefields of the SR C-scan data in the t-s and f-k domains, respectively. k, w and t are the wavenumber, angular frequency and time, respectively. x and y are the spatial cartesian coordinates.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 18 average pooling layer. Then, the convolution, ReLU layer and sigmoid function provide non-linearity between the channels and make multiple channel-wise features to be emphasized for a non-mutually exclusive relationship. Next, the upscaling step, as the third step, extends the feature data of the image to the SR resolution size. It consists of the deconvolution layer with 256 kernels of 3 × 3 size and the stride of 1, which increases the size of each pixel by four times in this network. Finally, the SR image is generated through a single convolution layer comprised of three kernels of 3 × 3 size and stride of 1 in the reconstruction step, as the fourth step, of Figure 1.

SR GPR Image-Based f-k Analysis
Once the SR GPR images are generated, the SR C-scan images can be obtained in the t-s domain for the subsequent f-k analysis, as shown in Figure 2. Note that the proposed f-k analysis can be easily extended to B-or D-scan images, although the sequential C-scan images along the depth direction are used in this study. Even though the SR C-scan images contain meaningful wave signals mixed with noise components, it is difficult to remove the undesired noises caused by arbitrary underground medium inhomogeneity and data measurement procedure in the t-s domain. On the other hand, the noise components can be effectively eliminated in the f-k domain. Moreover, wavefield decomposition in the f-k domain is able to precisely analyze wave propagation directivity by reconverting the decomposed wavefield data to the t-s domain data as shown in Figure 2. First, the SR C-scan data in the t-s domain are transformed into the f-k domain through 3D Fourier transform which is given by: where and denote the electromagnetic wavefields of the SR C-scan data in the t-s and f-k domains, respectively. , and are the wavenumber, angular frequency and time, respectively. and are the spatial cartesian coordinates. Subsequently, the tailored f-k filter is designed. First, a lowpass filter is applied in the f domain so that the measurement noises can be eliminated outside the excitation frequency range as shown in Figure 3. The filtering frequency bandwidth can be determined by considering the excitation frequency range. The k domain filter is then developed. Since the electromagnetic waves propagating along arbitrary underground media are partially and randomly reflected from the media's inhomogeneity and numerous small porosities, these reflection signals randomly appear and look like non-propagating wave components in the C-scan images of the t-s domain. It physically means that no spatially propagating wave can be observed if there is no certain object inside the underground media. Even if there are meaningful signals reflected from a certain object, such physical phenomenon can be clearly observed, as shown in Figure 3a. The wave energy is highly Subsequently, the tailored f-k filter is designed. First, a lowpass filter is applied in the f domain so that the measurement noises can be eliminated outside the excitation frequency range as shown in Figure 3. The filtering frequency bandwidth can be determined by considering the excitation frequency range. The k domain filter is then developed. Since the electromagnetic waves propagating along arbitrary underground media are partially and randomly reflected from the media's inhomogeneity and numerous small porosities, these reflection signals randomly appear and look like non-propagating wave components in the C-scan images of the t-s domain. It physically means that no spatially propagating wave can be observed if there is no certain object inside the underground media. Even if there are meaningful signals reflected from a certain object, such physical phenomenon can be clearly observed, as shown in Figure 3a. The wave energy is highly concentrated near zero k x value, which are undesired noise components to be eliminated. One more interesting thing to see here is that these randomly reflected signals in the f-k domain may have dominant energy due to their high repetition rate, although each reflection signal intrinsically has small amplitude. Based on the physical observation, the k domain filter is established using a Laplacian of Gaussian window (Φ k ): where σ is the standard deviation.
Once the f-k filter is designed, the filtered SR C-scan data (U f ) can be obtained in the f-k domain. Figure 3b shows that the measurement noises, as well as non-propagating components, are clearly filtered out, and meaningful wave components remain.
Furthermore, the electromagnetic wavefields can be decomposed in the f-k domain so that the wave propagation directivity can be precisely analyzed, which is useful to recognize underground objects' size and location as well as to classify the object type. For instance, a ±x directional window filter (Φ ±k x ) can be designed for decomposing U f to the +x or −x directional wavefield (U ±k x ) in the f-k domain, which is given by: In a similar fashion, it can be readily extended to the ±y directional filter.
concentrated near zero value, which are undesired noise components to be eliminated. One more interesting thing to see here is that these randomly reflected signals in the f-k domain may have dominant energy due to their high repetition rate, although each reflection signal intrinsically has small amplitude. Based on the physical observation, the k domain filter is established using a Laplacian of Gaussian window (Φ ): where is the standard deviation.
Once the f-k filter is designed, the filtered SR C-scan data ( ) can be obtained in the f-k domain. Figure 3b shows that the measurement noises, as well as non-propagating components, are clearly filtered out, and meaningful wave components remain.
Furthermore, the electromagnetic wavefields can be decomposed in the f-k domain so that the wave propagation directivity can be precisely analyzed, which is useful to recognize underground objects' size and location as well as to classify the object type. For instance, a ± directional window filter (Φ ± ) can be designed for decomposing to the + or − directional wavefield ( ± ) in the f-k domain, which is given by: In a similar fashion, it can be readily extended to the ± directional filter. Next, the resultant C-scan data ( ± ) in the t-s domain can be reconstructed using the following inverse 3D Fourier transform: As one of the representative examples, only the − directional wavefield ( ) remains, and it reveals much higher signal-to-noise ratio (SNR) than without pixel information loss and distortion, as shown in Figure 2. Note that the wavefield decomposition process is optional in the algorithm, thus ± can be replaced by in Equation (5), resulting in the filtered SR C-scan data ( ) in the t-s domain.
(a) (b) Figure 3. The representative f-plots at the center of (a) before and (b) after filtering. Figure 3. The representative f -k x plots at the center of k y (a) before and (b) after filtering.
Next, the resultant C-scan data (E ±k x ) in the t-s domain can be reconstructed using the following inverse 3D Fourier transform: U ±k x k x , k y , w e i(k x x+k y y+wt) dk x dk y dw Remote Sens. 2020, 12, 3056 6 of 18 As one of the representative examples, only the −x directional wavefield (E −k x ) remains, and it reveals much higher signal-to-noise ratio (SNR) than E without pixel information loss and distortion, as shown in Figure 2. Note that the wavefield decomposition process is optional in the algorithm, thus U ±k x can be replaced by U f in Equation (5), resulting in the filtered SR C-scan data (E f ) in the t-s domain.

Numerical and Experimental Validations
The proposed f-k analysis technique is numerically and experimentally validated through 3D GPR simulation using gprMax [43] and in-situ 3D GPR data obtained from urban roads in Seoul, Korea.

Numerical Validation
The target 3D model is comprised of 8 × 2.975 × 2.75 m 3 soil layer, 8 × 0.525 × 2.75 m 3 air layer and steel pipe with a diameter of 500 mm, as depicted in Figure 4. It was modelled so that the pipe was buried perpendicular to the GPR scanning direction inside the soil layer. Note that the pipe was intentionally selected in this study, because it is one of the representative wave scatters which can be clearly reflected in all GPR channels constituting the C-scan images. Here, the relative permittivity values of air, soil and pipe were set to 1, 5 and infinity, respectively. The transmitter (Tx) was 50 mm apart from the receiver (Rx), and the GPR data reflected from the pipe were acquired by moving the Tx and Rx antennas along the soil layer surface, as shown in Figure 4. A finite difference time domain method [44] was used to simulate electromagnetic wave propagation. To simulate the similar conditions with the real-world GPR scanning, the spatial discretization was set to 20 mm, which is equivalent to 20 km/h scanning speed with 20 GPR channels in the real-world application. Here, the 20 GPR channels are able to cover a road width of 1.5 m. The excitation electromagnetic wave was normalized by the second derivative of a Gaussian waveform with a center frequency of 1.8 GHz. In addition, Gaussian random noises, which are equivalent to 25% of magnitude of the maximum value of the GPR signal, were artificially added to simulate the arbitrary underground medium inhomogeneity and undesired measurement noises.

Numerical and Experimental Validations
The proposed f-k analysis technique is numerically and experimentally validated through 3D GPR simulation using gprMax [43] and in-situ 3D GPR data obtained from urban roads in Seoul, Korea.

Numerical Validation
The target 3D model is comprised of 8 × 2.975 × 2.75 m 3 soil layer, 8 × 0.525 × 2.75 m 3 air layer and steel pipe with a diameter of 500 mm, as depicted in Figure 4. It was modelled so that the pipe was buried perpendicular to the GPR scanning direction inside the soil layer. Note that the pipe was intentionally selected in this study, because it is one of the representative wave scatters which can be clearly reflected in all GPR channels constituting the C-scan images. Here, the relative permittivity values of air, soil and pipe were set to 1, 5 and infinity, respectively. The transmitter (Tx) was 50 mm apart from the receiver (Rx), and the GPR data reflected from the pipe were acquired by moving the Tx and Rx antennas along the soil layer surface, as shown in Figure 4. A finite difference time domain method [44] was used to simulate electromagnetic wave propagation. To simulate the similar conditions with the real-world GPR scanning, the spatial discretization was set to 20 mm, which is equivalent to 20 km/h scanning speed with 20 GPR channels in the real-world application. Here, the 20 GPR channels are able to cover a road width of 1.5 m. The excitation electromagnetic wave was normalized by the second derivative of a Gaussian waveform with a center frequency of 1.8 GHz. In addition, Gaussian random noises, which are equivalent to 25% of magnitude of the maximum value of the GPR signal, were artificially added to simulate the arbitrary underground medium inhomogeneity and undesired measurement noises.     Figure 6a shows the representativeplots at 300 MHz. As expected, the non-propagating components caused by incoherent noise components are highly concentrated on the origin of the plane. To remove the non-propagating components, the f-k filter is applied to shown in Figure 6a using Equation (2). The lowpass filter was designed by fitting an exponential function with a rate parameter of 0.05. As for the k domain filter, was set to 1 considering k to all excitation frequency ranges. After applying the f-k filter, the non-propagating components are remarkably reduced in , while meaningful wave components reflected from the pipe remain, as shown in Figure 6b. Subsequently, and are obtained by applying Φ ± using Equation (4) as shown in Figure 6c,d, respectively.  Figure 6a shows the representative k x -k y plots at 300 MHz. As expected, the non-propagating components caused by incoherent noise components are highly concentrated on the origin of the k x -k y plane. To remove the non-propagating components, the f-k filter is applied to U shown in Figure 6a using Equation (2). The lowpass filter was designed by fitting an exponential function with a rate parameter of 0.05. As for the k domain filter, σ was set to 1 considering k to all excitation frequency ranges. After applying the f-k filter, the non-propagating components are remarkably reduced in U f , while meaningful wave components reflected from the pipe remain, as shown in Figure 6b. Subsequently, U −k x and U +k x are obtained by applying Φ ±k x using Equation (4) as shown in Figure 6c,d, respectively.  Figure 6a shows the representativeplots at 300 MHz. As expected, the non-propagating components caused by incoherent noise components are highly concentrated on the origin of the plane. To remove the non-propagating components, the f-k filter is applied to shown in Figure 6a using Equation (2). The lowpass filter was designed by fitting an exponential function with a rate parameter of 0.05. As for the k domain filter, was set to 1 considering k to all excitation frequency ranges. After applying the f-k filter, the non-propagating components are remarkably reduced in , while meaningful wave components reflected from the pipe remain, as shown in Figure 6b. Subsequently, and are obtained by applying Φ ± using Equation (4) as shown in Figure 6c,d, respectively.  Figure 7 shows the resultant t-s domain images corresponding to Figure 6, which are reconstructed using Equation (5). Compared to Figure 7a, it is clearly observed that the incoherent and random noises are significantly eliminated in Figure 7b. To quantitatively estimate the results, the SNR values of the representative A-scan signals along the vertical white dash-dotted lines in Figure 7a,b were compared. Figure 8 shows the A-scan signals with the reference signals obtained by smoothing spline curve fitting. Figure 8a reveals that the A-scan signal of is quite different from the reference signal, resulting in SNR of 19.2 dB as summarized in Table 1. Once the f-k filter is applied, Figure 8b shows that the A-scan signal of is well matched with the reference signal, which has 54.1 dB SNR as shown in Table 1. It can be confirmed that the proposed f-k filter is very effective in removing incoherent and random noise components. In addition, Figure 7c,d, respectively, show that and are successfully decomposed along the − and + directions. Again, the wavefield decomposition is very powerful in identifying the underground object boundary and classifying the object type.  Figure 7 shows the resultant t-s domain images corresponding to Figure 6, which are reconstructed using Equation (5). Compared to Figure 7a, it is clearly observed that the incoherent and random noises are significantly eliminated in Figure 7b. To quantitatively estimate the results, the SNR values of the representative A-scan signals along the vertical white dash-dotted lines in Figure 7a,b were compared. Figure 8 shows the A-scan signals with the reference signals obtained by smoothing spline curve fitting. Figure 8a reveals that the A-scan signal of E is quite different from the reference signal, resulting in SNR of 19.2 dB as summarized in Table 1. Once the f-k filter is applied, Figure 8b shows that the A-scan signal of E f is well matched with the reference signal, which has 54.1 dB SNR as shown in Table 1. It can be confirmed that the proposed f-k filter is very effective in removing incoherent and random noise components. In addition, Figure 7c,d, respectively, show that E −k x and E +k x are successfully decomposed along the −x and +x directions. Again, the wavefield decomposition is very powerful in identifying the underground object boundary and classifying the object type.

Experimental Validation Using In-Situ 3D GPR Data
The proposed f-k analysis technique was also experimentally validated using 3D GPR data collected from complex urban roads in Seoul, Korea. Figure 9a shows the 3D GPR-mounted van for the field tests. The 3D GPR consisted of bow-tie monopole 20 Tx and Rx antennas generate a step frequency with wide frequency bandwidth ranging from 100 MHz to 3 GHz. Here, the 20-channeled GPR device has 1.5 m scanning width, which can typically cover a single road lane as shown in Figure  9a,b. The average scanning speed was approximately 20 km/h for avoiding traffic congestion in urban roads. The GEOSCOPE MK IV data acquisition system shown in Figure 9c, which has 3 GHz sampling rate and 250 ns time range, was used in the field tests.

Experimental Validation Using In-Situ 3D GPR Data
The proposed f-k analysis technique was also experimentally validated using 3D GPR data collected from complex urban roads in Seoul, Korea. Figure 9a shows the 3D GPR-mounted van for the field tests. The 3D GPR consisted of bow-tie monopole 20 Tx and Rx antennas generate a step frequency with wide frequency bandwidth ranging from 100 MHz to 3 GHz. Here, the 20-channeled GPR device has 1.5 m scanning width, which can typically cover a single road lane as shown in Figure 9a   , the Φ ± window filter is similarly applied using Equation (4), as shown in Figures  12c,d and 13c,d.  Figures 12 and 13 show the representative k x -k y plots at 500 MHz of the pipe cases 1 and 2. Similarly, the non-propagation components caused by the incoherent noises are concentrated on the origin of the k x -k y plane in U, as shown in Figures 12a and 13a. Then, U f 's of Figures 12b and 13b show that the undesired non-propagating components are remarkably reduced by applying the filtering parameters, i.e., rate parameter of 0.05 and σ of 1 in Equation (2). To decompose U f into U −k x and U +k x , the Φ ±k x window filter is similarly applied using Equation (4) , the Φ ± window filter is similarly applied using Equation (4), as shown in Figures  12c,d and 13c,d.    Figures 14a and 15a. E −k x and E +k x are also successfully decomposed using Equation (5), as shown in Figure 14c,d and Figure 15c, and are also successfully decomposed using Equation (5), as shown in Figures 14c,d and 15c,d, respectively.   Similar to the simulation one, the quantitative comparison results using SNR are summarized in Table 2. Both the pipe cases 1 and 2 show about 75% improvement after applying the proposed technique, which is consistent with the simulation results.

Discussion
The proposed f-k analysis technique based on SR GPR images was well validated via the numerical simulation and field tests. Note that the SNR improvement rate of simulation is higher than the experimental ones, because the incoherent noises were simply assumed using ideal Gaussian random noises in the simulation. One more interesting thing is that the pipe cases 1 and 2 show the similar SNR improvement rates, which means that the proposed f-k analysis technique is robust against the test environmental variation. In other words, the performance of the proposed technique would be consistent regardless of underground site conditions. Although the pipe cases, which is one of the most dominant features for clearer validation, are shown in the paper, the proposed technique can be easily extended to other types of underground objects. In addition, the cylindrical coordinate can be employed for the directivity analysis depending on the target objects' shapes, rather than the Cartesian coordinate.

Conclusions
This paper proposes a frequency-wavenumber (f-k) technique of 3D ground penetrating radar (GPR) data, which enables one to effectively eliminate incoherent noises and precisely analyze the electromagnetic wave propagation directivity. This technique is newly proposed and validated using super resolution (SR) artificially generated by the deep learning network. Three-dimensional GPR data collected using the existing GPR devices typically suffer from the lack of resolution problem, making it difficult to be analyzed in the f-k domain. To avoid the f-k analysis distortion, a deep learning-based SR GPR image enhancement network incorporated with the f-k analysis is developed. The proposed technique is able to effectively eliminate incoherent noises caused by arbitrary underground medium inhomogeneity and undesired measurement noises, which is one of the biggest technical conundrums in real-world GPR data interpretation. In addition, electromagnetic wave propagation directivity can be precisely analyzed through wavefield decomposition, which is another strong benefit of the f-k analysis. The proposed f-k analysis technique is successfully validated via 3D GPR simulation, as well as field tests, revealing the pretty consistent and outstanding performances. The proposed f-k analysis would be a promising tool for 3D GPR data interpretation especially obtained from complex urban roads.
As the follow-up work, it is warranted that wavefield decomposition-based underground object characterization is thoroughly studied using more GPR data obtained from various in-site roads. Moreover, the proposed f-k analysis can be combined with deep learning-based automated data classification, making it possible to outperform the existing deep learning networks. It is envisioned that this novel f-k analysis can be helpful for not only underground object identification but also concrete structure inspection using GPR.