Estimation of Moisture Content in Railway Subgrade by Ground Penetrating Radar

Abstract

China is strongly dependent on railway transportation, but the frost heaving of the subgrade in cold regions has seriously affected the safety and comfort of trains. Moisture content is an essential parameter in the subgrade frost heave. Non-destructive and efficient geophysical methods have great potential in measuring the moisture content of railway subgrade. In this paper, we use the common mid-point (CMP) measurement of ground penetrating radar (GPR) to estimate the propagation velocity of electromagnetic waves in a subgrade application. We establish a synthetic model to simulate the railway subgrade structure. The synthetic CMP gathers acquired from shallow and thin layers are seriously disturbed by multiple waves and refraction waves, which make the routine velocity analysis unable to provide accurate velocities. Through the analysis of numerical simulation results, it is found that the primary reflection waves, multiple waves, and refraction waves are dominant in different offset ranges of CMP gather. Therefore, we propose a solution of the optimal gather at a certain range of offset dominated by the primary reflection wave to calculate the velocity spectrum and extract the accurate velocities for the subgrade model. The relative dielectric constants of the corresponding layers are calculated after the stacking velocities are converted into the interval velocities. Then, the moisture content is obtained by the Topp formula, which expresses the relationship between dielectric constant and moisture content. Finally, we apply the optimal gather scheme and the above interpretation process to the GPR data acquired at the railway site, and we form a long moisture content profile of the railway subgrade. Compared with the polarizability measured by the induced polarization (IP) method, it is found that the regions with high moisture content correspond to polarizability anomalies with different strengths. The comparison shows the reliability of GPR results to some extent.

1. Introduction

China has a vast territory, a large population, and an unbalanced distribution of resources and industries. Therefore, China’s railway plays an important role in economic and social development. Subgrade frost heaving is a common geological phenomenon in northern China, and its treatment is the main task of railway maintenance in winter. Frost heaving of the subgrade seriously affects the safety and comfort of train operation. Temperature, soil properties, and moisture content are the three factors of the formation of subgrade frost heaving. The main cause of soil freezing is the moisture content in soil. Therefore, it is crucial to measure the moisture content of railway subgrade quickly and accurately for solving the problem of frost heaving.

Drilling samples can obtain the most accurate information about the subgrade moisture content [1]. However, this method is time-consuming and expensive for large-scale measurement. Ground penetrating radar (GPR) is a geophysical technique for the non-destructive and efficient detection of shallow layers by electromagnetic (EM) waves, which is widely used in the field of environment and engineering. Using GPR to measure the moisture content has become an important branch of GPR applications [2]. In recent years, a lot of papers using GPR to measure soil moisture content have been published [3,4,5,6,7,8,9,10,11,12]. Some of them focus on ground water table measurement. As an example, Tsoflias et al. observe that the radar signal response is correlated to changes in the water saturation of the fracture and provide spatial information about the saturation of the fracture [5]. Saintenoy and Hopmans study the water table detection sensitivity of GPR and state a power type relationship between the reflected signal amplitude and the slope of the soil retention curve [6]. Others are concerned about the measurement of soil moisture content. Two review papers are valuable. Huisman et al. (2003) presents the basic theoretical principles of GPR and how they can be used to investigate the spatiotemporal variation of soil water content [10]. In particular, they propose four categories of methods to determine soil water content with GPR. Klotzsche et al. (2018) provide an update on the review of recent advances in vadose zone applications of GPR with a particular focus on new possibilities, e.g., for multi-offset and borehole GPR measurements [11].

A large number of GPR applications in measuring soil moisture content has benefited from the research of Topp et al., which established the relationship between the relative dielectric constant and soil moisture content [13]. It is essential to estimate the velocity of EM waves, since it is directly related to the dielectric constant of medium. The common mid-point (CMP) or wide-angle reflection and refraction (WARR) sounding mode in GPR are primarily used to obtain an estimate of the radar wave velocity versus depth in the ground by varying the antenna spacing and measuring the change in the two-way travel time [2]. The velocity analysis of CMP gather, coming from seismic exploration, is a common and accurate way to obtain wave velocity. Some scholars have carried out research studies to get a better velocity estimation in GPR. Feng et al. improved the signal to noise ratio by CMP antenna array and data processing technology and obtained a good result in velocity analysis, which they successfully applied to landmine detection [14]. Liu et al. obtain the EM wave velocity through the envelope velocity analysis method, which can monitor the underground dynamic water level [15]. Lu described the method to obtain the soil moisture content and monitored the underground water level [16]. Liu et al. used GPR profiles and moisture content to estimated hydraulic conductivity parameters and accurately distinguished a slight change of groundwater level [17]. Pue et al. proposed a modified velocity analysis method to successfully estimate depth and propagation velocity with small offset air-coupled GPR configurations by accounting for the refraction at the surface [18]. Yi et al. proposed a high-resolution velocity analysis by applying the $l$-1 norm regularized least-squares for pavement inspection [19]. In recent years, many other related papers have been published to improve the performance of CMP velocity analysis in the estimation of EM wave velocity [20,21,22]. However, there are few research studies that deal with shallow and thin multi-layered structures such as the railway subgrade. In such a case, the reliability of the routine velocity analysis process needs to be retested.

In this paper, we introduce the basic theory used here firstly. Then, two synthetic models are established to simulate the railway subgrade structure and the problems of velocity analysis in a shallow thin layer are analyzed by numerical simulation. After that, the solution of optimal gather is proposed. The proposed scheme is applied to the measured data of railway subgrade, and a long moisture content profile of railway subgrade is obtained. Finally, the polarization features from induced polarization (IP) measurement show that the proposed scheme is reliable in a certain range.

2. Methods

2.1. Velocity Analysis of CMP Data

The common mid-point (CMP) measurement in a GPR survey is a multi-offset measurement that is primarily used to obtain an estimated velocity of the EM wave in the ground. By increasing the antenna preparation symmetrically, the path of the EM wave varies while keeping the point of reflection fixed. Thus, the difference in two-way travel time enables wave velocity to be estimated [2].

The velocity analysis by CMP gather assumes that the subsurface media is horizontally layered. The two-way travel time of reflected waves can be obtained by

$$t\left(x_i\right)=\sqrt{t_0^2+\frac{x_i^2}{v_{rms}^2}}$$

(1)

where $x_i$ is the antenna offset of the i-th channel, $t_0$ is the two-way travel time of zero offset, and $v_{rms}$ is the root mean square velocity.

The velocity spectrum can be obtained by transforming the data from the offset versus two-way time domain to the stacking velocity versus two-way zero-offset time domain [23]. The velocity spectrum shows the picks that correspond to the best coherency of the signal along a hyperbolic trajectory over the entire spread length of the CMP gather. The coherency can be computed in several ways [24]. However, in the velocity analysis of GPR data, the simple and effective measures are the average stacked amplitude ($A$) and the energy of the average stacked amplitude ($E$), which can suppress interference signals with large amplitude. The $A$ and $E$ can be defined by

$$A=\frac{1}{N}{\displaystyle \sum }_{j=1}^M\left|{\displaystyle \sum }_{i=1}^N{f}_{i,j+r_i}\right|$$

(2)

$$E={\displaystyle \sum }_{j=1}^M{\left(\frac{1}{N}{\displaystyle \sum }_{i=1}^N{f}_{i,j+r_i}\right)}^2$$

(3)

where $f_{i,j+r_i}$ is the amplitude value of the j+r-th point on the i-th trace in which r is a delay associated with normal moveout (NMO). $N$ is the number of traces in the CMP gather. $M$ is the number of sampling points in a stack time window. Each stacking velocity ($v_{rms}$) is scanned according to Equation (1). When the scanning velocity is equal to the $v_{rms}$, $E$ reaches the maximum value, which shows a focused energy cluster in the velocity spectrum. After the stacking velocity is picked up, the interval velocity can be calculated by Dix formula [25], and the expression is as follows:

$$v_{int,n}=\sqrt{\frac{t_{0,n}{v}_{rms,n}^2-t_{0,n-1}{v}_{rms,n-1}^2}{t_{0,n}-t_{0,n-1}}}$$

(4)

where $v_{int,n}$ is the interval velocity of the nth layer; $v_{rms,n}$ is the stacking velocity to the bottom of the nth layer; $t_{0,n}$ is the two-way travel time with zero offset to the bottom of the nth layer; $v_{rms,n-1}$ is the stacking velocity to the bottom of the n − 1th layer; and $t_{0,n-1}$ is the two-way travel time with zero offset to the bottom of the n − 1th layer.

After obtaining the interval velocities ($v_{int}$) of radar wave propagating in the subsurface, the relative dielectric constant ${\epsilon }_r$ in low loss medium can be obtained by

$$v_{int}=\frac{c}{\sqrt{{\epsilon }_r}}$$

(5)

where $c=0.3\mathrm{m}/\mathrm{ns}$ is the velocity of EM wave in air.

2.2. Optimal Gather in Velocity Analysis

As the offset increases, the primary reflections from the interface of each layer will be interfered by multiple reflections and refractions in CMP gathers. In the media composed of thin multi-layers, this interference is so serious that it makes the velocity spectra unable to give correct velocities of EM waves. In addition, the CMP gathers obtained in small offsets can be interfered by direct wave. Since the layers are shallow and thin, the reflected wave is close to the direct wave, whose strong energy makes the reflected wave difficult to be identified. Therefore, it is crucial to use high-quality primary reflections from layer interfaces for velocity analysis to obtain accurate velocities.

The optimal gathers include mainly the primary reflections who suffer little interference from multi-waves, refracted waves, and direct waves. These gathers are determined and used to calculate the velocity spectrum. The reflection angle can be used to explain the operation of optimal gather.

The reflection angle $\theta$ can be calculated using the sketch shown in Figure 1. $x_1$ represents a small offset, and the CMP gathers acquired in $x_1$ are supposed to be contaminated with direct waves. $x_2$ represents a large offset, and the CMP gathers acquired in $x_2$ are supposed to be contaminated with multi-waves and refracted waves. The blue arrows in Figure 1 show the corresponding antenna locations when effective reflection from the interface can be obtained, and the reflection angle $\theta$ ranges from ${\theta }_1$ to ${\theta }_2$. The optimal gather for the velocity analysis should be in this angle range.

In this paper, we can analyze the influence range of interference waves in CMP gathers through numerical simulation because the railway subgrade structure is known. It can help to determine the optimal gather in field CMP data.

2.3. Topp Formulae

A well-known empirical relationship between the relative dielectric constant and volumetric moisture content at frequencies between 1 MHz and 1 GHz was developed by Topp et al. [13]. This empirical relationship is independent of soil type, soil density, soil temperature, and soluble salt content, and it is determined by compiling data for many soils under varying moisture conditions:

$${\epsilon }_r=3.03+9.3\theta +146.0{\theta }^2-76.7{\theta }^3$$

(6)

where $\theta$ is the volumetric soil moisture content, which is equal to the product of the water saturation and porosity. A polynomial expression for moisture content as a function of the relative dielectric constant,${\epsilon }_r$, was reported as well:

$$\theta =-5.3\times {10}^{-2}+2.92\times {10}^{-2}{\epsilon }_r-5.5\times {10}^{-4}{\epsilon }_r^2+4.3\times {10}^{-6}{\epsilon }_r^3.$$

(7)

2.4. Induced Polarization

Induced polarization (IP) is a geophysical method that is mainly used in metal ore exploration and groundwater search. In recent years, with the development of sensitive instruments and methods to improve the signal-to-noise ratio, the IP method has been widely used in the field of environment and engineering [26,27,28]. IP may reflect the moisture content through the membrane effect, which is related to variations in the mobility of ions in fluids at grain scales in a porous medium under the influence of the electrical field. The polarizability generated by the water-bearing sand gravel layer is higher than that of non-water bearing structure.

In time domain IP, a primary current is injected in the ground for a period T using a pair of electrodes A and B. The difference of electrical potential is measured using several bipoles of electrodes M and N. The secondary voltage on a bipole MN decays over time after the shut-down of the primary current. The rate of this decay reflects the strength of the polarization effect. The simplest way to measure the IP effect with time-domain equipment is to compare the residual voltage $V\left(t\right)$ existing at a time t after the current is cut off with the steady voltage $V_c$ during the current flow interval. It is not possible to measure potential at the instant of cutoff because of large transients caused by breaking the current circuit. On the other hand, $V\left(t\right)$ must be measured before the residual has decayed to the noise level. Therefore, the polarizability, $\eta \left(t\right)$ can be obtained by

$$\eta \left(t\right)=\frac{V\left(t\right)}{V_c}\times 100\%.$$

(8)

3. Velocity Analysis of CMP Models

3.1. Railway Subgrade

Figure 2 shows a typical railway subgrade structure, which can be roughly divided into 4 layers. These are constituted with ballast, graded gravel, medium to coarse-grained sand, and A, B group filling from top to bottom. There is a sheet of geotextile in the middle of the 0.2 m sand layer (marked by a red line in Figure 2).

The horizontal layered structure of the railway subgrade satisfies the assumption of velocity analysis through CMP gather. However, due to the shallow depth and the thin multi-layers, the reflected waves from the interfaces of horizontal layers will be interfered by multiple waves and refracted waves. Once the hyperbolic trajectory of the reflected wave in CMP gathers is damaged, the velocity spectrum will not obtain the picks representing the best coherency and failure to pick up the accurate velocity.

3.2. Modeling

In this section, we carry out CMP forward simulation on a simple model and a railway subgrade model to study the interference caused by multiple waves and refraction waves and propose solutions. The GPR forward modeling uses an open source software gprMax that simulates EM wave propagation by the Finite-Difference Time-Domain (FDTD) method [29].

3.2.1. A Simple Model

A simple model (Figure 3) is used firstly to analyze the sources of interference. In Figure 3, the air layer is marked by a white color with ${\epsilon }_r=1$. The thickness of the blue layer is 0.7 m with ${\epsilon }_r=4$. The velocity of EM wave 0.15 m/ns in the layer can be obtained by Equation (5). The red layer is with the ${\epsilon }_r=10$. The synthetic CMP data are generated by a 100 MHz Ricker wavelet and the grid size of 0.05 m. The mid-point of CMP gather is fixed at coordinate (15,12), as shown in Figure 3. The increment of antenna separation is 0.2 m, and 64 traces are acquired in the CMP gather.

The synthetic CMP gather is shown in Figure 4. It can be seen from Figure 4a that there are several discontinuous events. In order to analyze the sources of waveforms in CMP gathers, we calculate the arrival times of reflected, multiples, and refracted waves according to the model parameters, and the calculated values are shown by curves with different color. In Figure 4b, the black dotted line marked with ① is the direct air wave that the wave directly propagates from the transmitting antenna to the receiving antenna. The black dotted lines marked with numbers ②–④ are the refracted waves propagating in the air after they are reflected twice or more by the layer interface. The red dot line in Figure 4b represents the primary reflected wave from the layer interface, which we need. However, the continuity of the event is obviously destroyed by multiple reflections. The green dotted lines in Figure 4b represent reflection multiples. Equation (1) is the signal received by the antenna after it is reflected twice by the interface, and (2)–(4) are reflected multiple times by the interface before it is received by the antenna. The above analysis shows that the existence of multiple reflections brings serious interference to velocity analysis.

3.2.2. A Railway Subgrade Model

The typical railway subgrade model is shown in Figure 5. The white layer is air. The thickness of the blue layer corresponding to ballast is 0.7 m with ${\epsilon }_r=4$. The next layer marked green is graded gravel that is 0.6 m thick, and the ${\epsilon }_r=5.5$. The thin yellow layer of 0.2 m thickness is medium to coarse-grained sand with ${\epsilon }_r=7$. The thickness of the orange layer is A, B group filling that is 1.1 m thick and ${\epsilon }_r=8.5$. The red part is a layer of raw soil with the ${\epsilon }_r$ of 10. The mid-point of CMP gather is fixed at coordinate (15,12), as shown in Figure 5. The increment of antenna separation is 0.2 m, and 64 traces are acquired in the CMP gather.

Figure 6a shows the synthetic CMP gather of the railway subgrade model. Similarly, we use the model parameters to calculate the arrival time of radar waves, and the calculated wave curves are shown with different colors in Figure 6b. The analysis for the simple model in Section 3.2.1 makes it easy to distinguish the various waves generated in this complex model.

In Figure 6b, the black dotted line marked with ① is the direct wave in air. The next four black dotted curves with the labels of ②–⑤ are refracted waves, and the amplitudes of these curves decrease in turn. There are two green dotted lines in Figure 6b that represent multiple reflections. The upper green line represents the multiple wave that crosses the first layer twice and the second layer once, and its energy is mainly obtained at the offset of 5–12.8 m. The lower green curve crosses the first layer three times and the second layer once, and its energy is mainly obtained at the offset of 9–12.8 m.

The three red dotted lines represent the primary reflections from the layer interfaces, which are the target reflections. The top red curve is reflected from the first interface, whose energy is mainly obtained at offsets less than 5 m. However, the third layer that is 0.2 m thick is too thin to be identified.

The generation of these marked reflections can be explained by the reflection angle shown in Figure 1. First, when the reflection angle is small, the reflected wave is close to the direct wave whose strong energy makes the reflected wave difficult to be identified. Second, with the increase of reflection angle, the energy of multiple reflections increase, while the primary reflections become weaker until they cannot be detected. Therefore, the reflected wave can be distinguished only when the reflection angle is within a certain range. Next, we deal with optimal gathers carefully to ensure that the velocity spectrum can provide reliable velocity.

We calculate the velocity spectra with different offsets for CMP gathers of the railway subgrade model. Since the initial offset in the field CMP data is 0.6 m (shown in Section 4), the same initial offset is used in the velocity analysis of the model. Figure 7a,c,e,g show the velocity spectra after normalization with the CMP gathers at offset 0.6–3 m, 0.6–4 m, 0.6–5 m, and 0.6–6 m, respectively. Figure 7b,d,f,h shows the curves of the average amplitude, corresponding to Figure 7a,c,e,g, which is the mean value of the amplitude of all velocities at each arrival time. The velocity can be picked up by the correspondence between the energy cluster in the velocity spectrum and the wave crest in the average amplitude curve.

In the eight subgraphs of Figure 7, the velocity of EM waves from each layer is represented by energy clusters in the velocity spectra and wave peaks in the average amplitude curves marked by ①, ②, and ③. Other energy clusters and peaks, such as those that appear at the arrival time greater than 50 ns, are from the contribution of multiple waves.

In Figure 7a, when the offset (0.6–3 m) is small and the number of gathers is less, the energy cluster is decentralized, which brings trouble to velocity picking. As the offset increases, the energy clusters become more and more focused in Figure 7c,e. Combined with Figure 7d,f, velocity picking becomes easier. Although Figure 7e shows that the energy clusters are more focused than that in Figure 7c, the velocity of EM waves of the first layer is clearer in Figure 7c. In Figure 7g,h, it is difficult to pick up the velocity of first layer due to the large offset, which is 0.6–6 m. Therefore, a smaller offset makes the energy decentralized, while a larger offset makes the velocity of the first layer impossible to pick up. The optimal gather in this railway subgrade model is ended at the offset between 4 and 5 m with the initial offset of 0.6 m.

Then, the velocities ($v_{rms}$) are picked up from the velocity spectrum in Figure 7e. The $v_{rms}$ are converted to the $v_{int}$ by the Dix formula of Equation (4), which is plotted by a green line in Figure 8. The black line shows the velocities that are calculated from the ${\epsilon }_r$ of the model. It shows that the estimated $v_{int}$ in the first and second layer are basically in accordance with the model velocities. The third layer in the model is a medium to coarse-grained sand that is 0.2 m thick. GPR with the central frequency of 100 MHz failed to recognize this thin layer. Therefore, the velocity spectrum is also not able to give the velocity of this layer. It leads to the difference between the model velocity and $v_{int}$ around 20 ns in Figure 8. The model velocity of the fourth layer is 0.103 m/ns, while the $v_{int}$ of that is 0.106 m/ns (

Table 1. The detailed information estimated from the velocity spectrum in Figure 7e.

Layer Sequence	Time (ns)	${\mathit{v}}_{\mathit{r}\mathit{m}\mathit{s}}$ (m/ns)	${\mathit{v}}_{\mathit{i}\mathit{n}\mathit{t}}$ (m/ns)	Layer Thickness (m)	Dielectric Constant	Moisture Content
1	9.4	0.15	0.15	0.705	4	5.5%
2,3	23.3	0.138	0.129	0.89	5.39	8.9%
4	44.7	0.124	0.106	1.13	8.1	15%

). The failure to pick up the $v_{rms}$ of the third layer is part of the reason for the error. After the dielectric constant (${\epsilon }_r$) is obtained by Equation (5), the moisture content can be calculated by Equation (7). Table 1 shows the details of the estimated information from the velocity spectrum in Figure 7e.

It can be seen from Table 1 that the estimation of the first layer is accurate. The errors occur when the second and third layers are estimated as one. It also affects the estimation of the fourth layer. Our ultimate goal is to get moisture content under the subgrade. The model moisture content of the four layers is 5.53%, 9.17%, 12.6%, and 15.8%, respectively. Compared with the estimated water content in the last column of Table 1, although the big error occurs in layers 2 and 3, the vertical variation trend of moisture content is still credible, which is essential to identify a meaningful distribution of moisture content. According to the above optimal gather method and processing flow, the moisture content estimation in the railway subgrade is conducted in the following case study.

4. Case Study

We measured the GPR data in a certain railway in Heilongjiang, which is the coldest area in China. In Figure 9, a SIR-20 GPR system that employed the shielded antenna with the central frequency of 100 MHz from GSSI (Geophysical Survey Systems, Inc.) is used in the field survey. The CMP gather is acquired with the initial offset of 0.6 m between the transmitter and receiver. The increment of antenna separation is 0.2 m and the maximum antenna transceiver distance is about 14 m. The parameter setting in acquisition is consistent with that in the model. CMP gathers from 56 mid-points are acquired, and the separation between two adjacent mid-points is 10 m. In order to avoid the influence of sleepers, CMP data are acquired along the outside of the rail (see the photo in Figure 9).

The preprocessing for the original CMP gathers includes zero-time correction, gain control, and frequency domain filtering. Figure 10a,b show the original CMP gather acquired at the first CMP point and its profile after preprocessing, respectively.

The CMP gathers acquired at the offset of 0.6 m to 5 m are used in velocity analysis, which is decided from the tests on the model in Section 3.2.2. Figure 11 shows the velocity analysis result for the CMP gathers in Figure 10b. In Figure 11a, the velocity spectrum shows five energy clusters which are marked by ①–③, ⓐ, and ⓑ. Only the energy clusters labeled ①–③ have corresponding peaks in the average amplitude curve, which is also marked by ①–③ in Figure 11b, while ⓐ and ⓑ cannot be identified. The hyperbolic trajectories of the signals are plotted according to the velocities picked up from the five energy clusters, as shown in Figure 12. In Figure 12a, the hyperbolic trajectories generated by velocities ①–③ correspond to relatively complete hyperbolic events, which are meaningful. On the other hand, the clusters ⓐ and ⓑ, in Figure 12b, are invalid because they produce trajectories without corresponding hyperbolic events. Therefore, the energy clusters ①–③ are the reflections from the interfaces of the railway subgrade, and we extract the $v_{rms}$ from them to calculate the $v_{int}$. The $v_{int}$ is shown in Figure 11c, and the moisture content is calculated as shown in Figure 11d.

As shown in Figure 11d, the vertical moisture content ranges from 2.24% to 7.1% with depth, and there is no abnormal trend in moisture content change.

Table 2. The detailed information estimated from the velocity analysis in Figure 11.

Layer Sequence 1	Time (ns)	${\mathit{v}}_{\mathit{r}\mathit{m}\mathit{s}}$ (m/ns)	${\mathit{v}}_{\mathit{i}\mathit{n}\mathit{t}}$ (m/ns)	Layer Thickness (m)	Dielectric Constant	Moisture Content
1	6.5	0.182	0.182	0.59	2.71	2.24%
2	19.5	0.166	0.157	1.02	3.63	4.6%
3	34.5	0.155	0.139	1.05	4.63	7.1%

¹ The layers here correspond to the marks ①–③ in Figure 11a,b.

shows the details of related information from the estimates. The estimated total thickness is consistent with that of the actual railway subgrade structure, but some errors exist in the thickness of each layer. We think that the main reason for these errors is the velocity pick-up errors. It is also related to the resolution provided by the 100 MHz central frequency of the radar antenna, which may have a limitation in detecting thin layers such as in this case.

Figure 13 shows the velocity analysis for another CMP gather acquired at the 27th mid-point. Four energy clusters are considered to be valid, which are labeled by ①–④ in Figure 13a, and the corresponding wave crests can be found in the average amplitude curve (Figure 13b). The velocities, 0.2 m/ns (8 ns), 0.186 m/ns (14.5 ns), 0.150 m/ns (34 ns), and 0.135 m/ns (46 ns), are picked up. These velocities ($v_{rms}$) are converted to the $v_{int}$, which is shown in Figure 13c, and four velocity layers are identified. The thicknesses of the first two layers, which are considered as layers of ballast and graded gravel, are accurately estimated. The third layer of medium to coarse-grained sand with the thickness of 0.2 m is not distinguished. A velocity increasing can be found in the fourth layer of the A, B group filling. In Figure 13d, the profile shows the moisture content of 11.8% in the A, B group filling, and it increases to 27.1% at the bottom of the A, B group filling at the depth of 2.47 m. We direct attention to the abnormal increase of moisture content occurring at the bottom of the A, B group filling layer.

The total of 56 CMPs are acquired with the separation of 10 m between two adjacent mid-points. The velocity spectrum of each CMP gather is calculated by the optimal gather with an offset of 0.6 to 5 m. The valid energy clusters in the velocity spectrum are determined by judging their correspondence with the peak of the average amplitude curves and the events in CMP gathers. The stacking velocities ($v_{rms}$) are picked up and converted to interval velocities ($v_{int}$) by the Dix formula. Then, the relative dielectric constant (${\epsilon }_r$) can be obtained by Equation (5)—that is, ${\epsilon }_r={(c/v_{int})}^2$. At last, the moisture content can be calculated by Equation (7) of the Topp formula. The moisture content curves from 56 mid-points are interpolated into a profile of moisture content, as shown in Figure 14a.

In Figure 14a, high moisture content areas are mainly located at around 130–140 m, 180–190 m, 260–270 m, 310 m, 370 m, 420–450 m, and 520–540 m. There are two known culverts with a size of 4 m × 4 m along this survey line. One is located between 150 and 160 m and the other is at 400–410 m with a depth of about 2 m. The positions of the culverts are marked by white boxes in Figure 14a. The culverts are filled with air, and their existence affects the moisture estimation of their surroundings. These regions feature large velocities of EM waves, small dielectric constants, and small moisture content.

IP measurements are also conducted along the same survey line. A four-electrode sounding method is used in the IP measurement. In the IP acquisition, current electrodes A and B are stainless steel electrodes, and Pb-PbCl nonpolarizable electrodes are used for the potential electrodes M and N to avoid the charge-up effect. The distance between the electrodes of M and N is fixed to 0.2 m, and the distance between electrodes A to B is gradually increased. These separations are 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, and 10.0 m, respectively. The IP measuring points are consistent with the mid-points in CMP measurement, and the separation of two adjacent measuring points is 10 m. The polarizability of each measuring point is calculated by Equation (8). These 1D polarizability data are used to create a polarizability profile by interpolation, as shown in Figure 14b.

In Figure 14b, there are some anomalies of polarizability whose values are higher than that of the surroundings. Compared with Figure 14a, it is found that the positions of polarizability anomalies are consistent with those of the abnormal area of moisture content, although the level of polarizability and moisture content are not consistent. For example, there is an obvious high-water content area around 250 m in Figure 14a, but it shows a polarizability anomaly at the same position, which is weak. In Figure 14b there is an obvious polarization anomaly at 100–150 m, which is stronger than that seen at 250 m. However, the moisture content evaluated at 100–150 m is lower than that evaluated at 250 m in Figure 14a.

These observations indicate that the relationship between polarizability and moisture content is complex. There is no doubt that an IP survey can find water. However, the dominant factor controlling the strength of polarizability should be the ion concentration in water, not only the amount of water. Therefore, we do not use the strength of polarization to show the moisture content, since we have no information about the ion concentration in soil solution. We compare the positions of abnormal polarizability with that of moisture content from GPR data. We state that no matter how strong or weak, as long as anomalies of the moisture content and polarizability appear in the same position, they can be considered as representing the same target. As a result, this comparison shows that the moisture content extracted from velocity analysis is reliable in the case of railway subgrade.

5. Conclusions

In this paper, we use GPR to detect the moisture content of railway subgrade. The flow of the processing and interpretation includes: first, to obtain the stacking velocity of each layer of subgrade by the velocity analysis of CMP gathers; second, to convert the stacking velocities into interval velocities by the Dix formula and then calculate the dielectric constant; finally, to obtain the moisture content distribution of subgrade by the Topp formula. In this flow, it is the most important thing to pick up the accurate stacking velocity.

The shallow and thin multi-layers of railway subgrade make CMP gather seriously disturbed by multiple waves and refraction waves, which makes the routine velocity analysis unable to provide accurate velocity. In the railway subgrade case, we propose a solution of the optimal gather from a certain offset range dominated by the primary reflection to calculate the velocity spectrum. The extraction of the valid velocities can be completed by judging their correspondence with the peak of the average amplitude curves and the events in CMP gathers. The proposed optimal gather and interpretation flow are applied to the CMP data acquired at a railway site, and a long profile of moisture content distribution is obtained. IP measurement along the same survey line shows the abnormal regions of polarizability whose positions are consistent with large water content areas extracted from GPR data.

In the field survey, a GPR system equipped with a shielded antenna with a central frequency of 100 MHz is used. This frequency of the antenna has a limitation in the detection of the shallow and thin layer structure of the subgrade. The shielded antenna with the center frequency of 200 MHz may give a better result which is known through numerical simulation (not shown in the text). Unfortunately, we do not have such equipment. However, it is gratifying to see that the results given by this paper show that the antenna with a center frequency of 100 MHz is feasible. In conclusion, the velocity analysis of CMP gathers in GPR has great potential in quantitatively estimating the moisture content of the railway subgrade in the construction and maintenance of railway in cold regions.