*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 ε*r* = 1. The thickness of the blue layer is 0.7 m with ε*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 ε*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.

**Figure 3.** A simple model.

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 -1 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 -2 –-4 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.

**Figure 4.** Synthetic common mid-point (CMP) gather for the simple model (**a**) with the direct wave, refracted waves, and reflected waves marked on (**b**).

#### 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 ε*r* = 4. The next layer marked green is graded

gravel that is 0.6 m thick, and the ε*r* = 5.5. The thin yellow layer of 0.2 m thickness is medium to coarse-grained sand with ε*r* = 7. The thickness of the orange layer is A, B group filling that is 1.1 m thick and ε*r* = 8.5. The red part is a layer of raw soil with the ε*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 5.** The typical railway subgrade model. The labels of I–IV correspond to the layers I–IV in Figure 2.

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.

**Figure 6.** CMP gather for the railway subgrade model (**a**) with the direct wave, refracted waves, and reflected waves marked on (**b**).

In Figure 6b, the black dotted line marked with -1 is the direct wave in air. The next four black dotted curves with the labels of -2 –-5 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.

**Figure 7.** Velocity spectra and curves of average amplitude for the railway subgrade model are calculated using different offsets. (**<sup>a</sup>**,**c**,**e**,**g**) are velocity spectra with offsets of 0.6–3 m, 0.6–4 m, 0.6–5 m, and 0.6–6 m, respectively; (**b**,**d**,**f**,**h**) are curves of average amplitude calculated from (**<sup>a</sup>**,**c**,**e**,**g**), respectively.

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 -1 , -2 , and -3 . 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 (*vrms*) are picked up from the velocity spectrum in Figure 7e. The *vrms* are converted to the *vint* 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 ε*r* of the model. It shows that the estimated *vint* 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 *vint* around 20 ns in Figure 8. The model velocity of the fourth layer is 0.103 m/ns, while the *vint* of that is 0.106 m/ns (Table 1). The failure to pick up the *vrms* of the third layer is part of the reason for the error. After the dielectric constant (<sup>ε</sup>*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.

**Figure 8.** Comparison between the interval velocity and model velocity for the railway subgrade model.


**Table 1.** The detailed information estimated 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 ge<sup>t</sup> 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).

**Figure 9.** A photograph of the experiment site showing the CMP measurement.

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.

**Figure 10.** Original CMP gather (**a**) and the profile after preprocessing (**b**).

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 -1 –-3 , a , and b . Only the energy clusters labeled -1 –-3 have corresponding peaks in the average amplitude curve, which is also marked by -1 –-3 in Figure 11b, while a and b 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 -1 –-3 correspond to relatively complete hyperbolic events, which are meaningful. On the other hand, the clusters a and b , in Figure 12b, are invalid because they produce trajectories without corresponding hyperbolic events. Therefore, the energy clusters -1 –-3 are the reflections from the interfaces of the railway subgrade, and we extract the *vrms* from them to calculate the *vint*. The *vint* is shown in Figure 11c, and the moisture content is calculated as shown in Figure 11d.

**Figure 11.** Velocity analysis of the CMP gather acquired at the first mid-point to obtain the moisture content of the railway subgrade. (**a**) Velocity spectrum with the offset of 0.6–5 m, (**b**) curve of average amplitude, (**c**) interval velocity, and (**d**) moisture content.

**Figure 12.** CMP gather with corresponding hyperbolic trajectory (red lines) which are produced by (**a**) the velocities picked up from the valid energy clusters -1 , -2 , and -3 , and (**b**) velocities picked up from the invalid energy clusters a and b .

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 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.

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


1 The layers here correspond to the marks -1 –-3 in Figure 11a,b.

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 -1 –-4 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 (*vrms*) are converted to the *vint*, 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.

**Figure 13.** Velocity analysis of the CMP gather acquired at the 27th mid-point to obtain the moisture content of the railway subgrade. (**a**) Velocity spectrum with the offset of 0.6–5 m, (**b**) curve of average amplitude, (**c**) interval velocity, and (**d**) moisture content.

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 (*vrms*) are picked up and converted to interval velocities (*vint*) by the Dix formula. Then, the relative dielectric constant (<sup>ε</sup>*r*) can be obtained by Equation (5)—that is, ε*r* = (*c*/*vint*)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.

**Figure 14.** Comparison between moisture content obtained by velocity analysis of CMP (**a**) and polarizability obtained by induced polarization (IP) (**b**). The white boxes are culverts. The black triangles represent the locations of the first and the 27th mid-points of CMP gathers, respectively.

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 a ffects 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 e ffect. 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.
