Three-Dimensional Quantitative Recognition of Filler Materials Ahead of a Tunnel Face via Time-Energy Density Analysis of Wavelet Transforms

: Advanced geological prediction of tunnels has become an indispensable task to ensure the safety and effectiveness of tunnel construction before excavation in karst areas. Geological disasters caused by unfavorable geological conditions, such as karst caves, faults, and broken zones ahead of a tunnel face, are highly sudden and destructive. Determining how to predict the spatial location and geometric size of unfavorable geological bodies accurately is a challenging problem. In order to facilitate a three-dimensional quantitative analysis of the ﬁller material ahead of the tunnel face, a biorthogonal wavelet with short support, linear phase, and highly matching waveform of ground penetrating radar (GPR) wavelet is constructed by lifting a simple and general initial ﬁlter on the basis of lifting wavelet theory. A method for a time-energy density analysis of wavelet transforms (TEDAWT) is proposed in accordance with the biorthogonal wavelet. Fifteen longitudinal and horizontal survey lines are used to detect void ﬁllers of different heights. Then, static correction, DC bias, gain, band-pass ﬁltering, and offset processing are performed in the original GPR proﬁle to enhance reﬂected signals and converge diffraction signals. A slice map of GPR proﬁle is generated in accordance with the relative position of longitudinal and horizontal survey lines in space. The wavelet transform analysis of a single-channel signal of each survey line is performed by adopting the TEDAWT method because of the similar rule of the single-channel signal of GPR on the waveform overlay and the ability of the constructed wavelet basis to highlight the time-frequency characteristics of GPR signals. The characteristic value points of the ﬁrst and second interfaces of the void ﬁllers can be clearly determined, and the three-dimensional spatial position and geometric sizes of different void ﬁllers can be obtained. Therefore, the three-dimensional visualization of GPR data is realized. Results show that the TEDAWT method has a good practical application effect in the quantitative identiﬁcation of void ﬁllers, which provides a basis for the interpretation of advanced geological prediction data of tunnels and for the construction decision.


Introduction
Infrastructure monitoring during construction and service has made significant progress in recent decades [1][2][3][4], benefitting from the rapid advances in sensing techniques [5][6][7] and signal-processing algorithms [8][9][10]. As important infrastructure, tunnel construction has attracted much attention in monitoring technology to mitigate geological challenges, such as those from the karst environment, which accounts for 10% of the environment in the world. The prediction of karst characteristics through advanced monitoring is crucial to reduce the risk of tunnel construction [11][12][13][14].
conditions of interpretation. Therefore, professional knowledge in digital signal processing and interpretation techniques is required during the construction of shape, size, and spatial model of the unfavorable geological body ahead of the tunnel face. The current wavelet transform method developed on the basis of Fourier transform has the characteristics of multiresolution, and it is suitable for analyzing and processing nonlinear and nonstationary signals. For example, wavelet-based methods, with great feature extraction capacity, are common data-processing tools in structural health monitoring [45] and damage detection [46] and have many applications, such as concrete crack detection [47], bolt looseness monitoring [48], corrosion detection [49], debonding monitoring [50], and void detection in concrete-filled steel tubes [51]. Thus, wavelet time-frequency localization analysis enables critical components (e.g., signal threshold denoising, feature extraction, spectrum analysis, waveform prediction, and image recognition) toward good signal interpretation capability for GPR [52][53][54]. Nevertheless, when the classical wavelet is used to analyze the time-frequency localization of GPR signals, the choice of wavelet bases is not standardized and targeted, and the results of different wavelet bases may change greatly [55,56]. Hence, determining how to use the latest signal-processing methods to construct a wavelet basis that is suitable for the analysis of the time-frequency characteristics of GPR signals and for processing the GPR signal by employing wavelet theory is difficult.
In this study, we construct a biorthogonal wavelet basis with short support, linear phase, and highly matching waveform of GPR wavelet on the basis of the lifting wavelet theory framework and then add it to the wavelet analysis toolbox. A method for the time-energy density analysis of wavelet transforms (TEDAWT) is presented using the biorthogonal wavelet. The GPR profiles obtained from 15 longitudinal and horizontal survey lines determine the approximate location of void fillers. Subsequently, the singlechannel signal reflecting the characteristics of the void fillers can be analyzed by using the TEDAWT method combined with the good time-frequency localization properties of wavelet transform. The starting and end points of the void fillers are determined on the basis of the mutation or abnormal position on time-energy density curves. Lastly, the spatial positions of the void fillers are generated in accordance with the correlation of spatial positions of survey lines and the calculation results of mutation positions with different single-channel signals. The three-dimensional visualization of the spatial position and geometric size of the void fillers is realized, which provides a reference for the interpretation of advanced geological prediction data of tunnel engineering and is highly advantageous for the safety and efficiency of onsite construction.
The remaining sections of this paper are proposed in the following manner. The theoretical background of the construction of biorthogonal wavelet basis and the TEDAWT method are introduced in Section 2. Section 3 elucidates the measurement methodology adopted for the experiment. The results and discussion addressed in Section 4 are to prove the validity and reliability of the proposed method. Lastly, the conclusions are shown in Section 5.

Construction of Biorthogonal Wavelet Basis
Assuming that an initial biorthogonal filter bank is h, g, h, g , the vanishing moments of the initial filter function at the decomposition and reconstruction ends are N and N, respectively. From the definition of vanishing moment, the following equations exist [57][58][59].
g(z) = (z − 1) N q(z) g(z) = (z − 1) N q(z). (1) After further lifting, new filters h new and g new can be generated. Setting the lifting goal to increase the vanishing moment from N to N, lifting can be regarded as designing a suitable lifting operator s(z). Hence, the new filter g new can be divided by (z − 1) N , as shown as follows.
(z − 1) N q(z) + h(z)s(z 2 ) Given that h(1) = 0, the lifting operator s(z) must contain the factor (z − 1) N . Therefore, the lifting operator s(z) can be defined as follows.
We then obtain the following.
In accordance with the lifting target requirement of vanishing moment, the above equation is satisfied as follows.
Suppose the vanishing moments of the main wavelet and its dual wavelet are both three after lifting. The initial filter is a Haar wavelet with a first-order vanishing moment. The scale and dual-scale filter coefficients of the new biorthogonal wavelet filter bank are obtained by using dual and primal lifting of the initial filter and combining with the particle swarm algorithm for optimization search, as shown as follows.
Meanwhile, the coefficients of wavelet and dual-wavelet filters satisfy the following.
On the basis of the scale filter h 1 (z) and dual-scale filter h 1 (z) obtained above, a new biorthogonal wavelet basis can be constructed.

TEDAWT Method
Assuming a signal function ψ(t) ∈ L 2 (R) (L 2 (R) represents a square integrable space), its continuous Fourier transform is ψ(ω). When ψ(ω) satisfies the wavelet admissibility condition, the following is obtained: where the function ψ(t) is defined as the mother wavelet. If the mother wavelet ψ(t) is stretched and translated, a wavelet basis sequence is obtained. For the real numbers a, b, and a = 0, the wavelet basis sequence is as follows: where a is the stretch factor, and b is the translation factor.
For any signal f (t) ∈ L 2 (R), its continuous wavelet transform is as follows.
In the equation, < f (t), ψ a,b (t) > represents the inner product of f (t) and ψ a,b (t), and . If the biorthogonal wavelet basis constructed above satisfies the admissibility condition of Equation (9), the continuous wavelet transform of the signal is complete, and energy conservation is maintained during the wavelet transform. In accordance with a Moyal inner product theorem, the following equation is given.
Equation (12) shows that the integral of the square amplitude of the wavelet transform coefficients is proportional to the energy of the analyzed signal. When analyzing nonstationary signals, the instantaneous energy density at a certain point in the time-frequency phase space cannot be known due to the Heisenberg uncertainty principle. Therefore, the instantaneous frequency cannot be known at any given moment. However, in Equatio (12), W f (a, b) 2 /C ψ a 2 can be regarded as the energy density function on the time-scale plane, and W f (a, b) 2 ∆a∆b/C ψ a 2 can be regarded as the signal energy at stretch factor a and translation factor b at the center, with stretch interval ∆a and translation interval ∆b. On the basis of the conception of energy density, Equation (12) can be rewritten as follows: where the following is obtained.
In a sense, the stretch factor a in wavelet transform indirectly corresponds to the frequency. Therefore, Equation (14) represents the energy distribution of all frequency bands of the signal over time and is called the time-energy density function.

GPR Methodology
GPR is a nondestructive geophysical detection method without the need for drilling or digging and mainly includes three components, namely, control unit, transmitter, and receiver. The GPR control unit transmits high-frequency electromagnetic waves to the targeted object by using a transmitter. A part of electromagnetic waves becomes reflected when they encounter the targeted object; they then follow the reflection law of electromagnetic waves and transmit to the receiver on the ground due to the differences in relative dielectric constant between the targeted object and the surrounding rock. The control unit images the targeted object with different electromagnetic characteristics, and the recorded images are shown on the display of the control unit. Another part of electromagnetic waves penetrates through the targeted objected and spreads underground. Therefore, the detection principle of GPR is to emit electromagnetic waves into the stratum along the object surface via moving the antenna and then measuring the two-way travel time, strength, and energy of the signal reflected from the targeted object. A schematic of GPR detection and imaging is shown in Figure 1, and a schematic of GPR detection on the tunnel face is shown in Figure 2.

GPR Methodology
GPR is a nondestructive geophysical detection method without the need for drilling or digging and mainly includes three components, namely, control unit, transmitter, and receiver. The GPR control unit transmits high-frequency electromagnetic waves to the targeted object by using a transmitter. A part of electromagnetic waves becomes reflected when they encounter the targeted object; they then follow the reflection law of electromagnetic waves and transmit to the receiver on the ground due to the differences in relative dielectric constant between the targeted object and the surrounding rock. The control unit images the targeted object with different electromagnetic characteristics, and the recorded images are shown on the display of the control unit. Another part of electromagnetic waves penetrates through the targeted objected and spreads underground. Therefore, the detection principle of GPR is to emit electromagnetic waves into the stratum along the object surface via moving the antenna and then measuring the two-way travel time, strength, and energy of the signal reflected from the targeted object. A schematic of GPR detection and imaging is shown in Figure 1, and a schematic of GPR detection on the tunnel face is shown in Figure 2.   GPR produces reflection and refraction on the targeted object by using the characteristics of electromagnetic waves through impedance discontinuity, which conforms to the law of reflection and refraction. The energy generated by electromagnetic waves depends on the reflection coefficient R and refraction coefficient T as follows: where 1 ε and 2 ε are the relative dielectric constants of the medium above and below the interface, respectively. Equations (15) and (16) present that the greater the contrast in GPR produces reflection and refraction on the targeted object by using the characteristics of electromagnetic waves through impedance discontinuity, which conforms to the law of reflection and refraction. The energy generated by electromagnetic waves depends on the reflection coefficient R and refraction coefficient T as follows: where ε 1 and ε 2 are the relative dielectric constants of the medium above and below the interface, respectively. Equations (15) and (16) present that the greater the contrast in the relative dielectric constants between the upper and lower media at the interface, the greater the reflection coefficient of electromagnetic waves, and the greater the reflected energy of electromagnetic waves gain. After receiving the GPR signal via the receiving antenna, the control unit of GPR converts the detection signal into a digital signal. The GPR control unit synthesizes the recorded digital signal to a time profile to determine the depth of the targeted object. Subsequently, in accordance with the recorded two-way travel time and the calculated speed, the depth of the targeted object can be calculated as follows: where v = c/ √ ε r , c is the propagation speed of electromagnetic waves in the air (0.3 m/ns), and ε r is the relative dielectric constant of the medium.
From Equation (17), the depth of the targeted object is given by the following.

Experimental Setup
An indoor model experiment is one of the most intuitive and reliable research methods. In the model experiment of void fillers, the model size and material are primary elements that affect the results of model experiment. The similarity ratio of the model size should be controlled within 15 to 30 where possible. On the basis of the dimension of tunnel caves, namely, small caves (height less than 1.5 m), medium caves (height from 1.5 m to 3 m), and large caves (height more than 3 m), the size of the simulated karst cave is set from 5 cm to 20 cm; four different heights, 5, 10, 15, and 20 cm, are considered; and the length and width are 20 cm. A large amount of yellow mud, wet sand, and other filling media are poured into the tunnel karst caves for real application. In this experiment, empty boxes of different sizes are filled with 25% moisture content of loess mixture, which is used to simulate void fillers. Figure 3 shows the physical drawings of void fillers with different sizes.

Experimental Setup
An indoor model experiment is one of the most intuitive and reliable research methods. In the model experiment of void fillers, the model size and material are primary elements that affect the results of model experiment. The similarity ratio of the model size should be controlled within 15 to 30 where possible. On the basis of the dimension of tunnel caves, namely, small caves (height less than 1.5 m), medium caves (height from 1.5 m to 3 m), and large caves (height more than 3 m), the size of the simulated karst cave is set from 5 cm to 20 cm; four different heights, 5, 10, 15, and 20 cm, are considered; and the length and width are 20 cm. A large amount of yellow mud, wet sand, and other filling media are poured into the tunnel karst caves for real application. In this experiment, empty boxes of different sizes are filled with 25% moisture content of loess mixture, which is used to simulate void fillers. Figure 3 shows the physical drawings of void fillers with different sizes. The location accuracy of the karst caves ahead of the tunnel face recognized by GPR depends on the propagation speed of electromagnetic waves, that is, the dielectric constant of the surrounding rock. Hence, a material with a similar relative dielectric constant as limestone is chosen to serve as the surrounding medium in the tunnel model experiment. The relative dielectric constant of dry sand is from 4 to 6 and close to that of limestone. Dry sand is also selected in the experiment due to its reusability and low cost. To facilitate the molding of the sand box, the sand box is surrounded by concrete. The length, width, and height are set to 4, 3, and 1.5 m, respectively, to minimize the interference of the concrete boundary on the signal collected by GPR.
For easy understanding, the horizontal direction is set to the x-axis, the y-axis is perpendicular to the x-axis, and the z-axis is vertically downward. The four model boxes are buried along the y-axis direction every 0.6 m from low to high, with a depth of 10 cm in the sand box. The GPR antenna passes over the model box and marks each detection situation when gathering data. Three longitudinal survey lines, which are denoted as x 1 , The location accuracy of the karst caves ahead of the tunnel face recognized by GPR depends on the propagation speed of electromagnetic waves, that is, the dielectric constant of the surrounding rock. Hence, a material with a similar relative dielectric constant as limestone is chosen to serve as the surrounding medium in the tunnel model experiment.
The relative dielectric constant of dry sand is from 4 to 6 and close to that of limestone. Dry sand is also selected in the experiment due to its reusability and low cost. To facilitate the molding of the sand box, the sand box is surrounded by concrete. The length, width, and height are set to 4, 3, and 1.5 m, respectively, to minimize the interference of the concrete boundary on the signal collected by GPR.
For easy understanding, the horizontal direction is set to the x-axis, the y-axis is perpendicular to the x-axis, and the z-axis is vertically downward. The four model boxes are buried along the y-axis direction every 0.6 m from low to high, with a depth of 10 cm in the sand box. The GPR antenna passes over the model box and marks each detection situation when gathering data. Three longitudinal survey lines, which are denoted as x 1 , x 2 , and x 3 , are arranged perpendicularly to the x-axis to describe the three-dimensional size information of void fillers; x 2 is arranged along the center of the model box; and x 1 and x 3 are along the model box's edge layout. Twelve horizontal survey lines are arranged perpendicularly to the y-axis, followed by y 1 , y 2 , . . . , y 12 , and the GPR tests of four model boxes in the laboratory are shown in Figure 4.

GPR Data Collection
The detection depth of GPR is generally related to the antenna frequency. The lower the frequency of the GPR receiving-transmitting antenna, the greater the detection depth and the lower the resolution; otherwise, the shallower the detection depth, the higher the resolution. Given that the height of the model box is 5, 10, 15, and 20 cm and from previous research experience, a shielded GPR antenna with 1600 MHz center frequency can be used to collect signals. The Ground Probing Radar Systems produced by Ingegneria Dei Sistemi, Italy is selected in the experiment, and each signal collects 384 data points within a 15 ns sampling time window. Considering that the specific position of void fillers in the sand box is obtained, during the detection process, the signals start to be collected from the sand box boundary at least 0.5 m to reduce the impact of boundary conditions on the GPR profile. In the horizontal sampling process, a signal is collected every 1 cm; 300 sig-

GPR Data Collection
The detection depth of GPR is generally related to the antenna frequency. The lower the frequency of the GPR receiving-transmitting antenna, the greater the detection depth and the lower the resolution; otherwise, the shallower the detection depth, the higher the resolution. Given that the height of the model box is 5, 10, 15, and 20 cm and from previous research experience, a shielded GPR antenna with 1600 MHz center frequency can be used to collect signals. The Ground Probing Radar Systems produced by Ingegneria Dei Sistemi, Italy is selected in the experiment, and each signal collects 384 data points within a 15 ns sampling time window. Considering that the specific position of void fillers in the sand box is obtained, during the detection process, the signals start to be collected from the sand box boundary at least 0.5 m to reduce the impact of boundary conditions on the GPR profile. In the horizontal sampling process, a signal is collected every 1 cm; 300 signals are collected in the longitudinal directions x 1 , x 2 , and x 3 ; and 200 signals are collected in the horizontal directions y 1 , y 2 , . . . , y 12 . Only three longitudinal survey lines and four horizontal survey lines are listed in this paper in consideration of the space limitations. The results are separately shown in Figures 5 and 6. Among them, the red dotted line boxes represent the location of the target body on the GPR time profile. The horizontal direction is the true distance along the GPR survey line in meters; the vertical direction is the two-way travel time of GPR, measured in nanoseconds.   The energy emitted using the GPR transmitting antenna is not limited to vertical propagation into the underground; therefore, it can also be transmitted laterally through the air or the interface between the air and the ground. The energy received by the GPR receiving antenna is mainly interfered by the direct wave between the transmitting and receiving antennas. Thus, the direct wave can be directly used as a zero mark to calculate the depth of object. In addition, the energy received by the GPR receiving antenna comes from the target reflection caused by the impedance discontinuity from shallow to deep. Figures 5 and 6 show that the direct wave reflection of GPR on the surface has strong amplitudes and continuous phases, causing the electromagnetic wave energy to decay rapidly and the deep effective signal to weaken. Hence, the direct wave must be removed to improve the reflection of targeted object via GPR detection.
Each model box occupies a certain width in the direction of the survey line. The GPR antenna collects 200 or 300 single-channel samplings along each survey line from the beginning to the end and then forms a GPR profile. Therefore, the first interface of model boxes in different sizes can be clearly displayed in the GPR profile, and the main feature is a local and discontinuous flat hyperbola. With the increase in the dimension of the model box, the attenuation of electromagnetic waves strengthens, whereas the reflection at the second interface of void fillers weakens. The GPR receiver usually has an ultrawide bandwidth, and the GPR system can easily collect various noises and interferences, resulting in relatively weak reflection signals caused by impedance discontinuities in the stratum. Filtering and migration techniques should be employed to perform signal filtering and enhancement on the original GPR profile to improve the resolution and accuracy of GPR images.

Basic Signal Processing
The original signals shown in Figures 5 and 6 are subjected to static correction, DC bias, gain, band-pass filtering, and offset processing by using the standard analysis software of GPR. The GPR profiles of void fillers are shown in Figures 7 and 8. The energy emitted using the GPR transmitting antenna is not limited to vertical propagation into the underground; therefore, it can also be transmitted laterally through the air or the interface between the air and the ground. The energy received by the GPR receiving antenna is mainly interfered by the direct wave between the transmitting and receiving antennas. Thus, the direct wave can be directly used as a zero mark to calculate the depth of object. In addition, the energy received by the GPR receiving antenna comes from the target reflection caused by the impedance discontinuity from shallow to deep. Figures 5 and 6 show that the direct wave reflection of GPR on the surface has strong amplitudes and continuous phases, causing the electromagnetic wave energy to decay rapidly and the deep effective signal to weaken. Hence, the direct wave must be removed to improve the reflection of targeted object via GPR detection.
Each model box occupies a certain width in the direction of the survey line. The GPR antenna collects 200 or 300 single-channel samplings along each survey line from the beginning to the end and then forms a GPR profile. Therefore, the first interface of model boxes in different sizes can be clearly displayed in the GPR profile, and the main feature is a local and discontinuous flat hyperbola. With the increase in the dimension of the model box, the attenuation of electromagnetic waves strengthens, whereas the reflection at the second interface of void fillers weakens. The GPR receiver usually has an ultrawide bandwidth, and the GPR system can easily collect various noises and interferences, resulting in relatively weak reflection signals caused by impedance discontinuities in the stratum. Filtering and migration techniques should be employed to perform signal filtering and enhancement on the original GPR profile to improve the resolution and accuracy of GPR images.

Basic Signal Processing
The original signals shown in Figures 5 and 6 are subjected to static correction, DC bias, gain, band-pass filtering, and offset processing by using the standard analysis software of GPR. The GPR profiles of void fillers are shown in Figures 7 and 8.   Figures 7 and 8 show that after being processed using the standard software of GPR, the reflection energy of the GPR signal of void fillers is greatly enhanced, the image resolution is significantly improved, and the diffraction waves converge. The reflection of the second interface of void fillers can be determined clearly by comparing with Figures 5 and  6. Similarly, the same operation process can be performed on eight other GPR data collected using the horizontal survey lines, and the standard-processed GPR profiles can also be obtained. The GPR profiles obtained from three longitudinal survey lines and 12 horizontal survey lines are sliced, and the results are shown in Figures 9 and 10.   Figures 7 and 8 show that after being processed using the standard software of GPR, the reflection energy of the GPR signal of void fillers is greatly enhanced, the image resolution is significantly improved, and the diffraction waves converge. The reflection of the second interface of void fillers can be determined clearly by comparing with Figures 5 and 6. Similarly, the same operation process can be performed on eight other GPR data collected using the horizontal survey lines, and the standard-processed GPR profiles can also be obtained. The GPR profiles obtained from three longitudinal survey lines and 12 horizontal survey lines are sliced, and the results are shown in Figures 9 and 10.  Figures 7 and 8 show that after being processed using the standard software of GPR, the reflection energy of the GPR signal of void fillers is greatly enhanced, the image resolution is significantly improved, and the diffraction waves converge. The reflection of the second interface of void fillers can be determined clearly by comparing with Figures 5 and  6. Similarly, the same operation process can be performed on eight other GPR data collected using the horizontal survey lines, and the standard-processed GPR profiles can also be obtained. The GPR profiles obtained from three longitudinal survey lines and 12 horizontal survey lines are sliced, and the results are shown in Figures 9 and 10.

Time-Energy Density Analysis
In consideration of the space limitations in this paper, the longitudinal x2 line and the horizontal y 2 , y 5 , y 8

Time-Energy Density Analysis
In consideration of the space limitations in this paper, the longitudinal x 2 line and the horizontal y 2 , y 5 , y 8 , and y 11 lines are used as examples for explanation. The typical single-channel signals reflecting the characteristics of void fillers are extracted from the GPR profile on the longitudinal survey line x 2 and the horizontal survey lines y 2 , y 5 , y 8 , and y 11 . The results are shown in Figures 11 and 12.

Time-Energy Density Analysis
In consideration of the space limitations in this paper, the longitudinal x2 line and the horizontal y 2 , y 5 , y 8 Figures 11 and 12 show that the longitudinal and horizontal single-channel signals have a similar attenuation rule on the time axis. The GPR reflection signals of the upper and lower interfaces of the large void filler are clearly separated, while the GPR reflection signals of the upper and lower interfaces of the small void filler are superimposed on each other. That is, as the size of the void filler increases, the GPR reflection signals of the upper and lower interfaces of the void fillers are more clearly separated. Therefore, a single-channel signal can be used as an example for analysis. The biorthogonal wavelet basis constructed above is used as the wavelet basis for wavelet transform. The timeenergy density analysis program for wavelet transform is written in MATLAB, and then the different single-channel signals are analyzed using the TEDAWT method. The results are shown in Figures 13 and 14. Figure 13 shows that the two prominent extreme points on the TEDAWT curves of single-channel GPR signals on the longitudinal survey line x 2 are two characteristics at the first and second interfaces of void fillers, that is, characteristic points t 1 and t 2 . The time differences between the two characteristic points are 0.585, 1.375, 2.107, and 3.015 ns. Given that the reflection time of GPR detection is a two-way travel time and in accordance with Equation (18), the height values of void fillers are calculated as 4.01, 9.41, 14.43, and 20.64 cm. The relative errors of using this method to identify void fillers are −19.90%, −5.86%, −3.83%, and 3.21%. From the analysis results of Figure 14, the heights of different void fillers are 3.81, 9.02, 14.03, and 20.64, that is, the relative errors are −23.87%, −9.83%, −6.48%, and 3.21%, respectively.   Figure 14, the heights of different void fillers are 3.81, 9.02, 14.03, and 20.64, that is, the relative errors are −23.87%, −9.83%, −6.48%, and 3.21%, respectively.

Comparative Analysis
Extensive analysis and research show that the characteristics of wavelet bases of the Daubechies series are good, compact, and supportive and show regularity and smoothness; thus, they are widely used in GPR signal analysis and processing. The following uses db2, db4, db6, and db8 wavelet bases to perform wavelet transform on the same signal y 2 to obtain the modulus curves of wavelet transform coefficients. The results are shown in Figure 15. It can be observed from Figure 15 that, by using different wavelet bases to analyze the same GPR signal, the modulus curves of wavelet transform coefficients obtained are not exactly the same. The time of the first local maximum point is 0.819, 0.907, 0.819, and 0.761, respectively. The main reason is that different wavelet bases have different time-frequency localization characteristics. By comparing and analyzing the modulus curves of Figures 15 and 14a, it was observed that the resolution of the modulus curve obtained by using the Daubechies series wavelet bases is poor, and it is difficult to accurately obtain the two modulus maximum points of a GPR signal.

Comparative Analysis
Extensive analysis and research show that the characteristics of wavelet bases of the Daubechies series are good, compact, and supportive and show regularity and smoothness; thus, they are widely used in GPR signal analysis and processing. The following uses db2, db4, db6, and db8 wavelet bases to perform wavelet transform on the same signal y 2 to obtain the modulus curves of wavelet transform coefficients. The results are shown in Figure 15. It can be observed from Figure 15 that, by using different wavelet bases to analyze the same GPR signal, the modulus curves of wavelet transform coefficients obtained are not exactly the same. The time of the first local maximum point is 0.819, 0.907, 0.819, and 0.761, respectively. The main reason is that different wavelet bases have different time-frequency localization characteristics. By comparing and analyzing the modulus curves of Figures 14a and 15, it was observed that the resolution of the modulus curve obtained by using the Daubechies series wavelet bases is poor, and it is difficult to accurately obtain the two modulus maximum points of a GPR signal.
The db5 wavelet basis is employed to perform wavelet transform on the GPR signal on the survey line y 2 , y 5 , y 8 , and y 11 , and the modulus curves of wavelet transform coefficient are obtained. The results are shown in Figure 16. Comparing the modulus curves in Figure 16 with the modulus curves in Figure 14, the resolution of the modulus curves obtained by using the db5 wavelet basis is relatively poor; thus, it is difficult to distinguish the first local point in the modulus curves with the second local point (such as Figure 16c,d), which is not beneficial to automatically identify feature points in the modulus curve.

Comparative Analysis
Extensive analysis and research show that the characteristics of wavelet bases of the Daubechies series are good, compact, and supportive and show regularity and smoothness; thus, they are widely used in GPR signal analysis and processing. The following uses db2, db4, db6, and db8 wavelet bases to perform wavelet transform on the same signal y 2 to obtain the modulus curves of wavelet transform coefficients. The results are shown in Figure 15. It can be observed from Figure 15 that, by using different wavelet bases to analyze the same GPR signal, the modulus curves of wavelet transform coefficients obtained are not exactly the same. The time of the first local maximum point is 0.819, 0.907, 0.819, and 0.761, respectively. The main reason is that different wavelet bases have different time-frequency localization characteristics. By comparing and analyzing the modulus curves of Figures 15 and 14a, it was observed that the resolution of the modulus curve obtained by using the Daubechies series wavelet bases is poor, and it is difficult to accurately obtain the two modulus maximum points of a GPR signal.  Figure 16. Comparing the modulus curves in Figure 16 with the modulus curves in Figure 14, the resolution of the modulus curves obtained by using the db5 wavelet basis is relatively poor; thus, it is difficult to distinguish the first local point in the modulus curves with the second local point (such as Figure 16c,d), which is not beneficial to automatically identify feature points in the modulus curve.  The db5 wavelet basis is employed to perform wavelet transform on the GPR signal on the survey line y 2 , y 5 , y 8 , and y 11 , and the modulus curves of wavelet transform coefficient are obtained. The results are shown in Figure 16. Comparing the modulus curves in Figure 16 with the modulus curves in Figure 14, the resolution of the modulus curves obtained by using the db5 wavelet basis is relatively poor; thus, it is difficult to distinguish the first local point in the modulus curves with the second local point (such as Figure 16c,d), which is not beneficial to automatically identify feature points in the modulus curve.  When the wavelet transform modulus maximum method is used to analyze the singularity characteristics of a GPR signal, the stretch factor has a great influence on the result of wavelet transform. When the stretch factors are 2, 4, 6, and 8, the db5 wavelet basis is used to perform wavelet transform on the GPR signal on the survey line y 2 , and the modulus curves of wavelet transform coefficients are obtained. The results are shown in Figure 17. It can be observed from Figure 17 that under different stretch scales, the obtained wavelet transform coefficient modulus curves are not completely the same, and the first local maximum points obtained are 0.937, 0.995, 1.024, and 0.878, respectively. Therefore, when performing wavelet transform on the GPR signal, the selected stretch factor is different, and the singularity of GPR signal obtained is also different. Finally, compared the modulus curves of Figure 17 with the modulus curves of Figure 14a, the resolution of the curve obtained by the TEDAWT method proposed is higher than the wavelet transform modulus maximum method. 2 ulus curves of wavelet transform coefficients are obtained. The results are shown in Figure  17. It can be observed from Figure 17 that under different stretch scales, the obtained wavelet transform coefficient modulus curves are not completely the same, and the first local maximum points obtained are 0.937, 0.995, 1.024, and 0.878, respectively. Therefore, when performing wavelet transform on the GPR signal, the selected stretch factor is different, and the singularity of GPR signal obtained is also different. Finally, compared the modulus curves of Figure 17 with the modulus curves of Figure 14a, the resolution of the curve obtained by the TEDAWT method proposed is higher than the wavelet transform modulus maximum method.

Quantitative Recognition
The height values of different void fillers obtained from 480 single-channel signals on each GPR survey line are determined to obtain detailed statistical laws, and the recognition results are plotted in Figure 18.

Quantitative Recognition
The height values of different void fillers obtained from 480 single-channel signals on each GPR survey line are determined to obtain detailed statistical laws, and the recognition results are plotted in Figure 18.    20, 9.42, 14.17, and 20.79 and variances of 0.11, 0.24, 0.18, and 0.21, respectively. The identification errors of void fillers are −16.09%, −5.79%, −5.54%, and 3.94%. With a continuous increase in the size of the void filler, the identification error is smaller, thus indicating the accuracy of the method proposed in this paper is higher. The extremum points obtained from the first and second interfaces of different void fillers are connected, and the results are plotted in Figure 19.   20, 9.42, 14.17, and 20.79 and vari ances of 0.11, 0.24, 0.18, and 0.21, respectively. The identification errors of void fillers are −16.09%, −5.79%, −5.54%, and 3.94%. With a continuous increase in the size of the void filler, the identification error is smaller, thus indicating the accuracy of the method pro posed in this paper is higher. The extremum points obtained from the first and second interfaces of different void fillers are connected, and the results are plotted in Figure 19.  Figure 19 shows that the positions and geometric sizes of void fillers can be obtained on the basis of the relative position of each survey line passing through the void filler and the height value between the first and second interfaces. Subsequently, connecting the characteristic value points of each filling with a curved surface forms a three-dimen sional curved surface diagram. The results are shown in Figure 20, from which the specifi positions and geometric sizes of each filling are obtained. Therefore, the TEDAWT method Depth (m) Figure 19. Three-dimensional scatter diagrams formed from the characteristic values of void fillers. Figure 19 shows that the positions and geometric sizes of void fillers can be obtained on the basis of the relative position of each survey line passing through the void fillers and the height value between the first and second interfaces. Subsequently, connecting the characteristic value points of each filling with a curved surface forms a three-dimensional curved surface diagram. The results are shown in Figure 20, from which the specific positions and geometric sizes of each filling are obtained. Therefore, the TEDAWT method can successfully identify the positions of the first and second interfaces of void fillers at the singular point of GPR signals and then accurately reveal the spatial position and three-dimensional geometry of void fillers in the sand box. It provides a basis for the quantitative identification of unfavorable geological bodies in the advanced geological prediction of tunnels, which is of great significance to the safe and efficient operation of engineering projects.
can successfully identify the positions of the first and second interfaces of void fillers at the singular point of GPR signals and then accurately reveal the spatial position and threedimensional geometry of void fillers in the sand box. It provides a basis for the quantitative identification of unfavorable geological bodies in the advanced geological prediction of tunnels, which is of great significance to the safe and efficient operation of engineering projects.

Conclusions
In this study, a biorthogonal wavelet basis that is short support, linear phase, and highly matching waveform of GPR wavelet was constructed and added to the wavelet analysis toolbox by using the lifting wavelet theory framework. Then, the TEDAWT method was proposed and applied to analyze the mutation or abnormal position of singlechannel signals reflecting the characteristics of void fillers. The experimental results demonstrated that the proposed method can effectively determine the spatial position and geometric size of void fillers. The main contributions of this paper are summarized as follows: (1) Static correction, DC bias, gain, band-pass filtering, and offset processing are performed on the original GPR profile on the basis of the principle of GPR detection to converge the reflected waves of the GPR signal and return the diffracted waves to their positions. The reflections of the second interface can be displayed clearly in the GPR profile. The approximate location and rough shape of void fillers on each survey line can be determined visually by forming the GPR profiles obtained from the longitudinal and horizontal survey lines. (2) In accordance with a similar rule of single-channel signal of GPR on the waveform overlay, the single-channel signal in the void fillers can be analyzed by using the TEDAWT method, suppressing the interference of random noise. The curve resolution is high, and the characteristic value points of the first and second interfaces of void fillers are determined clearly. In accordance with the characteristic value points obtained using longitudinal and horizontal survey lines, the height values and space sizes of different void fillers are formed, and a three-dimensional visualization of void fillers is realized. As the sizes of void fillers increase, the accuracy of identification by using the TEDAWT method increases. This condition means that the

Conclusions
In this study, a biorthogonal wavelet basis that is short support, linear phase, and highly matching waveform of GPR wavelet was constructed and added to the wavelet analysis toolbox by using the lifting wavelet theory framework. Then, the TEDAWT method was proposed and applied to analyze the mutation or abnormal position of single-channel signals reflecting the characteristics of void fillers. The experimental results demonstrated that the proposed method can effectively determine the spatial position and geometric size of void fillers. The main contributions of this paper are summarized as follows: (1) Static correction, DC bias, gain, band-pass filtering, and offset processing are performed on the original GPR profile on the basis of the principle of GPR detection to converge the reflected waves of the GPR signal and return the diffracted waves to their positions. The reflections of the second interface can be displayed clearly in the GPR profile. The approximate location and rough shape of void fillers on each survey line can be determined visually by forming the GPR profiles obtained from the longitudinal and horizontal survey lines. (2) In accordance with a similar rule of single-channel signal of GPR on the waveform overlay, the single-channel signal in the void fillers can be analyzed by using the TEDAWT method, suppressing the interference of random noise. The curve resolution is high, and the characteristic value points of the first and second interfaces of void fillers are determined clearly. In accordance with the characteristic value points obtained using longitudinal and horizontal survey lines, the height values and space sizes of different void fillers are formed, and a three-dimensional visualization of void fillers is realized. As the sizes of void fillers increase, the accuracy of identification by using the TEDAWT method increases. This condition means that the TEDAWT method has good superiority in the detection and analysis of GPR signal singularity. (3) Forecasting the unfavorable geological conditions, such as faults, karst caves, and weak broken zones ahead of the tunnel face, is helpful to prepare in advance during tunnel excavation. Geological disasters, such as water and mud inrush and collapse caused by blind excavation, can also be avoided. Compared with the two-dimensional profile analysis method of GPR, the three-dimensional visualization method can make a more intuitive judgment on the scale, geometric size, and distribution shape