A Co-offset GPR Data Inversion Based on Ray Theory

We use Monte Carlo method to establish a randomly undulating homogeneous layered medium model. Then the co-offset GPR data for the built geological model is simulated by FDTD. The ray-tracing based inversion is applied for interval velocity estimation. During the inversion, the offset between antennas, the velocity value of the first layer, the picked amplitudes values of each reflection layer and reference amplitude are the input data. The thickness and velocity of each layer are calculated by this recursive method. It is proved that this inversion method is feasible for the layered geological model with random undulation, and the fact that undulation, i.e. RMS height and correlation length, influences the inversion results is discussed.


Introduction
GPR is a high-resolution geophysical technology based on the propagation of electromagnetic wave in the frequency range between 1MHz and 3GHz [1]. In seismic research, velocity models are usually obtained by analyzing data collected at multiple source-receiver offsets [2]. However, most commercial GPR systems are equipped with a single receiver antenna, so the acquisition of multiple offset data sets is very demanding [3]. GPR data collection in common offset is simple, fast, and can be used for largescale survey along the gird or line. The obtained data show the change of the physical properties of the underground medium within a certain range. Therefore, we can use the co-offset GPR data to estimate velocity.

Monte Carlo method for stochastic rough surface modeling
The basic idea of Monte Carlo method is to filter it with power spectrum in frequency domain, do inverse fast Fourier transform (IFFT) to get the height fluctuation of rough surface. Figure 1 and Figure 2 show the samples of random rough surface under different RMS height and correlation length. The RMS height and the correlation length related to the scale of the rough surface in the vertical and the horizontal direction respectively.

Forward method based on FDTD
FDTD method differentiates Maxwell equation in time domain and space domain. Using leap-frog---alternating calculation of electric and magnetic fields in space domain, and the change of electromagnetic field is simulated by updating in the time domain to achieve the purpose of numerical calculation [4]. FDTD method is relatively simple in concept, accurate for any complex model [5][6]. This paper uses MATLAB code based on FDTD and PML absorbing boundary provided by Irving and Knight in 2006.

Inversion method based on geometric ray theory
The velocity analysis method we used was proposed by Forte et al. in 2014. We assume that the propagation signal is a plane EM wave. Near each trace location, the underground medium is assumed to be homogeneous, isotropic, non-magnetic ( = 1), non-conductive (σ = 0), and non-dispersible.
The changes of antenna coupling, inherent attenuation and scattering effect are ignored. Normal GPR surveys are performed in transverse electric (TE) broadside configuration. Therefore, we only consider the TE mode (Forte et al., 2014). According to these assumptions, the proposed method requires offset(x), velocity of the first layer medium ( 1 ), amplitude of direct wave ( 1 ), amplitude of reflected wave ( ) and travel time ( ). Given the first layer velocity( 1 ), offset(x), and travel time of the first layer ( 1 ), we can get the thickness of the first layer(ℎ 1 ) by equation (1): Then the angle of incidence is obtained by equation (2): Using the amplitude values picked up ( 1 1 ), the reflection coefficient of the first layer 1 can be obtained by equation (3).
Then Snell equation gives the velocity of GPR signal in the second layer: Where 2 is obtained by rearranging the Fresnel equation of TE mode [1]: If we know the thickness of the first − 1 layers, the GPR signal velocities of the first n layers and the of the nth interface reflected wave, then the thickness of the nth layer (ℎ ) is the only positive solution of the following third-degree equation [1]: In horizontally layered media, the incident angle of the th interface is equal to the transmitted angle of the ( − 1)th interface. Such angle ( ) is related to the horizontal projection (∆ ) of the travel path in the kth layer (Fig.2) through the following equation: Considering the small value of and the Snell's equation, we obtain: According to geometric conditions: Then, we can use the Fresnel equation for TE antenna configuration to calculate the first − 1 reflection and transmission coefficients relative to the nth reflected wave: with = 1,2, ⋯ , − 1.

Figure 2.
Schematic diagram of electromagnetic wave ray path in horizontally layered media (take three-layer medium as an example).S represents the transmitter, R represents the receiver,x is the offset,ℎ and are respectively the thicknesses and EM velocities of the layers, and ∆ are respectively the incident angles and the horizontal projections of the travel path.
Considering the symmetry of the travel path (Fig.2), and using the above coefficients and amplitudes we picked ( 1 ), we can obtain the incident and reflected amplitudes of the nth interface by the follow equation: Then, the reflection coefficient of the th interface is given by: = The velocity in the ( + 1)th layer is given by the Snell's equation: with +1 = Arctan ( By iterating this method for all the interpreted reflections in a GPR trace, we can obtain the thicknesses and velocities of all the imaged layers.

Numerical examples
First, we established three undulating stratum models of size 10m × 5m, with the same correlation length and different RMS heights. The specific parameters are shown in Table 1. Then, numerical simulation is performed on the established model with the following parameters: After obtaining the forward results, we pick up the reflection amplitude, perform the above iterative inversion on all traces of data, and perform linear interpolation on the space between the two tracks. The figures of numerical models and calculated results are as follows:   The top blue layer of the diagram of models is air. The transmitting and receiving antennas are located at 0m on the y-axis.
Due to the position setting of transmitting antenna and the receiving antenna, we can only calculate the dielectric constant of the medium between 0.75m − 8.75m.
It can be seen from above examples that this method can be used for the inversion of co-offset GPR data. It can perfectly distinguish the shape and position of the reflection interface, but as the RMS increases, the calculation error of the permittivity of the medium is also increasing.
Then, we established three undulating stratum models with the same RMS heights and different correlation length. The specific parameters are shown in Table 3. The figures of numerical models and calculated results are as follows:   It can be seen from above examples, as the correlation length increases, the calculation error of the permittivity of the medium is also decreasing.

Conclusions
It is proved that the ray-tracing based inversion is feasible in the layered geological model with random undulation, as the RMS increases, that is, the degree of undulation of the reflection interface increases, the calculation error of the permittivity of the medium is also increasing; as the correlation length increases, that is, the degree of undulation of the reflection interface decreases, the calculation error of the permittivity of the medium is also decreasing.