Abstract
How to obtain dynamic parameters of rock masses quickly and precisely is a popular and difficult problem, which plays a very important part in engineering design or construction. Currently, the methods used to obtain these parameters are in situ testing method, empirical formula, and so on. However, these methods have some shortcomings, such as large investment and long construction period, which cannot obtain the dynamic parameters precisely and quickly in the engineering scale. In this study, a new method of estimating the rock parameters based on the measured field blasting vibration signals is proposed according to theory of elastic stress wave. In addition, an improved method for S-wave identification used in engineering scale was proposed and then the numerical simulation is given to verify the feasibility. Comparison of the numerical identification results and theoretical results clearly show that the improved method is available in S-wave identification with errors less than 2%. By identifying the arrival times of P and S waves, the propagation velocities of P and S waves are calculated and the parameters of rock mass can be obtained at last. Through analyzing the measured field blasting vibration signals in Fengning pumped-storage power station, the dynamic elastic modulus of rock mass inversed by vibration signals is about 2.2~2.9 times of its static elastic modulus, while the inversed dynamic Poisson's ratio is 0.9~0.975 times of the static.
1. Introduction
With the rapid development of society and economy in China, more and more hydropower stations are being built in the southwest and northwest area. Bench blasting is still the main excavation method in hydraulic engineering construction, and this construction method may lead to several safety problems of hydropower projects, such as the deformation and stability problems under the action of blast vibration dynamic load or seismic load [1–3]. The determination of dynamic parameters of rock mass such as dynamic elastic modulus and dynamic Poisson's ratio in large scale is prerequisite of solving these dynamic problems. At the same time, the design organization often provides the parameters of elastostatic, not suitable to dynamical calculation. Therefore, obtaining reliable rock dynamic parameters quickly and efficiently has become a hotspot of study in the field of rock dynamics.
Because of the importance of the rock mass dynamic parameters, various studies on the subject with different approaches have been conducted by many researchers. Some scholars studied the dynamic parameters in different loading rates by laboratory experiment [4]. However, because of the scale effect of rock mass influenced by sample sizes, sample disturbance, and sampling conditions, the value of dynamic parameters obtained by laboratory tests cannot reflect the properties of rock mass. To obtain the parameters accurately, in situ testing was proposed by some scholars. However, there were still many problems because of the restrictions in complex operation, high costs and time-consuming [5–7]. To overcome the above shortcomings, many researchers have introduced ultrasonic wave test to estimate the mechanical parameters of rock mass [8–10]. Jiang et al. [11] applied ultrasonic wave testing to investigate the quality of rock specimens and researched the relationship between static and dynamic parameters. Asef and Najibi [12] were to study the ratio between static and dynamic parameters under different confining pressures. Their study results inferred that it was feasible to predict these parameters based on seismic velocities at different confining pressures. By using the ultrasonic technique, the dynamic elastic modulus of ten carbonate rocks was predicted accurately. A new formula was proposed to calculate static Young’s modulus effectively. On the other hand, studies have shown that mechanical parameters of rock mass could be well described by some empirical equations [13, 14]. Shen et al. [15] deduced the relationships between P wave modulus and mechanical parameters of rock mass by using the field test data. Their studies also pointed that calculating these parameters by using P wave modulus was more accurate and reliable than P wave velocity. Barton [16] deduced the empirical equation to predict deformation modulus through the relationship among the rock mass quality Q-value, P wave velocity and deformation modulus. Song et al. [17] summarized the existed empirical methods for estimating the rock modulus, then discussed the applicable conditions and predicted results. Last, they developed a more practical prediction formula using the data monitored from Maerdang dam project.
The purpose of this paper is to study the method to calculate dynamic parameters of rock mass in engineering scale efficiently. In the following discussion, the relationships between the dynamic parameters and the body wave velocities were explained from the theory of seismic wave propagation. Then the S-wave identification based on field blasting vibration signals is discussed and the results of S-wave identification are validated through numerical simulation. Subsequently, back analysis of rock dynamic parameters is discussed. Finally, the reasonableness of these parameters is analyzed. The present study is of interest to geotechnical engineering due to the results potentially providing guidance for the practical engineering design and construction.
2. Relationship between Velocities of Seismic Waves and Dynamic Parameters
According to the basic theory on transmission of elastic stress wave in continuous medium [18, 19], the equations of particle motion are deduced and can be expressed as
where is the bulk density; and are the Lame constants; , , and are the displacement in , , and three directions; is the volumetric strain; is the Laplacian and this Laplacian can be expressed as .
As shown in (1), the right side of each equation can be divided into three components, which represent the tension or compression deformation, pure shear deformation, and force density contributing to the movement of rigid body respectively. To solve (1), partial differentiation on , , and is performed to the three equations in (1). And by adding these equations, the equation on volumetric strain can be established as follows:
Partial differentiation on and is performed to the last two equations in (1), respectively. Then volumetric strain could be eliminated by subtracting them from each other, and the equation can be given as
Equation (3) can be further rewritten as
where is the angular velocity which is rotating on the X-axis and could be described as .
In the same way, applying the operation curl to both sides of (1), the following are obtained:
where and are the rotational components which are symmetrical to the y and z axis.
Equations (4)-(6) are the vector equations and can be rewritten as
From both (2) and (7), wave motions in the elastic medium have different propagation velocity and vibration mode, which can be divided into two propagation forms: propagation in dilatational disturbance and propagation in rotational disturbance. Equation (2) describes that a dilatational disturbance could be transmitted through the body and the propagation velocity can be obtained by (2), which the value is . The properties of rotational disturbance propagation are depicted by (7). By solving the equation, the velocity of the rotational disturbance propagation can be given and the value is.
For the propagation of dilatational disturbance, the particle vibration direction is consistent with the wave propagation direction, so it is called longitudinal waves. And because of the measured firstly, it is also called primary wave, or P wave for short. For the propagation in rotational disturbance, because of the particle vibration direction perpendicular to wave propagation direction, we call this type of wave transverse wave. On the other hand, the propagation velocity is next only to P wave, so it can be called secondary wave, or S wave in short.
In summary, the expressions of the P and S wave velocities are given as follows:
where and are the dynamic elastic modulus and dynamic Poisson ratio, respectively.
Obviously, the velocities of wave propagation are related to bulk density as well as the elastic parameters of the medium and this is the basis and foundation of the theory of studying the mechanical characteristics of medium by studying the wave propagation in rock mass.
The bulk density is the basic physical parameter for the rock mass. The main reason for the different densities of rock mass is the difference of composition and texture of rock. For the bulk density research, to quantify the density changes, the empirical formulas of the bulk density and P wave velocity have been established by some scholars [21, 22]. Considering the calculation accuracy and convenience, the bulk density of rock mass is calculated by the empirical formulas (10) proposed by Gardner [21]. where is in units of and is in units of .
3. Method of S-Wave Identification
From the preceding discussion, it is clear that to obtain the dynamic parameters of rock mass, the prerequisite is to get the propagation velocities of P and S waves by identifying the onset time of P and S waves. Many scholars of various countries in recent years have already done a lot of research on the field of P and S wave recognition. Fortunately, for P wave, the first arrival time can be easily identified because of the fastest propagation in the rock mass. In this paper, a key assumption is that the P wave first arrival has been identified precisely by other methods and the method proposed by Bear and Kradolfer [23] is recommended. However, due to the affect by the P wave coda, S wave phase identification could be difficult for measured vibration signals. So the S wave identification is mainly discussed in this part.
So far, the methods and techniques on the S wave identification in engineering scales have come directly from studies of the natural earthquakes. In earthquakes, Roberts et al [24] adopt the method of time-shift polarization analysis to identify the S wave phase depending on the different particle motion directions of P and S waves. Gentili et al [25] proposed a method to automatic identifying the P and S wave onset time by using the artificial neural networks and the algorithm of STA/ LTA. In the S wave identification studies, these methods focus on the characteristics of the S waves different from the P waves, such as the linear polarization and the vibration directions [26–28]. Although there are many research results in earthquakes, some problems in S wave identification still exist in engineering scales, such as, the bad separation of P and S waves affected by the anisotropy of rock mass, the shorter duration of blast vibration signals compared to the natural seismic signals, etc. So studying the S wave phase identification in blasting vibration signals is important meaning. Basic thought of this paper is to identify the S wave phase according to the differences in polarization, vibration direction and the energy of P and S wave based on the work of P wave phase identification already done.
3.1. Determination of the Length of Calculation Window
During the analysis of blasting vibration signals, the calculation window’s length is a key parameter and the precise study of it is needed. However, the method is difficult to determinate the window’s length because of the variety of the vibration signals in blasting process. For convenience, the length of window could be calculated through the following proposed by Cichowicz [29]:where and are the sampling rate and the P wave predominant frequency, respectively.
Based on Nyquist theory, data recorded in 2 ms since P wave first arrival are used to calculate and analyze. According to some scholars’ findings [30, 31], the predominant frequency in (11) can be calculated by the following:In (12), and are the displacement power spectrum and the velocity power spectrum, respectively.
In order to identify and analyze the S wave in the blasting vibration signals, we should rotate the measured vibration data coordinate from the ground motion monitoring system defined by the vertical (Z), horizontal (X), and horizontal tangential (Y) directions into the wave particle motion system defined by L, Q, and T, as shown in Figure 1, where the L component represents the direction of P wave particle motion and the Q and T components represent the S wave (SV and SH wave) particle motion directions. The covariance matrix with a set of data in length from the start of P wave arrival of the three components can be calculated as follows:In (13), the covariance between two components is expressed as follows:
Based on Kanasewich’s research [32], the P wave particle motion direction is identical with the eigenvector corresponding to the maximum eigenvalue of the auto covariance matrix. The coordination system in the form of X, Y, and Z could be rotated into the L, Q, and T coordinate system by the three eigenvectors of (13). The rotation expression of the coordinate system is shown as follows:where are the cosine value between the eigenvector corresponding to the th eigenvalue to the X, Y and Z axes.
3.2. Identification of S-Wave Phase
The polarization properties of blast-induced seismic waves can be analyzed by the eigenvalues and corresponding eigenvectors of the covariance matrix. In practice, it is effective in S wave phase identification by using three parameters, such as the deflection angle D(t), the degree of polarization P(t), and the ratio between transversal and total energy related to the covariance matrix E(t). The deflection angle D(t) as the first parameter is defined as the value of angle normalization between the direction of P wave particle motion(L) and the direction of the eigenvector corresponding the largest eigenvalue. The deflection angle D(t) is given by where is the angle between L and the eigenvector associated to the largest eigenvalue at t time with the P wave arrival time as the starting point. Because of the mutually vertical of the polarization directions, P and S waves, the value of D(t) is changed from 0 to 1, which represented P wave arrival and S wave arrival. To calculate the second parameter, polarization degree, the equation proposed by Samson [33] is given aswhere , , and are the eigenvalues of the covariance matrix at t time. The value of P(t) is 1 for the P and S waves first arrival, and the value of P(t) is in the range of 0 to 1 in P wave coda region. The third parameter is the ratio between transversal and total energy related to the covariance matrix E(t) in L, Q, and T coordinate system, which can be calculated by the equation (18) as follows:In L, Q, and T coordinate system, the value of this parameter E(t) can reach 1, which means the first arrival of S waves and reach 0 for the P waves first arrival. The third parameter E(t) can expend the discrepancy between P and S waves by squaring the amplitude of P and S wave and makes the S waves more easily identifiable. Finally, using the characteristic function with three parameters D(t), P(t), and E(t), S waves first arrival can be identified rapidly. The characteristic function can be expressed as
The characteristic function can be used to improve the degree of the linear polarization of vibration signals. In order to make a sharp distinction between the P and S waves, each element of the characteristic function CF(t) is squared. During the propagation of seismic waves, the S wave arrival can make the value of characteristic function CF(t) increase markedly and, in practical application, we can consider the point with the dramatic increase in the value of CF(t) corresponding time as the S wave arrival time in the P wave coda area. In general, it cannot be always identified as the S wave arrivals because of the equivocal time for dramatic increase, as shown in Figure 2. In practice, to make the problem simple, a weight identification function is built by introducing an optimizing weight parameter . As the weight parameter, it can be expressed as So the weight identification function can be described byThe S wave arrival may lead to a rapid rise in characteristic function and weight parameter amplitude. Then the time signal corresponding to the maximum of is the time signal demanded when S wave is arrived. For good visual presentation, the value of is normalized as identification mark.
(a)
(b)
(c)
3.3. Numerical Simulation Experiment
In order to verify the S waves identification method correctness, a numerical study of the single hole blasting was implemented. As the symmetry of the numerical simulation model, only 1/4 of the model was established, as shown in Figure 3. The length, width, and height of this model are 20 m, 10 m, and 9 m, respectively. The total numbers of either nodes and elements of this model are 260012 and 240304, respectively. The depth of borehole is 3 m and at the upper 0.9 m is stemming. The diameter of the hole is 90 mm and the charge diameter is 32 mm. The symmetry boundary was employed in the normal direction of the plane; meanwhile, for the other faces of this model, expect for the free surface, the nonreflecting boundary was applied to prevent the waves reflection at the model surface and the imaginary monitoring points were arranged on the top surface of this model, as shown in Figure 3. In this simulation, the parameters of the rock and explosive are given in Tables 1 and 2 [20], respectively. The S wave identification results on the two monitoring points are given in Figure 4 and Table 3.
(a) R=10 m
(b) R=15 m
Based on the rock parameters and (8) and (9), the P and S wave velocities can be calculated with their respective values: 4248 m/s and 2545 m/s. Then the theoretical value of the S wave arrival time for the monitoring points can be calculated, as shown in Table 3. Comparing the values obtained by the S waves identification method with the theoretical values, the deviation between them was less than 2%. The simulation results indicated that this identification method has high accuracy and can meet the requirements of the engineering.
4. Field Experimental Studies on Back Analysis and Calculation of Rock Dynamic Parameters
4.1. Field Experiment Background and Data Sources
The field blasting tests were carried out in Fengning pumped storage power station to calculate the parameters value of rock mass by monitoring the blasting vibration signals. The Fengning pumped storage power station, with a capacity of 3600 MW, is located in Fengning manchu autonomous county in Hebei, China. The project construction will be divided into two phases. During underground power house blasting excavation, the second-phase construction would cause inevitable effects on the first stage engineering. To guarantee the reliability and security of the underground power house of second stage excavation, it is important to obtain dynamic mechanical parameters of rock mass quickly and efficiently. Taking the excavation blasting of underground powerhouse in the second stage as the background, the site blasting experiments have been conducted. Based on the topographic maps, the experiment sites could be considered flat. The rocks in the test area are mainly intact granite with high uniaxial compression strength and mechanical intensity.
There are arranged six vertical boreholes in the blasting test zone. Hole-by-hole blasting was applied to the site test by using the half-second delay detonator. The layout of blasting design is shown in Figure 5. The No. 2 emulsion explosive was used in the site blasting experiments. The parameters of these boreholes are listed in Table 4.
Vibration monitoring produced by field blasting was carried based on the blasting safety regulation and code for blasting safety monitoring of hydropower and water resources engineering. As shown in Figure 5, there are 6 monitoring points set up along the measuring line on the ground surface. The monitoring instrument system has several components: three axis velocity detector, signal gather recorder equipment, and data processing system. In order to monitor the blasting vibration, the device of TC-4850 is given in Figure 6. The measurement range of amplitude is 0.001~35.4 cm/s, and the measure precision is 0.001 cm/s. Meanwhile the system is of wide detecting vibration frequency band (0~1000 Hz) and an advanced power against interference. Due to the space limited of test area, the distance from the explosion of monitoring point positions was ranged from 10~150 m. In light of Aldas’s work [34], the vibration signals induced by blasting can be recorded by the device of TC-4850 accurately in these measure points. Figure 7 shows the original horizontal tangential velocity recorded by monitoring point (10#).
4.2. Rock Dynamic Parameters Calculation Results and Analysis
By the methods of P and S waves identification proposed above, the field blasting vibration signals carried on geology exploration tunnel in Fengning pumped storage power station were analyzed. The identification results are given in Table 5. In the Table 5, data in bracket are the values of S waves arrival time. Figure 8 shows application of S wave identification method to the typical monitoring point (10#).
(a) vibration time-history of borehole 2
(b) vibration time-history of borehole 3
(c) vibration time-history of borehole 4
(d) vibration time-history of borehole 5
Then putting the data in Table 5 into the (22), the velocities of P and S waves in the region between the blasting source and each monitoring points were obtained and these values are given in Table 6. According to (8) and (9), the P and S waves propagation velocities depend largely on the density, module of elasticity, and Poisson ratio of rockmass, which can be used to reflect the rockmass dynamic parameters. So the parameters such as density, module of elasticity and Poisson ratio of rockmass can be calculated and the calculation results are also given in Table 6. On the other hand, the uniaxial compression experiment and on-site tests were performed to investigate the static mechanics parameters in the experiment zone by Beijing Engineering Corporation Limited, Power Construction Corporation of China, during the design phase. And the suggestion values given are shown in Table 7.where is the distance between the explosion source and monitoring point, in units of m; is the time difference between P wave arrival and S wave arrival, in units of ms.
Based on the contrast of the bulk density from Tables 6 and 7, the difference between the values obtained by back analysis and the suggestive values is relatively good with errors less than 0.15%, which indicates that it is a convenient and effective way to use the P wave velocity to estimate the bulk density. It also provides a good guidance on forecasting the bulk density of rock mass in construction site.
In addition, from the point of view of elastic modulus, the dynamic elastic modulus () obtained by backanalysis from the vibration signals locates in the range of 53 to 70 GPa and the suggestive value of the static elastic modulus () in test zone is 24 GPa. It can be found that is much greater than and ratio between them is calculated as 2.2~2.9. From the Poisson ratio view, the dynamic Poisson ratio() ranges from 0.216 to 0.234 and the suggestive value of static Poisson ratio () is 0.24. So it is easy to see that is slightly less than with the ratio range 0.9 to 0.975. Meanwhile studying the dynamic and static mechanics parameters of rockmass for the Xishan Rock Cliff Statue, Jiang et al. [11] found that the ratio between the and was 2.2~3.2. The results were in line with the results of this paper.
5. Conclusions
By numerical simulation analysis and analyzing the measured vibration signals of Fengning pumped storage power station, the following conclusions can be drawn from this paper:
The S wave identification method in this study is based on the differences in polarization, vibration directions, and energy of P and S waves. The advantages of adopting the weight identification function, , are that they can include some properties of the S wave and can be distinguished from the P wave. In addition, the indicator can enlarge the energy difference between the P and S waves. Then by the numerical simulation, the S wave identification method is effective with errors less than 2%.
Measured vibration signals were used to back analyze the dynamic parameters of rock mass. The results demonstrated that the ratio between the dynamic elastic modulus and the suggestive value of static modulus is range 2.2 to 2.9. The ratio of dynamic and static Poisson ratio is calculated as 0.9~0.975. Meanwhile, the bulk density computed from P wave velocity is basically equal to the recommended value by Beijing Engineering Corporation Limited. These inversion results were consistent with previous studies, which indicated the applicability of the inversion method of measured vibration signals.
It is generally known that the velocities of P and S waves mainly relate to the initial defects of rock mass, including holes, microcracks, and joints. So for this inversion method, the rock mass dynamic parameters obtained by measured vibration signals are mainly about the average properties of rock mass. It is necessary to strengthen the analysis the influence of these defects to the parameters of rock mass in further studies.
The method proposed in this paper possesses many merits, for example, short calculating times, simple operation, wide detecting range, etc. It can provide large-scale mean dynamic mechanical parameters of rock mass in engineering construction. It has provided a new idea and method to calculate the dynamic parameters of rock mass.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work is supported by the National Science Fund for Distinguished Young Scholars (51125037) and the Natural Science Foundation of China (51779190). The authors wish to express their thanks to all supports.