Experimental Research on Evaluation of Soil Water Content Using Ground Penetrating Radar and Wavelet Packet-Based Energy Analysis

: Soil water content is one of the most important factors affecting the safety and stability of buildings or structures, especially in roadbeds, slopes, earth dams and foundations. Accurate assessments of soil water content can ensure the quality of construction, reduce construction costs and prevent accidents, among other beneﬁts. In this study, ground penetrating radar (GPR) was used to detect and evaluate changes in soil water content. The GPR signal is usually nonstationary and nonlinear; however, traditional Fourier theory is typically suitable for periodic stationary signals, and cannot reﬂect the law of the frequency and energy of the GPR signal changing with time. Wavelet transform has good time-frequency localization characteristics, and therefore represents a new method for analyzing and processing GPR signals. According to the time-frequency characteristics of GPR signals, in this paper, a new biorthogonal wavelet basis which was highly matched with the GPR waveform was constructed using the lifting framework of wavelet theory. Subsequently, an evaluation method, namely, the wavelet packet-based energy analysis (WPEA) method, was proposed. The method was utilized to calculate the wavelet packet-based energy indexes (WPEI) of the GPR single-channel signals for clay samples with water contents ranging from 10% to 24%. The research results showed that there was a highly correlated linear relationship between the WPEI and the soil water contents, and the relationship between the two was ﬁtted with a linear ﬁtting function. The feasibility of the method was veriﬁed by comparing our results with those obtained using classical wavelet bases to perform the wavelet packet transform. The large-area, continuous scanning measurement method of GPR was shown to be suitable for evaluations of soil water contents in roadbeds, slopes, earth dams, and foundations.


Introduction
Soil is an essential element for human survival, production and life. Detailed knowledge of the physical and mechanical parameters of soil is critical for scientific research in fields such as geophysics, ecology, agricultural engineering, water conservation, and civil engineering [1,2]. In the civil engineering, the physical properties (i.e., density, saturation, porosity, etc.) and mechanical properties (i.e., stress, compression coefficient, compression modulus, swelling, shrinkage, etc.) of soil are deeply affected by water content (i.e., moisture content) [3,4]. Changes in the soil water content will not only cause huge changes in soil strength, deformation, frost heave, and infiltration, but will also have a significant impact on the stability of the foundations of buildings and structures, the subsidence of In recent years, with the rapid development of digital signal processing technology, wavelet transform has been widely used in image compression, edge detection, speech synthesis, digital watermarking, signal filtering, geophysical prospecting, medical detection, and engineering calculations [44]. The procedure has good time-frequency characteristics, and as such, it opens up a new method for GPR signals analysis and processing. Wavelet analysis technology has been used in a large amount of research on GPR signal threshold denoising, attribute analysis, feature extraction and waveform predictions [45][46][47]. However, when performing wavelet transform on a GPR signal, the question of which wavelet basis is used completely depends on the experience of the relevant personnel [20,48]. Moreover, the selected wavelet bases do not have good similarity with the GPR signals being analyzed, and the inherent properties of the GPR signals cannot be accurately extracted, which hinders the in-depth study of wavelet transform in the analysis and processing of GPR signals [49,50]. Therefore, the question of how to use modern signal processing methods to construct a wavelet basis that matches the characteristics of the GPR signals must be urgently addressed if wavelet transform is to be applied to GPR signal processing.
The aim of this paper was to construct a symmetrical and tightly supported biorthogonal wavelet basis that was highly matched to the GPR waveform. Based upon this, a wavelet packet-based energy analysis (WPEA) method was proposed. The feasibility of the method was verified by setting up soils with different water contents in a model box in order to carry out simulation experiments. Then, the WPEA method was used to evaluate the soil water content in the model experiment. In the laboratory, eight kinds of clay samples with different water contents were used, and the GPR method was applied to carry out four repetitive experiments on the samples. The WPEA method was used to calculate the wavelet packet-based energy indexes (WPEI) of the GPR signals under different soil water contents, and the relationship curves between the WPEI and the soil water content were obtained. The effectiveness of the method proposed in this paper was verified by comparing our results with those obtained using classical wavelet bases for wavelet packet transform.
The rest of this paper is organized as follows. Section 2 describes the construction of a biorthogonal wavelet basis, wavelet packet decomposition, the WPEA method, and the GPR detection principle. The forward simulation experiment verifies the feasibility of using GPR to test the soil water content. Section 3 describes the GPR detection procedure and the data collection process in the case of soils with different water contents in the indoor experiment. The response law of the GPR signals and the WPEI of soils with different water contents were analyzed and obtained, and the effectiveness of the wavelet basis proposed in this paper is verified in Section 4. Finally, conclusions are summarized in Section 5.

Construction of Biorthogonal Wavelet Basis
In discrete wavelet or stationary wavelet transform, the Mallat algorithm is mainly used in wavelet analysis. Usually, the signal is decomposed and reconstructed by a bank composed of four filters h, g, h, g , and the signal is divided into low-frequency approximation and high-frequency detail information. The filters g and h respectively represent the dual wavelet filter and dual scale filter at the decomposition end, also known as the decomposition high pass filter and decomposition low pass filter. Filters g and h denote the wavelet filter and scale filter at the reconstruction end, also known as the reconstruction high pass filter and the reconstruction low pass filter. Suppose the vanishing moments of the filter bank at the decomposition and reconstruction ends are N and N, the wavelet filter function has the following relationship based on the definition of vanishing moments [50,51]: (2), there must be a factor ( 1) N z − in the lifting operator ( ) s z . Therefore, the lifting operator ( ) s z can be defined as [52] ( where U z ( ) is Laurent polynomial to be solved. Substituting Equation (3) into Equation (2) gives: − , N′ is the vanishing moment of the wavelet filter after primal lifting, N is the vanishing moment of the initial wavelet filter.
where P z ( ) is the Laurent polynomial about z. According to the lifting target requirement of vanishing moment of the new wavelet function, the Equation (5) satisfies the following relational expression [52,54]: Assuming that the lifting goal is to increase the vanishing moment of the wavelet filter g from N to N , then the lifting process can be regarded as designing a suitable lifting operator s(z) so that the wavelet filter g new after being lifted can be divisible by (z−1) N [52,53] ( where Q(z) is the Laurent polynomial about z.
Since h(1) = 0 in Equation (2), there must be a factor (z − 1) N in the lifting operator s(z). Therefore, the lifting operator s(z) can be defined as [52] where U(z) is Laurent polynomial to be solved. Substituting Equation (3) into Equation (2) gives: where ∆N = N − N, N is the vanishing moment of the wavelet filter after primal lifting, N is the vanishing moment of the initial wavelet filter. Given where P(z) is the Laurent polynomial about z. According to the lifting target requirement of vanishing moment of the new wavelet function, the Equation (5) satisfies the following relational expression [52,54]: The expressions of Haar wavelet filters at the decomposition and reconstruction end are h(z) = 1 2 (1 + z −1 ) and h(z) = 1 2 (1 + z −1 ), and the expressions of wavelet filters at the decomposition and reconstruction end are g(z) = 1 2 (z−1) · z −1 and g(z) = 1 2 (z−1) · z −1 . Assume that the vanishing moments of the wavelet filter function after lifting are 3 at the decomposition and reconstruction ends, and the initial filter is a Haar wavelet filter Remote Sens. 2021, 13, 5047 5 of 21 bank with 1st-order vanishing moments. By performing dual lifting and original lifting operations on the initial filter h, g, h, g , combined with the particle swarm algorithm for optimization search, the scale filter and dual scale filter coefficients of the new biorthogonal wavelet filter bank are finally obtained, namely At the same time, the coefficients of the wavelet filter and the dual wavelet filter are as follows: After calculating the scale filter coefficients ( h 1 and h 1 ) at the decomposition end and reconstruction end, the scale function ( ϕ new and ϕ new ) and wavelet function ( ψ new and ψ new ) forms corresponding to this wavelet can be obtained through an iterative algorithm based on the two-scale equation (see Formula (9)).
The function results ϕ new , ψ new , ϕ new , ψ new corresponding to the biorthogonal wavelet basis can be acquired by adding the scale filter coefficients ( h 1 and h 1 ) into the MATLAB wavelet toolbox.

Decomposition of Wavelet Packet Transform
Recently, various signal processing methods have been widely used for GPR signal analysis and processing. Fourier transform, fractal analysis, wavelet/wavelet packet transform and Hilbert-Huang transform are the most popular methods [49,[55][56][57]. In these methods, wavelet packet analysis can not only decompose the low-frequency part of the GPR signal, but can also further decompose the high-frequency part that is not refined by the wavelet analysis [58]. In addition, wavelet packet analysis can perform adaptive frequency band decomposition according to the inherent frequency characteristics of the GPR signal, which greatly improves the time-frequency resolution. Therefore, wavelet packet analysis has broad range of application prospects. In the decomposition process of the wavelet packet, the original GPR signals can be decomposed into several frequency band subsignals in a relatively short time window, and each subsignal includes an approximate part ("A", i.e., the low-frequency subsignal) and a detailed part ("D", i.e., the high-frequency subsignal). Figure 2 shows the decomposition process of a three-level wavelet packet for GPR signals. Component A1 is the low frequency of the first layer, and component D1 is the high-frequency of the first layer. Component AA2 is the low-frequency of the low-frequency part of the second-layer, and component DA2 is the high-frequency of the low-frequency part of the second-layer; component AD2 is the low-frequency of the high-frequency part of the second-layer, and component DD2 is the high-frequency of the high-frequency part of the second-layer. The other components can be deduced according to this logic. D", i.e., the high-frequency subsignal). Figure 2 shows the decomposition process of a three-level wavelet packet for GPR signals. Component A1 is the low frequency of the first layer, and component D1 is the high-frequency of the first layer. Component AA2 is the low-frequency of the low-frequency part of the second-layer, and component DA2 is the high-frequency of the low-frequency part of the second-layer; component AD2 is the low-frequency of the high-frequency part of the second-layer, and component DD2 is the high-frequency of the high-frequency part of the second-layer. The other components can be deduced according to this logic. If the original GPR single-channel signal ( ) X t is decomposed, a total of 2 i wavelet packet decomposition coefficients of the i -th layer are obtained, and the mathematical expression is given by [59,60]:

The WPEA Method
The GPR single-channel signal is subjected to wavelet packet transform. According to the sampling theorem, if the sampling frequency of the GPR single-channel signal is s f , where f is the abbreviation of frequency, and s is the abbreviation of sampling, then the Nyquist frequency (i.e., minimum sampling frequency) of the signal is 2 / s f [61,62]. The GPR single-channel signal is decomposed to the i -th layer, and the corre- If the original GPR single-channel signal X(t) is decomposed, a total of 2 i wavelet packet decomposition coefficients of the i-th layer are obtained, and the mathematical expression is given by [59,60]: where j = 1, 2, · · ·, 2 i − 1, 2 i .

The WPEA Method
The GPR single-channel signal is subjected to wavelet packet transform. According to the sampling theorem, if the sampling frequency of the GPR single-channel signal is f s , where f is the abbreviation of frequency, and s is the abbreviation of sampling, then the Nyquist frequency (i.e., minimum sampling frequency) of the signal is f s /2 [61,62]. The GPR single-channel signal is decomposed to the i-th layer, and the corresponding lowest frequency band is 0~f s /2 i+1 . Assuming that the energy corresponding to the reconstructed subsignal X i,j (t) in each frequency band is E i,j (t), the following formula may be applied [63] where x j,k (t)(k = 1, 2, · · ·, m) (m is the number of sampling points) represents the discrete point amplitude of the reconstructed subsignal.
Assuming that the total energy of a GPR single-channel signal is E total (t), the following relationship exists: The energy can reflect the propagation law of the GPR signals, but can also reflect the inherent properties of the medium [64,65]. In this study, the WPEA method was adapted to calculate the WPEI of the GPR single-channel signals under different soil water contents.

Detection Principle
The GPR is a geophysical detection method that uses media with different dielectric constants to produce different responses to electromagnetic wave reflection signals in order to determine the characteristics of underground media. It primarily includes a control unit, a transmitter and a receiver. The control unit provides instructions to the transmitter, which, in turn, emits high-frequency electromagnetic waves to the underground medium. During the propagation of electromagnetic waves into underground media, due to the difference in the relative dielectric constant of different media, some of the electromagnetic waves are reflected when they encounter the interface. The reflected waves follow the law of reflection, and are received by the antenna. The other part of the electromagnetic wave will pass through the interface and propagate to the deep part of the ground, where it is absorbed by the underground medium. As the GPR antenna moves, the control unit records the electromagnetic wave reflection signals at each sampling point. By analyzing the characteristics of the amplitude, phase, and frequency of the reflected signal, the stratum distribution and properties of the underground medium can be determined. A schematic diagram of GPR detection imaging is presented in Figure 3.

Forward Simulation
The finite-difference time-domain method can be used to realize the forward simulation process of GPR [66,67]. Assuming that the stratum is a uniform and continuous medium in a semi-infinite space, the reflection and refraction of electromagnetic waves used for forward simulation occur on a two-dimensional plane. To study the reflection law of electromagnetic waves in the stratum with different water contents, the stratum model shown in Figure 4a was designed. The simulation area of the GPR was a 1 m × 1 m square area; the upper left corner is the origin of the coordinates, the abscissa is the detection distance, and the ordinate is the detection depth. The detection target was a rectangular cavity with a length of 20 cm, a height of 10 cm, and a buried depth of 10 cm. Since dry sand could be reused, it was convenient to bury various model boxes. The background medium of the stratum was dry sand, the relative dielectric constant was 3, and the relative dielectric constant of air was set at 1.
GprMax is open source software that uses the finite-difference time-domain method to simulate electromagnetic wave propagation; as such, it can be used for GPR forward simulations [66,67]. The central-frequency of the antenna was 1600 MHz in the forward simulation of the GPR. The boundary absorption condition was a perfectly matched layer, and the excitation source employed was the Ricker wavelet. The spatial grid step size and the sampling step size were both 2.5 mm. The number of sampling channels was 365 and the sampling time was 30 ns. The GPR forward simulation software was adapted to perform forward simulations on the cavity model shown in Figure 4a. The obtained results are shown in Figure 4b.

Forward Simulation
The finite-difference time-domain method can be used to realize the forward simulation process of GPR [66,67]. Assuming that the stratum is a uniform and continuous medium in a semi-infinite space, the reflection and refraction of electromagnetic waves used for forward simulation occur on a two-dimensional plane. To study the reflection law of electromagnetic waves in the stratum with different water contents, the stratum model shown in Figure 4a was designed. The simulation area of the GPR was a 1 m × 1 m square area; the upper left corner is the origin of the coordinates, the abscissa is the detection distance, and the ordinate is the detection depth. The detection target was a rectangular cavity with a length of 20 cm, a height of 10 cm, and a buried depth of 10 cm. Since dry sand could be reused, it was convenient to bury various model boxes. The background medium of the stratum was dry sand, the relative dielectric constant was 3, and the relative dielectric constant of air was set at 1.  The influence of the change of soil water content on the electromagnetic w propagation characteristics is mainly reflected in differences in the relative dielec constant and the electrical conductivity of soil. To study the effect of varying these rameters in different soils on the electromagnetic wave reflection signal, the cavity shown in Figure 4a, was filled with soils with different water contents.
According to the formulas of Topp and Liu, the relative dielectric constant a electrical conductivity of soil with different water contents were calculated [68,69]. T two parameters were selected as variables to distinguish different media in the forw simulation of GPR. The relationships between the relative dielectric constant, electr conductivity and water content of the medium are as follows [68,69]: .
where w is soil water content (%), ε is relative dielectric constant, σ is electr conductivity (S/cm), and v is electromagnetic wave velocity (m/ns). Table 1 presents the physical properties of soil with increasing moisture content. T GPR forward simulation software was used to perform forward simulations on models with different water contents. Based upon 5087 iterative calculations, the res of the GPR forward simulation of soils with different water contents are shown in   GprMax is open source software that uses the finite-difference time-domain method to simulate electromagnetic wave propagation; as such, it can be used for GPR forward simulations [66,67]. The central-frequency of the antenna was 1600 MHz in the forward simulation of the GPR. The boundary absorption condition was a perfectly matched layer, and the excitation source employed was the Ricker wavelet. The spatial grid step size and the sampling step size were both 2.5 mm. The number of sampling channels was 365 and the sampling time was 30 ns. The GPR forward simulation software was adapted to perform forward simulations on the cavity model shown in Figure 4a. The obtained results are shown in Figure 4b.
The influence of the change of soil water content on the electromagnetic wave propagation characteristics is mainly reflected in differences in the relative dielectric constant and the electrical conductivity of soil. To study the effect of varying these parameters in different soils on the electromagnetic wave reflection signal, the cavity, as shown in Figure 4a, was filled with soils with different water contents.
According to the formulas of Topp and Liu, the relative dielectric constant and electrical conductivity of soil with different water contents were calculated [68,69]. The two parameters were selected as variables to distinguish different media in the forward simulation of GPR. The relationships between the relative dielectric constant, electrical conductivity and water content of the medium are as follows [68,69]: where w is soil water content (%), ε is relative dielectric constant, σ is electrical conductivity (S/cm), and v is electromagnetic wave velocity (m/ns). Table 1 presents the physical properties of soil with increasing moisture content. The GPR forward simulation software was used to perform forward simulations on soil models with different water contents. Based upon 5087 iterative calculations, the results of the GPR forward simulation of soils with different water contents are shown in Figure 5.  The radagram presented in Figure 5 indicates that the greater the difference in physical properties between dry sand and cavity-filled soil, the stronger the reflectio electromagnetic waves at the interface of dry sand and cavity-filled soil. From the G time profile shown in Figure 5a-h, the middle single-channel signals were extracted. The radagram presented in Figure 5 indicates that the greater the difference in the physical properties between dry sand and cavity-filled soil, the stronger the reflection of electromagnetic waves at the interface of dry sand and cavity-filled soil. From the GPR time profile shown in Figure 5a-h, the middle single-channel signals were extracted. The WPEA method proposed in Section 2.3 was employed to calculate the four single-channel signals, and the WPEI of the GPR signals in the case of soils with different water contents were obtained. The results are presented in Figure 6. WPEA method proposed in Section 2.3 was employed to calculate the four single-channel signals, and the WPEI of the GPR signals in the case of soils with different water contents were obtained. The results are presented in Figure 6. It can be seen from Figure 6 that when the differences in physical properties, i.e., relative dielectric constant and electrical conductivity between the cavity-filled soil and the dry sand, are small, the electromagnetic wave reflects weakly at the interface of the two media, the signal amplitude and energy are low, and the WPEI is small. As the water content increases, the differences in the physical properties between the soil and the dry sand increase; additionally, the amplitude and energy of the GPR reflected signal gradually increase, as does the WPEI of the obtained GPR signals. When the soil water content exceeds 18%, the GPR reflection at the interface between the soil and the dry sand is strong, and the WPEI of the obtained GPR signals is large, but the range of increase is small.
To determine the influence of the relative dielectric constant and electrical conductivity of the soil on the WPEI of the GPR signals, the following two working conditions were established. In the first case, when the soil water content was 18%, the corresponding electrical conductivity was 368.670 μS/cm; this value remained unchanged. When the soil water content changed from 10% to 24%, the corresponding relative dielectric constants were 5.343, 6.115, 6.983, 7.941, 8.987, 10.116, 11.325, 12.611, respectively. In the second case, when the soil water content was 18%, the corresponding relative dielectric constant was 8.897; this value also remained unchanged. When the soil water content changed from 10% to 24%, the corresponding electrical conductivities were 88.681 μS/cm, 126.085 μS/cm, 180.261 μS/cm, 258.134 μS/cm, 368.670 μS/cm, 523.129 μS/cm, 735.272 μS/cm, 1021.533 μS/cm, respectively.
A GPR forward simulation was performed on the geoelectric models in the two cases, and then the GPR single-channel signals were analyzed by the WPEA method. The analysis results are shown in Figure 7.
(a) (b) Figure 6. WPEI of cavity-filled soil with increasing water contents.
It can be seen from Figure 6 that when the differences in physical properties, i.e., relative dielectric constant and electrical conductivity between the cavity-filled soil and the dry sand, are small, the electromagnetic wave reflects weakly at the interface of the two media, the signal amplitude and energy are low, and the WPEI is small. As the water content increases, the differences in the physical properties between the soil and the dry sand increase; additionally, the amplitude and energy of the GPR reflected signal gradually increase, as does the WPEI of the obtained GPR signals. When the soil water content exceeds 18%, the GPR reflection at the interface between the soil and the dry sand is strong, and the WPEI of the obtained GPR signals is large, but the range of increase is small.
To determine the influence of the relative dielectric constant and electrical conductivity of the soil on the WPEI of the GPR signals, the following two working conditions were established. In the first case, when the soil water content was 18%, the corresponding electrical conductivity was 368.670 µS/cm; this value remained unchanged. When the soil water content changed from 10% to 24%, the corresponding relative dielectric constants were 5.343, 6.115, 6.983, 7.941, 8.987, 10.116, 11.325, 12.611, respectively. In the second case, when the soil water content was 18%, the corresponding relative dielectric constant was 8.897; this value also remained unchanged. When the soil water content changed from 10% to 24%, the corresponding electrical conductivities were 88.681 µS/cm, 126.085 µS/cm, 180.261 µS/cm, 258.134 µS/cm, 368.670 µS/cm, 523.129 µS/cm, 735.272 µS/cm, 1021.533 µS/cm, respectively.
A GPR forward simulation was performed on the geoelectric models in the two cases, and then the GPR single-channel signals were analyzed by the WPEA method. The analysis results are shown in Figure 7. Figure 7a shows that when the electrical conductivity of the soil remained unchanged, the WPEI of the GPR signals continued to increase as the relative dielectric constant of the soil increased. Figure 7b illustrates that when the relative dielectric constant of the soil was constant, as the electrical conductivity of the soil increased, the WPEI of the GPR signals decreased. constant was 8.897; this value also remained unchanged. When the soil water con changed from 10% to 24%, the corresponding electrical conductivities were 88.681 μS 126.085 μS/cm, 180.261 μS/cm, 258.134 μS/cm, 368.670 μS/cm, 523.129 μS/cm, 735 μS/cm, 1021.533 μS/cm, respectively.
A GPR forward simulation was performed on the geoelectric models in the cases, and then the GPR single-channel signals were analyzed by the WPEA method. analysis results are shown in Figure 7.

Experimental Setup
To ensure the reliability of this GPR test, four model boxes were set up. The heights of the model boxes were set to 5 cm, 10 cm, 15 cm and 20 cm, and the length and width were each 20 cm. In this experiment, the clay water content was set to 10%, 12%, 14%, 16%, 18%, 20%, 22%, and 24%. Using to the clay with different water contents, a total of eight GPR detection experiments were carried out. In the same test, the clay with the same water content was poured into four model boxes for the experiment at the same time. Figure 8 shows clay samples with different water contents. Among them, Figure 8a is clay with 10% water content, Figure 8b with 12% water content, Figure 8c with 14% water content, Figure 8d with 16% water content, Figure 8e with 18% water content, Figure 8f with 20% water content, Figure 8g with 22% water content, and Figure 8h with 24% water content.  Figure 7a shows that when the electrical conductivity of the soil remained unchanged, the WPEI of the GPR signals continued to increase as the relative dielectric constant of the soil increased. Figure 7b illustrates that when the relative dielectric constant of the soil was constant, as the electrical conductivity of the soil increased, the WPEI of the GPR signals decreased.

Experimental Setup
To ensure the reliability of this GPR test, four model boxes were set up. The heights of the model boxes were set to 5 cm, 10 cm, 15 cm and 20 cm, and the length and width were each 20 cm. In this experiment, the clay water content was set to 10%, 12%, 14%, 16%, 18%, 20%, 22%, and 24%. Using to the clay with different water contents, a total of eight GPR detection experiments were carried out. In the same test, the clay with the same water content was poured into four model boxes for the experiment at the same time. Figure 8 shows clay samples with different water contents. Among them, Figure 8a is clay with 10% water content, Figure 8b with 12% water content, Figure 8c with 14% water content, Figure 8d with 16% water content, Figure 8e with 18% water content, Figure 8f with 20% water content, Figure 8g with 22% water content, and Figure 8h with 24% water content. To facilitate molding, the sand box was made of reinforced concrete with a length, width and height of 4 m, 3 m and 1.5 m respectively. The horizontal direction of the box was set as the x-axis, the vertical inward was the y-axis, and the vertical downward was the z-axis. The four model boxes were buried in order from low to high every 0.6 m along the y-axis, and the embedding depth was 10 cm. The specific layout of the four model boxes is shown in Figure 9. To facilitate molding, the sand box was made of reinforced concrete with a length, width and height of 4 m, 3 m and 1.5 m respectively. The horizontal direction of the box was set as the x-axis, the vertical inward was the y-axis, and the vertical downward was the z-axis. The four model boxes were buried in order from low to high every 0.6 m along the y-axis, and the embedding depth was 10 cm. The specific layout of the four model boxes is shown in Figure 9.

GPR Data Collection
The detection depth of GPR is related to the antenna frequency. When the antenna frequency of GPR is higher, the attenuation speed of the electromagnetic wave in the formation process is faster, the resolution is higher, and the detection depth is lower. The heights of the model boxes in the sand box was 5 cm, 10 cm, 15 cm and 20 cm, and the buried depth was 10 cm. In this experiment, a ground probing radar system produced by Ingegneria Dei Sistemi, Italy, was used, and a 1600 MHz central-frequency shielded antenna was used to collect GPR signals from the model boxes. The clay sample shown in Figure 8a was filled with the four model boxes and then buried in the sand box. Since the model boxes were buried in specific positions, the GPR signals were collected at least 0.5 m from the sand box boundary during the GPR detection process to reduce interference from the boundary environment. During the sampling process, the GPR antenna passed above the model boxes, collecting one signal per 1 cm and a total of 300 single-channel signals. A total of 457 data points were collected for each single-channel signal within the 18 ns of the sampling time window. The GPR imaging results of the clay in the four model boxes are shown in Figure 10, where the positions of the four red boxes represented the reflections of the four model boxes in the GPR image.

GPR Data Collection
The detection depth of GPR is related to the antenna frequency. When the antenna frequency of GPR is higher, the attenuation speed of the electromagnetic wave in the formation process is faster, the resolution is higher, and the detection depth is lower. The heights of the model boxes in the sand box was 5 cm, 10 cm, 15 cm and 20 cm, and the buried depth was 10 cm. In this experiment, a ground probing radar system produced by Ingegneria Dei Sistemi, Italy, was used, and a 1600 MHz central-frequency shielded antenna was used to collect GPR signals from the model boxes. The clay sample shown in Figure 8a was filled with the four model boxes and then buried in the sand box. Since the model boxes were buried in specific positions, the GPR signals were collected at least 0.5 m from the sand box boundary during the GPR detection process to reduce interference from the boundary environment. During the sampling process, the GPR antenna passed above the model boxes, collecting one signal per 1 cm and a total of 300 single-channel signals. A total of 457 data points were collected for each single-channel signal within the 18 ns of the sampling time window. The GPR imaging results of the clay in the four model boxes are shown in Figure 10, where the positions of the four red boxes represented the reflections of the four model boxes in the GPR image.
In Figure 10, the horizontal direction is the true detection distance in meters, and the vertical direction is the two-way travel time in nanoseconds. Figure 10 shows that the direct wave signals have a large amplitude and high energy compared to the electromagnetic wave reflection signals; these characteristics seriously interfere with the reflection signals of the model boxes and reduce the quality and resolution of GPR images. In addition, more hyperbolic reflections and diffractions appeared in the GPR images as the clay moisture content in the model boxes increased, that is, the greater the contrast of the relative dielectric constant between the clay and the dry sand, the more local diffraction in the GPR images was observed. m from the sand box boundary during the GPR detection process to reduce interference from the boundary environment. During the sampling process, the GPR antenna passed above the model boxes, collecting one signal per 1 cm and a total of 300 single-channel signals. A total of 457 data points were collected for each single-channel signal within the 18 ns of the sampling time window. The GPR imaging results of the clay in the four model boxes are shown in Figure 10 In Figure 10, the horizontal direction is the true detection distance in meters, and the vertical direction is the two-way travel time in nanoseconds. Figure 10 shows that the direct wave signals have a large amplitude and high energy compared to the electromagnetic wave reflection signals; these characteristics seriously interfere with the reflection signals of the model boxes and reduce the quality and resolution of GPR images. In addition, more hyperbolic reflections and diffractions appeared in the GPR images as the clay moisture content in the model boxes increased, that is, the greater the contrast of the relative dielectric constant between the clay and the dry sand, the more local diffraction in the GPR images was observed.

Preprocessing of GPR Signals
The GPR postprocessing analysis software was used to perform static correction, DC bias, gain, and band-pass filtering on the GPR signals shown in Figure 10. The results are shown in Figure 11. The results show that the resolution of the GPR images in Figure 11 is high compared with that in Figure 10 when the direct wave and noise are removed; as such, the reflected signals of the model boxes can be displayed more clearly and intuitively. The results show that the resolution of the GPR images in Figure 11 is high compared with that in Figure 10 when the direct wave and noise are removed; as such, the reflected signals of the model boxes can be displayed more clearly and intuitively.

Time-Domain Signals Response
The clays shown in Figure 8b-h were filled separately into the model boxes, and then buried in sand boxes one by one. GPR was employed to detect the model boxes in the sand box, and the GPR time profiles of eight kinds of clays with different water contents were obtained. Then, from each GPR time profile, a typical single-channel signal that could reflect a model box was extracted; the results are shown in Figure 12. Figure 12a-d show the obtained GPR single-channel signals in the first to fourth model boxes. Each model box included eight kinds of soil with different water contents.
Remote Sens. 2021, 13, x FOR PEER REVIEW 16 The clays shown in Figure 8b-h were filled separately into the model boxes, then buried in sand boxes one by one. GPR was employed to detect the model box the sand box, and the GPR time profiles of eight kinds of clays with different water tents were obtained. Then, from each GPR time profile, a typical single-channel si that could reflect a model box was extracted; the results are shown in Figure 12.

WPEI of GPR Signals
According to the wavelet basis described in Section 2.1, the GPR single-cha signals shown in Figure 12 was subjected to wavelet packet decomposition. Then WPEA method described in Section 2.3 was used to calculate the wavelet packet co cients of the GPR single-channel signals, and the WPEI of the four model boxes different water contents was obtained. The results are shown in Figure 13. Figure 13 shows the WPEI of the GPR single-channel signals in the first to fourth model boxes.

WPEI of GPR Signals
According to the wavelet basis described in Section 2.1, the GPR single-channel signals shown in Figure 12 was subjected to wavelet packet decomposition. Then, the WPEA method described in Section 2.3 was used to calculate the wavelet packet coefficients of the GPR single-channel signals, and the WPEI of the four model boxes with different water contents was obtained. The results are shown in Figure 13. Figure 13a-d shows the WPEI of the GPR single-channel signals in the first to fourth model boxes.
These results suggest that the WPEI of the GPR signals continued to increase as the water contents of clay increased, with a strong correlation being observed between these two parameters. However, in the case of the same water content in the four model boxes, the WPEI values of the GPR signals were slightly different. That is to say, the slopes of the four fitted curves are not exactly the same. One possible reason is that the size of the model box and the unevenness of the clay filling had a certain impact on the attenuation law of the electromagnetic waves inside the model boxes.

WPEI of GPR Signals
According to the wavelet basis described in Section 2.1, the GPR single-channel signals shown in Figure 12 was subjected to wavelet packet decomposition. Then, the WPEA method described in Section 2.3 was used to calculate the wavelet packet coefficients of the GPR single-channel signals, and the WPEI of the four model boxes with different water contents was obtained. The results are shown in Figure 13. Figure 13a-d shows the WPEI of the GPR single-channel signals in the first to fourth model boxes. Four model experiments were repeated in four model boxes, and the obtained laws were consistent with each other, thereby verifying the reliability of the method. A linear function was used to curve-fit the WPEI of the clay with different water contents in the four model boxes. The obtained linear equations were: 14  These results suggest that the WPEI of the GPR signals continued to increase as the water contents of clay increased, with a strong correlation being observed between these two parameters. However, in the case of the same water content in the four model boxes, the WPEI values of the GPR signals were slightly different. That is to say, the slopes of the four fitted curves are not exactly the same. One possible reason is that the size of the model box and the unevenness of the clay filling had a certain impact on the attenuation law of the electromagnetic waves inside the model boxes.

Comparative Analysis
To verify the feasibility of the wavelet basis defined in this paper in the wavelet packet analysis of the GPR signals, classical wavelet bases were also used to perform wavelet packet transform on the GPR single-channel signals, and the analysis results of the two were compared. Many researchers have shown that the wavelet bases of the

Comparative Analysis
To verify the feasibility of the wavelet basis defined in this paper in the wavelet packet analysis of the GPR signals, classical wavelet bases were also used to perform wavelet packet transform on the GPR single-channel signals, and the analysis results of the two were compared. Many researchers have shown that the wavelet bases of the Daubechies series have the characteristics of regularity, smoothness and tight support. As such, they have been widely used for the analysis and processing of GPR signals [70,71]. The db6 wavelet basis belongs to the wavelet of order 6 in the Daubechies wavelet series, and the db6 wavelet basis is an orthogonal wavelet which is tightly supported and asymmetric. The db6 wavelet basis (wavelet of order 6) was adapted to perform wavelet packet transform on the GPR single-channel signals, as shown in Figure 12. The wavelet packet coefficients were calculated using the WPEA method, and the WPEI with different water contents in the four model boxes were obtained. The results are shown in Figure 14.  Figure 14 demonstrates that as the water content of the clay increases, the WPEI of the GPR signal shows a gradually increasing trend, even though this is not always consistent. In Figure 14a, the WPEI of the GPR signals gradually increases before the soil moisture content is 22%, but the WPEI decreases when the soil moisture content is 24%. In Figure 14b, before the soil moisture content reaches 18%, the WPEI of the GPR signals gradually increases. When the soil moisture content is 18%, the value of the WPEI reaches the maximum. After that, the value of the WPEI drops. In Figure 14c, the WPEI of the GPR signals basically shows an increasing trend before the soil moisture content reaches 20%. The value of the WPEI reaches the maximum when the soil moisture content is 20%, and then decreases. In Figure 14d, except for the fact that the soil water content is 20%, the WPEI of the GPR signals basically increases with the increase of the soil water content. As such, a linear function was used to curve-fit the WPEI of the clay with different water contents in the four model boxes. The obtained linear equations were: 14

Discussion
Previous studies have mostly used the oven-drying method to determine the soil moisture content, which causes serious damage to the soil structure [10][11][12][13]. In addition, a limited number of on-site samples hinders the application range of this method [14]. In this study, the spatial distribution information of soil water content can be obtained using the GPR method. The use of GPR nondestructive testing methods to assess soil moisture  Figure 14 demonstrates that as the water content of the clay increases, the WPEI of the GPR signal shows a gradually increasing trend, even though this is not always consistent. In Figure 14a, the WPEI of the GPR signals gradually increases before the soil moisture content is 22%, but the WPEI decreases when the soil moisture content is 24%. In Figure 14b, before the soil moisture content reaches 18%, the WPEI of the GPR signals gradually increases. When the soil moisture content is 18%, the value of the WPEI reaches the maximum. After that, the value of the WPEI drops. In Figure 14c, the WPEI of the GPR signals basically shows an increasing trend before the soil moisture content reaches 20%. The value of the WPEI reaches the maximum when the soil moisture content is 20%, and then decreases. In Figure 14d, except for the fact that the soil water content is 20%, the WPEI of the GPR signals basically increases with the increase of the soil water content. As such, a linear function was used to curve-fit the WPEI of the clay with different water contents in the four model boxes. The obtained linear equations were: y 1 = 1.357 × 10 14 x 1 − 1.015 × 10 13 , y 2 = 5.205 × 10 13 x 2 − 3.761 × 10 12 , y 3 = 4.249 × 10 13 x 3 − 2.759×10 12 , and y 4 = 6.399 × 10 13 x 4 − 3.368 × 10 12 , and the correlation coefficients were 0.85, 0.74, 0.85 and 0.92, respectively. This shows that the use of db6 wavelet as the basis for wavelet packet transform does reveal that the WPEI gradually increases in clays with different water contents.

Discussion
Previous studies have mostly used the oven-drying method to determine the soil moisture content, which causes serious damage to the soil structure [10][11][12][13]. In addition, a limited number of on-site samples hinders the application range of this method [14]. In this study, the spatial distribution information of soil water content can be obtained using the GPR method. The use of GPR nondestructive testing methods to assess soil moisture content has an obvious advantage, that is, GPR can detect the moisture content over large areas in a short time [1,6].
The results of the forward simulation and indoor model experiments in this paper showed that the changes of the electromagnetic wave characteristics in GPR measurements were related to the measured values of moisture content in each soil sample, as indicated in the literature. Increasing the soil water content in the model box will change the relative permittivity and electrical conductivity in the soil, and affect the propagation speed of electromagnetic waves in the soil, thus impacting the GPR signal, which helps us to learn more about the physical properties of soil by analyzing the GPR signal [42,65].
By using wavelet packet transform with good time-frequency localization characteristics, the nonlinear and nonstationary GPR signal features may be extracted effectively [44][45][46][47]. However, there is no good similarity between the conventional wavelet bases and the GPR signals [49][50][51]. According to the time-frequency characteristics of GPR signals, this paper constructed the most suitable basis for wavelet packet analysis of GPR signals. Using this basis to perform wavelet packet transform of GPR signals helped us to better understand the propagation law of electromagnetic waves under different soil moisture conditions, and conveniently provided a way for us to study the relationship between the time-frequency characteristics of GPR signals and soil moisture content.
The accuracy of the GPR method used to estimate the soil moisture content was determined based on precise and appropriate empirical relationships [32,68]. In order to evaluate the moisture content of the soil, a fitting function between the WPEI of the GPR signals and the soil moisture content was established. In this way, the fitting function could be employed to evaluate the soil moisture content of roadbeds, slopes, earth dams and foundations over large areas. It is worth noting that the correlation coefficient between the WPEI of the GPR signal and the soil moisture content exceeded 0.85, which is better than or comparable to the results reported in related studies [20,39].
As shown in Figure 13, the WPEI of the GPR signals had a highly correlated linear relationship with the soil water content, but the result was obtained under perfect experimental conditions. It is also worth noting that the propagation characteristics of GPR electromagnetic waves in soil were easily affected by the physical and mechanical properties of the soil, including mineral composition, particle size, shape and gradation, density, compactness, porosity, degree of cementation, and stress history. In addition, the influences of groundwater level and seepage on the imaging results of GPR were not considered in this paper.

Conclusions
In this paper, a method for detecting soil water content using GPR technology is proposed. With the increase of soil water contents, the WPEI of the GPR signals also gradually increase. The WPEA method can be used to calculate the WPEI of the GPR signals to quantitatively evaluate the soil water content. The test results show that there is a highly correlated linear relationship between the WPEI and the soil water content, which can expand the application of GPR in slope monitoring, ground disaster prevention, subgrade, earth dam, and foundation compaction. In future work, the team will conduct research on the following three aspects: First, the soil water content detection experiment will be further refined to ensure that the soil density, compactness and porosity are basically constant when the soil water content changes. Second, a variety of soil samples will be used in model tests, including gravelly soil, sandy soil and cohesive soil, as well as soil samples of different particle sizes, shapes, and gradations. Third, according to the relationship between the WPEI of the GPR signals and the soil water content, a suitable model will be established to quantitatively evaluate the soil water content.
Author Contributions: S.Z., T.L. and G.F. contributed to the design and implementation of the research, S.Z., L.Z. and Y.G. contributed to the analysis, processing and interpretation of the data, S.Z., L.Z. and T.L. wrote the article with input from all authors. All authors have read and agreed to the published version of the manuscript.