Numerical Simulation of Polarity Characteristics of Vector Elastic Wave of Advanced Detection in the Roadway

The high-order staggering grid Finite-Difference (FD) scheme based on first-order velocity-stress elastic wave equation has been deduced. The calculation method of PML boundary condition and stability condition established in this study can be used for numerical simulation of advanced detection of elastic wave in roadway, with the obtaining of high-precision seismogram. Then we systematically analyze the polarity of vector wave field in post-source observation system. The results indicate that the relationship between the vector wave field and the polarity of direct wave is related to reflection coefficient on the interface, while the polarity relationship between horizontal and vertical components of vector wave field is related to vertical position of the interface. During data processing for advanced detection of elastic waves, the sign of the reflection coefficient on the interface ahead can be determined based on the polarity relationship between reflected wave and direct wave from the seismograms; the soft and hard rock and other geological information on both sides of the interface is thus be determined. In addition, the direction of source wave depends on polarity relationship between horizontal and vertical components of reflected wave and is used to achieve the separation of up going and down going waves.


INTRODUCTION
Advanced detection of elastic wave is a commonly used method for predicting the abnormal geological structure in the roadway ahead and at the working surface. Compared with ground seismic prospecting, advanced detection is featured by high resolution and good application prospect in safe production of the coal mines. So far, the advanced detection of elastic wave has found extensive application in roadway tunneling (Zeng et al., 2001).
In recent years, researchers have performed many studies on the numerical simulation of advanced detection of elastic wave in mines. Liu (2008) carried out numerical simulation on elastic wave advance detection based on ray theory. Zhu (2008) performed numerical simulation on reverberation characteristics at roadway roof and bottom boundary caused by the earthquake focus at the roadway bottom. Wu (2009) studied the methods for advanced detection synthetic seismogram based on wave theory. Yang (2010) researched the characteristics of diffracted wave at the tectonic interface ahead of the roadway using numerical simulation. Ji (2012) investigated the generation mechanism of channel wave from actual seismogram using 3D numerical simulation and analyzed its frequency dispersion characteristics. The forward modeling described above has been carried out in the purpose of identifying various types of reflected waves generated by excited elastic waves upon encounter with geological abnormal bodies in coal roadway. But much fewer works have dealt with the attributes of the vector wave field and in particular, the polarity characteristics of vector wave field, which is the subject of this study. The post-source observation system is used to acquire the polarity characteristics of vector wave field and the results obtained are of certain guidance and application values for the advanced detection of elastic wave in roadway.

PRINCIPLE OF HIGH-ORDER STAGGERED-GRID NUMERICAL FD METHOD
Elastic wave equation: Under two-dimensional isotropic condition, the first-order velocity-stress elastic wave equation has the following form (Virieux, 1986): where, V x & V z : The horizontal and vertical components of mass point vibration, respectively λ & µ : Lame coefficients, related to the propagation speed of P and S waves in the medium and the density of the medium ρ : Medium density x & z : Represent two directions of spatial coordinate system t : Time σ xx , σ zz & σ xz : The three components of stress Equation (1) could be solved with staggered-grid finitedifference method. Each elastic wave field component, elastic parameter space grid position and time grid position was shown in the Table 1.
Thus, high-order staggered-grid FD method is adopted for difference discretization of Eq. (1). The FD scheme for first-order velocity-stress elastic wave function in isotropic medium with the precision of O (∆t 2 , ∆t 2N ) can be written as: 1 2, 1 2 1 2, 2 1 2 1 2, 2 1 2 1 The rest differential format of V x , σ zz and σ xz could be obtained similarly. C n as differential coefficient, which can be solved by the method described by Mu and Pei (2005).

Absorbing boundary conditions:
In numerical simulation of elastic wave propagation in a limited area, artificial boundary reflection must be considered. The commonly used boundary condition is the absorbing boundary condition based on paraxial approximation theory for scalar wave equation and many new algorithms for boundary condition developed on this basis. The advantages of these algorithms include less calculation and complete absorption of reflected wave at normal incidence. But they are only effective within a certain range of incidence angle and frequency. Berenger (1994a, b) first proposed the PML absorbing boundary condition and applied it to numerical simulation of electromagnetic wave propagation, with quite good results. Collino and Tsogka (2001) applied PML to heterogeneous anisotropic media and achieved good results as well. We have applied PML boundary condition in numerical simulation of advanced detection of elastic wave in roadway, also with good results. The application of PML boundary condition at the edge of computational domain will prevent all reflected waves from entering PML. The seismic waves entering into the PML as a damping medium would decay rapidly. Figure 1 is the diagram of the PML absorbing boundary. The elastic wave field components are decomposed into two parts during numerical simulation and in the case of VX, one along the x direction and the other along z direction. Figure 2 shows the schematic diagram of PML. The attenuation coefficient is introduced on this basis, to derive the wave equation of VX component containing the damping factor: where, d(x) and d(z) are respectively the attenuation factor along x and z direction. On PML: region B, d (x) = 0, d (z) >0; region C, d (x) >0, d (z) = 0; region D, d (x) >0, d (z) >0. All of them are expressed with the attenuation model deduced by Collino and Tsogka (2001): where, R : The theoretical reflection coefficient : The thickness of the PML x : The distance of entering into the PML region of the elastic wave The attenuation coefficient curve along the x direction was shown as below with vp = 3000 m/s, δ = 10 m, R = 10 -8 .
As Fig. 3 that the attenuation coefficient increases rapidly with the increasing distance of wave field entering into PML, exhibiting a "parabolic-like" distribution. This indicates that elastic wave entering into the PML absorbing boundary rapidly decays with the increasing distance due to damping effect in the absorbing layer and ultimately tends towards the ideal boundary reflection coefficient (Collino and Tsogka, 2001). In the case of VX component, the high-order FD scheme within the PML absorbing boundary is expressed as: ( The differential format of the rest of V z , σ xx , σ xz and σ zz in the PML absorbing layer with the accuracy of O (∆t 2 , ∆t 2N ) could be obtained similarly.

Source and stability conditions:
Numerical stability is another concern that must be addressed in FD scheme. According to Lax equivalence theorem, the convergence of difference algorithm also needed to be ensured by stability. The commonly used stability analysis method of FD method is Fourier spectral analysis method proposed by Von Neumann. Grid ratio is the major parameter affecting the stability of algorithm in numerical simulation. The stability condition adopted for FD method by Dong (2000) is: = The spatial order of FD algorithm ∆t = Time step = Staggered-grid FD coefficient In order to ensure computational precision and frequency dispersion, the size of spatial grid is first determined based on the speed and dominant frequency of the source. Next the time step is determined according to the grid ratio specified in by stability condition.
The hidden small-scale structures in the roadway or the working surface are the research objects in advanced detection of elastic wave in roadway. Thus, a higher resolution of seismic exploration is required. Taking the actual conditions of roadway exploration into account, Ricker wavelet with dominant frequency of 300 Hz is used as the excitation source for numerical simulation.

NUMERICAL EXAMPLES
To study the polarity characteristics of vector wave field in advanced detection of elastic wave, four models of single dipping interface are designed as shown in    Fig. 3a the interface is located at the fore upward of roadway and the interface is located at the anterior inferior of roadway in Fig. 3b. Detailed parameters of models are shown in Table 1. In the four models the position of source is located after the geophone array, at (75, 100). The pure P wave source is used to generate 300 Hz Rick sub wave. A total of 20 channels of geophones are arranged, with the geophone placed between (90, 99) and (120,99) in model 1 and model 2 and the geophones are between (90, 101) and (120, 101) in model 3 and model 4. The array length is 30 m, the nearest offset 15 m and the channel spacing 1.5 m. The spatial sampling interval for numerical simulation is ∆x = ∆z = 0.5 m and the temporal sampling interval ∆t = 0.05 ms. The time for seismograph is T = 200 ms. A wave field snapshot is taken every 10 ms. PML absorbing boundary has a thickness of 20 grids. Forward simulation is performed on model, second order in time and eighth order in space. Figure 4 single shot record of lower interface in front of the roadway. The wave field snapshots in Fig. 5 and 6 shows that different types of waves are generated from the energy conversion taking place upon the incidence of elastic waves onto the interface ahead: reflected PP wave, TS wave, transmitted PP wave and converted PS wave. As shown by Fig. 5a and b, when the reflection coefficient of the interface in front of the post-source observation system is less than 0, the PP wave on the wave field snapshots has negative polarity; the upper half of the PS wave has positive polarity and the second half negative. PP wave shows positive polarity; the polarity of PS wave is the same as that of horizontal component. From single-shot records in Fig. 7a and b the direct wave on horizontal component, PP wave and PS wave have negative polarity; the direct wave on vertical component and PP wave have positive polarity. Different from direct wave, PS wave has negative polarity. The polarity characteristics of all kinds of waves in single-shot records for model 2, 3 and 4 are derived in Fig. 8.

CONCLUSION
• The high-order staggered-grid FD method based on the first-order velocity-stress wave equation is deduced, with some discussion on PML absorbing boundary, stability condition of the algorithm and the selection of high-frequency source. After that, forward modeling of seismic wave field in advanced detection of elastic wave in roadway has been carried out and seismograms of high resolution are obtained. • By advanced detection of elastic wave in roadway using post-source observation system, the relationship between reflection coefficient on the interface and polarity of vector wave field is studied. The results indicate: when R<0, the polarity of horizontal components of PP wave and PS wave are the same as that of direct wave; when R>0, the polarity of horizontal components of PP wave and PS wave is opposite to that of direct wave; the polarity of vertical component of PS wave is opposite to that of horizontal component. Thus, the sign of reflection coefficient on interface ahead can be determined by the polarity information of the vector wave field from the seismograms and hence the litho logical information on both sides of the interface in front of the roadway. • The relationship between the location of the interface and the polarity of horizontal and vertical components of vector wave field is also examined.
The results indicate: the polarity of the horizontal and vertical components of reflected PP wave from the interface located in front of and above the roadway is the same; while the polarity of horizontal and vertical components of PS wave is the opposite. The polarity characteristics of vector wave field located in front of and below the interface are altogether different. The reflected wave can be determined as from either the upper or lower interface by the polarity relationship between vertical or the horizontal components of the vector wave field. It is very important for effective separation of up going and down going waves