# Study on the Attenuation Characteristics of Seismic Wave Energy Induced by Underwater Drilling and Blasting

Authors: Yaxiong Peng; Ying Su; Li Wu; Chunhui Chen
Journal: Shock and Vibration
 (2019)

---

## Abstract

Seismic wave induced by underwater drilling and blasting may exert adverse effects to the surrounding structures. Peak particle velocity (PPV) has been widely used to reflect the intensity of seismic wave, while the method fails to take other influential parameters into consideration. Synthesizing PPV and duration, seismic wave energy (SWE) is employed to quantify the intensity, and the method of time-frequency analysis is adopted to study the characteristics of vibration frequency. Besides, this paper deduces the prediction formula for SWE via dimension analysis. Based on the practical engineering, the attenuation characteristics of SWE induced by underwater drilling and blasting is analysed by the wavelet transform and adaptive optimal kernel (AOK) time-frequency analysis method. Results show that the dominant frequencies at the directions of horizontal tangential, horizontal radial, and vertical are different and the energies in high-frequency bands are extremely low. Moreover, the lower the frequency is, the slower the attenuation of SWE is. Comparing with other prediction formulas of PPV, fitting the SWE with the help of the prediction formula in this study would achieve more accurate prediction results.


---

## Body

## 1. Introduction

Underwater drilling and blasting has become an important rock excavation method in construction projects such as harbours, wharfs, and channels. With the development of waterway transport, such projects have increased sharply and projects with large rock mass volume and large explosive charge have been the trend. In underwater blasting engineering, only 20% to 30% energy released by explosive explosion is spent to break rock mass, while the rest exerts adverse effects to the surrounding environment, such as seismic wave, shock wave, dust and fumes, and fly rocks. [1–3]. Among which, seismic wave is the most destructive one. That is to say, the most effective method to reduce the adverse effects is to control the intensity of seismic waves.Three key elements of seismic wave include peak particle velocity (PPV), frequency, and duration. PPV is frequently employed to quantify the intensity of seismic wave, and the value is deemed to be influenced by the quantity of explosive charge, the distance from the measured point to the blast-face, and the geological conditions of the excavation area. In a certain site, the geological conditions and the distance from the measured point to the blast-face are uniquely determined. Thus, field tests are applied to predict PPV and the quantity of explosive charges are adjusted to lower the value to the safe range. Through numbers of experiments, the prediction formulas of PPV have been derived by numerous scholars [4–7]. And, the elevation difference factor β was introduced to modify the prediction formula [8, 9]. The vibration attenuation law and elevation effect of underwater drilling and blasting were studied by nonlinear regression and robust regression [10, 11]. Gu et al. [12, 13] deduced a revised formula for water depth and obtained PPVs in water bottom and shore. In order to improve the prediction accuracy, the genetic algorithm has been widely used to fit PPV [13–15]. Khandelwal and Singh [16–18] applied various artificial neural networks to predict PPV, and all achieved good results. Liu et al. [19] studied PPV induced by underwater blasting based on the combination of grey relational analysis and genetic neural network.When the seismic wave frequency is close to the structure natural frequency, the structure may be damaged by resonance. Due to the continuous blasting vibration, the allowable vibration velocity of the structure may decrease dramatically. Therefore, the frequency and duration of seismic waves also have large influence on the nearby structure. The research team of Professor Lu studied the distribution characteristics and attenuation law of blasting frequency under different engineering conditions [20, 21]. Zhong et al. [22] discussed the influential parameters of duration and established the prediction formula with the help of dimensional analysis. In terms of underwater blasting engineering, Peng et al. [23] deduced the functional relationship among the vibration frequency and other relevant influential parameters. Lia et al. [24] proposed a novel model for predicting PPV, which embodied the effects of geological factor, rock properties, and engineering category.The above researches are based on the single factor for the intensity of seismic waves, which fail to take all the factors into consideration. Thus, numerous scholars have studied the seismic waves from the aspect of seismic wave energy (SWE). Sanchidrián et al. [25] divided the elements of explosive energy into 3 parts, the fragmentation energy, the seismic energy, and the kinetic energy, and deduced the corresponding calculation formulas, respectively. Li et al. [26] studied the distribution characteristics of SWE under different forms of blasting and obtained the prediction formula of energy attenuation. Ling and Li [27] and Yang et al. [28] studied the distribution characteristics and the safety control method of SWE with the help of wavelet transform and time-frequency analyses. Researches above indicate that the parameter of SWE can perfectly reflect the characteristics of seismic wave.Above researches are mainly based on land blasting, while underwater drilling and blasting is quite different from land blasting due to the presence of water media. Thus, the paper deduces the attenuation formula for the total energy of seismic wave based on the signal analysis of particle energy. Many field measured data are collected, and the distribution characteristics and attenuation law of SWE in different frequency bands are studied with the help of wavelet transform and AOK time-frequency analysis. Furthermore, the attenuation formula is adopted to predict the SWE in practical engineering.


## 2. Theoretical Study

### 2.1. Prediction Formula of PPV

It has been the consensus that the particle vibration velocity can be employed to describe the intensity of seismic wave. By analysing a large number of measured data, researchers deduced various prediction formulas of PPV, part of which is shown in Table1 [4–7]. In underwater drilling and blasting engineering, monitoring points can only be set on the riverbank slope. Under the circumstance of deep water depth, the influence that the elevation difference has on the propagation of seismic wave will also increase substantially. Chinese Construction Technical Specifications (CCTS) [9] proposed a prediction formula of PPV involving the effects of elevation, which is shown in Table 1. In the table, Q is the maximal charge per delay, R is the distance from the measured point to the blast-face, H is the altitude from the measured point to the blast-face, and k, α, and β are the parameters related to the transmission medium.Table 1
Different vibration prediction formulas [4–7, 9].
Prediction formula
Equation
USBM
v
=
k
R
/
Q
−
α
Langefors–Kihlstrom
v
=
k
Q
/
R
2
/
3
α
Bureau of Indian Standard
v
=
k
Q
/
R
2
/
3
α
Sadaovsk formula
v
=
k
Q
1
/
3
/
R
α
CCTS
v
=
k
Q
1
/
3
/
R
α
Q
1
/
3
/
H
β


### 2.2. Prediction Formula of SWE

When studying the influence that blasting seismic waves exert on the surrounding environment, the kinetic energy is the primary factor to be concerned. At timet, the SWE for an elementary mass Δm in space can be expressed as follows:(1)et=12Δm⋅v2t,where et is the SWE at time t, Δm is the elementary mass, and vt is the velocity at time t.To leave out the effect of point mass, the SWE per mass unitemt can be deduced as follows:(2)emt=etΔm=12v2t.Thus, the energy per mass unitE for the entire process of blasting vibration can be calculated by(3)E=∫0temtdt=12∫0tv2tdt.Blasting seismic wave is a random signal with transient characteristics. As the formation and propagation of blasting seismic wave are affected by numerous factors, it is too hard to establish the prediction formula which takes all the factors into consideration. Thus, merely several major influential parameters are selected for studies. A dimensional analysis (DA) technique has been performed on major influential parameters to propose a new formula [29].In general, PPV (v), seismic wave energy (E), seismic wave frequency (f), and seismic wave duration (t) are the primary parameters to describe the effect of blasting seismic wave [13, 30]. In underwater blasting engineering, maximal charge per delay and distance from the measured point to the blast-face have remarkable influence on SWE, and geological conditions and rock properties also have major effect on it [26, 31]. As a result, maximal charge per delay (Q), distance from the measured point to the blast-face (R), altitude from the measured point to the blast-face (H), hole depth (h), delay interval (Td), density of rock mass (ρ), and medium velocity of seismic waves (c) are chosen as the main physical quantities. Therefore, it is generally agreed that the propagation process of seismic wave induced by underwater drilling and blasting involves 11 parameters, as shown in Table 2.Table 2
Major influential parameters.
Type
Symbol
Meaning
Dimension
Dependent variables
v
Peak particle velocity
LT−1
E
Seismic wave energy
L2T−2
f
Seismic wave frequency
T−1
t
Seismic wave duration
T
Independent variables
Q
Maximal charge per delay
M
R
Distance from the measured point to the blast-face
L
H
Altitude from the measured point to the blast-face
L
h
Hole depth
L
T
d
Delay interval
T
ρ
Density of rock mass
ML−3
c
Medium velocity of seismic wave
LT−1As seen from the table, totally 11 influential parameters are adopted to quantify the intensity of seismic wave induced by underwater drilling and blasting. According to theπ theorem, Q, R, and c are selected to be the basic physical variables [13]:(4)FQ,R,c,v,E,f,t,H,h,Td,ρ=0.Suppose thatπ is denoted as the dimensionless group, the value of π can be solved via the following composite formula:(5)π=QxRyczX,where X is the variable except the above basic physical variables in Table 2; x, y, and z are undetermined coefficients.According to the dimensionless homogeneous theorem, the variablev can be calculated via the following equation:(6)dimv=LT−1=MxLyLT−1z.Asx1 = 0, y1 = 0, and z1 = −1, the simplified formula is depicted as follows:(7)π1=c−1⋅v.For the convenience of analysis, other simplified formulas are given:(8)π2=c−2⋅E,π3=Rc−1⋅f,π4=R−1c⋅t,π5=Q1/3H−1⋅ρ−1/3,π6=R−1⋅h,π7=R−1c⋅Td,π8=Q1/3R−1⋅ρ−1/3.Combining equations (4), (7), and (8), π can be expressed as follows:(9)π=Φc−1v,c−2E,Rc−1f,R−1ct,Q1/3H−1ρ−1/3,R−1h,R−1cTd,Q1/3R−1ρ−1/3.Since the product of several dimensionless values is still a dimensionless value, the equation of similar dimensionless group related to SWE can be depicted as follows:(10)π2=π5⋅π8.According to equations (9) and (10), there is a certain functional relationship among c−2E, Q1/3H−1ρ−1/3, and Q1/3R−1ρ−1/3. As the parameters c and ρ for a certain site can be treated as constants [8], the function can be expressed as follows:(11)lnE=k1+γ1lnQ1/3R−k2−γ2lnQ1/3H.Referring to the Sadaovsk Formula, assume thatkE=−k1k2, αE=γ1, and βE=γ2, and the prediction formula of SWE can be derived as follows:(12)E=−k1k2Q1/3Rγ1⋅Q1/3Hγ2=kEQ1/3RαE⋅Q1/3HβE,where E is the SWE, cm2/s2; kE, αE, and βE are the parameters related to the transmission medium.


## 2.1. Prediction Formula of PPV

It has been the consensus that the particle vibration velocity can be employed to describe the intensity of seismic wave. By analysing a large number of measured data, researchers deduced various prediction formulas of PPV, part of which is shown in Table1 [4–7]. In underwater drilling and blasting engineering, monitoring points can only be set on the riverbank slope. Under the circumstance of deep water depth, the influence that the elevation difference has on the propagation of seismic wave will also increase substantially. Chinese Construction Technical Specifications (CCTS) [9] proposed a prediction formula of PPV involving the effects of elevation, which is shown in Table 1. In the table, Q is the maximal charge per delay, R is the distance from the measured point to the blast-face, H is the altitude from the measured point to the blast-face, and k, α, and β are the parameters related to the transmission medium.Table 1
Different vibration prediction formulas [4–7, 9].
Prediction formula
Equation
USBM
v
=
k
R
/
Q
−
α
Langefors–Kihlstrom
v
=
k
Q
/
R
2
/
3
α
Bureau of Indian Standard
v
=
k
Q
/
R
2
/
3
α
Sadaovsk formula
v
=
k
Q
1
/
3
/
R
α
CCTS
v
=
k
Q
1
/
3
/
R
α
Q
1
/
3
/
H
β


## 2.2. Prediction Formula of SWE

When studying the influence that blasting seismic waves exert on the surrounding environment, the kinetic energy is the primary factor to be concerned. At timet, the SWE for an elementary mass Δm in space can be expressed as follows:(1)et=12Δm⋅v2t,where et is the SWE at time t, Δm is the elementary mass, and vt is the velocity at time t.To leave out the effect of point mass, the SWE per mass unitemt can be deduced as follows:(2)emt=etΔm=12v2t.Thus, the energy per mass unitE for the entire process of blasting vibration can be calculated by(3)E=∫0temtdt=12∫0tv2tdt.Blasting seismic wave is a random signal with transient characteristics. As the formation and propagation of blasting seismic wave are affected by numerous factors, it is too hard to establish the prediction formula which takes all the factors into consideration. Thus, merely several major influential parameters are selected for studies. A dimensional analysis (DA) technique has been performed on major influential parameters to propose a new formula [29].In general, PPV (v), seismic wave energy (E), seismic wave frequency (f), and seismic wave duration (t) are the primary parameters to describe the effect of blasting seismic wave [13, 30]. In underwater blasting engineering, maximal charge per delay and distance from the measured point to the blast-face have remarkable influence on SWE, and geological conditions and rock properties also have major effect on it [26, 31]. As a result, maximal charge per delay (Q), distance from the measured point to the blast-face (R), altitude from the measured point to the blast-face (H), hole depth (h), delay interval (Td), density of rock mass (ρ), and medium velocity of seismic waves (c) are chosen as the main physical quantities. Therefore, it is generally agreed that the propagation process of seismic wave induced by underwater drilling and blasting involves 11 parameters, as shown in Table 2.Table 2
Major influential parameters.
Type
Symbol
Meaning
Dimension
Dependent variables
v
Peak particle velocity
LT−1
E
Seismic wave energy
L2T−2
f
Seismic wave frequency
T−1
t
Seismic wave duration
T
Independent variables
Q
Maximal charge per delay
M
R
Distance from the measured point to the blast-face
L
H
Altitude from the measured point to the blast-face
L
h
Hole depth
L
T
d
Delay interval
T
ρ
Density of rock mass
ML−3
c
Medium velocity of seismic wave
LT−1As seen from the table, totally 11 influential parameters are adopted to quantify the intensity of seismic wave induced by underwater drilling and blasting. According to theπ theorem, Q, R, and c are selected to be the basic physical variables [13]:(4)FQ,R,c,v,E,f,t,H,h,Td,ρ=0.Suppose thatπ is denoted as the dimensionless group, the value of π can be solved via the following composite formula:(5)π=QxRyczX,where X is the variable except the above basic physical variables in Table 2; x, y, and z are undetermined coefficients.According to the dimensionless homogeneous theorem, the variablev can be calculated via the following equation:(6)dimv=LT−1=MxLyLT−1z.Asx1 = 0, y1 = 0, and z1 = −1, the simplified formula is depicted as follows:(7)π1=c−1⋅v.For the convenience of analysis, other simplified formulas are given:(8)π2=c−2⋅E,π3=Rc−1⋅f,π4=R−1c⋅t,π5=Q1/3H−1⋅ρ−1/3,π6=R−1⋅h,π7=R−1c⋅Td,π8=Q1/3R−1⋅ρ−1/3.Combining equations (4), (7), and (8), π can be expressed as follows:(9)π=Φc−1v,c−2E,Rc−1f,R−1ct,Q1/3H−1ρ−1/3,R−1h,R−1cTd,Q1/3R−1ρ−1/3.Since the product of several dimensionless values is still a dimensionless value, the equation of similar dimensionless group related to SWE can be depicted as follows:(10)π2=π5⋅π8.According to equations (9) and (10), there is a certain functional relationship among c−2E, Q1/3H−1ρ−1/3, and Q1/3R−1ρ−1/3. As the parameters c and ρ for a certain site can be treated as constants [8], the function can be expressed as follows:(11)lnE=k1+γ1lnQ1/3R−k2−γ2lnQ1/3H.Referring to the Sadaovsk Formula, assume thatkE=−k1k2, αE=γ1, and βE=γ2, and the prediction formula of SWE can be derived as follows:(12)E=−k1k2Q1/3Rγ1⋅Q1/3Hγ2=kEQ1/3RαE⋅Q1/3HβE,where E is the SWE, cm2/s2; kE, αE, and βE are the parameters related to the transmission medium.


## 3. Case Study

### 3.1. Project Profile

The construction project of Jiulongpo to Chaotianmen reach in the upper Yangtze River is located in the central downtown area of Chongqing, China. The length of the reach is around 22 km. As the designed water level of the Three Gorges is up to 175 m and the reach is located in the upper fluctuating backwater area, the maintenance workload will increase with years and may even form new silt shoals. There are 5 reefs to be exploded in the dredging project, and the total amount of the reef is about 154,000 m3. The measured data of this research are collected from the Zuanzhaozi Reef, which is saturated with moderate and strong weathering sandstones. The length, width, and height of the reef are, respectively, 210 m, 90 m, and 171 m. The construction period lasts for one year. And, due to seasonal changes of the water level of the Yangtze River, the water depth of construction ranges from 8 to 14 m. The location and the site conditions are depicted in Figure 1.Figure 1
Project location and site photos.The underwater holes were drilled with the help of a drill-burst ship. The single-row and multihole blasting was employed, and the number of holes was selected in line with the cross-sectional shapes of reefs. In order to adapt to underwater blasting environment,2# emulsion explosives were adopted, and the detonation velocity was about 3500 m/s. Besides, the electric network was selected to detonate the explosives, and electric detonators with millisecond delay time were used to control the delay interval. The detailed design parameters are shown as Table 3.Table 3
Detailed blasting design parameters.
Hole diameter (mm)
Hole depth (m)
Extra-deep (m)
Hole space (m)
Array pitch (m)
Number of hole
Delay time (ms)
110
7.2∼13.7
1.5
2
2.5
9∼12
25∼125


### 3.2. Data Monitoring

The data of blasting seismic wave were collected by TC-4850 blasting vibration meters, and the monitoring points were set on riverbank. The vibration meter has three channels and can collect seismic signals from three directions, the horizontal radial, the horizontal tangential, and the vertical direction. The received frequency is set at 4000 sps, and the sampling duration is 1 s. 30 groups of data are collected in this study, and the distance from the measured point to the blast-face is measured by GPS. The statistics for monitoring data are shown in Table4, and the typical signal of seismic wave is shown in Figure 2.Table 4
Statistics for monitoring data.
Parameter
Unit
Min.
Max.
Mean
Maximal charge per delay
kg
32
56
46.4
Total charge
kg
320
560
440
Distance from the blast-face to the monitoring point
m
65
260
87.5
Altitude from the blast-face to the monitoring point
m
20
40
29.2
Depth of blasthole
m
8.6
14.8
12.2
Peak particle velocity
cm/s
0.250
1.963
1.027Figure 2
Typical signal of seismic wave.


## 3.1. Project Profile

The construction project of Jiulongpo to Chaotianmen reach in the upper Yangtze River is located in the central downtown area of Chongqing, China. The length of the reach is around 22 km. As the designed water level of the Three Gorges is up to 175 m and the reach is located in the upper fluctuating backwater area, the maintenance workload will increase with years and may even form new silt shoals. There are 5 reefs to be exploded in the dredging project, and the total amount of the reef is about 154,000 m3. The measured data of this research are collected from the Zuanzhaozi Reef, which is saturated with moderate and strong weathering sandstones. The length, width, and height of the reef are, respectively, 210 m, 90 m, and 171 m. The construction period lasts for one year. And, due to seasonal changes of the water level of the Yangtze River, the water depth of construction ranges from 8 to 14 m. The location and the site conditions are depicted in Figure 1.Figure 1
Project location and site photos.The underwater holes were drilled with the help of a drill-burst ship. The single-row and multihole blasting was employed, and the number of holes was selected in line with the cross-sectional shapes of reefs. In order to adapt to underwater blasting environment,2# emulsion explosives were adopted, and the detonation velocity was about 3500 m/s. Besides, the electric network was selected to detonate the explosives, and electric detonators with millisecond delay time were used to control the delay interval. The detailed design parameters are shown as Table 3.Table 3
Detailed blasting design parameters.
Hole diameter (mm)
Hole depth (m)
Extra-deep (m)
Hole space (m)
Array pitch (m)
Number of hole
Delay time (ms)
110
7.2∼13.7
1.5
2
2.5
9∼12
25∼125


## 3.2. Data Monitoring

The data of blasting seismic wave were collected by TC-4850 blasting vibration meters, and the monitoring points were set on riverbank. The vibration meter has three channels and can collect seismic signals from three directions, the horizontal radial, the horizontal tangential, and the vertical direction. The received frequency is set at 4000 sps, and the sampling duration is 1 s. 30 groups of data are collected in this study, and the distance from the measured point to the blast-face is measured by GPS. The statistics for monitoring data are shown in Table4, and the typical signal of seismic wave is shown in Figure 2.Table 4
Statistics for monitoring data.
Parameter
Unit
Min.
Max.
Mean
Maximal charge per delay
kg
32
56
46.4
Total charge
kg
320
560
440
Distance from the blast-face to the monitoring point
m
65
260
87.5
Altitude from the blast-face to the monitoring point
m
20
40
29.2
Depth of blasthole
m
8.6
14.8
12.2
Peak particle velocity
cm/s
0.250
1.963
1.027Figure 2
Typical signal of seismic wave.


## 4. Characteristics of SWE

### 4.1. Time-Frequency Analysis

#### 4.1.1. Wavelet Transform

Mallat [32] proposed the concept of multiresolution analysis when constructing orthogonal wavelet basis, illustrated the multiresolution characteristics of wavelet in geometric space, and deduced the fast algorithm for the decomposition and reconstruction of wavelet transform. The basic idea is to transform the signal ft into two orthogonal subsignals, ajf and djf, with scale j. Among which, ajf, the approximation signal, is corresponding to the low-frequency component, while djf, the detail signal, is for the high one. Next, transform the subsignal ajf into signals aj+1f and dj+1f with scale j+1. Thus, the signal is decomposed and reconstructed with multiresolution and multiscale analysis. The decomposition and reconstruction of wavelet with scale 3 are shown in Figure 3.Figure 3
Decomposition and reconstruction of wavelet with scale 3.


#### 4.1.2. AOK Time-Frequency Analysis

Baraniuk and Jones [33] proposed the idea of adaptive optimal kernel (AOK) time-frequency analysis. The main idea is to select appropriate short-time ambiguity function and adaptive time-varying kernel function. Thus, more details of subsignals can be expressed. AOK time-frequency distribution is defined as follows:(13)Pt,f=12π∫−∞∞∫−∞∞At;θ;τΦoptt;θ;τe−jθτ−jftdθdτ,where Φoptt;θ;τ is the optimal kernel function; At;θ;τ is the short-time ambiguity function; θ is the angle function in polar coordinates; τ is the time interval; and f is the frequency.The definition of the short-time ambiguity functionAt;θ;τ is as follows:(14)At;θ;τ=∫−∞∞hu−τ2wu−t−τ2hu+τ2fu−t+τ2ejθudu,where wu is the symmetric window function and t is the center of wu.Provided thatu>T,wu=0, the kernel function of the signal just can be calculated out in the range of t−T,t+T, while the short-time ambiguity function can express the details accurately. The value of θoptt;θ;τ can be readily calculated out with the help of the short-time ambiguity function. As At;θ;τ changes with time, θoptt;θ;τ will be alterable too. θoptt;θ;τ can be obtained by solving the following optimization function:(15)maxΦ∫02π∫0+∞At;t;φΦt;t;φ2rdrdφ.The constraint conditions are shown as follows:(16)Φt;t;φ=exp−r22σ2φ,12π∫02π∫0+∞Φt;t;φ2rdrdφ=12π∫02πσ2φdφ≤α,where σφ is the extension function; φ is the angle between the radial and horizontal direction, φ=arctanτ/θ and r=θ2+τ21/2; and α is the energy volume of radial gaussian kernel. Some self-components of the kernel functions will be filtered out with too small α, while the cross components cannot be effectively removed with overlarge α, and the rational value of α is [1, 5] [34].


#### 4.1.3. Calculation of SWE

When the seismic signals are decompose, the number of layers should be determined by the specific signals and the frequency bands of the blasting vibration meters. Assume that the signal is decomposed to the layern, namely, there are total n + 1 sets of frequency bands. The signals with random frequency bands are analysed by the AOK time-frequency analysis method, and Pt,f in certain moment and certain frequency can be obtained. Thus, the SWE can be solved by the following equation:(17)En=∫−∞∞∫−∞∞Pt,fdtdf=∑t=1n∑f=1met,f,where En is the SWE corresponding to the nth frequency band.The proportion of SWE corresponding to each band can be expressed asrn, as follows:(18)rn=EnE×100%.


### 4.2. Characteristics of SWE

In terms of the sampling features of the measured seismic wave, the sampling interval is set at 4000sps. According to the Nyquist theorem, the Nyquist frequency of the measured signal is 2000 Hz. The measured signals are decomposed into 8 layers with the help of wavelet transform, and the bandwidth of the eighth low-frequency wavelets (a8) is 0∼7.8125 Hz. The bandwidth of each band is depicted in Table5.Table 5
Bandwidth distribution of measured seismic wave.
Band
a8
d8
d7
d6
d5
d4
d3
d2
d1
Bandwidth (Hz)
0∼7.8125
7.8125∼15.625
15.625∼31.25
31.25∼62.5
62.5∼125
125∼250
250∼500
500∼1000
1000∼2000When adopting db8 wavelet base function to decompose and reconstruct the measured seismic wave, the relative error between the reconstructed signal and the measured signal is the minimal one. Thus, db8 wavelet base function is selected to decompose and reconstruct all the measured signals in this study, and the relative error is about 1.02 × 10−11∼2.07 × 10−11, which can fully meet the requirements of the practical project. The reconstructed signal of each band in vertical direction is shown in Figure 4.Figure 4
Reconstructed signal of each band in vertical direction.The time-frequency characteristics and energy density of each reconstructed signal can be obtained with the help of AOK time-frequency analysis, as shown in Figure5. As depicted in the time-frequency spectrum, the major centralized area can be clearly distinguished. The centralized time is about 0.1 s, and the frequency mainly ranges from 0 to 125 Hz. Comparing all the time-frequency characteristics and energy density blow, it is found that the energies of low frequency bands (d5∼d8 and a8) are huge, while those of high-frequency bands (d1∼d4) are tiny. The durations of d1 and d2 occupy more than half of time, and the energies are extremely low, which denotes that the bands are mainly signal noise in the whole process of blasting.Figure 5
Time-frequency spectrum for each reconstructed signal: (a) as1, (b) a8, (c) d8, (d) d7, (e) d6, (f) d5, (g) d4, (h) d3, (i) d2, and (j) d1.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)According to equations (13)∼(16), the energies for reconstructed signals in three directions, the total energy, and the energy proportions can be calculated out (Table 6). The energies of reconstructed signals in each band are depicted in Figure 6.Table 6
Energy components, total energy, and energy proportions of reconstructed signals.
Frequency band
Horizontal tangential
Horizontal radial
Vertical direction
E
n (cm2/s2)
P
n (%)
E
n (cm2/s2)
P
n (%)
E
n (cm2/s2)
P
n (%)
a8
407.193
40.275
1.9167
1.566
9.520
6.333
d8
358.930
35.501
5.555
4.539
38.802
25.814
d7
146.421
14.482
31.778
25.964
41.292
27.471
d6
57.142
5.652
38.028
31.070
36.239
24.108
d5
36.035
3.564
31.078
25.392
22.797
15.166
d4
4.244
0.420
8.192
6.693
1.340
0.892
d3
0.788
0.0780
3.889
3.178
0.321
0.213
d2
0.176
0.0170
1.536
1.255
0.00242
0.00147
d1
0.111
0.0110
0.419
0.343
0.00251
0.00153
E (cm2/s2)
1011.039
122.392
150.317
E
t (cm2/s2)
1283.748Figure 6
Energies of reconstructed signals in each band.As seen from Table5 and Figure 6, the energy proportions of seismic wave in the bands of d5∼d8 and a8 are totally 85%, while energies in other bands are extremely low, which indicates that the dominant frequency ranges from 0 to 125 Hz. As the high-frequency bands of d1 and d2 are mainly noise signals, which can be left out when studying SWE. The bands of d3 and d4 are in the range of 125∼500 Hz, which are higher than the natural vibration frequencies of nearby structures.The frequency bands in the range of 0∼125 Hz are the key to study SWE. However, there are significant differences in the SWE of different directions. The value of horizontal tangential energy is the maximum, and horizontal radial energy is the minimum. The range of dominant frequency of horizontal tangential, horizontal radial, and vertical direction are, respectively, 0∼31.25 Hz, 15.625∼125 Hz, and 7.8125∼62.5 Hz.


### 4.3. Attenuation Law of SWE

5 sets of monitoring data with different distances from the blast-face to the monitoring points (65 m, 75 m, 85 m, 105 m, and 125 m) in the same blasting are selected to analyse the characteristics of SWE. The energy components of the 5 monitoring points in 3 directions, horizontal tangential, horizontal radial, and vertical direction, as well as the total energy are calculated, as depicted in Table7. The attenuation laws of energy components along with distance can be obtained from Figure 7.Table 7
SWE in different distances from the blast-face to the monitoring points.
Distance from the measured point to the blast-face
Energy componentsE (cm2/s2)
Total energyEt (cm2/s2)
Horizontal tangential
Horizontal radial
Vertical direction
65
1011.039
122.392
150.317
1283.748
75
302.307
50.321
58.712
411.34
85
115.844
25.273
27.135
168.252
105
19.398
8.5198
9.275
37.192
125
10.485
7.047
6.425
23.957Figure 7
Attenuation laws of energy components along with distance.As seen from Table6 and Figure 7, the attenuation laws of SWE along with distance in each direction are roughly the same. At the same distance from the measured point to the blast-face, horizontal tangential energy is the minimal one and attenuates rapidly, while horizontal radial energy is the largest one and changes relatively slowly. All the seismic wave energies present typical forms of exponential decay, that is to say, when near the blast-face, the energies attenuate rapidly, and with the increase of distance, the attenuation rates are obviously reduced. When it comes to zones far away from the blast-face, the energies attenuate extremely slowly. Fit the total energy of seismic wave with the help of exponential function, and the attenuation formula is shown in the following equation:(19)Et=1.651×1017R−7.783,r2=0.994.Five major bands, with frequency of 0∼125 Hz, are employed to study the attenuation law of SWE. The attenuation laws of the total energy and energy pro portion of each band are, respectively, depicted in Figures8 and 9.Figure 8
Attenuation laws of the total energies along with distance.Figure 9
Attenuation laws of the energy proportions along with distance.As seen from the above figures, the total SWE of each band all experiences exponential attenuation with the increase of distance. Among which, the total energies of the bands a8 and d8 present the slowest attenuation, which indicates that the lower the frequency is, the more gentle the energy attenuation is. With the increase of distance, the energy proportions of the bands d5∼d7 constantly decrease, while those of the bands a8 and d8 increase gradually. As the distance from the measured point to the blast-face is 125 m, the total energy proportion of the bands in the frequency of 0∼12.625 Hz is up to 95.393%.


## 4.1. Time-Frequency Analysis

### 4.1.1. Wavelet Transform

Mallat [32] proposed the concept of multiresolution analysis when constructing orthogonal wavelet basis, illustrated the multiresolution characteristics of wavelet in geometric space, and deduced the fast algorithm for the decomposition and reconstruction of wavelet transform. The basic idea is to transform the signal ft into two orthogonal subsignals, ajf and djf, with scale j. Among which, ajf, the approximation signal, is corresponding to the low-frequency component, while djf, the detail signal, is for the high one. Next, transform the subsignal ajf into signals aj+1f and dj+1f with scale j+1. Thus, the signal is decomposed and reconstructed with multiresolution and multiscale analysis. The decomposition and reconstruction of wavelet with scale 3 are shown in Figure 3.Figure 3
Decomposition and reconstruction of wavelet with scale 3.


### 4.1.2. AOK Time-Frequency Analysis

Baraniuk and Jones [33] proposed the idea of adaptive optimal kernel (AOK) time-frequency analysis. The main idea is to select appropriate short-time ambiguity function and adaptive time-varying kernel function. Thus, more details of subsignals can be expressed. AOK time-frequency distribution is defined as follows:(13)Pt,f=12π∫−∞∞∫−∞∞At;θ;τΦoptt;θ;τe−jθτ−jftdθdτ,where Φoptt;θ;τ is the optimal kernel function; At;θ;τ is the short-time ambiguity function; θ is the angle function in polar coordinates; τ is the time interval; and f is the frequency.The definition of the short-time ambiguity functionAt;θ;τ is as follows:(14)At;θ;τ=∫−∞∞hu−τ2wu−t−τ2hu+τ2fu−t+τ2ejθudu,where wu is the symmetric window function and t is the center of wu.Provided thatu>T,wu=0, the kernel function of the signal just can be calculated out in the range of t−T,t+T, while the short-time ambiguity function can express the details accurately. The value of θoptt;θ;τ can be readily calculated out with the help of the short-time ambiguity function. As At;θ;τ changes with time, θoptt;θ;τ will be alterable too. θoptt;θ;τ can be obtained by solving the following optimization function:(15)maxΦ∫02π∫0+∞At;t;φΦt;t;φ2rdrdφ.The constraint conditions are shown as follows:(16)Φt;t;φ=exp−r22σ2φ,12π∫02π∫0+∞Φt;t;φ2rdrdφ=12π∫02πσ2φdφ≤α,where σφ is the extension function; φ is the angle between the radial and horizontal direction, φ=arctanτ/θ and r=θ2+τ21/2; and α is the energy volume of radial gaussian kernel. Some self-components of the kernel functions will be filtered out with too small α, while the cross components cannot be effectively removed with overlarge α, and the rational value of α is [1, 5] [34].


### 4.1.3. Calculation of SWE

When the seismic signals are decompose, the number of layers should be determined by the specific signals and the frequency bands of the blasting vibration meters. Assume that the signal is decomposed to the layern, namely, there are total n + 1 sets of frequency bands. The signals with random frequency bands are analysed by the AOK time-frequency analysis method, and Pt,f in certain moment and certain frequency can be obtained. Thus, the SWE can be solved by the following equation:(17)En=∫−∞∞∫−∞∞Pt,fdtdf=∑t=1n∑f=1met,f,where En is the SWE corresponding to the nth frequency band.The proportion of SWE corresponding to each band can be expressed asrn, as follows:(18)rn=EnE×100%.


## 4.1.1. Wavelet Transform

Mallat [32] proposed the concept of multiresolution analysis when constructing orthogonal wavelet basis, illustrated the multiresolution characteristics of wavelet in geometric space, and deduced the fast algorithm for the decomposition and reconstruction of wavelet transform. The basic idea is to transform the signal ft into two orthogonal subsignals, ajf and djf, with scale j. Among which, ajf, the approximation signal, is corresponding to the low-frequency component, while djf, the detail signal, is for the high one. Next, transform the subsignal ajf into signals aj+1f and dj+1f with scale j+1. Thus, the signal is decomposed and reconstructed with multiresolution and multiscale analysis. The decomposition and reconstruction of wavelet with scale 3 are shown in Figure 3.Figure 3
Decomposition and reconstruction of wavelet with scale 3.


## 4.1.2. AOK Time-Frequency Analysis

Baraniuk and Jones [33] proposed the idea of adaptive optimal kernel (AOK) time-frequency analysis. The main idea is to select appropriate short-time ambiguity function and adaptive time-varying kernel function. Thus, more details of subsignals can be expressed. AOK time-frequency distribution is defined as follows:(13)Pt,f=12π∫−∞∞∫−∞∞At;θ;τΦoptt;θ;τe−jθτ−jftdθdτ,where Φoptt;θ;τ is the optimal kernel function; At;θ;τ is the short-time ambiguity function; θ is the angle function in polar coordinates; τ is the time interval; and f is the frequency.The definition of the short-time ambiguity functionAt;θ;τ is as follows:(14)At;θ;τ=∫−∞∞hu−τ2wu−t−τ2hu+τ2fu−t+τ2ejθudu,where wu is the symmetric window function and t is the center of wu.Provided thatu>T,wu=0, the kernel function of the signal just can be calculated out in the range of t−T,t+T, while the short-time ambiguity function can express the details accurately. The value of θoptt;θ;τ can be readily calculated out with the help of the short-time ambiguity function. As At;θ;τ changes with time, θoptt;θ;τ will be alterable too. θoptt;θ;τ can be obtained by solving the following optimization function:(15)maxΦ∫02π∫0+∞At;t;φΦt;t;φ2rdrdφ.The constraint conditions are shown as follows:(16)Φt;t;φ=exp−r22σ2φ,12π∫02π∫0+∞Φt;t;φ2rdrdφ=12π∫02πσ2φdφ≤α,where σφ is the extension function; φ is the angle between the radial and horizontal direction, φ=arctanτ/θ and r=θ2+τ21/2; and α is the energy volume of radial gaussian kernel. Some self-components of the kernel functions will be filtered out with too small α, while the cross components cannot be effectively removed with overlarge α, and the rational value of α is [1, 5] [34].


## 4.1.3. Calculation of SWE

When the seismic signals are decompose, the number of layers should be determined by the specific signals and the frequency bands of the blasting vibration meters. Assume that the signal is decomposed to the layern, namely, there are total n + 1 sets of frequency bands. The signals with random frequency bands are analysed by the AOK time-frequency analysis method, and Pt,f in certain moment and certain frequency can be obtained. Thus, the SWE can be solved by the following equation:(17)En=∫−∞∞∫−∞∞Pt,fdtdf=∑t=1n∑f=1met,f,where En is the SWE corresponding to the nth frequency band.The proportion of SWE corresponding to each band can be expressed asrn, as follows:(18)rn=EnE×100%.


## 4.2. Characteristics of SWE

In terms of the sampling features of the measured seismic wave, the sampling interval is set at 4000sps. According to the Nyquist theorem, the Nyquist frequency of the measured signal is 2000 Hz. The measured signals are decomposed into 8 layers with the help of wavelet transform, and the bandwidth of the eighth low-frequency wavelets (a8) is 0∼7.8125 Hz. The bandwidth of each band is depicted in Table5.Table 5
Bandwidth distribution of measured seismic wave.
Band
a8
d8
d7
d6
d5
d4
d3
d2
d1
Bandwidth (Hz)
0∼7.8125
7.8125∼15.625
15.625∼31.25
31.25∼62.5
62.5∼125
125∼250
250∼500
500∼1000
1000∼2000When adopting db8 wavelet base function to decompose and reconstruct the measured seismic wave, the relative error between the reconstructed signal and the measured signal is the minimal one. Thus, db8 wavelet base function is selected to decompose and reconstruct all the measured signals in this study, and the relative error is about 1.02 × 10−11∼2.07 × 10−11, which can fully meet the requirements of the practical project. The reconstructed signal of each band in vertical direction is shown in Figure 4.Figure 4
Reconstructed signal of each band in vertical direction.The time-frequency characteristics and energy density of each reconstructed signal can be obtained with the help of AOK time-frequency analysis, as shown in Figure5. As depicted in the time-frequency spectrum, the major centralized area can be clearly distinguished. The centralized time is about 0.1 s, and the frequency mainly ranges from 0 to 125 Hz. Comparing all the time-frequency characteristics and energy density blow, it is found that the energies of low frequency bands (d5∼d8 and a8) are huge, while those of high-frequency bands (d1∼d4) are tiny. The durations of d1 and d2 occupy more than half of time, and the energies are extremely low, which denotes that the bands are mainly signal noise in the whole process of blasting.Figure 5
Time-frequency spectrum for each reconstructed signal: (a) as1, (b) a8, (c) d8, (d) d7, (e) d6, (f) d5, (g) d4, (h) d3, (i) d2, and (j) d1.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)According to equations (13)∼(16), the energies for reconstructed signals in three directions, the total energy, and the energy proportions can be calculated out (Table 6). The energies of reconstructed signals in each band are depicted in Figure 6.Table 6
Energy components, total energy, and energy proportions of reconstructed signals.
Frequency band
Horizontal tangential
Horizontal radial
Vertical direction
E
n (cm2/s2)
P
n (%)
E
n (cm2/s2)
P
n (%)
E
n (cm2/s2)
P
n (%)
a8
407.193
40.275
1.9167
1.566
9.520
6.333
d8
358.930
35.501
5.555
4.539
38.802
25.814
d7
146.421
14.482
31.778
25.964
41.292
27.471
d6
57.142
5.652
38.028
31.070
36.239
24.108
d5
36.035
3.564
31.078
25.392
22.797
15.166
d4
4.244
0.420
8.192
6.693
1.340
0.892
d3
0.788
0.0780
3.889
3.178
0.321
0.213
d2
0.176
0.0170
1.536
1.255
0.00242
0.00147
d1
0.111
0.0110
0.419
0.343
0.00251
0.00153
E (cm2/s2)
1011.039
122.392
150.317
E
t (cm2/s2)
1283.748Figure 6
Energies of reconstructed signals in each band.As seen from Table5 and Figure 6, the energy proportions of seismic wave in the bands of d5∼d8 and a8 are totally 85%, while energies in other bands are extremely low, which indicates that the dominant frequency ranges from 0 to 125 Hz. As the high-frequency bands of d1 and d2 are mainly noise signals, which can be left out when studying SWE. The bands of d3 and d4 are in the range of 125∼500 Hz, which are higher than the natural vibration frequencies of nearby structures.The frequency bands in the range of 0∼125 Hz are the key to study SWE. However, there are significant differences in the SWE of different directions. The value of horizontal tangential energy is the maximum, and horizontal radial energy is the minimum. The range of dominant frequency of horizontal tangential, horizontal radial, and vertical direction are, respectively, 0∼31.25 Hz, 15.625∼125 Hz, and 7.8125∼62.5 Hz.


## 4.3. Attenuation Law of SWE

5 sets of monitoring data with different distances from the blast-face to the monitoring points (65 m, 75 m, 85 m, 105 m, and 125 m) in the same blasting are selected to analyse the characteristics of SWE. The energy components of the 5 monitoring points in 3 directions, horizontal tangential, horizontal radial, and vertical direction, as well as the total energy are calculated, as depicted in Table7. The attenuation laws of energy components along with distance can be obtained from Figure 7.Table 7
SWE in different distances from the blast-face to the monitoring points.
Distance from the measured point to the blast-face
Energy componentsE (cm2/s2)
Total energyEt (cm2/s2)
Horizontal tangential
Horizontal radial
Vertical direction
65
1011.039
122.392
150.317
1283.748
75
302.307
50.321
58.712
411.34
85
115.844
25.273
27.135
168.252
105
19.398
8.5198
9.275
37.192
125
10.485
7.047
6.425
23.957Figure 7
Attenuation laws of energy components along with distance.As seen from Table6 and Figure 7, the attenuation laws of SWE along with distance in each direction are roughly the same. At the same distance from the measured point to the blast-face, horizontal tangential energy is the minimal one and attenuates rapidly, while horizontal radial energy is the largest one and changes relatively slowly. All the seismic wave energies present typical forms of exponential decay, that is to say, when near the blast-face, the energies attenuate rapidly, and with the increase of distance, the attenuation rates are obviously reduced. When it comes to zones far away from the blast-face, the energies attenuate extremely slowly. Fit the total energy of seismic wave with the help of exponential function, and the attenuation formula is shown in the following equation:(19)Et=1.651×1017R−7.783,r2=0.994.Five major bands, with frequency of 0∼125 Hz, are employed to study the attenuation law of SWE. The attenuation laws of the total energy and energy pro portion of each band are, respectively, depicted in Figures8 and 9.Figure 8
Attenuation laws of the total energies along with distance.Figure 9
Attenuation laws of the energy proportions along with distance.As seen from the above figures, the total SWE of each band all experiences exponential attenuation with the increase of distance. Among which, the total energies of the bands a8 and d8 present the slowest attenuation, which indicates that the lower the frequency is, the more gentle the energy attenuation is. With the increase of distance, the energy proportions of the bands d5∼d7 constantly decrease, while those of the bands a8 and d8 increase gradually. As the distance from the measured point to the blast-face is 125 m, the total energy proportion of the bands in the frequency of 0∼12.625 Hz is up to 95.393%.


## 5. Prediction of SWE and Discussion

Thirty groups of measured signals are employed to analyse the characteristics of SWE, and the total energy of each bandEt is fitted by equation (12). And, the prediction formulas are adopted to fit the measured PPVs. The results are shown in Table 8.Table 8
Fitting results of different prediction formulas.
Prediction formula
Parameters
R
2
k
α
β
USBM
4.870
1.077
—
0.828
Langefors–Kihlstrom
0.413
2.413
—
0.692
Bureau of Indian Standard
0.413
1.207
—
0.702
Sadaovsk formula
123.016
1.467
—
0.842
CCTS
107.06
1.572
−0.223
0.875
Equation (12)
6.148
5.951
0.211
0.911The relationship between the measured total energyEt and the predicted one is depicted in Figure 10. Figures 11–15 demonstrate the deviations of the measured PPVs when compared with the predicted PPVs fitted by USBM, Langefors–Kihlstrom, Bureau of Indian Standard, Sadaovsk formula, and CCTS. To check the capacity performance of the prediction models, values of mean absolute percent error (MAPE) are solved via the following equation:(20)MAPE=1n∑i=1ny−y′y×100,where y and y′ are the measured and predicted values and n is the total number of the data. When the MAPE value is zero, the model’s performance is prefect.Figure 10
Measured energy against predicted energy fitted by energy prediction formula.Figure 11
Measured PPVs against predicted PPVs fitted by USBM.Figure 12
Measured PPVs against predicted PPVs fitted by Langefors–Kihlstrom.Figure 13
Measured PPVs against predicted PPVs fitted by the Bureau of Indian Standard.Figure 14
Measured PPVs against predicted PPVs fitted by Sadaovsk formula.Figure 15
Measured PPVs against predicted PPVs fitted by CCTS.The performance indices of all the mentioned models are shown in Table9. In the prediction of PPVs, the MAPE values obtained by USBM, Langefors–Kihlstrom, Bureau of Indian Standard, Sadaovsk formula, and CCTS [4–7, 9] are, respectively, 18.459%, 27.426%, 27.435%, 15.897%, and 13.566%, which indicates that the accuracy of CCTS in predicting PPV is better than that of other formulas.Table 9
Performance indices of all prediction models.
Predictive model
R
2
MAPE (%)
Energy prediction formula
0.96485
11.154
USBM [4]
0.83419
18.459
Langefors and Kihlstrom [5]
0.70283
27.426
Bureau of Indian Standard [6]
0.70281
27.435
Sadaovsk formula [7]
0.88291
15.897
CCTS [9]
0.91139
13.566The performance indicesR2 and MAPE of the energy prediction formula deduced in this study are, respectively, 0.96485 and 11.154%, which are superior to the results of other formulas. That is to say, the energy prediction formula possesses excellent performance capacity and favorable prediction accuracy. It can be concluded that adopting SWE as the judgment standard of seismic intensity and fitting the energy with the help of the prediction formula in this study can achieve more accurate prediction results, which provides solid basis for the safety control to the harmful effects of seismic waves.


## 6. Conclusion

In the blasting engineering, PPV has been the major index to describe the intensity of seismic waves. Nevertheless, many countries’ safety criterion ignored the effects of duration time and vibration frequency. The paper regards the seismic wave energy as the index to quantify seismic intensity, puts forward the prediction formula of seismic wave energy, and obtains the distribution characteristics and attenuation law of seismic wave energy. Conclusions can be drawn as follows:(1)
The energy proportions of seismic wave induced by underwater drilling and blasting in the bands of d5∼d8 and a8 are totally 85%, while energies in high-frequency bands are extremely low. The range of dominant frequency of horizontal tangential, horizontal radial, and vertical direction are, respectively, 0∼31.25 Hz, 15.625∼125 Hz, and 7.8125∼62.5 Hz.(2)
All the seismic wave energies present typical forms of exponential decay, that is to say, when near the blast-face, the energies attenuate rapidly, and with the increase of distance, the attenuation rates are obviously reduced. When it comes to zones far away from the blast-face, the energies attenuate extremely slowly.(3)
Comparing with other prediction formulas of PPV, adopting seismic wave energy as the judgment standard of seismic intensity and fitting the energy with the help of the prediction formula in this study can achieve more accurate prediction results, which provides solid basis for the safety control to the harmful effects of seismic waves.



---

*Source: 4367698-2019-09-26.xml*