2D Wavelet Decomposition and F-K Migration for Identifying Fractured Rock Areas Using Ground Penetrating Radar

Abstract

The quality of the surrounding rock is crucial to the stability of underground caverns, thereby requiring an effective monitoring technology. Ground-penetrating radar (GPR) can reconstruct the subterranean profile by electromagnetic waves, but two significant issues, called clutter and hyperbola tails, affect the signal quality. We propose an approach to identify fractured rocks using 2D Wavelet transform (WT) and F-K migration. F-K migration can handle the hyperbola using Fourier analysis. WT can mitigate clutter, distinguish signal discontinuity, and provide signals with a good time-frequency resolution for F-K migration. In the simulation, the migration result from horizontal detail coefficients highlight the crack locations and reduce the scattering signals. Noise has been separated by 2D WT. Hyperbola tails are decomposed to vertical and diagonal detail coefficients. Similar promising results have been achieved in the field measurement. Therefore, the proposed approach can process GPR signals for identifying fractured rock areas.

1. Introduction

Underground water-sealed caverns are of strategic significance for energy resource storage (e.g., gas [1], oil [2]) owing to their advantages, including considerable capacities, reliable structural safety and limited land resource occupation. Surrounding rock quality is a primary factor determining the structural stability of these storage caverns [3], thereby requiring monitoring strategies and rock reinforcement approaches. Although rock quality evaluation has been extensively developed [4,5] and mature numerical models for analyzing fractured rock deformation has been proposed [6,7,8], field detection of rock formations is still limited. The classic opening and pit sampling [5] measure is wasteful and time-consuming, accompanied by the potential risks of breaking intact rocks or ignoring some fractured regions in this discontinuous procedure. Force and acceleration transducer arrays have been applied for monitoring deformation, stress and dynamic responses (e.g., [9,10]) but limited in shallow ground. Seismic methods use high-energy elastic waves to monitor underground caverns (e.g., [11,12]), but the signals are produced by explosions or rock mass failure. Electrical resistivity tomography is suitable for detecting underground terrain (e.g., [13]) but has resolution limitations. Therefore, it is significant to develop effective non-destructive technology for monitoring the surrounding rock quality of underground caverns.

Ground-penetrating radar (GPR) is a non-destructive detection technology to explore the unseen subsurface world, developed based on the theory of electromagnetic wave propagation in materials. It behaves like an echo listener recording the wave reflections from the underground, and this approach has been proven effective in many geological surveys (e.g., [14,15]). For the fractured surrounding rock, the interfaces of numerous fissures become significant reflection and diffraction sources. Nonetheless, the resolution scale is yet the primary issue obstructing the radargram interpretation. The sizes of the joint fissures are small, and densely distributed cracks can generate complicated signals. Although high-frequency electromagnetic waves have a high spatial resolution, the wave energy attenuates quickly within a short propagation distance [16]. Some researches can quantify small voids, objects and throughout cracks in the centimeter or millimeter level (e.g., [17,18,19]), but the complicated situation with joint fissures deeply beneath the ground remains unrevealed in the literature. Another issue hiding crack locations arises from the hyperbola signal of any scatter in the B-scan radargram [20]. The hyperbola tails of crack groups cover the entire B-scan image, which makes the actual locations indistinguishable. Besides, instrumental and environmental clutter increases the difficulty in signal processing [21]. Clutter suppression and hyperbola elimination for fractured rock identification are subjected to researches.

Approaches for eliminating the hyperbola effect has evolved during the last decades. The direct method is to calculate the hyperbola and then either sum the signals to the focusing point [22] or mitigate the hyperbola [23]. This method has low resolution. The two mainstream approaches are migration and inversion by the forward model, orig-inally proposed for processing seismic signals [24]. Migration converts the hyperbola signal to the focused position [25]. Kirchhoff’s wave equation [26] and Fourier transform [27] based migration is extensively applied. Recent researches concern improving the migration approach (e.g., increasing accuracy of velocity estimation [28]) according to different application requirements and combining it with deep learning approaches [29]. Full-wave inversion establishes the forward model and adjusts the parameters to achieve the detected profile [30]. This method applies in simple models with fewer parameters (e.g., [31,32]), limited by high computational complexity. Recent researches concern accelerating inversion algorithms by deep learning approaches [33] and Bayesian inference [34]. In our research, we utilize the migration approach owing to the complicated situations of fractured rocks.

Wavelet transform (WT) is a time-frequency analysis approach to interpret signals [35]. It can distinguish signal discontinuities and decompose the signal into different mode functions containing distinct features. For the fractured surrounding rock identification, signal discontinuities appear at crack interfaces, and WT is expected to capture these characteristics. WT has demonstrated good de-clutter practices in processing GPR signals (e.g., [36,37,38]) but can only eliminate part of the hyperbola tails (e.g., [39,40]). Migration can handle hyperbola tails, but different migration approaches present different performance in terms of resolution and target intensity and sometimes preserve undesired scattering signals [20]. Besides, clutter also affects the migration results [41]. Therefore, some researchers complemented strategies to enhance the migration results (e.g., [42,43]). A typical migration approach is F-K migration, which operates frequency and wavenumber in the Fourier domain and then inverses to the time domain [44]. WT is functional-similar in signal analysis to the adopted Fourier transform and can provide good frequency resolution. Therefore, WT can assist F-K migration for clutter suppression and hyperbola elimination, and enhance the migration results.

In this paper, we propose a novel approach based on two-dimensional (2D) WT and F-K migration for identifying fractured rock areas using GPR. The 2D WT can decompose the radargram to different frequency levels and then feed the signals into F-K migration. The remainder of this paper is organized as follows: Section 2 describes the approach for processing GPR signals; Section 3 analyzes the results of numerical simulation and field measurements; Section 4 and Section 5 discuss and conclude this paper.

2. Materials and Methods

2.1. Signal Pre-Processing

Step 1: dewow filtering

Dewow filtering refers to the procedure of zero offset removal, which is primarily designed to eliminate the effect of direct current bias and the low-frequency signal trend contained in the data. The ‘wow’ phenomenon partially arises from the swamping or saturation of early arrivals and inductive coupling effect, thereby generating a distortion of the mean amplitude away from zero [25]. This can affect the subsequently processed results, particularly, the offset will be amplified to an indistinguishable extent in step 2. Efficient dewow filtering approach can be found in [25].

Step 2: Amplitude gain

Amplifying signal amplitude is required to enhance the appearance of later arrivals due to the attenuation from wave propagation and geometry spreading loss. Since the GPR antenna launches spherical waves, not plane waves, the signal spreads out, accelerating the attenuation. According to the electric field solution of Maxwell’s equations, the electromagnetic waves attenuate exponentially with propagation distance in a low-loss uniform medium [25]. Heavy scattering and refraction of the electromagnetic wave at defects constitute the geometry spreading loss. Based on the assumption of evenly distributed fissures, the wave energy decreases uniformly with the propagation distance. In this research, we combine the common gain functions, called exponential and linear functions [45,46], to compensate for the signal attenuation.

2.2. 2D Wavelet Decomposition

WT is a signal decomposition and time-frequency representation approach [47]. It utilizes function groups as band-pass filters to decompose the signals into different frequency levels. Since the wavelet functions are localized in the time domain, the decomposed results present the local modes in the time series and provide good time-frequency resolution. Identifying distinct local modes is significant for GPR radargram interpretation, as the distinct signals are reflected from the subsurface interfaces with dielectric contrasts. WT can decompose signals over oscillatory waveforms that reveal signal characteristics and provide sparse representations of regular signals that may include transients and singularities. Therefore, it has the potential to distinguish between clutter and targeted signals, as well as between hyperbola tails and focusing signals. The wavelet basis consisting of the scaling function ${\phi }_{0,0}$ and wavelets ${\psi }_{j,n}$ (Equation (1)) determines the decomposition paths. It can represent the signals by the wavelet inner-product coefficients and recover the signals by summing the coefficients (Equation (2)) [47].

$${\psi }_{j,n}\left(t\right)=\frac{1}{\sqrt{2^j}}\psi \left(\frac{t-2^{j}n}{2^j}\right),(j,n)\in {\mathbf{Z}}^{\mathbf{2}}$$

(1)

$$\begin{array}{c}\hfill \langle Y,{\psi }_{j,n}\rangle ={\int }_{-\infty }^{\infty }Y{\psi }_{j,n}\left(t\right)dt\\ \hfill Y=\sum _{j=-\infty }^{\infty }\sum _{n=-\infty }^{\infty }\langle Y,{\psi }_{j,n}\rangle {\psi }_{j,n}\end{array}$$

(2)

where wavelets ${\psi }_{j,n}$ that dilate and translate from the mother wavelet $\psi$ constitute the orthonormal basis, Y is the GPR signal function, $\langle -,-\rangle$ represents the inner product. An orthonormal basis is a complete orthonormal system for the Hilbert space, representing signals with no redundant information. Although different scaling functions and wavelet groups can have similar bandwidths, they will influence the local mode reconstruction, e.g., sharp curves or smooth oscillation. In this research, the fourth-order Daubechies wavelets are utilized to decompose GPR signals, and the scaling function and the mother wavelet are illustrated in Figure 1.

The down-sampling procedure after band-pass filtering at each wavelet decomposition stage reduces the signal resolution, which ignores minor signal responses and affects the identification of the interface locations from radargrams. Stationary wavelet transform (SWT) was designed to complement the translation-invariance, achieved by removing the down-samplers and up-samplers in the discrete wavelet transform and up-sampling the filter coefficients [48]. The signal size remains constant in the decomposition procedure. Keeping the signal size is significant, as the minimum recognizable scale is determined by the signal intervals. 2D SWT is appropriate for analyzing 2D datasets or pictures, as these signals have continuity in different orientations. The subsurface formations have both transverse and perpendicular continuities, and utilizing 2D SWT to interpret radargrams can reveal the horizontal relationships among diverse A-scan traces. The B-scan radargram is decomposed into multiple horizontal detail coefficients $D{H}_i$, vertical detail coefficients $D{V}_i$, diagonal detail coefficients $D{D}_i$ and approximation coefficients $A_i$, $i=1,2,\dots ,N$ (Equation (3)).

$$Y(x,t)=\sum _{i=1}^N(D{H}_i+D{V}_i+D{D}_i)+A_N,$$

(3)

where $Y(x,t)$ is the B-scan radargram and N is the decomposition number. In this research, we use the Matlab function ‘swt2’ for 2D SWT decomposition, and the sole parameter required is the fourth-order Daubechies wavelet.

2.3. F-K Migration

Frequency-wavenumber migration, or F-K migration, is an approach to convert hyperbola signals to object locations, firstly applied to process seismic signals [27]. It utilizes the exploding source model, where the scattered signal field is originated by the explosion at the object locations [49], to solve the wave equation (Equation (4)). F-K migration calculates the wave-field $\phi$ at $t=0$ when the explosion happens and the waves still locate at sources. The essence of F-K migration is the Fourier transform, which is derived from the general summation expression of a wave function (Equation (5)) [50].

$$(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial z^2}-\frac{1}{v^2}\frac{{\partial }^2}{\partial t^2})\phi (x,z,t)=0$$

(4)

$$\phi (x,z,t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\varphi (k_x,\omega )e^{-j(k_{x}x+k_{z}z-\omega t)}d{k}_{x}d\omega$$

(5)

where v is the propagation velocity (assumed constant in K-F migration) of electromagnetic waves, $k_x$ and $k_z$ are wave-numbers in the x and z directions, $\omega =v\sqrt{k_x^2+k_z^2}$ represents the frequency, and $\varphi$ is the Fourier transform from the surface field $\phi (x,0,t)$. $\phi (x,0,t)$ equals to the measured signal $Y(x,t)$ as GPR acquires the waves propagated to the surface. The target image is estimated by the initial wave-field $I(x,z)$ at $t=0$ (Equation (6)).

$$I(x,z)=\phi (x,z,0)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\varphi (k_x,\omega )e^{-j(k_{x}x+k_{z}z)}d{k}_{x}d\omega$$

(6)

By resampling $\varphi (k_x,\omega )$ to the $k_x-k_z$ domain, the migration results can be calculated by the inverse fast Fourier transform ${\mathcal{F}}^{-1}$ (Equation (7)). In this research, F-K migration is programmed by Matlab, and the parameter of propagation velocity can be calculated according to [25].

$$\begin{array}{c}\hfill I(x,z)={\mathcal{F}}_{k_x}^{-1}{\mathcal{F}}_{k_z}^{-1}\Phi (k_x,k_z)\\ \hfill \Phi (k_x,k_z)=\frac{v{k}_z}{\sqrt{k_x^2+k_z^2}}\varphi (k_x,v\sqrt{k_x^2+k_z^2})\end{array}$$

(7)

F-K migration is a linear function since the Fourier transform satisfies the linear superposition principle. The migration results of radargram can be regarded as the summation of migration from all SWT coefficients (Equation (8)).

$$M\left(Y(x,t)\right)=\sum _{i=1}^N(M\left(D{H}_i\right)+M\left(D{V}_i\right)+M\left(D{D}_i\right))+M\left(A_N\right),$$

(8)

where $M(·)$ represents the F-K migration function. The SWT coefficients have a good time-frequency resolution, and each component can contain clutter, targets, or hyperbola interference. Effective SWT signals containing the targets can enhance the migration results, while those occupied by clutter or hyperbola interference should be discarded. Thus the final target profile is reconstructed by selected migration components. Figure 2 presents the procedures of our proposed approach.

2.4. Test Site

Field measurements were carried out in an underground cavern in Huangdao, China, before grouting reinforcement, with the measurement scheme shown in Figure 3a. The underground caverns were built inside the layer of granite for structural stability. We use the GPR instrument from ‘MALA-Geoscience’, moving along the cavern surface and transmitting 250 MHz electromagnetic waves into surrounding rock sections. The wave frequency is 250 MHz to reach the depth around 7 m. The measurement length of an example section is $0\le x\le 20$ m, and the time window is 130 ns. Figure 3b presents the original radargram. Since opening the rocks and obtaining the actual rock distribution are difficult, a finely processed radargram by expertise is used as ‘ground truth’ (Figure 3c).

3. Results

3.1. Simulated Signal Analysis

Scenario 1: identifying sparse rock cracks

Since visibly obtaining the actual fissure distribution inside the surrounding rocks is difficult, the numerically generated radargram is analyzed to compare the reconstructed profile with the preset geometric model. The open-source software, ‘gprMax’ [51], can simulate electromagnetic wave propagation inside the subterranean sections and have been extensively applied in evaluating GPR signal processing approaches (e.g., [52,53]). We assume that the signal reflections are all generated by rock fissures, not strata interfaces, since the underground caverns are built inside the single layer of granite. In this paper, we simulated the GPR signals from a subterranean area with the measurement length $0\le x\le$ 2 m and depth $0\le z\le$ 1 m. The subterranean material is granite with relative permittivity of 7 [54] and conductivity of 0.012 [55]. The model includes six horizontal 5 mm-thick cracks (the direction is limited by ‘gprMax’, vertically or horizontally) filled by air, as shown in Figure 4a. The simulated grid step is $▵x=▵z=0.005$ m, and the time window is 20 ns. GPR transmits an 800 MHz ricker wave and receives A-scan traces with 0.01 m intervals. The simulated radargram after pre-processing is presented in Figure 4b. The hyperbola tails in the radargram affect the visibility of crack sizes and locations.

Firstly, the radargram is decomposed by 2D SWT, and Figure 5 illustrates the 2nd layer coefficients for an example. The approximation coefficient is not considered in migration since it preserves both target signals and hyperbola tails similar to the pre-processed radargram. Vertical and diagonal detail coefficients only contain hyperbola tails, while hyperbola centers corresponding to cracks are invisible. These two coefficients are excluded for migration. The horizontal detail coefficient preserves the target signal and discards part of the hyperbola tails, suitable for further F-K migration. Therefore, only horizontal detail coefficients $D{H}_i$ in Equation (8) are considered in the migration step.

The radargram is decomposed into four layers by 2D SWT. F-K migration processes the four horizontal detail coefficients and the pre-processed radargram for comparison, as shown in Figure 6. Without SWT decomposition, the migrated radargram (Figure 6e) highlights the crack locations, but preserves scattering signals, which are produced by hyperbola tails. Two-dimensional SWT has reduced the scatter in the migration results especially from the 3rd and 4th-layer horizontal coefficients (Figure 6c,d). However, when decomposing to deep SWT layers (e.g., 4th-layer in Figure 6d), the highlighted crack sizes increase. This is because the deep-layer signals have low frequency but large wavelengths. Therefore, $D{H}_3$ is more promising for avoiding scattering signals in migration and highlighting the crack locations and sizes. The radargrams contain artifacts marked in the dashed circles, which are unavoidable in all migration results. These artefacts are generated by multiple reflections between two upper cracks (center locations (0.52 m, 0.25 m) and (1.14 m, 0.45 m)).

Scenario 2: identifying sparse rock cracks from noisy environments

Noise is a significant issue affecting GPR signals. We add Gaussian white noise to the radargram to simulate the random noise. The same model in Figure 4a is utilized, and the polluted radargram (the signal-to-noise ratio is $-5$ db) is shown in Figure 7. Noisy speckles occupy the radargram, affecting the visibility of target signals and hyperbola tails.

The polluted radargram is decomposed into four layers by 2D SWT. Similar to scenario 1, only horizontal detail coefficients preserve target signals and mitigate part of hyperbola tails, as shown in Figure 8. The radargram is not decomposed into deeper layers that increase the crack sizes in the migration result. The coefficients $D{H}_1$ and $D{H}_2$ are occupied by noise, which means 2D SWT can separate noise out. $D{H}_3$ is a mixture of noise and target signals, unsuitable for crack identification. Noise nearly disappears in $D{H}_4$, and target signals become visible. Noise, target signals, and hyperbola tails have different time-frequency features, and therefore they can be separated by 2D SWT.

The F-K approach migrates the polluted radargram in Figure 7 and the 4th-layer $D{H}_4$, as illustrated in Figure 9. Although highlighting the crack locations, Figure 9a contains a similar noise level to Figure 7. Besides, scattering signals appear surrounding the focused locations in Figure 9a. The noise intensity decreases dramatically in Figure 9b since 2D SWT has separated the noisy layers. The migration result from $D{H}_4$ has highlighted the target locations and eliminated noise and hyperbola tails. However, the identified crack sizes increase due to the large wavelength in deep SWT layers.

Scenario 3: identifying fractured rock areas with dense fissures

Fissure groups other than the single throughout crack are common in fractured rock areas. A subterranean area with the measurement length $0\le x\le$ 2 m and depth $0\le z\le$ 1 m is simulated by ‘gprMax’. There are three predefined fractured rock portions with rock fissures randomly distributed inside (fissure widths of 5 mm, lengths of 1–5 cm, and intervals of 1–5 cm), as shown in Figure 10a. Other model parameters are the same as scenario 1. Figure 10b illustrates the corresponding radargram after pre-processing. Dense hyperbola tails have occupied large subsurface areas, affecting the identification of fractured rock areas.

The radargram is firstly decomposed by 2D SWT, and Figure 11 illustrates the 2nd-layer coefficients for an example. Similar to scenario 1, the approximation coefficient contain both target signals and hyperbola tails. Vertical and diagonal detail coefficients only contain hyperbola tails, while hyperbola centers corresponding to fractured rocks are invisible. These coefficients are not considered for migration. The horizontal detail coefficient preserves target signals and mitigates part of the hyperbola tails. Therefore, only horizontal detail coefficients $DH$ are suitable for migration to identify the fractured rocks.

The radargram is decomposed into four layers by 2D SWT. F-K migration processes the four horizontal detail coefficients and the pre-processed radargram for comparison, as shown in Figure 12. Without SWT, scattering signals appear surrounding the fractured rocks highlighted in the migration result (Figure 12e). Migration from horizontal detail coefficients, especially of the deep layers, can mitigate the scattering signals. Artefacts (marked in the dashed circles) arising from multiple reflections also exist in the migration results, but SWT decomposition demonstrates potential to reduce this effect in $D{H}_3$ and $D{H}_4$. Different from scenario 1, the wavelength effect on identifying crack sizes is less significant when handling fractured rock areas. Limited by ‘gprMax’ software, the simulated cracks are horizontal, and $D{H}_3$ and $D{H}_4$ have revealed the horizontal directions of the cracks. Both $D{H}_3$ and $D{H}_4$ present promising results in highlighting the fractured areas and eliminating scattering signals.

Scenario 4: identifying fractured rock areas with dense fissures from noisy environment

Similarly, we consider the noise effect by adding Gaussian white noise to the radargram (Figure 10b), with the polluted profile (signal-to-noise ratio −5 db) shown in Figure 13. Noisy speckles occupy the radargram, affecting the visibility of target signals and hyperbola tails.

The polluted radargram is decomposed into four layers by 2D SWT. Similar to previous scenarios, horizontal detail coefficients that preserve target signals and mitigate part of hyperbola tails are illustrated in Figure 14. The 2D SWT has separated noise to $D{H}_1$ and $D{H}_2$, and the noise nearly disappear in $D{H}_4$. Therefore, 2D SWT can distinguish noise, targets, and hyperbola tails.

The F-K approach migrates the polluted radargram in Figure 13 and the 4th-layer $D{H}_4$, as illustrated in Figure 15. Although highlighting the crack locations, Figure 15a contains noise speckles polluting the crack signals and scattering signals appear surrounding the focused locations. The noise intensity decreases dramatically in Figure 15b since 2D SWT has separated the noisy layers. The migration result from $D{H}_4$ has highlighted the target locations and eliminated noise and hyperbola tails. Similar to scenario 3, the wavelength effect on crack lengths is little when handling fractured rock areas.

3.2. Field Signal Analysis

Figure 16 presents the pre-processed radargram from field measurement. Signals weaker than 2.5% of the maximum amplitude are excluded for better showing the fractured areas. These weak signals are mostly noise, which will not affect fractured rock identification but block the visible presentation from figures. The entire section is covered by intensive signals except for the upper-right portion (with the length $12\le x\le 17$ m and time $20\le t\le 65$ ns), which makes the fractured rocks indistinguishable.

Firstly, the pre-processed radargram is decomposed by 2D SWT, with the 2nd-layer coefficients shown in Figure 17. Similar to the simulation results, hyperbola tails have been decomposed to $D{V}_2$ and $D{D}_2$. The approximation coefficient presents similar distribution to the pre-processed radargram. Only $D{H}_2$ discards part of the noise and hyperbola tails, and preserves the target signals. Therefore, the horizontal coefficients are considered for further migration.

The F-K approach migrates $D{H}_1$, $D{H}_2$, $D{H}_3$, and the pre-processed radargram for comparison, as illustrated in Figure 18. The signals inside region U and the lower part of region V are weak in Figure 3c, but direct F-K migration preserves intensive signals (Figure 18d). This is affected by clutter. The 2D SWT can mitigate the clutter and thus the signals inside region U and the lower part of region V in Figure 18c become less intensive. Region W is identified as severely fractured rock areas. 2D SWT has demonstrated good practice in enhancing the migration results.

4. Discussion

In this paper, a novel approach based on 2D WT and F-K migration is proposed to identify fractured rock areas from GPR radargrams. WT is theoretically similar to Fourier transform, the foundation of F-K migration, but alters the orthogonal basis. Therefore in signal processing, WT usually provides better time-frequency resolution [47]. By applying the proposed approach, three significant development has been achieved: (1) Noise, hy- perbola tails, and target signals are separated into different WT coefficients; (2) Scattered signals surrounding cracks in F-K migration results have been eliminated; (3) Noise has been mitigated in the focused radargram. The remained scattering energy in migrated profiles does not affect the identification of sparse objects (e.g., [44,56] and Figure 6e) and strata (e.g., [57]), but increases the identified crack group areas (Figure 12b) significantly. Therefore, eliminating the remained scattered signal is important for our application.

One phenomenon that arises from 2D WT is the wavelength effect. When identifying sparse cracks (Figure 6d and Figure 9b), the imaged crack length increases. The reason is that the wavelength increases as WT decomposes signals to deeper low-frequency components. Normally, the imaged crack size on an A-scan trace is influenced by the antenna frequency [58] (imaged size = wave-number × wavelength) and thus does not equal the actual size vertically. The 2D WT decomposes the radargram vertically and horizontally, and thereby the two-direction wavelengths vary in different decomposition layers. This results in the horizontal deformation of imaged cracks. To accurately calculate the object sizes, the frequency of processed radargram needs to compare with the antenna frequency vertically and with the original signal frequency horizontally. Nonetheless, the wavelength effect is little when identifying crack groups (Figure 12d) and thus affect little on our applications.

In our research, we apply 800 MHz waves in simulation and 250 MHz waves in field measurements. The choice of central frequency is determined by the required detection depth and resolution. The high-frequency electromagnetic waves have high resolution, but the energy attenuates quickly [25]. If using 250 MHz waves in simulation, the major difference is that the imaged cracks become thicker due to large wavelengths.

Our future researches will concern practical investigations of oblique and non-rectilinear rock fractures. The crack directions are limited by grids in numerical simulation and thus cannot cover all the crack conditions. In literature, long and throughout cracks (e.g., [59,60]) in different directions has been investigated in practical cases, but the conditions are not same as the fractured rocks. Although our approach can process radargrams from simulating other objects, including ‘cylinders’ and ‘spheres’, simulation of oblique and non-rectilinear rock fractures are impossible. Therefore, further researches should include more physical experiments to cover all the rock conditions.

The physical experiments are also important for validation since we cannot open the entire rock section in field measurements. The cracks or the fractured rock areas with different rock qualities can be preset, and thus our approach can be further evaluated. On this basis, we can establish the relationship between GPR radargrams and rock quality values in criterion.

5. Conclusions

In this paper, a novel approach based on 2D SWT and F-K migration is proposed to identify fractured rock areas from GPR radargrams. The 2D SWT can mitigate clutter, distinguish signal discontinuity, and provide signals with a good time-frequency resolution. It decomposes the radargram to different frequency levels and then feeds the signals into F-K migration. F-K migration can focus the scattered hyperbola signals back to the target locations. In practice, the detail coefficients achieved by 2D SWT can contain clutter, hyperbola tails, and target signals separately, and thus they are selected for the migration step.

The simulation and measurement results demonstrate that: comparing with direct F-K migration results, migration from the horizontal detail coefficients by 2D SWT can mitigate scattering signals and eliminate noise. 2D SWT decomposes the hyperbola tails to the vertical and diagonal detail coefficients, and separates the noise to the first few horizontal detail coefficients. Therefore, the corresponding migration results of later $DH$ highlight the crack locations and directions with little signal interference. The wavelength effect affects the identification of single cracks but not the crack groups. When decomposing to deep SWT layers, the frequency decreases and the wavelength increases. This results in the increased scale corresponding to the crack sizes. Nonetheless, the promising results indicate the effectiveness in identifying fractured rock areas. Future research will concern more practical cases since the simulations of oblique and non-rectilinear rock fractures are impossible.