Inverse Analysis of Rock Mass Dynamic Parameters from Blasting Vibration Signals

The accurate assessment of rock mass mechanical parameters is quite critical for the effective planning and implementation of blasting operations. Inverse analysis methodologies have emerged as a pivotal tool in the computation of dynamic parameters of rock masses in blasting zones, contributing significantly to the refinement of blasting techniques, cost efficiency, and precision in execution.

Presently, the predominant methods for inverse analysis of rock mass mechanical parameters encompass experimental approaches, application of empirical formulas, and ultrasonic testing techniques. Experimental approaches are bifurcated into laboratory tests [1], [2], where rock samples from project sites are examined to ascertain mechanical parameters such as tensile and compressive strength through static experiments. However, this method may not adequately represent the rock structure in blasting areas. Conversely, in-situ tests aim to surmount the limitations of laboratory tests by objectively mirroring the mechanical characteristics of rocks in situ [3], albeit at a higher cost and complexity.

In the realm of theoretical and empirical research, scholars globally have made strides in estimating mechanical parameters of rock masses in blasting zones. Ivars et al. [4] introduced the synthetic rock mass (SRM) technology, facilitating the simulation of jointed rock masses and enabling the creation of models that encapsulate the spatial distribution of rock mass joints for numerical calculation of mechanical parameters. Chun et al. [5] leveraged the rock mass rating (RMR) system to analyze field data for deriving the mechanical parameters of the site-specific rock mass. Furthermore, Zhang [6] formulated equations that establish a correlation between rock quality designation (RQD) and rock mass mechanical parameters, drawing on extensive field data.

Ultrasonic testing methods have been identified as providing rapid inverse calculation of mechanical parameters in blasting areas, with advanced processing of blasting vibration waveforms enhancing signal fidelity. It was found by Huang [7] that acoustic waves, during their propagation, carry substantial information, the alteration of which can reflect the physical and mechanical properties of rock masses. The variational mode decomposition (VMD) algorithm, introduced by Dragomiretskiy and Zosso [8], achieves adaptive signal frequency domain separation, effectively mitigating mode mixing and splitting phenomena, and has gained widespread application in signal decomposition. In pursuit of enhanced prediction accuracy for peak blasting vibration velocity, Li et al. [9] devised a sparrow search algorithm (SSA) optimized back propagation (BP) neural network model, addressing the limitations of traditional BP neural network models in avoiding local optima. Similarly, Tian et al. [10] proposed a signal denoising and smoothing model based on complementary ensemble empirical mode decomposition (CEEMD), establishing an objective function predicated on filtering algorithm similarity and smoothness for optimal solution calculation.

Sarkar et al. [11] conducted measurements of various rock properties, such as density, slake durability index, uniaxial compressive strength (UCS), and P-wave velocity. These measurements underpinned the establishment of a correlation between P-wave velocity and rock mechanical properties, thus facilitating the determination of rock mechanical parameters via ultrasonic testing. Huang et al. [12] explored acoustic emission (AE) effects and failure characteristics of sandstone under varying confining pressures through uniaxial and triaxial compression acoustic emission tests. Their findings indicated a decrease in fractal dimension with increased confining pressure, suggesting a higher orderliness in rock failure. Xue et al. [13] employed a cube-cluster model to effectively characterize the spatial correlation of scattered AE events in granite samples during the rock failure process under triaxial compression conditions.

Further advancements were made by Lawal et al. [14], who proposed predictive models utilizing gene expression programming (GEP), adaptive neuro-fuzzy inference system (ANFIS), and sine cosine algorithm optimized artificial neural network (SCA-ANN) for ground vibrations from blasting in granite quarries. Additionally, Zhou et al. [15] developed a hybrid denoising method combining local mean decomposition (LMD), multiscale fuzzy entropy (MFE), and singular value decomposition (SVD), significantly enhancing denoising efficiency.

Wang et al. [16] employed a novel signal smoothing denoising model, integrating optimal variational mode decomposition (OVMD) with multiscale permutation entropy (MPE). Through the analysis of both simulated and actual engineering signals, it was demonstrated that this model is proficient in accurately identifying signal frequency information. The primary two components were shown to effectively represent the significant elements of the superimposed signal, proving particularly adept for high-precision data sequence analysis, feature extraction, and denoising, with a notable 22.05\% improvement in denoising error over the ensemble empirical mode decomposition (EEMD)-MPE model. Sun et al. [17] introduced the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) algorithm, which incorporates adaptive noise signals into empirical mode decomposition (EMD). This innovation was found to effectively counter the inherent mode mixing issue inherent in EMD, enhancing the time-frequency analysis precision of seismic wave signals resulting from blast-induced ground vibrations. Moreover, Sun et al. [18] proposed a frequency slice wavelet transform-Fourier decomposition method (FSWT-FDM) for blasting vibration signals. This method, which combines FDM with FSWT, reduces ineffective signal components during decomposition, offering new perspectives for analyzing vibration signals in open-pit mining blasting. Liu et al. [19], in their analysis of open-pit slope blasting vibration tests, utilized the principal empirical mode decomposition (PEMD) method to eliminate noise in vibration signals, effectively revealing the frequency spectrum characteristics of slope blasting vibration signals at various elevations.

In summation, the predominant focus of current research is on ultrasonic testing and empirical formulae for the inversion of rock mass mechanical parameters and the processing of blasting vibration signals. However, there remains a paucity of studies that integrate blasting vibration signals with the mechanical properties of mine rocks for parameter inversion. Most existing methods rely on sample collection for mechanical parameter testing, which may not comprehensively represent regional rock mass characteristics. This study, employing initial arrival recognition methods for P-waves and S-waves within vibration signal energy and integrating elastic wave theory, deduces expressions linking P-wave velocity and S-wave velocity to rock mass dynamic parameters. Consequently, it introduces an inverse analysis method for these parameters based on blasting vibration signals. Applied in the Liaoning Province's Guanbaoshan Open-pit Iron Mine, this method has been proven to yield rational and reliable rock mass dynamic parameters, significantly elevating blasting efficiency and reducing associated costs.

The identification of the initial arrival of P-waves in blasting vibration signals is a critical aspect of rock mass dynamic parameter analysis. Commonly implemented techniques, such as the short-term average/long-term average (STA/LTA) method, spectral analysis, and pattern recognition, primarily focus on specific signal aspects. Nonetheless, their efficacy is limited when applied to complex, non-stationary signals like those generated by blasting vibrations. A notable characteristic of these signals is that the energy parameter variations around the P-wave initial arrival are more distinct than amplitude fluctuations, thus providing a more reliable identification marker.

In response to these challenges, a novel approach for P-wave initial arrival recognition, based on the energy of vibration signals, has been developed. This method incorporates the processing of trend items in blasting vibration signals and introduces the envelope derivative operator (EDO). EDO effectively amalgamates frequency domain information with the time-varying envelope of the signal, demonstrating superior capacity in tracking instantaneous changes in the signal. This characteristic renders it particularly effective for accurately demodulating signals, even in environments with substantial noise, thereby making it apt for processing non-stationary blasting vibration signals. The computational steps for the EDO involve the following:

In this equation, $H[x(t)]$ denotes the Hilbert transform of $x(t)$, focusing solely on the instantaneous amplitude without considering the instantaneous frequency. The EDO, akin to the Teager-Kaiser (TK) energy operator method, employs a derivative function as a weighted filter, leveraging the properties of the Fourier transform:

The EDO denoted as $\Gamma x(t)$ for a continuous signal $x(t)$ is defined as:

where, $\dot{x}(t)=\frac{dx(t)}{dt}$, and $H[x(t)]$ represents the Hilbert transform of $x(t)$. For discrete signals $x(m)$, the TK energy operator typically employs a forward difference approach, while the EDO utilizes the central difference method. The discrete formulation of the EDO is presented as follows:

where, $h(m)=H(x(m))h$ represents the $m$-th sampling point within the data window. Subgraph (a) of Figure 1 illustrates the demodulation results of both the TK energy operator and EDO applied to a given signal. Conversely, Subgraph (b) of Figure 1 displays a signal comprising two components, which can be mathematically represented as:

As depicted in Figure 1, both the TK operator and EDO exhibit proficient tracking capabilities for energy changes within the signal. However, at certain extremal points, the TK operator yields negative values due to its modulation term, leading to singular values in subsequent energy separation processes. In contrast, the EDO demonstrates advantageous non-negative characteristics. The total energy of the blasting vibration wave can be quantified as:

where, $x_i$, $y_i$, and $z_i$ denote the amplitudes of the recorded blasting vibrations. The time history curves for the signal components $x_i$, $y_i$, and $z_i$, and the total energy, depicted in Figure 2, reveal a pronounced increase in energy coinciding with the initial arrival of the P-wave.

Adaptations of certain principles from the STA/LTA algorithm [20] were made to accentuate changes in energy, thus aiding the identification of the P-wave's initial arrival. The STA/LTA method involves computing the energy of the vibration within designated long and short calculation windows to pinpoint the P-wave's onset. The formula for the STA/LTA calculation is as follows:

The chosen characteristic function, sensitive to amplitude and frequency changes, is defined as:

where, $X(i)$ denotes the value at moment $i$. To enhance accuracy, a weighting coefficient, based on Eq. (8), is introduced:

This coefficient aptly mirrors the signal’s variations in the frequency domain. From Eqs. (7) to (9), it becomes evident that the characteristic function, essential for identifying the moment of P-wave initial arrival, is computed as:

The analysis reveals that upon the P-wave's signal arrival, there is an abrupt shift in energy, characterized by a sharp increase. Consequently, the moment corresponding to the maximum value of Eq. (10) is determined to be the time of the P-wave's initial arrival.

Following the precise identification of the P-wave's initial arrival, the recognition of S-waves, particularly within an engineering context, often encounters challenges due to the interference of P-wave tails and ambient noise, factors influenced by propagation distance and site conditions. This research advances the understanding of S-wave recognition by comparing the physical properties of P-waves and S-waves to ascertain the initial arrival time of S-waves in measured blasting vibration waves. A critical factor in analyzing P-waves and S-waves is the calculation window width. The relationship between this window width and the dominant frequency of P-waves is generally expressed as:

where, $l$ denotes the window width length, $f_p$ represents the dominant frequency of P-waves, and $\Delta t$ denotes the instrument's sampling time. The calculation of the angle between the polarization directions of S-waves and P-waves is determined by the following formula:

where, $M$ symbolizes the covariance matrix of P-waves over the length $l$, with the angle $\alpha$ between its maximum eigenvalue and the polarization direction of P-waves as the deflection angle. The normalization of $\alpha$ is achieved through:

Furthermore, the formula for calculating the polarization degree of P-waves and S-waves is given as:

where, $\lambda_1, \lambda_2$ and $\lambda_3$ represent the eigenvalues of $M$ corresponding to each respective window. The transformation of the original coordinate system of the vibration signal, specifically radial $(X)$, tangential $(Y)$, and vertical $(Z)$ directions, into the $LQT$ coordinate system is executed using the following formula:

where, $u_{i j}\,\,(j=1,2,3)$ denotes the cosine of the angle between the $i$-th principal direction and the radial $(X)$, tangential $(Y)$, and vertical $(Z)$ directions. The ratio of the transverse energy of P-waves and S-waves to the total energy is calculated as per the subsequent formula:

The method for identifying the initial arrival of S-waves is expressed through:

where, $\sqrt{{{Q}^{2}}+{{T}^{2}}}$ represents the transverse energy. Analysis of prior research, alongside Eqs. (11) to (17), leads to the conclusion that the moment yielding the maximum value of Eq. (17) is designated as the initial arrival time of the S-wave. The normalization process for Eq. (17) is formulated as:

The propagation of blasting vibration waves through rock mass media imparts key insights into the dynamic parameters of the rock, such as vibration velocity, frequency, duration, and waveform. These parameters are intricately influenced by the rock's mechanical properties. Additionally, the propagation process itself carries vital physical information pertinent to these mechanical characteristics. An analysis of the measured blasting vibration signals further elucidates the rock mass's elastic modulus, Poisson's ratio, and density. According to elastic wave theory, the correlation between P-wave and S-wave velocities and the dynamic parameters of the rock mass is established as follows:

This equation encompasses longitudinal wave velocity $V_P$ (km/s), transverse wave velocity $V_S$ (km/s), the dynamic elastic modulus $E_d$ (GPa) of the rock mass, dynamic Poisson's ratio ${\mu }_{\text{d}}$, and the density $\rho$ (g/cm³) of the rock mass.

Eq. (19) demonstrates that the dynamic parameters of the rock mass can be deduced from the propagation velocities of blasting vibration waves, as articulated in the subsequent formulas:

Rock mass density is recognized as a fundamental physical property. Extensive research by scholars globally [21], [22], [23] has delved into the correlation between density and P-wave velocity. It has been observed that different types of rock masses exhibit diverse functional relationships between density and P-wave velocity, attributable to their distinct properties. Commonly, Gardner's empirical formula, introduced in 1974, is employed in these assessments:

The calculation of P-wave and S-wave velocities within the rock mass at various measuring points and the blast source is formulated as:

where, l signifies the distance (m) between the blast source and the measuring point, and $\Delta t$ represents the time difference $(ms)$ between the initial arrivals of P-waves and S-waves.

From Eqs. (19) to (23), it is inferred that in practical applications, the analysis of the waveform of the measured blasting vibration signal at the measuring points, combined with the recognition of the initial arrivals of P-waves and S-waves, and the computation of P-wave and S-wave velocities between the points, facilitates the determination of the rock mass's elastic modulus, Poisson's ratio, and density in the area between the points. The inversion process of these parameters is delineated in Figure 3.

To authenticate the precision of the method employed for identifying the initial arrival times of P-waves and S-waves based on vibration signal energy, the LS-DYNA dynamic analysis program was utilized to simulate the propagation of blasting vibrations in a rock mass within a predetermined area. The rock mass was represented using 3D Solid 164 elements, featuring dimensions of 20m×10m×15m, a borehole depth of 4m, stemming of 1m, a borehole diameter of 1m, and a charge diameter of 0.8m. Vibration monitoring points were strategically positioned above the model, as depicted in Figure 4.

Figure 5 demonstrates the efficacy of P-wave and S-wave recognition in the blasting vibration waveform at a distance of 15m from the blast center.

The rock's physical and mechanical parameters, essential for the analysis, are summarized in Table 1. Utilizing Eq. (23), the P-wave velocity was computed to be 5113m/s, and the S-wave velocity was determined to be 3064m/s.

The proposed method was applied to identify the initial arrivals of P-waves and S-waves. These results were then juxtaposed with the theoretical initial arrival times. Table 2 details the accurately calculated initial arrival times of P-waves and S-waves. Upon comparison with the theoretical initial arrival times, it was found that the differences in the identified initial arrival times of P-waves and S-waves were nominal, with errors remaining below 3%. This outcome validates the efficacy of the recognition method and its compatibility with the requirements of field engineering.

A field experiment was conducted at the 285m level in the Guanbaoshan project, Liaoning Province. Area A, the focus of the blasting operation, was a rock region comprising 81 blastholes. The holes were spaced 6.0 m apart, with a row spacing of 4.2 m. A delay time of 17 ms was implemented, with hole depths varying between 14.0 m and 15.6 m. The maximum charge per hole was set at 202 kg, resulting in a total designed charge of 15044 kg. The objective of the study was to deduce the challenges associated with blasting in the subsequent area B, through the analysis of vibration conditions in area A. Eight measurement points were strategically placed on-site, distanced 8-10 m apart, along the central line of blasting area B, as illustrated in Figure 6.

For the experiment, the Topu integrated iSV-420 and the TC-4850 blasting vibrometer, developed by the Chinese Academy of Sciences, were employed. These instruments are capable of measuring three-dimensional velocity, acceleration, and other vibration parameters, making them apt for the on-site measurement of engineering blasting vibrations. Figure 7 shows the layout of the measurement points, while Figure 8 displays the arrangement of the on-site vibrometers.

The signal measured at point 1\# during the experiment is presented in Figure 9. The blasting vibration waveform exhibited a baseline deviation, necessitating correction before further analysis. The adjusted blasting vibration signal is depicted in Figure 10.

The methods outlined earlier for identifying the initial arrivals of P-waves and S-waves were applied to the vibration waveform data collected. Table 3 displays the identification results for the initial arrival times of P-waves and S-waves at each measuring point. Additionally, this paper presents the characteristic function curves for the recognition of P-waves and S-waves in the blasting vibration waveform at a representative measuring point, point 4\#, as illustrated in Figure 11.

Upon substituting the initial arrival times of P-waves and S-waves, as identified in Table 3, into Eq. (23), the average velocities of P-waves and S-waves between various measuring points in blasting area B were calculated. The results of these calculations are summarized in Table 4.

To corroborate the reliability and accuracy of the inversion method, ultrasonic testing was performed in blast zone B, with an ultrasonic gauge (Figure 12) measuring the velocities of P-waves and S-waves. The results of these ultrasonic measurements are summarized in Table 5.

The P-wave velocity in blast area B varied from 4.856 km/s to 5.442 km/s, and the S-wave velocity ranged from 2.927 km/s to 3.273 km/s, as indicated in Table 4. Ultrasonic measurements (Table 5) showed P-wave velocities varying from 5.075 km/s to 5.639 km/s, and S-wave velocities from 3.081 km/s to 3.443 km/s. The average errors between the P-wave and S-wave velocities obtained from the measured blasting vibration signal waveforms and those measured by ultrasonics were 4.67% and 4.74%, respectively, as shown in Table 6. This level of accuracy confirms that the method is suitable for practical engineering applications.

The velocities of P-waves and S-waves derived from Table 4 were incorporated into Eqs. (20) to (22) for the calculation of the rock mass's mechanical parameters in the experimental area: dynamic elastic modulus, dynamic Poisson's ratio, and density. Table 7 encapsulates these results. Furthermore, Table 8 presents the static parameters of the rock mass recommended by the Anshan Iron and Steel Group Mining Design Institute for blast area B.

Figure 13 depicts a comparison of the rock mass parameters. The dynamic Poisson's ratios, obtained from the inversion of blasting vibration waveforms, varied from 0.186 to 0.226, compared to the institute's provided ratio of 0.22, resulting in a comparative ratio of 0.85 to 1.03. The density derived from the waveform inversion, relative to the density value provided by the Anshan Iron and Steel Group Mining Design Institute, exhibited a comparative ratio of 0.98 to 1.01. The dynamic elastic modulus of the rock mass, ascertained from waveform inversion, ranged from 53 to 69 GPa, in contrast to the institute's recommended static elastic modulus of 26 GPa, yielding a comparative ratio of 2.0 to 2.6. These inversion results fall within the anticipated dynamic-static modulus ratio range of 2.2 to 2.9, substantiating the viability of deducing rock mass mechanical parameters from blasting vibration signal waveforms.

This study has addressed the gap in studies focusing on the inversion of mechanical parameters of rock masses within blasting areas using field blasting vibration data. Previously, methods predominantly revolved around sample collection for mechanical parameter testing without an in-depth analysis of the rock mass properties in the measured area. This paper has investigated the disparities in frequency, deflection angle, polarization, and energy between P-waves and S-waves in blasting vibration signals. Expressions for recognizing the initial arrivals of P-waves and S-waves have been derived, and based on elastic wave theory, the dynamic interplay between P-wave velocity, S-wave velocity, and the dynamic parameters of the rock mass has been analyzed. This approach offers an improved representation of the variations in frequency, energy, and deflection angle around the initial arrival time of S-waves. The key findings and prospects of this study are as follows:

Firstly, the average relative error between the P-wave and S-wave velocities, as recognized from the blasting vibration signal waveforms and those measured by ultrasonics in the field, is 4.67\% and 4.74\% respectively. These values fall within the acceptable range for engineering applications, underscoring the method's practical feasibility.

Secondly, by incorporating elastic wave theory, a relationship has been established between P-wave and S-wave velocities from blasting vibrations and the dynamic parameters of the rock mass. The proposed method, premised on the recognition of initial arrivals of P-waves and S-waves, inverts the dynamic parameters of rock from blasting vibration signals. Engineering practice demonstrates that the inverted density values approximate those provided by the research institute, with a comparative ratio of 0.98 to 1.01; the dynamic-static Poisson's ratio comparative ratio is 0.85 to 1.03; and the dynamic-static modulus of elasticity comparative ratio is 2.0 to 2.6. These inversion results align with prior research, validating the reliability and applicability of the dynamic parameters of the rock mass obtained through this method.

Finally, the influence of surface waves on the accuracy of P-wave and S-wave recognition was not examined in this paper, potentially affecting the results. The signal pickup process may inadvertently capture noise, and classifying the initial motion of P-waves and S-waves could be subject to inaccuracies. Additionally, the assumption that the blasting area represents an isotropic continuous medium, without considering the presence of joints, fractures, and other conditions within the rock mass, may introduce errors in the study of inverting the dynamic parameters using P-wave and S-wave velocities.