Identification of Unstable Subsurface Rock Structure Using Ground Penetrating Radar: An EEMD-Based Processing Method

Featured Application

The approach for extracting features from GPR radargrams was developed based on EEMD. The proposed method was utilized to eliminate the interference signals and reconstruct subsurface profiles. This radargram interpretation approach was applied to identify surrounding cracked rocks of underground caverns.

Abstract

Surrounding rock quality of underground caverns is crucial to structural safety and stability in geological engineering. Classic measures for rock quality investigation are destructive and time consuming, and therefore technology evolution for efficiently evaluating rock quality is significantly required. In this paper, the non-destructive technology ground penetrating radar (GPR) assisted by an ensemble empirical mode decomposition (EEMD)-based signal processing approach is investigated for identifying unstable subsurface rock structures. By decomposing the pre-processed GPR signals into multiple intrinsic mode functions (IMFs) and residues, one typical IMF can preserve the distinct local modes and is considered to reconstruct the subterranean profile. Promising results have been achieved in simple scenarios and filed measurements. The reconstructed profiles can accurately illustrate the subsurface interfaces and eliminate the interference signals. Unstable rock structures have been identified in further field applications. Therefore, the developed approach is efficient in unstable rock structure identification.

1. Introduction

With the rapid development of underground spaces for energy storage [1], transport [2] and mining [3], the surrounding rock stability is increasingly significant for the safe operation of underground caverns. Rock formation can fundamentally determine the geological strength [4], as cracked surrounding rocks risk collapse and require rock reinforcement approaches. To visibly evaluate the surrounding rock quality, opening and pit sampling approaches (e.g., [5,6]) are extensively utilized and can accurately quantify the attribute parameters. However, intermittent sampling fails to detect the entire rock section continuously, while mechanically opening the rocks is destructive and time-consuming. Some researches (e.g., [7,8]) established sensor arrays for stress-strain monitoring, but the detection positions were discontinuously distributed inside the shallow ground. Therefore, developing a field-efficacious non-destructive detection technology to scan the surrounding rock profiles continuously is significant for rock structural stability.

Ground penetrating radar (GPR) is a remote sensing technology originally designed to detect the ice thickness while primarily applied in discovering the unseen objects beneath the earth’s surface, e.g., with applications in mapping subsurface stratigraphy [9], archaeological prospection [10] and detecting pavement defects [11]. Since GPR works as an echo-listener recording electromagnetic waves reflected from interfaces with dielectric differences, this non-destructive technology has the potential to distinguish fresh and weathered rock masses. However, signal interference is an unavoidable issue affecting radargram interpretation. Inversion and migration approaches (e.g., [12,13]) were developed based on Maxwell’s equations to eliminate signal interference, but excessive derivative and integral calculations are time-consuming. Machine learning (e.g., [14,15]) provides automatic approaches for target identification from radargrams, but required considerable datasets with similar underground attributes are difficult to acquire from field measurement. Jin and Duan [16] developed a discrete grid method to eliminate hyperbolic interference in rock detection, but with many parameters determined by experience. Since electromagnetic waves would be reflected and interfered frequently in cracked rock sections, a novel approach is required to improve efficiency in handling signal interference.

Modal analysis is a common signal analysis approach in structural health monitoring. Generally, it characterizes the vibration signals in frequency and time domains, which provide attribute and stability information of the tested structure. The modal difference exists between intact and cracked rock masses, and GPR can capture this subterranean structural feature. Different from classic approaches that concern the physics of wave propagation for GPR signal processing, modal analysis methods consider the physics of signal strength and frequencies. Some researches have attempted to interpret radargrams with these methods, e.g., wavelet transform [17,18] and empirical mode decomposition [19,20], but their applications are either at the de-noise stage or in simple object identification. To our best knowledge, no research has modally analyzed GPR radargrams with complicated targets like cracked rock masses. In the cracked subsurface area, numerous rock fissures are nearly continuously distributed, therefore instantaneous changes of GPR signals becoming significant in modal analysis. Empirical mode decomposition (EMD) was originally proposed by Huang et al. [21] for nonlinear and non-stationary signal analysis and then extensively applied to extract modes from signals (e.g., [22,23]). It can decompose signals into multiple intrinsic mode functions (IMFs) which focus on the local time scale and contain instantaneous frequency information, thereby having the potential to analyze instantaneous changes in radargrams. The decomposition results correspond to the natural mode determined by the signal itself, different from wavelet transform that acquires modes inside certain filter bands. Mode mixing sometimes appears and reduces the decomposition effectiveness, and in order to overcome this disadvantage, ensemble empirical mode decomposition (EEMD) was developed utilizing a white-noise-assisted method [24].

In this paper, an EEMD-based approach is developed to analyze GPR data for the novel application in identifying unstable subsurface rock structures. The remainder of this paper is organized as follows: Section 2 introduces the developed signal processing approach for interpreting GPR radargrams; Section 3 presents the results of both numerically simulated and field-measured GPR signals; Section 4 discusses and concludes this paper.

2. Materials and Methods

2.1. Pre-Processing of GPR Data

To prepare GPR signals well for further EEMD analysis, three significant steps are utilized to pre-process original GPR data: time-depth conversion, ‘de-wow’ and amplitude gain.

Firstly, the time-depth conversion is designed to extract location information from the recorded time series. According to propagation velocity invariant assumption and two-way travelling time (transmission and reflection) of electromagnetic waves in the subterranean area, the recorded depth of interfaces can be determined [25]:

$$Z=\frac{vt}{2}=\frac{ct}{2\sqrt{{\mu }_r{\epsilon }_r}}$$

(1)

where, v is the propagation velocity of electromagnetic waves, t is the two-way travelling time, c represents the propagation velocity in the vacuum, ${\mu }_r$ is the relative magnetic permeability, and ${\epsilon }_r$ is the relative permittivity.

Secondly, ‘de-wow’ processing is utilized to eliminate low-frequency distortions resulting from the inductive coupling effect or early saturation [25]. An effective approach with good practice to handle the ‘wow’ issue is simply removing the average amplitude from the signals [26]:

$$A(Z)=a(Z)-\frac{1}{N}\sum _{k=1}^{N}a(Z_k)$$

(2)

where, N is the number of discrete signal points, $a(Z_k)$ represents the original signal amplitude, and $A(Z)$ represents the processed signal.

Thirdly, amplitude gain is utilized to amplify the signals against energy loss including propagation attenuation and geological spreading loss [25]. Considering both exponential attenuation and evenly distributed rock fissures, radar signals can be recovered to realistic levels by multiplying Equation (3) [16].

$$F(Z)=e^{2\alpha (Z-Z_0)}+2\phi (Z-Z_0)$$

(3)

where, $\alpha$ is the exponential attenuation parameter, $\phi$ is the geological spreading loss parameter, and $Z_0$ is the starting depth for amplitude gain.

2.2. Ensemble Empirical Mode Decomposition

In 1998, Huang et al. [21] proposed EMD for Hilbert transform with physical meanings to analyze nonlinear and non-stationary signals in multiple resolutions. It offers a self-adaptive approach to decompose any complicated signal series into a finite number of IMFs, which focus on the local time scale and involve only one oscillation mode in each waveform. The instantaneous frequency defined in an IMF corresponds well with the local physics of the studied system, which would have an expected performance in interpreting short-distance features from GPR radargrams.

According to the definition, an arbitrary signal $X(Z)$ and its Hilbert transform $Y(Z)$ form the complex conjugate pair to obtain the analytical signal $P(Z)$.

$$P(Z)=X(Z)+iY(Z)=a(Z)e^{i\theta (Z)}$$

(4)

where, i represents the imaginary unit, $a(Z)=\sqrt{X^2(Z)+Y^2(Z)}$ is the signal amplitude, and $\theta (Z)=arctan(\frac{Y(Z)}{X(Z)})$ is the signal phase. Therefore, the instantaneous frequency is defined as:

$$\omega =\frac{d\theta (Z)}{dZ}$$

(5)

In order to handle the issue that the instantaneous frequency has no precise calculation approach, signal data are limited by narrow bands to provide physical meanings for the instantaneous frequency, and these narrow bands can be obtained by EMD procedure.

Processing a single signal series, the EMD method applies the following procedures:

• Firstly record the initial data, $r_1=A(Z)$, and initialize from $i=1$, $j=1$;
• Find the upper and lower envelopes ($U_j,L_j$) of the data by cubic spline interpolation and then calculate their means, $m_j=(U_j+L_j)/2$;
• Setting $h_{i,1}=r_i$, obtain the rest part by removing the means $m_j$, $h_{i,j+1}=h_{i,j}-m_j$, $j=j+1$;
• Repeat the procedure ii. and iii. until $h_{i,j+1}$ satisfies the decomposition threshold, and we can obtain the $ith$ IMF as $C_i=h_{i,j+1}$. In Huang’s method, the threshold parameter $SD$ should be between 0.2 and 0.3;

  $$SD=\sum _{k=1}^j\frac{m_k^2}{h_{i,j}^2}$$

  (6)
• $i=i+1$, calculate the rest part of signals $r_i=r_{i-1}-C_{i-1}$. Repeat the above steps ii. to iv. to obtain diverse IMFs and the final residue data $r_n$.

  $$A(Z)=\sum _{i=1}^{n-1}{C}_i+r_n$$

  (7)

After the EMD algorithm, the analytical signal for Hilbert transform with physical meanings is expressed as:

$$P(Z)=\sum _{i=1}^{n-1}{C}_i{e}^{i\int {\omega }_i(Z)dZ}$$

(8)

However, the significant drawback of EMD, mode mixing appearing frequently, has depressed its performance and restricted the applications. This effect either mingles signals of widely disparate scales in an IMF, or residents similar-scale signals in different IMFs. A noise assisted method, EEMD, was then developed by Wu and Huang [24] to overcome the problem of mode mixing. In the EEMD approach, the true IMFs are defined as the ensemble mean of multiple decomposition trials, each with initialized data of the original signals plus a random white noise with narrow amplitude, and this thereby improves mode separating as other researches demonstrate [27,28].

The EEMD algorithm can be summarized as follows:

• Initializing from $i=1$, generate a random white noise and add it to pre-processed signals, $A_i(Z)=A(Z)+nois{e}_i$. The standard deviation of noise amplitudes is suggested to be 1/5 that of the signals;
• Apply EMD to decompose $A_i(Z)$ into different IMFs $C_{i,1},C_{i,2},\dots ,C_{i,n-1}$;
• $i=i+1$, repeat the above two steps for N times, with the minimum value of N should satisfy Equation (9), where ${\epsilon }_n$ is the final standard deviation of error arising from added noises, and $\epsilon$ is the standard deviation of noise amplitude;

  $$ln{\epsilon }_n+\frac{\epsilon }{2}lnN=0$$

  (9)
• Obtain the true decomposition results by calculating the ensemble means of all corresponding IMFs.

  $$C_k=\frac{1}{N}\sum _{i=1}^N{C}_{i,k},k=1,2,\dots ,n-1$$

  (10)

In our signal processing procedures, EEMD is employed to extract the distinct mode which should arise from fissures or other shattered rock interfaces. Figure 1 shows representative GPR signals from two measurement traces with and without buried objects beneath respectively and their fast Fourier transform results. The distinct mode appears visibly around the depth of buried objects (0.6–0.8 m), which would be evident local modal features for IMFs to capture.

The algorithm for processing GPR data utilizing EEMD is developed:

• Pre-process GPR profiles using time-depth conversion, ‘de-wow’ and amplitude gain;
• Initializing from $j=1$, obtain the vertical signal series from a single measurement trace $A(x_j,Z)$, where $x_j$ is the abscissa of the measurement point;
• Apply EEMD to decompose $A(x_j,Z)$ into different IMFs, and preserve the IMF that can distinguish the mode of abnormal objects;
• $j=j+1$, repeat the above two steps to process the vertical signals from all measurement traces, and then we can reconstruct the final analytical subterranean profile by all preserved IMFs.

3. Results

3.1. Investigation in Simple Scenarios

The basic approach for identifying the unstable rock structures is to identify the rock interfaces, which is firstly investigated in this paper. Two simple scenarios, a practical radargram with buried objects and a numerically simulated radargram with predefined cracked rock areas, are under investigation to evaluate the accuracy of the developed signal processing approach in locating object interfaces.

Firstly, GPR acquired the subterranean profile with two pre-buried pipelines (32 cm diameters), which are sparsely distributed with around 2.5 m interval. The pipelines are discernible enough and will not interfere each other’s signals because of the large burying distance. The central frequency of GPR attenna is 700 MHz. Figure 2 illustrates the pre-processed radargram, with the measurement length L ($0\le L\le 6$ m) and detection depth Z ($0\le Z\le 2.5$ m). The pre-processing parameters are $\alpha =0.56$, $\phi =0.1$ and $Z_0$ = 0.25 m. Two red circles are manually supplemented to the illustration corresponding to the locations and sizes of the pipelines. Although hyperbolic signals appear at pipeline locations, their two ‘pseudo’ tails (marked as dashed boxes in Figure 2) have extended outside the pipeline interfaces, which will result in misidentifying normal areas as our targets. These hyperbolic interference signals are generated by the wide-angle transmission of electromagnetic waves [16]. Intensive signals in the near-surface portions (with depth less than 0.75 m) are polluted by noises and without consideration in this research.

Two representative signal traces (seen in Figure 3) on Lines 1 and 2 with and without buried pipelines respectively are compared to analyze the mode difference. Distinct local modes appear around $Z=1$ m due to the existence of pipelines. Figure 4 shows the decomposition results by EEMD on both traces. IMF1s are noise-polluted mode functions with limited physical meanings because white noises are added to signals in EEMD procedures. IMF2s, IMF5s and IMF6s present little evidence of mode difference between Lines 1 and 2, thereby not suitable for further object identification. This missing characteristic reappears in IMF3s and IMF4s, as visible oscillations exist around Z = 1 m on Line 1 (corresponding positions marked in Figure 4). Therefore, IMF3s and IMF4s on all traces of the pre-processed radargram are chosen to reconstruct the subsurface profile.

Figure 5 illustrates the reconstructed profiles by IMF3s and IMF4s respectively, which are decomposed from all signal traces. It is reasonable to eliminate the IMF signals weaker than the added white noises (amplitudes less than 500) to avoid artificial noise issues, and this would affect little on the actual modes (strength 20–40 times). In the IMF3s profile (Figure 5a), the two hyperbolic ‘tails’ of each pipeline have been eliminated. Meanwhile, the pipelines have been revealed from the noisy subterranean environment with highlighted locations and sizes. In contrast, these commendable results have not appeared as expected in the IMF4s profile (Figure 5b). Interference signals remain outside the pipeline interfaces, and near-surface signals have extended to deeper areas (from 0.75 m in Figure 2 to 1 m in Figure 5b) thereby polluting the pipeline signals. Although signals inside the pipeline portions are intensive, this reconstructed IMF4s profile is insufficient to identify the objects. Therefore, only IMF3s acquired by the developed approach can reconstruct a promising underground profile. However, unexpected visible noises appear in the deepest 10% portion, which may arise from the end effects of EEMD [29], and GPR time windows should be at least 10% larger than targeted depths in further applications.

Secondly, a numerically simulated radargram of cracked rock areas was generated by the software ‘gprMax’ [30]. ‘gprMax’ is an open source software to simulate GPR signals based on propagation physics of electromagnetic waves and the Finite-Difference Time-Domain method, widely applied in investigating signal processing approaches (e.g., [31,32]). Numerical simulation can determine the isotropy and homogeneity of the underground medium, thus ensuring object signals are not interfered by other unknown structures. Figure 6a illustrates the simulated rock profile, with the scanning length L ($0\le L\le 2$ m) and detection depth Z ($0\le Z\le 1$ m). The entire section is divided into mesh grids with resolution $0.005\times 0.0025$ m. Certain attributes have been assigned to the rock section, with the relative permittivity of 8, the conductivity of 0.02 S/m, the relative permeability of 1 and the magnetic loss of 0. There are three predefined cracked portions (sizes seen in the caption of Figure 6), each containing multiple randomly distributed cracks with widths of 2.5 mm, lengths and intervals ranging from 1 cm to 5 cm. These cracks are simulated as perfect electric conductors. The GPR antenna is simulated by the Ricker waveform with the central frequency of 800 MHz, while the detection time window is $4\times {10}^{-8}$ s. Figure 6b presents the corresponding pre-processed radargram. The pre-processing parameters are $\alpha =0.13$, $\phi =0.2$ and $Z_0=0.1$ m. The predefined cracked areas generate not only intensive responses at their locations but also interference signals affecting other locations. Over half of the profile is occupied by GPR signals, which increases difficulties in cracked rock identification.

Two representative GPR traces on Lines 1 and 2 with cracks at different depths are compared to analyze the mode difference, as shown in Figure 7. The shallow cracked portion on Line 1 generates distinct local modes between 0.1 and 0.4 m depths (marked in the dotted circle). EEMD results of both traces are illustrated in Figure 8. Only mode functions after IMF5s are presented as IMF1s to IMF4s are entirely occupied by added white noises. IMF6s and IMF7s have captured the modal characteristics of both signal traces, as marked in dotted circles. However, the oscillations in IMF7s have extended beyond the actual crack depths. Therefore, IMF6s are chosen to reconstruct the GPR profile.

Figure 9 presents the reconstructed profile by IMF6s. Signals with amplitude weaker than 3 are eliminated since noises would affect our identification results. In the reconstructed radargram, responses are concentrated inside the cracked portions (marked as dashed rectangles), and interference signals have been almost eliminated. However, some weak signals surrounding the cracked portions are preserved (inside portions B, C, D and E) since eliminating all interference signals is difficult. Nonetheless, the intact portion F, which was covered by intensive signals in Figure 6b, has been revealed by the proposed approach. These results indicate the applicability of EEMD in distinguishing the intact and cracked rock areas. End effects sometimes occur in our decomposition procedure thus generating noises in the portion G, which indicates that GPR time windows should be at least 10% larger than targeted depths to avoid these effects.

3.2. Investigation in Field Measurement

Filed measurement was conducted in an underground cavern in Huangdao, China to evaluate the surrounding rock stability before grouting reinforcement. Figure 10 illustrates a brief scheme of detecting surrounding rock sections. A GPR instrument is attached on the cavern wall to transmit electromagnetic waves into the surrounding rocks. Once encountering the rock cracks, the electromagnetic waves will be reflected back to the receiving antenna (RA) and the response signals are recorded on the vertical trace. Owing to the wide transmission angle of the transmitting antenna (TA), reflections from non-vertical directions will generate ‘pseudo’ signals that affect the identification of those from vertical directions. Details of basic detection principles and ‘pseudo’ signals can be seen in [25]. During the investigation of the entire surrounding rock section, the GPR instrument will continuously move along the measurement direction while acquiring numerous vertical traces. In our research, a MALA 2-dimensional GPR instrument assisted by the data acquisition software ‘Groundvision’ is adopted. The central frequency of the GPR antenna is 250 MHz, and the horizontal resolution (intervals between neighbouring vertical traces) is 3 cm.

A representative subterranean section is firstly investigated to analyze the detectability of shattered rock masses by the developed approach, which is compared with the classic migration approaches (e.g., [33,34,35]) utilizing the commercial software ‘Reflexw’. The original radargram is presented in Figure 11a. Although strong reflections occur in the near-surface layers (amplitudes over 20,000), signals beneath 2 m are too weak (amplitudes below 100) for post-processing. The visible distribution in Figure 11a will differ significantly from the actual distribution due to electromagnetic attenuation, and later arrivals beneath 6 m nearly disappear in this situation. Nonetheless, the surrounding rock conditions are complicated as intensive responses have been received from most areas. Figure 11b shows the pre-processed radargram of the surrounding rock section, with the measurement length L ($0\le L\le 20$ m) and detection depth Z ($0\le Z\le 6.86$ m). The pre-processing parameters are $\alpha =0.56$, $\phi =0.1$ and $Z_0=0.69$ m. The depth-variant amplitude gain function has reconstructed the signals to normal levels (amplitudes around 10,000), while signals weaker than 200 (2% of normal amplitudes) have been eliminated to avoid noises. Although the right-middle area H is conveniently recognized as intact rock masses, other areas are almost covered by intensive signals, which increases difficulties in identifying the unstable structures with fragmented rocks.

Two representative measurement traces from Figure 11b, Line 1 with intermittent signals and Line 2 with later heavy signals, are decomposed to compare the local mode difference, as shown in Figure 12 and Figure 13. Besides the added white noises, IMF1s include the local mode in the near-surface layers, but this aberrant strong reflection area is not considered in the identification. Similarly, the compositions after IMF4s are also mostly occupied by the near-surface signals, unable to recognize the joint fissures at deep locations. In contrast, the IMF2 of Line 2 has captured the discernible fluctuation after 3.5 m, and this composition is chosen for further analysis.

Utilizing EEMD to decompose all vertical measurement traces, the IMF2s are preserved to reconstruct the amplitude profile of the surrounding rock section, seen in Figure 14. like the previous procedure, the IMF2 signals weaker than the added white noise (amplitudes less than 200) has been removed. For comparison, Figure 15 presents the imaging results finely processed by the commercial GPR software. The distribution of signals in areas U, W, X, Y (Figure 14) is consistent with the corresponding areas $U1$, $W1$, $X1$, $Y1$ (Figure 15), while signals in lower part of area V are not as intensive as the area $V1$. It is noticeable that comparing with the pre-processed profile (Figure 11), the noisy signals in area Y nearly disappear, which implies the high surrounding rock quality in this area. Although W and X are concentrated areas of fragmented rocks, the lower part of W presents higher rock quality after EEMD processing and the area X is not as cracked as shown in Figure 11. These results indicate that the developed approach has the potential to discover the subterranean joint fissures to further identify the surrounding rock stability.

3.3. Identifying Unstable Rock Structures

Two surrounding rock sections were investigated for further applying the developed approach in identifying unstable rock structures.

Figure 16 shows the original and pre-processed radargrams of a surrounding rock section, with the measurement length L ($0\le L\le 16$ m) and detection depth Z ($0\le Z\le 6.86$ m), and the reconstructed profile by IMF2s is illustrated in Figure 17. The original radargram is insufficient for post-processing as the signals beneath 2 m are over 200 times weaker than those in the near-surface layer (amplitudes over 20,000). Later arrivals beneath 6 m are nearly invisible owing to the electromagnetic attenuation. This issue has been handled by the amplitude gain function, with pre-processed results shown in Figure 16b. However, rock identification is still difficult since intensive interference signals have nearly occupied the entire section. Contrarily, the cracked and intact rock masses are distinguishable in the reconstructed profile (Figure 17). The intact portions I, J and K that present evidence in Figure 16b have been entirely revealed, while the intact portions M and N covered by intensive signals in Figure 16b also become visible. Other areas are identified as unstable rock structures with different cracked situations, thus requiring further reinforcement measures. The deepest 10% area marked in the red rectangle is difficult to identify since the end effects of EEMD would affect this area.

Figure 18 shows the original and pre-processed radargrams of another surrounding rock section, with the measurement length L ($0\le L\le 18$ m) and detection depth Z ($0\le Z\le 6.26$ m), and the reconstructed profile by IMF2s is illustrated in Figure 19. Similar to the previous results, the signals beneath 2 m in the original radargram are over 200 times weaker than those from the near-surface layer (amplitudes over 20,000), which results from the electromagnetic attenuation. This issue has been handled by the amplitude gain function, with pre-processed results shown in Figure 18b. However, intensive signals have nearly occupied the entire section, making our identification difficult. Stable rock structures have been revealed after reconstructing the radargram. Different from the previous rock section in Figure 18, the deep portion R around 4–6 m of this section is almost cracked, which means the rock structure is unstable. The intact portions P and Q that show evidence in Figure 18b have been entirely revealed, while the intact portion O covered by intensive signals in Figure 18b also become visible. Therefore, our developed approach shows good practice in unstable rock structure identification. The deepest signals in the red rectangle are affected by the end effects of EEMD, and therefore the GPR time window should be 10% larger than the targeted depth in further applications.

4. Discussion and Conclusions

Two significant research developments are achieved in this paper: successfully developing mode analysis approach for complicated radargram interpretation, and successfully applying the non-destructive technology GPR for unstable rock structure identification.

Since surrounding rock quality attracted attention in underground engineering, standards to evaluate rock quality has continuously evolved, but opening and pit sampling has remained as the field detection strategy yet [36]. Although seismic velocity measurements have been utilized to indirectly estimate rock quality [37], this energy required technology is not suitable in underground caverns. GPR has been applied in object detection, stratum identification and underground water detection, but rarely in such a complicated condition with intensive wave reflection from intensive rock interfaces. Modal analysis has been applied to interpret simple GPR radargrams before, e.g., wavelet transform for de-noising [17] or pipeline identification [18] and EMD for object detection [20], but not in complicated conditions like cracked subsurface rock masses.

In this paper, an EEMD-based signal processing approach including three steps, pre-processing, signal decomposition and profile reconstruction, is investigated for interpreting GPR radargrams for unstable subsurface rock structure identification. Amplitude gain function has restored the signals to an analyzable level, against the electromagnetic attenuation which reduces the energy of responses beneath 2 m to 200 times weaker than that in the near-surface layer. The pre-processed radargrams present some different distributions with the original ones, as well as the evidence of intact rock portions (e.g., portion H in Figure 11b). By decomposing the pre-processed GPR signals into multiple IMFs and residues, one typical IMF can preserve the distinct local modes, which arise from instantaneous changes of fissure interfaces, and is chosen to reconstruct the subterranean profile. Different IMFs can present different signal features, including noises, near-surface-layer signals and object responses, while object responses can appear in different IMFs. For example, pipeline responses are preserved in the IMF3 and IMF4 in Figure 4, but only IMF3s profile can accurately show pipeline areas; simulated crack responses are preserved in the IMF6 and IMF7 in Figure 8, but oscillations in the IMF7 have extended outside the cracked portions. Choosing appropriate IMFs to reconstruct the profile is significant for further unstable rock identification.

After considering appropriate IMFs to reconstruct the radargram, the intact and cracked rock portions become distinguishable. The developed approach can eliminate the hyperbolic interference signals, but some weak signals remain surrounding the cracked portions in the numerical simulation (Figure 9). Despite these small biases, the identified cracked areas correspond well with the pre-defined ones, and the intact portion G covered by intensive signals become visible in Figure 9. In the field measurements, the intact areas that have shown evidence in pre-processed radargrams are entirely revealed (e.g., I, J and K in Figure 17, P and Q in Figure 19), while those covered by intensive signals in pre-processed radargrams also become visible (e.g., Y in Figure 14, M and N in Figure 17, O in Figure 19). Other portions containing strong responses in the reconstructed radargrams are identified as cracked areas. Our results indicate the ability of GPR assisted by the developed signal processing approach to identify shattered rock masses, which require reinforcement strategy to increase stability, and this non-destructive technology can monitoring the surrounding rock conditions during cavern operation. The end effects of EEMD would affect deepest signals, marked as red rectangles in Figure 9, Figure 17 and Figure 19, and therefore the GPR time window should be at least 10% larger than the targeted depth in further applications.

We also want to mention that this research only identifies the fissure interfaces inside the single-type rock formation (granite). Different types of rock formations can determine the rock strength that also influences the surrounding rock stability. Further research is planned in following aspects: Firstly, we will discover the quantity relationship between rock quality and GPR radargrams; secondly, we will investigate the signal processing approach to identify the rock stability with different types of rock formations; thirdly, we will investigate the quantity relationship between rock strengths with different types of rock formations and GPR signal features.