Digital Processing of Seismic Data from Open-Pit Mining Blasts

This article describes an approach of mathematical processing of signals (seismograms) from five blasthole charges from experimental blasting, each 3 m deep, with equal explosive weight (1 kg), and equidistant (3 m) from one other. The seismic explosive waves were measured at a 13 to 25 m distance. This article provides spectral analysis, wavelet analysis, and fractal analysis results. It defines the dependence of dominant frequency and amplitude on the distance to the blast center. According to the experimental data, the dominant frequency is calculated as y = 1.0262x0.2622 and the amplitude dependency as y = 18.139x−2.276. Furthermore, the analysis shows that 80% of the entire signal is concentrated in half the area of frequency range, i.e., the low frequency zone is of the most interest. This research defines the dependence of distance on the energy value of signal wavelet analysis. It is demonstrated that, according to the experimental data, the 12th frequency range is closely correlated with the distance values. This article gives the definitions of entropy, correlation dimension, and predictability time. This experiment shows that entropy and correlation dimension decrease but predictability time increases when the distance to the blast center increases. This article also describes the method for determining optimal drilling and blasting parameters, and concludes with the possibility of applying the analytical results to predicting and enhancing drilling and blasting operations.


Introduction
Determining the parameters of drilling and blasting operations is currently a relevant and very important task [1][2][3]. The physical nature of blast waves, namely, their seismic characteristics, is especially interesting for drilling and blasting operations. Blast waves, in the form of elastic waves, propagate through the rock in all directions with decreasing amplitude. If the vibration intensity remains high enough, nearby objects could be seriously damaged [4][5][6]. Practically, for the determination of drilling and blasting parameters, few tasks should be needed. Direct tasks include the determination of seismic characteristics of blast waves (amplitude, velocity) in order to estimate the impact on nearby objects (possible damages estimation). Inverse tasks, such as blast parameter selection, can achieve particular amplitude at a certain point, and factor estimation can affect the blast wave amplitude [7,8]. The mined rock blasting control is quite complicated, since it is difficult to determine the blast parameters and rock properties impacting the operations. These values are very often considered to be disturbances, as their effect is changeable and sometimes unpredictable. There are several approaches for predicting drilling and blasting parameters-one of them is modeling [9,10], which includes a few methods. The first method, the simplest and most popular, is the empirically fitted equation linking the basic target parameters and various factors affecting these parameters, such as rock properties, explosive weight, well depth, etc. The process of this dependence determination is quite diverse and varies from searching for the equation in the process of statistical data processing to strictly formalized formulas that apply variable empirical coefficients. An example of the latter is M.A. Sadovsky's approach [11], which is applicable in Russia. The rest of the world uses approaches formulated by USBM and others [12][13][14][15]. Another more modern approach involving ECM (Error Correction Model) focuses on the application of forecast models [16][17][18]. These models are based on artificial neural networks (ANNs) and appear to be more precise, easily adjustable, and more predictable, and they avoid strictly formalized mathematical relations. These models consider factors and relationships, the influence of which could be missed in other approaches. These ANN-based models have a range of features and disadvantages. Regarding the features, forecasting only after a certain period of time showed that the stored amount of data is enough for the training and validation of the neural network. While the ANN is being trained, the object remains uncontrollable, without any possibility to obtain the model data. So, practically, it is easier to use empirical equations. The models based on empirical data are especially relevant, but this statement is true only for a certain forecast horizon. Usually, the forecast horizon is determined strictly by the area where the empirical equation coefficients were found. In this regard, it is a relevant task to search for mathematical methods enabling us to find the empirical equation coefficients applicable to blasting process forecast, and to expand the forecast horizon within a short period of time without loss of accuracy.
This paper describes how the mathematical analysis of seismograms helps to find the empirical equations and dependences used for calculations. A principal hypothesis of this research is that the application of mathematical methods for signal processing to seismograms can provide detailed information on drilling and blasting operations, applicable to the determination of empirical coefficients in characteristic equations within the time period but not exceeding the system computation cycle. This hypothesis is formulated considering a certain level of drilling and blasting automation, when seismogram values are automatically entered into a control system, and the mathematical methods described below are used for data processing algorithms that aim to optimize the operations and provide real-time recommendations on effective parameters of drilling and blasting operations. Therefore, the "system computation cycle" means the period, defined as the real time in the automation system. Spectral, wavelet, and fractal methods of analysis are selected as mathematical methods in this research.
The purpose of this research is to define mathematical methods enabling prompt seismic blast wave amplitude prediction. In order to forecast the propagation of blast waves, it is necessary to consider two principal values: peak particle velocity (PPV) and frequency. Thus, the principal target of research devoted to the blast wave effect forecast involves the determination of these values [5,6,14]. In most cases, the values are predicted based on two or three blasting and drilling parameters (mainly charge mass and well depth), and other parameters (soil characteristics, weather conditions, etc.) are not considered.
However, there are a few research papers where the PPV is forecasted based on a large scope of data, for example, the work in [19], which involves 150 parameters based on ANN. Moreover, some research has studied the application of fuzzy logic and genetic algorithms. The most important examples of such research are listed below.
The classification and regression tree (CART) method for PPV prediction is proposed in [20]. It uses 51 datasets as the initial ones, and maximum charge per delay and distance from the blast face as the input variables.
A dependence relationship between the simple fractal dimension and complexity and the blast wave frequency structure is discussed in [21]. It describes the relationships between the blast wave parameters and multifractal parameters.
The process of devising a model for PPV forecasting is described in [22]. It describes an empirical equation and ANN, enabling the authors to forecast the explosion-induced PPV. The results are based on the datasets of 95 values collected from a granite quarry site in Malaysia (95 blasting works were precisely monitored in a granite quarry site in Malaysia).
A forecast model for soil vibration assessment based on gene expression programming (GEP) is presented in [23]. The results are based on the datasets of 102 values collected from a granite quarry in Malaysia.
The process of ANN modeling for PPV prediction using five input parameters, namely, load, interval, maximum charge per deceleration stage, distance from the face to the observation point, and rock properties (burden, spacing, maximum charge per delay, distance from blast face to monitoring point and rock quality designation) is described in [24]. The results are based on the datasets of 112 values collected from limestone mines in Iran.
The wavelet analysis application for bearings troubleshooting is discussed in [25]. The analysis shows how the signal is decomposed into different frequencies which are analyzed in order to identify various faults and assess the severity level.
The spectra analysis, energy wavelet and fractal characteristics of the blast wave obtained experimentally at different distances from the explosion center are discussed in [26]. This research reveals the dependence of dominant frequency, wavelet transform energy, and their nature in low frequency and high frequency areas, as well as correlation dimension on the distance to the explosion center.

Method
Blast waves have different characteristics that can be extracted using Fourier transform [27], wavelet transform, and fractal analysis [28][29][30]. The wavelet wave signal decomposition is used to derive the wavelet coefficients. These coefficients are linked with blast wave characteristics. S(t) decomposition is computed as [31]: where f i,j t j is the reconstructed signal of the wavelet packet at the point (i, j), j = 1, 2, . . . , 2 i − 1 (i = 1, 2, 3, 4, 5). According to wavelet theory [32], the primary function for spectral analysis of wavelet packet energy should be selected. This was the 'dB8 wavelet. The sampling frequency of the equipment is 1024 Hz, so the corresponding Nyquist frequency is 512 Hz. The wave is decomposed into five layers, with S 5,j energy denoted by E 5,j : where x j,k (j = 0, 1, 2, . . . , 2 5 − 1; k = 1, 2, . . . , m) is the amplitude of discrete points of the reconstructed signal S 5,j and m is the number of discrete sampling points. The energy percentage in each frequency band is computed as: The wavelet analysis results were estimated by means of correlation analysis. Within the fractal analysis, it is necessary to determine several characteristics, including Delay period, Delay time for the gap, Embedding dimension, Entropy, Predictability period, and Correlation dimensionality [33]. For the first time, the series analysis was described by Takens [34]. Based on this approach, it is possible to predict the system dynamics by its time series if the following hold: (1) The system is indeed dynamic (has a finite degree of freedom) (2) The forecast horizon of the system is clear (by the entropy of the time series) (3) The forecast accuracy can be determined.
The first characteristic of the dynamical system description is the possibility of finitedimensional description, namely, the determination of factors influencing the system dynamics. For this, the embedding value is determined, that is, the minimum number of dynamic variables that uniquely describes the system behavior. The C-C algorithm was chosen [35] for delay time calculation.
According to the C-C algorithm, the correlation integral for the embedded time series is computed as: where H is the Heaviside step function, N is the size of the dataset, t is the index lag, and M = N-(m-1)t is the number of embedded points in m-dimensional space. The mathematical sign of . . . represents the sup-norm, and C(m, N, r, t) measures the fraction of the pairs of points { x i } , i = 1,2, . . . , M, whose sup-norm separation is no greater than the radius of neighborhood r.
If the stochastic process {x i } is independent identically distributed, so The Brock-Dechert-Scheinkman (BDS) statistic corresponding to equation According to [26], the best embedding dimension m is computed as: The entropy was calculated using the following standard formula: where q is the attractor dimensional, and p is the probability of finding an element of the attractor value in ith hypercube with edge ε covering part of the n-dimensional Euclidean space.

Experiments
To prove the stated hypothesis, the experiment was carried out. Blast wave seismograms were obtained and analyzed using the signal processing methods described above.
The experiment included five trial blasts with the same amount of explosive (1 kg), simulating major rock blasting. The experimental blasts (single) had stable parameters (charge mass, rock characteristics, etc.) and only one was different: the distance to the blast center (from 13 to 25 m). The experiments were carried out near the existing gas pipeline. Interestingly, the distance from the explosive charges to the gas pipeline was 18 m. Explosive mass was 1 kg Nitronit PAS-60 (a 0.7 kg TNT (Trinitrotoluene equivalent)). The exploded rocks were composed of granite with a rock hardness ratio of up to f = 12 according to the scale of Protodyakonov, the strength value at natural humidity ranged from 86 MPa to 119 MPa, and bulk density was 2.7 g/cm 3 .
The experiment scheme is illustrated in Figure 1. For the experimental measurements, the Blastmate III seismic station was used. The farthest charge was exploded first, followed by all the others. The blasts were carried out in a granite massif with Nitronit PAS-60 (ПAC-60) explosives. For more exact results, the soil was excavated (1 m depth) up to the granite rocks. The well depth was 3 m. To exclude the possible scattering of the blasted rock mass, the blasting area was covered with a sand layer (1.5 m).
The experiment pictures are presented in Figure 2.

Results and Discussion
The experiment results are reflected by seismograms showing the ground flow velocity change with respect to three coordinates-radial, tangential and vertical. Three components of the ground velocity at a distance of 13 m are shown in Figure 3. For the experimental measurements, the Blastmate III seismic station was used. The farthest charge was exploded first, followed by all the others. The blasts were carried out in a granite massif with Nitronit PAS-60 (ΠAC-60) explosives. For more exact results, the soil was excavated (1 m depth) up to the granite rocks. The well depth was 3 m. To exclude the possible scattering of the blasted rock mass, the blasting area was covered with a sand layer (1.5 m).
The experiment pictures are presented in Figure 2. For the experimental measurements, the Blastmate III seismic station was used. The farthest charge was exploded first, followed by all the others. The blasts were carried out in a granite massif with Nitronit PAS-60 (ПAC-60) explosives. For more exact results, the soil was excavated (1 m depth) up to the granite rocks. The well depth was 3 m. To exclude the possible scattering of the blasted rock mass, the blasting area was covered with a sand layer (1.5 m).
The experiment pictures are presented in Figure 2.

Results and Discussion
The experiment results are reflected by seismograms showing the ground flow velocity change with respect to three coordinates-radial, tangential and vertical. Three components of the ground velocity at a distance of 13 m are shown in Figure 3.

Results and Discussion
The experiment results are reflected by seismograms showing the ground flow velocity change with respect to three coordinates-radial, tangential and vertical. Three components of the ground velocity at a distance of 13 m are shown in Figure 3. The seismograms were processed in several ways, as described below. The Fourier transform provides the spectra as shown in Figure 4. Table 1     The seismograms were processed in several ways, as described below.
The Fourier transform provides the spectra as shown in Figure 4. Table 1 presents the dominant frequency and amplitude for every component of the signal.  The seismograms were processed in several ways, as described below.
The Fourier transform provides the spectra as shown in Figure 4. Table 1 presents the dominant frequency and amplitude for every component of the signal.     [36,37]. High frequency vibrations generally have no effect and should be classified just as white noise. The dependence of the dominant frequency on the distance is described in Figure 5.
From Table 1, it is seen that at low frequencies (up to 30 Hz) the dominant frequency occurs only for Long signal (radial component of velocity). For blast wave oscillation analysis, it is only necessary to consider the low frequency oscillations [36,37]. High frequency vibrations generally have no effect and should be classified just as white noise. The dependence of the dominant frequency on the distance is described in Figure 5. The dependence of the amplitude on the distance is described in Figure 6. As seen, the dominant frequency increases as the distance to the explosion center increases, as a power function, which can be described by the following equation: y = 1.0262x 0.2622 . The amplitude decreases as the distance to the explosion center increases, as a power function, which can be described by the equation: y = 18.139x −2.276 . Thus, it can be argued that there is a power-law dependence between the dominant frequency and the distance to the explosion center, as well as between the amplitude and the distance to the explosion center. It is worth noting that the dependence of dominant frequency increasing Distance to the blast center, м The dependence of the amplitude on the distance is described in Figure 6.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 7 of 14 From Table 1, it is seen that at low frequencies (up to 30 Hz) the dominant frequency occurs only for Long signal (radial component of velocity). For blast wave oscillation analysis, it is only necessary to consider the low frequency oscillations [36,37]. High frequency vibrations generally have no effect and should be classified just as white noise. The dependence of the dominant frequency on the distance is described in Figure 5. The dependence of the amplitude on the distance is described in Figure 6. As seen, the dominant frequency increases as the distance to the explosion center increases, as a power function, which can be described by the following equation: y = 1.0262x 0.2622 . The amplitude decreases as the distance to the explosion center increases, as a power function, which can be described by the equation: y = 18.139x −2.276 . Thus, it can be argued that there is a power-law dependence between the dominant frequency and the distance to the explosion center, as well as between the amplitude and the distance to the explosion center. It is worth noting that the dependence of dominant frequency increasing Distance to the blast center, м As seen, the dominant frequency increases as the distance to the explosion center increases, as a power function, which can be described by the following equation: y = 1.0262x 0.2622 . The amplitude decreases as the distance to the explosion center increases, as a power function, which can be described by the equation: y = 18.139x −2.276 . Thus, it can be argued that there is a power-law dependence between the dominant frequency and the distance to the explosion center, as well as between the amplitude and the distance to the explosion center. It is worth noting that the dependence of dominant frequency increasing with distance contradicts the results in [26], which states that the dominant frequency decreases with distance. This contradiction is due to some differences in the experiment procedure. The experiments described in [26] were carried out underground with long distance (from 30 to 110 m) measurements and heterogeneous rock. The current experimental blasts were carried out on the surface of hard rock with preliminary soft soils excavation and short distance (up to 30 m) measurements close to the existing gas pipeline. The dominant increase can be explained by the soil vibration heterogeneity, as well as the close gas pipeline vibration impact. We assume the following reasons for the research results. It is possible that vibrational amplitude at higher frequencies could be underestimated due to the influence of the gas pipeline natural vibration occurring at the same frequencies, but with phase mismatch. At the same time, when the oscillation phases coincide, the amplitude rate may, however, increase. In all cases, any source of vibration close to the measurement point influences the resulting amplitude rate. The research results could be more comprehensive if, in addition to the measurements at certain points, the instruments were also applied directly on the pipeline to record its natural oscillations and other parameters showing the resonance clearly. However, the experiment lacked the opportunity to do this, and currently a detailed explanation of the results does not seem possible. It is worth mentioning that basic seismic monitoring approaches do not involve natural vibration data from nearby objects, regardless of their effect on total vibration. For example, in [26,29] the measurements were carried out the same way, but the instruments were located at selected points without direct recording of any nearby object's natural vibration. However, additional research on the dominant frequency changes when possible resonance is considered, which may improve the accuracy of the drilling and blasting parameters. The natural vibration of every nearby object can cause a positive or negative effect on the dominant frequency, depending on the phase term. There is no doubt that the resonance probability is quite low, but this requires additional assessment and should be considered in drilling and blasting parameter variations. To sum up, according to the authors' opinion, for enhanced seismic safety of the guarded objects, it is important to estimate the natural vibration of the guarded objects and the vibration of the objects located close to the explosion site. The findings need detailed investigations and perhaps additional experiments with the simultaneous measurement of close objects' natural vibration and resonance survey. Furthermore, the obtained results prove the practical relevance of preliminary experiments, since pretested conditions of drilling and blasting operations are helpful for parameter optimization.
Using this dependence, it is possible to calculate the optimal drilling and blasting parameters for any specific task. For example, using this dependence and Sadovsky's approach [11], it is possible to calculate the required mass of an explosive substance for a certain dominant frequency.
The result of the signal wavelet analysis is illustrated in Figure 7. As one can see in the graph, the percentage distributions of energy at different distances are closely related to each other. To confirm this dependence, it was necessary to estimate the energy distribution correlation (in percent) based on the signal wavelet analysis obtained at different distances from the explosion center. This wavelet analysis shows that, on average,~80% of the entire signal is concentrated in the area up to the 16th frequency range (half of the frequency spectrum). That is because the low frequency zone is of the most interest. The range of 80 percent is suggestive of the Pareto principle [38]. The correlation matrix for the energy distribution (in percent) based on the signal wavelet analysis at different distances from the explosion center is shown in Table 2.  The correlation matrix for the energy distribution (in percent) based on the signal wavelet analysis at different distances from the explosion center is shown in Table 2.  Table 2 it is seen that the correlation ratio between the signals is very high. This result can be interpreted as proof of the experiment's validity. This tight correlation indicates the data's homogeneity, which means that all the necessary terms of the experiment are complied with.
Within the experiment, it was hypothesized that each parameter that influences blast wave propagation can be detected at a certain signal frequency range. So, the signal decomposition into its frequency components can show how energy values correlate to other parameters influencing the blast wave nature at every frequency range (for example, strong correlation of energy and explosive mass values at the fifth range probably means that for optimal explosive mass determination it is necessary to analyze the signal frequencies of the fifth range, etc.). To prove this hypothesis, the research shows one more correlation assessment of the distance to the blast site and the energy values for each frequency range. This becomes possible since all the other parameters (substance mass, rock nature) influencing the blast wave propagation used in all five experiments were identical. The assessment of the results is presented in the form of a graph in Figure 8. From Table 2 it is seen that the correlation ratio between the signals is very high. This result can be interpreted as proof of the experiment's validity. This tight correlation indicates the data's homogeneity, which means that all the necessary terms of the experiment are complied with.
Within the experiment, it was hypothesized that each parameter that influences blast wave propagation can be detected at a certain signal frequency range. So, the signal decomposition into its frequency components can show how energy values correlate to other parameters influencing the blast wave nature at every frequency range (for example, strong correlation of energy and explosive mass values at the fifth range probably means that for optimal explosive mass determination it is necessary to analyze the signal frequencies of the fifth range, etc.). To prove this hypothesis, the research shows one more correlation assessment of the distance to the blast site and the energy values for each frequency range. This becomes possible since all the other parameters (substance mass, rock nature) influencing the blast wave propagation used in all five experiments were identical. The assessment of the results is presented in the form of a graph in Figure 8. As illustrated, the correlation coefficient is more than ± 0.9, namely, −0.96 is estimated on the 12th frequency range. It can be assumed that, in order to predict the distance values, it is necessary to analyze the value of the wavelet energy transform at the 12th frequency range. Additionally, a close correlation is observed on the 2, 4-8, and 11-16 ranges. This evidence proves that, for distance estimation, it is necessary to analyze several frequency As illustrated, the correlation coefficient is more than ± 0.9, namely, −0.96 is estimated on the 12th frequency range. It can be assumed that, in order to predict the distance values, it is necessary to analyze the value of the wavelet energy transform at the 12th frequency range. Additionally, a close correlation is observed on the 2, 4-8, and 11-16 ranges. This evidence proves that, for distance estimation, it is necessary to analyze several frequency ranges at once. However, to verify this assumption, we need additional experimental data.
The fractal analysis provided additional information on the seismograms [39]. Primarily, it showed the chaotic system delay time at different distances, as well as the embedding dimension determined from these data. It was also possible to determine a correlation dimension, entropy and predictability time for each signal. The results are presented in Table 3. The final values of the Entropy exponent show the chaotic nature of the system under consideration and exclude its randomness (since the exponent is a final number, which does not tend to infinity, and entropy greater than zero indicates the existence of chaos). According to the data in Table 3, entropy and correlation dimension values decrease as the distance from the blast center increases. This proves the decreasing influence of various parameters on the blast wave propagation. The entropy value decreases by a factor of 8.5, and the correlation dimension decreases by a factor of 14.8. Thus, chaos decreases 1.74 times slower than the number of parameters affecting the system. Both characteristics tend to zero, which means that the forecast horizon for all parameters is sufficient and the signal under investigation is not a random value at the entire experimental area. Predictability time was determined for each time series. The low-level forecast (for minimal distance to the blast) is possible for 5 points (values), and the max-level one is possible (for max distance) for 50 points. Predictability time shows that a chaotic time series can be continued for a limited number of steps when the system itself is unclear.

Motivation or Practical Implications
The seismic blast wave parameters forecasting at different points are especially relevant when drilling and blasting operations are carried out in the immediate vicinity of the protected objects (for example, gas pipeline construction). The blast wave peak velocity direct determination is impossible, as factors such as soil vibration and close objects' natural vibrations wield a major influence on the physical nature of the wave propagation and the close object's safety is vulnerable. Soil or object vibrations occurring at certain frequencies probably cause resonance phenomena. Thus, the averaged (or theoretical) values of the application of neglected amendments of formulas or empirical coefficients entail unpredictable consequences. Thereby, the practical value of experimentally determined methods for empirical coefficients or particles' peak vibration velocity forecasting during drilling and blasting operations, especially close to the protected objects, is of great relevance.
The practical result of this research, based on the experimental data, is the formulated method for determining the optimal parameters for drilling and blasting, including the following stages: 1 Obtain the values of experimental seismograms for specific conditions. This involves a series of experimental blasts with an equal small mass of explosive at different distances. Carry out fractal analysis. Define entropy, correlation dimension, and predictability time. Define the forecast horizon based on predictability time. Make sure that the obtained characteristics are within the area of the forecast horizon.
It is worth pointing out that additional experiments with a fixed distance but different masses of explosives use the same principles to determine the max mass of explosives in order to ensure the safety of protected objects during drilling and blasting operations.
This algorithm shows that this research has possible practical applications. In fact, the application area is much wider and, obviously, requires further research.

Conclusions
The principal hypothesis of this research is that the application of mathematical methods and signal processing to seismograms can provide extended data on drilling and blasting operations applicable to the determination of empirical coefficients in characteristic equations within the time period that do not exceed the system computation cycle.
That is, the analysis results of the seismograms associated with the blast wave characteristics propagating through the rock enable us to create an algorithm. This programmed algorithm can provide us with prompt and certain recommendations for optimizing the parameters of drilling and blasting. This is especially relevant considering the mineral sector's current transition to the digital economy [40]. The principal target of this research is to see the dependencies between the spectral, wavelet and fractal analyses results and the blast wave characteristics.
The research conclusions are as follows.
(1) For the dependence between the dominant frequency and amplitude and the distance to the blast center: the dominant frequency dependency is determined by the equation y = 1.0262x 0.2622 , amplitude dependency by y = 18.139x −2.276 (2) Of the entire signal, 80% is concentrated on half of the entire frequency range, i.e., the low frequency zone is the most interesting for the analysis. (3) The distance effect could be estimated by certain frequency range analysis. This research has identified that the 12th range has a close correlation ratio (more than 0.9), as well as 2, 4-8, 11, 12, and 14-16 ranges with a less close one, but also quite a strong dependency (more than 0.7). (4) Entropy and correlation dimension decrease and predictability time increases when the distance to the blast center increases.