Attenuation model of tunnel blast vibration velocity based on the influence of free surface

In tunnel blasting excavation, it is important to clarify the attenuation law of blast wave propagation and predict the blast vibration velocity effectively to ensure safe tunnel construction and protection design. The effects of the free surface area its quantity on the blast vibration velocity are considered, and free surface parameters are introduced to improve the existing blast vibration velocity prediction formula. Based on the Tianhuan railway Daqianshiling tunnel project, field blast vibration monitoring tests are performed to determine changes in the peak blasting vibration velocity based on the blast distance and free surface area. LS-DYNA is used to establish tunnel blasting excavation models under three operating conditions; subsequently, the attenuation law of blast vibration velocity and changes in the vibration response spectrum are analysed. Results show that the free surface area and number of free surfaces enable the blast vibration velocity to be predicted under various operating conditions: a smaller free surface area results in a narrower frequency band range, whereas more free surfaces result in a narrower frequency band range. The improved blast vibration velocity prediction formula is validated using field and numerical test data. It is indicated that the improved formula is applicable to various tunnelling conditions.


Attenuation model of tunnel blast vibration velocity based on the influence of free surface Baoxin Jia * , Linli Zhou, Jiaojiao Cui & Hao Chen
In tunnel blasting excavation, it is important to clarify the attenuation law of blast wave propagation and predict the blast vibration velocity effectively to ensure safe tunnel construction and protection design. The effects of the free surface area its quantity on the blast vibration velocity are considered, and free surface parameters are introduced to improve the existing blast vibration velocity prediction formula. Based on the Tianhuan railway Daqianshiling tunnel project, field blast vibration monitoring tests are performed to determine changes in the peak blasting vibration velocity based on the blast distance and free surface area. LS-DYNA is used to establish tunnel blasting excavation models under three operating conditions; subsequently, the attenuation law of blast vibration velocity and changes in the vibration response spectrum are analysed. Results show that the free surface area and number of free surfaces enable the blast vibration velocity to be predicted under various operating conditions: a smaller free surface area results in a narrower frequency band range, whereas more free surfaces result in a narrower frequency band range. The improved blast vibration velocity prediction formula is validated using field and numerical test data. It is indicated that the improved formula is applicable to various tunnelling conditions.
Blasting technology is vital to transportation, construction, water conservancy, and mining industries; however, in tunnel blasting excavation, the dynamic effects caused by blast vibration can be detrimental to the surrounding rock environment and lining structure. To avoid the damage and effect caused by blast vibration, the investigation into the blast vibration attenuation mechanism has received significant attention, among which the prediction and control of blast vibration intensity are crucial. The peak particle vibration velocity is representative of the blast vibration intensity, and the best criterion for estimating the vibration damage level of structures is the blast vibration velocity. Therefore, the peak particle vibration velocity is used as the control criterion for blast vibration 1 . Accurately predicting the blast-induced particle vibration velocity to ensure the safe and smooth operation of blasting projects is the main method for controlling blast vibration hazards [2][3][4][5] .
Studies pertaining to the blast vibration velocity as an object have been conducted extensively. For example, Cao et al. 4 analysed the transmission propagation law of blast stress waves and solved the safety threshold of blast vibration velocity based on the ultimate strength theory. Sun et al. 6 proposed a value of blast vibration velocity for shallow rock, and Singh et al. 7 determined the safety level of the ground structure under blast vibration. Both discussed the problem of determining the safety value of blast vibration from the blast vibration velocity. Zhuang et al. 8 monitored the blast vibration law of semi-coal rock laneway, obtained the blasting main vibration frequency from 50 to 250 Hz, and performed Sadoff formula regression fitting to calculate the safe allowable distance for blasting in different rock layers. Kim et al. 9 proposed geometric damping coefficients; Tian et al. 10 considered different paths of vibration propagation under laminated rock conditions and investigated the blast vibration propagation and attenuation laws. Yu et al. 11 conducted field tests to investigate the propagation law of blast vibration in the surface and blast-face side of shallow buried small clearance tunnels. Zhu et al. 12 preferred a blast prediction model that considered the blasting charge and blast centre distance, as well as one that considered the blasting charge and height difference; subsequently, they established a segmental prediction model. Wu et al. 13 elucidated the amplification effect of step terrains on the blast vibration velocity and established a vibration velocity prediction model for step terrains. Zhu et al. 14 established a multi-degree-of-freedom model of blast vibration and its differential equations of motion to propose a method for predicting the blast vibration velocity of a multihole cut blast in layered rock masses. Applying computer theory technology, Feng 15 used empirical formulas and a grey model to predict the peak vibration velocity. Zou et al. 16 used asymmetric triangular fuzzy numbers, which were used to determine the parameters of a prediction formula, whereas Ma et al. 17  www.nature.com/scientificreports/ Singh et al. 18 , and Wei et al. 19 used neural network algorithms, artificial intelligence, and finite element simulation techniques to establish prediction models, respectively. Yue et al. 20 combined the particle swarm algorithm and least-squares support vector machine model to predict the peak blast vibration velocity. Cai et al. 5 improved the Elman neural network model by considering the free surface area to predict the blast vibration velocity using the beetle antenna search algorithm. Whereas the above mentioned studies considered many factors such as total charge, blast distance, elevation, hole diameter, number of holes, and blasting duration to establish prediction models including the neural network model and least squares support vector machine, less research has been conducted regarding the effect of free surface parameters on blast vibration attenuation. In addition, in tunnel blasting excavation, the cutting holes create a free surface, and the area and number of free surfaces affect blasting excavation and tunnel blast vibration attenuation. In this study, based on the Sadoff empirical formula and the vibration velocity attenuation formula of Lu et al. 21 , free surface parameters are introduced into a formula to improve the prediction of the peak blast vibration velocity.

Tunnel blast vibration velocity prediction formula
In the blasting industry, blast vibration hazard effects are evaluated by predicting the peak particle vibration velocity, and blast vibration control measures are formulated; however, the prediction accuracy is a longstanding issue in the blasting industry. Based on extensive engineering experience and research, scholars have proposed different empirical formulas for the blast vibration velocity. Currently, the empirical formula of the peak blast vibration velocity proposed by Sadoff, which is shown in Eq. (1), is the most widely used in China's blasting industry: where V is the peak blast vibration velocity (cm/s); Q is the maximum single detonation charge (kg); R is the source distance (m); K and q are field coefficients.
The formula above does not reflect the effects of explosive properties, explosive type, charging method, hole size, and surrounding rock parameters on the blast vibration velocity. Lu et al. obtained an equation expressing particle blast vibration velocity using the fluctuation method to reflect the above mentioned factors. Based on the wavelet of the long column charge and column surface wave theory, the Helen solution for the stress wave field excited by the column charge package was analysed to derive an equation for the peak blast vibration velocity 21 , which is shown in Eq. (2).
where k is the improvement factor associated with the blasting conditions, p 0 the initial pressure of the blast-generated gas in the hole (Pa), b the radius of the hole (mm), ρ the rock density (kg/m 3 ), R the distance of the detonation point from the particle (m), c p the rock longitudinal wave velocity (m/s), and α the blast vibration attenuation index. When the charge is coupled, p 0 = p e , where p e = ρ e D 2 /2(γ + 1). When the charge is uncoupled, if the uncoupled factor b/a is smaller, then p 0 = p e (a/b) 2γ ; if b/a is larger, then p 0 = [ρ e D 2 /2(γ + 1)] γ 0 /γ p k (γ −γ 0 ) (a/b) 2γ 0 , where ρ e is the explosive density, D the explosive detonation velocity, γ the adiabatic index, p k the critical pressure of the explosive, p e the average explosive detonation pressure, and a the radius of the charge.
Generally, the surface of a blasted rock in contact with air is referred to as the free surface. The larger the free surface area and the greater the number of free surfaces, the more prominent is the blasting effect. The free surface is one of the important factors affecting the blasting effect. Hence, the prediction accuracy can be improved by investigating the effects of free surface factors on the blast vibration velocity and the influence law of the free surface into the prediction formula.
The free surface parameters were introduced into the vibration velocity calculation formula to improve the peak vibration velocity prediction formula based on the Sadoff formula and the peak vibration decay formula by Lu et al. This is because the larger the free surface area, the greater is the number of free surfaces, and the smaller is the surrounding blast vibration 22 . Therefore, the free surface area S, number of free surfaces m, and free surface index β were selected as the free surface parameters and then incorporated into Sadoff 's formula [as shown in Eq. where S is the free surface area, which is the surface area where the hole arrangement is located; m is the number of free surfaces, which is the number of blasting rock mass surfaces exposed to air; β is the free surface index. The remaining variables are the same as those in Eqs. (1) and (2). www.nature.com/scientificreports/

Tunnel blasting field tests
Field overview. Based on the Tianhuan railway Daqianshiling tunnel as an example, the total length of the tunnel was 2460 m, which reflects a medium-sized tunnel in a railway tunnel. The tunnel inlet and outlet sections DK69 + 225-DK69 + 495 and DK71 + 180-DK71 + 685 contained level IV-V surrounding rocks, excavated via the overrunning small conduit pre-supporting step method, and the middle section DK69 + 495-DK71 + 180 contained level II-III surrounding rocks, excavated via the full section method; both excavation methods are the blasting excavation method. Based on the geological survey data of the project and indoor rock mechanics test, it was discovered that the surrounding rock of the tunnel was quartz sandstone, which is a hard rock, and the surrounding rock parameters are listed in Table 1. The tunnel cycle feed was set as 1.6 m. The length of the steps was 10 m, which is considered short. The tunnel cross-section was a horseshoe-shaped four-centred circular arch tunnel, and the tunnel cross-section dimensions are shown in Fig. 1. The blasting excavation was performed using rock emulsion explosives with a TNT equivalent value of 0.73 23 . The explosive density was 1200 kg/m 3 , the explosive detonation velocity was 4500 m/s, and the adiabatic index was 3. In the blasting design, 153 holes were arranged in the outer contour of the tunnel, including two bum holes, seven cutting holes, 87 reliever holes, and 57 rim holes; the hole arrangement is shown in Fig. 2. The diameter of the hole was 35 mm, and the charging method was coupled. The total amount of charge in the entire section was 152.16 kg, 72.48 kg in the upper step and 79.68 kg in the lower step.
Tunnel blast vibration monitoring. The data acquisition system developed by the Institute of Geology, China Earthquake Administration, was used for the blast vibration monitoring test. It comprised five parts: sensors, cables, data hubs, converters, and signal collectors (Fig. 3). The signal acquisition frequency of the blast test was set to 50 kHz, and the sensor was a non-directional velocity sensor.
Seven monitoring points were established for blast vibration monitoring, and five blasting tests were conducted. Owing to the complexity of the field environment and the difference in the blasting effect of the rim holes, the free surface area of each blast changed, and the free surface areas of the fourth and fifth blasts increased because of the addition of an emergency avoidance hole in those areas. The monitoring points were fixed to the right side of the interior of the tunnel. The monitoring points were named using two digits, where the first digit represents the serial number of the blast, and the second digit represents the location of the monitoring point. For example, 3-2 indicates the third blast-monitoring point, No. 2. The peak vibration velocity of the extracted monitoring signal, the area of the free surface of each blast, and the distance of each monitoring point from the blast free surface are listed in Table 2.
Comparison of blast velocity prediction formulas. The burst source distance R, peak vibration velocity V, and free surface area S listed in Table 2 were used to establish the elements individually in the correspond- www.nature.com/scientificreports/ ing column vector, based on the prediction formula shown in Eqs. (1) to (4). The MATLAB Curve Fitting Tool was used for fitting, and the fitting curve is shown in Fig. 4, whereas the fitting parameters are listed in Table 3. It was clear that the Sadoff formula and Lu et al. formula for predicting the peak vibration velocity considering the blast vibration parameters were different. However, the relationship between the change in the peak vibration velocity and the propagation distance was the same; hence, the effects of the two fittings were the same.   www.nature.com/scientificreports/ The summed variance (SSE), standard deviation (RMSE), coefficient of determination (R-square) and adjusted coefficient of determination (adjusted R-square) in Table 3 were calculated and generated using the MATLAB Curve Fitting Tool. They were calculated as shown in Eqs. (5)-(8), respectively, as follows: In the equations above, w i is the data sample weight, y i the original data, ŷ i the predicted data, y i the original data mean, n the number of samples, and p the number of features.
The closer the SSE and RMSE are to 0, the closer are the R-square and adjusted R-square to 1. These values reflects a well fit model, indicating that the independent variables can better reflect the dependent variable. Therefore, these four fitting parameters can be used to compare the strengths and weaknesses of the fit of the prediction formulas. By comparing the data shown in Table 3, it can be concluded that the blast vibration velocity prediction formula with the addition of the free surface parameters is better in terms of fitting. The Fig. 4 shows that the distribution of scattered points on the curve is more scattered than the distribution of scattered points on the surface. This indicates that the effects of the free surface area and number of free surfaces on the peak blast vibration velocity are non-negligible. Therefore, free surface parameters must be incorporated into the base prediction formula.  www.nature.com/scientificreports/

Tunnel blasting numerical tests
Owing to certain limitations of the field tests, e.g. small variations in the free surface area, vibration velocity monitoring was performed only during the full section stage of tunnel excavation. Therefore, the effects of the free surface parameters of tunnel blasting on the blast vibration velocity is discussed based on numerical calculations. Two numerical models were established: a full section model (as shown in Fig. 5a) and a step model (as shown in Fig. 5b). Each model was segregated into two modules, one for the excavated section, and the other for the unexcavated section. Furthermore, the same physical and mechanical parameters were set for the surrounding rock in both sections. In the full section model, the length of the excavated section was 50 m, i.e. the depth of the tunnel was 50 m. In the step model, the length of the step was 10 m, and the remaining dimensions were the same as those of the full section model. The model grid size was 0.25 m. The hole was located in the centre of the palm face, occupying a cell with a depth of 3 m, and 12 cells were removed. Because the model boundary was limited, i.e. much smaller than the field, when the fixed boundary conditions are used, the reflection of the boundary on the blast wave is evident; therefore, the boundary conditions were set to the no-reflection boundary, 'BOUNDARY_ NON_REFLECTING' , to reduce the reflected wave generated by the model boundary to the model perturbation. The tunnel calculation model was segregated into three regions as follows: (1) Full section blasting-a full section area measuring 48.9515 m 2 , a free surface is present, and the source is located in the centre of the full section;

LS
(2) upper step blasting-an upper step area measuring 22.8886 m 2 , a free surface is present, and the source is located in the centre of the upper step section; (3) lower step blasting-a lower step area measuring 26.0629 m 2 , two free surfaces are present, and the source is located in the centre of the lower step section.
The Mohr-Coulomb material keyword, 'MAT_MOHR_COULOMB' , was selected for the surrounding rock constitutive model. The blast load was set using the 'LOAD_BLAST_ENHANCED' function 24 , which simulates the action of the blast wave with the geotechnical mass based on 'LOAD_BLAST_SEGMENT_SET' . This blast setting method disregards the shape of the explosive charge and the propagation of the blast shock wave in air; as such, the charge and air need not be modelled, which reduces computational time. Zhang et al. 25 and Gilson et al. 26 applied the method above and verified its feasibility.
In this method of applying blast loads, the TNT equivalent of the explosive must be defined; additionally, the blast action plane and blast source coordinates must be selected. The blast source coordinates were set as the centre coordinates of the hole (see Table 4), and the blast action plane was selected as the inner side of the hole. Because the cell measured 0.25 m and the actual hole diameter was 0.035 m, the area of one cell corresponded to 65 actual holes after conversion. Because the average charge of a hole was 1 kg (based on the engineering blasting design), 65 kg of explosives were set in the holes of all three types of models; furthermore, it was ensured  www.nature.com/scientificreports/ that the explosives would not exceed the maximum charge under various operating conditions. The emulsified explosive charge was converted to a TNT equivalent of 47.45 kg. The blasting parameters of the three types of models were unified and satisfied the requirements of the comparison test.
To verify that the model above can reflect the field test results, the length of the excavated section of the model was extended to 150 m, which eases comparison with field monitoring data. Because the field data were only available for the full section excavation stage, model validation was performed only for the full section model. Disregarding the change in the free surface area in the field test, the monitoring point vibration velocity of the first three blasting full sections was selected for the calibration data. The numerical model explosive quantity was adjusted to 111.08 kg TNT equivalent for calculation. A comparison between the numerical simulation and monitoring results is shown in Fig. 6, which shows that the numerical model established in this study can simulate the field test accurately.
Particle vibration velocity monitoring. The monitoring points were set at the side walls of the upper and lower step separation interfaces and arranged horizontally along the tunnel axis. The monitoring points for the three types of model were in the same position. Eight monitoring points were set in the full section model and the lower step model, beginning from the palm face; the first six monitoring points were spaced 5 m apart (20 units), and the last two monitoring points were spaced 7.5 m apart (30 units). In the upper step model, the eight monitoring points were retained, and two monitoring points were added on the same horizontal line on the step with an interval of 5 m (20 units) (see Fig. 7). The coordinates of each monitoring point are listed in Table 5, and the monitoring points were numbered based on the node number. The resultant velocity of the monitoring points was extracted via numerical model calculations. The effects of the free surface parameters on the blast vibration velocity were compared and analysed.
Each monitoring point blast distance can be determined using the blast source coordinates and monitoring point coordinates, and the blast distances are shown in Table 6. The velocity time course curves are shown in   www.nature.com/scientificreports/ The peak particle vibration velocity at each monitoring point was determined, as shown in Table 7, and the V-R scatter diagram is plotted in Fig. 9.
Analysis of peak particle vibration velocity variation. Attenuation of peak particle vibration velocity. The MATLAB Curve Fitting Tool was used to fit the scatter data of the three models using Eq. (1) while the full scatter was fitted, and the statistics of the fitted parameters are listed in Table 8.
As shown in Table 8, the fitting effect was favourable when fitting the blast vibration decay curves individually; however, the differences among the fitted curves were significant. When all the data were fitted together, the fit was inferior. This indicates that the change in the blasting free surface parameters significantly affected the tunnel blast vibration velocity decay. The results of the numerical calculations were fitted using Eqs. (1) and (4). The fitted curves are shown in Fig. 10, and the fitted parameters are listed in Table 9.
The numerical test fitting indexes indicate that the blast vibration velocity prediction formula considering the free surface parameters significantly improved the fitting indexes and yielded a better fitting effect than the base formula. This indicates that Eqs. (3) and (4) can be more widely used for predicting the peak vibration velocity in various tunnel blasting sections.
Vibration response spectrum characteristics. A spectral analysis of the vibration velocity waveforms at each monitoring point was performed, and the spectrum plots of each monitoring point of the three models were obtained using Fourier transform, as shown in Fig. 11.
The vertical axis amplitude was increased by 100 times, and the frequency range with an amplitude greater than 0.5 was regarded as the main frequency range of vibration. As shown in Fig. 11, the numerical model blast vibration was dominated by low-frequency vibrations. Near the free surface monitoring point, the vibration band was the widest, and the vibration response was the most significant. As the propagation distance increased, the vibration response weakened, and the amplitude decay and band range reduction stabilised. Figure 11a shows that the vibration frequency range of the full section model was primarily in the range of 0-400 Hz, and that the high-frequency components above 150 Hz decayed the fastest. Figure 11b shows that the vibration frequency range of the upper step model was primarily in the range of 0-350 Hz, and the frequency components above 90 Hz decayed rapidly. Figure 11c shows that the overall vibration amplitude of the lower step model was significantly lower than that of the first two models, and that the vibration frequency range reduced significantly reduced, primarily in the range of 0-280 Hz, with the frequency components above 90 Hz decaying rapidly. As shown in Fig. 11, the amplitude of the vibration speed waveform decreased gradually from low frequency to high www.nature.com/scientificreports/ frequency. The inflexion point of the amplitude of the vibration speed waveform with frequency was at 100 Hz. The amplitude decay rate was high, and the curve in the 0-100 Hz band was steep. After 100 Hz, the amplitude decay rate was low, and the curve was flat. Therefore, with 100 Hz as the cutoff point, the waveform amplitude percentage was calculated for 0-100 Hz, as shown in Table 10. As shown, as the propagation distance increased, the low-frequency component of the waveform at 0-100 Hz increased proportionally and dominated. All the results show that the vibration response from blasting contained both high-and low-frequency components; however, as the blast wave propagated, the high-frequency components decayed rapidly, and the vibration response retained only the low-frequency components and remained in a stable frequency band range.
Comparing the three numerical models, it was discovered that the widths of the initial frequency bands differed for different free surface areas and numbers. It was observed that a smaller free surface resulted in a      www.nature.com/scientificreports/ narrower frequency band range, whereas more free surfaces resulted in a narrower frequency band range. The step model filtered out more high-frequency components and retained fewer low-frequency components.

Discussion and conclusions
In previous studies, the effects of free surface parameters on the blast vibration velocity decay were rarely considered, although the free surface parameters affected the blasting effect. Therefore, based on the blast vibration decay prediction formula proposed previously, the blast vibration prediction formula was improved by considering the effects of two parameters, namely, the free surface area and the number of free surfaces, as shown in Eqs. (3) and (4), respectively. The progress of the improved formula was verified using field blast monitoring data. A numerical simulation was performed using LS-DYNA to analyse the decay of the blast vibration velocity under three different operating conditions of tunnel excavation, and the results verified the superiority of the improved formula. The conclusions of this study are as follows:  www.nature.com/scientificreports/ 1. The blast vibration velocity prediction formulas by Sadoff and Lu et al. were improved by introducing two influencing parameters, i.e. the free surface area S and the number of surfaces m, by incorporating a factor term, (S/mR 2 ) β . The original blast vibration velocity prediction formula yielded better prediction only under a single operating condition, and the prediction formula can be applied more widely after adding the free surface parameters. The fitting index indicated that the improved prediction formula demonstrated better model selection and fitting. 2. The initial vibration response of the blast vibration was significant; in addition, the amplitude decayed rapidly as the propagation distance increased and then stabilised in the form of an exponential decay. 3. The main frequency range of blast vibration in the initial stage indicated a wide range of frequency bands with clear high-frequency components; as the propagation distance increased, the high-frequency components decayed rapidly, and the frequency band narrowed. The final vibration response was dominated by low frequencies, primarily in the 0-100 Hz range.