Analysis and Research on Blasting Network Delay of Deep-Buried Diversion Tunnel Crossing Fault Zone Based on EP-CEEMDAN-INHT

The Empirical Mode Decomposition (EMD) of blasting seismic wave monitoring signal with noise can get IMFs with serious modal confusion and endpoint divergence. The Hilbert transform is constrained by the Bedrosian theorem, when it dealing with such IMFs will get negative instantaneous frequency, which leads to serious error in the identification of non-electric millisecond detonator initiation delay. EP-CEEMDAN-INHT is proposed and applied to the delay analysis of blasting network of deep buried diversion tunnel crossing fault zone. Comparing EP-CEEMDAN-INHT with EMD-HHT, it is found that EP-CEEMDAN-INHT can clearly display the time–frequency information contained in the measured blasting vibration signal, and EMD modal confusion and endpoint effect are well suppressed. The actual millisecond time interval obtained by EP-CEEMDAN-INHT can judge whether the detonator is in normal service. At the same time, the blasting millisecond interval with the best damping effect is 54.51–59.75 ms, which can realize the optimization of blasting network and has important practical significance for blasting safety control.


Introduction
The time delay error caused by the detonation of nonelectric millisecond detonator will gradually accumulate with the development of the blasting process, resulting in a big discrepancy between the actual delay and the theoretical design delay (Gong et al. 2015;Xu et al. 2013). The millisecond delay initiation network with a reasonable design can partially offset interference of the blasting seismic waves, thereby reducing the damage of the surrounding rocks (Tian et al. 2012). Therefore, it is very important to set a reasonable time interval for deep tunnels with complex geological conditions and easy to cause water and mud gushing during construction (Peng et al. 2020(Peng et al. , 2021. The identification of delay time interval of non-electric millisecond detonator in actual construction is helpful to test the reliability of detonator used in blasting and the rationality of initiation network design, and to optimize blasting design parameters (Wang 2007;Zhou 2014).
At present, a large number of scholars have researched identification method of blasting network delay, and the main researches are concentrated in two aspects. Wavelet analysis and traditional EMD-HHT (Huang 1998), in which wavelet analysis overly relies on prior basis functions (Li et al. 2005). The traditional EMD-HHT has problems such as endpoint effect and modal confusion , which affects the analysis accuracy to a certain extent. Its improved algorithm EEMD-HHT can suppress EMD modal confusion to a certain extent, and cannot prove that the artificially introduced white noise has been completely eliminated, the completeness and authenticity of signal is affected (Cheng et al. 2018).
Aiming at this phenomenon, this paper introduces the EP-CEEMDAN-INHT algorithm, which can simultaneously suppress the endpoint effect and modal confusion existing in EMD and improve the accuracy of blasting network delay recognition. EP-CEEM-DAN-INHT is completed in three steps, The first step is to perform endpoint processing(EP) on the blasting seismic wave monitoring signals obtained from actual construction; The second step performs the Complete Ensemble Empirical Mode Decomposition with Adaptive Noise(CEEMDAN) on the signal of the ''first step'';The third step is to obtain the normalized intrinsic mode function(IMF) through the Improved Normalized Hilbert Transform(INHT). The envelope of IMF with the largest energy proportion is extracted to obtain the high-precision delay identification result of blasting network. According to the identification result, the signal separation technology is used to obtain the optimal millisecond time interval of damping effect, so as to realize the delay optimization of blasting network.

The Principle of Endpoint Processing
The principle of Intrinsic Mode Function (IMF) generation is to calculate the local mean value of the signal based on the upper and lower envelopes. The upper (lower) envelope is obtained from the local maximum (small) value of the signal through cubic spline interpolation. However, the endpoints of the signal cannot be at the maximum or minimum at the same time, and there are situations where the endpoints are not extreme points. In these two cases, cubic spline interpolation cannot be performed, so the envelope will diverge at the endpoints (Wu and Riemenschneider 2014). The authenticity of the IMF is affected.
Through waveform matching, a wavelet with the highest matching degree with the endpoint data is found in the signal, and the wavelet is moved to the endpoint for decomposition, which can realize endpoint processing. Take the processing of the left end of the signal as an example. First, find the time when all the maximum points of the signal occur, denoted as t max1 , t max2 ,…, t maxi (i = 1,2,3…M), and the corresponding amplitude is denoted as x max1 , x max2 ,…, x maxi (i = 1,2,3…M), in the same way, find the moment when all the minimum points of the signal occur, denoted as t min1 , t min2 ,…, t mini (i = 1,2,3…N), the corresponding amplitudes are denoted as y min1 , y min2 ,…, y mini (i = 1,2,3…N), the specific steps are as follows.
Step 1: Take a time series including (t max1 , x max1 ) and (t min1 , y min1 ) and record it as l 0 , where l 0 contains k sampling points; Step 2: Divide the original signal into n time series with k sampling points, and record them as l 1 l 2 l 3 …l i …l n (1 B i B n); Step 3: According to Eq. (1), choose the l i corresponding to the minimum value of n, and the matching degree of l i and l 0 is the highest; Step 4: Shift the time series of 2 k sampling points composed of l i and l i-1 to the position of l 0 to realize the left Endpoint Processing (EP).

The Principle of CEEMDAN
Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) adds finite times of adaptive white noise in each stage of decomposition, which can achieve almost zero reconstruction error with less average times (Maria et al. 2011;Zhang et al. 2020;Peng et al. 2021). The specific steps are as follows: Step 1: Add the adaptive white noise B j (t) to the signal S(t) after the endpoint processing in the second section, where j is the number of times to add the noise, generally 50-100 times (Sun et al. 2020), 60 times in this paper. Then the signal of the jth order can be expressed as S(t) = S(t) ? a j B j (t)(j = 1,2,3…60), Where a is the standard deviation of the jth white noise addition. For the IMF1 obtained by CEEMDAN, see Eq.
Step 2: Constructing a new signal S(t) = S(t) ? a j-B j (t) and decomposing for 60 times. IMF2 is obtained, and the remainder is R 2 (t) = S(t) -IMF2.
Step 3: Repeat ''Step1'' and ''Step2'' until the end of the program. S(t) gets c IMF and the unique remainder R through CEEMDAN. The decomposition result is shown in Eq. (3).

The Principle of INHT
The traditional Hilbert transform is constrained by the Bedrosian theorem and can get negative instantaneous frequency when dealing with non-stationary signals.
In view of this phenomenon, Huang et al. 2009 proposed Normalized Hilbert Transform (NHT) to improve the calculation accuracy of instantaneous frequency. In this paper, we make some improvements on the basis of Huang's NHT. Through reference (Huang et al. 2009), we can find that Huang's IMF for NHT comes from EMD or EEMD, but EMD and EEMD inevitably produce the problems of modal confusion and endpoint effect in signal decomposition. In this paper, the first step of IMF is through the endpoint processing, and the second step is obtained through the EMD improved algorithm CEEMDAN. The IMF obtained through the above two steps can suppress the endpoint effect and modal confusion at the same time. The specific operation is as follows: Step 1: Take the absolute value of IMF obtained by EP-CEEMDAN, and find out all the maximum values in |IMFi|.
Step 2: Get the spline envelope of the maximum point of |IMFi| by cubic spline interpolation, and record it as x (t); Step 3: Normalization processing, take IMF1 as an example, record the spline envelope of all maximum points in |IMF1| as x 1 (t), and calculate f 1 (t) = IMF1/ x 1 (t); Step 4: if all | f 1 (t)| satisfy | f 1 (t)|B 1, stop. On the contrary, re assign the value of IMF1, make IMF1 = f 1 (t), repeat '' Step 2'' to get the spline envelope of all the maximum points in | f 1 (t)| and record it as x 2 (t), repeat ''the third step'' to get f 2 (t) = f 1 (t)/ x 2 (t), and check whether | f 2 (t)| meets | f 2 (t)|B 1, see Eq. (4) for details, where d is the number of repetitions, and generally 2-3 times can meet the demand; In Eq. (4), f d (t) is the frequency modulation part of IMF1, and its amplitude modulation part w d (t) can be expressed by Eq. (5).
Therefore, the normalized IMF1 can be expressed by Eq. (6). It is not difficult to find that the essence of normalization is to separate the frequency modulation and amplitude modulation components of IMF.
The phase function h(t) of IMF1 is obtained by using arctangent function. See Eq. (8).
The instantaneous frequency x(t) can be obtained by deriving h(t). The expression of x(t) is shown in Eq. (9).
Through five steps, the Improved Normalized Hilbert Transform (INHT) is realized. At this time, the instantaneous frequency is obtained by the frequency modulation component of IMF. After normalization, the Hilbert transform is no longer limited by Bedrosian theorem, and the obtained instantaneous frequency is more authentic.

Application of EP-CEEMDAN-INHT in Delay Analysis of Blasting Network for Deep Buried Diversion Tunnel Crossing Fault Zone
Taking the blasting excavation project of a diversion tunnel in Fujian Province as the research object, the total length of the diversion tunnel is 13.842 km, and the main geological structure of the area is fault. According to the survey report, there are more than 20 faults along the tunnel, mainly tensile faults, with a width of 2-6 m and a wide fracture zone on both sides. The geological conditions are complex and changeable, and there are several major technical problems in the construction process, such as water gushing, geothermal, rock burst and so on. Among them, large-scale continuous water inrush is the biggest technical difficulty. The diversion tunnel is formed by smooth millisecond blasting with non-electric millisecond detonator. The actual delay time of non-electric millisecond detonator blasting is calculated by EP-CEEMDAN-INHT to judge whether the performance of the corresponding batch of detonators meets the design requirements. According to the analysis results, the blasting network is optimized, and the millisecond time with the best damping effect is calculated according to the principle of waveform decomposition and superposition, so as to minimize the damage of blasting seismic wave to surrounding rock and reduce the risk of water gushing and mud bursting due to blasting construction. Figure 1 is the diagram of the blast holes layout of the diversion tunnel. By observing Fig. 1, it can be found that the blasting is carried out in five sections. The TC-4850 intelligent blasting vibration meter is used to arrange the measuring points along the tunnel axis to avoid the damage of the instrument by flying stones. The No.1 measuring point is set 50 m away from the blasting source, and the remaining four points are separated by 5 m, 10 m, 20 m and 40 m in turn. The typical blasting seismic wave signals monitored by TC-4850 intelligent blasting vibration meter at the blasting site are shown in Fig. 2.
In the actual construction, multiple detonations will be designed according to the specific working conditions. When each detonator detonates, it is bound to produce a certain range of mutation on its time history curve, which can be considered as the superposition of energy. The envelope of typical IMF components is solved by INHT, and the time interval corresponding to each mutation peak is calculated, which is the actual network delay time parameter.
Considering that the signal decomposed by EP-CEEMDAN algorithm will get multiple IMF components, each IMF carries certain time-frequency energy information of blasting seismic wave signal. The IMF component with the highest energy proportion can reflect the detailed information of time-frequency energy contained in the blasting seismic wave monitoring signal to the greatest extent. The amplitude envelope curve of the component is extracted by INHT. The time node corresponding to the envelope peak point represents the superposition of the energy of each section of the blasting network, and also Therefore, the time history curve in Fig. 2 is decomposed by EP-CEEMDAN algorithm to obtain the IMF component as shown in Fig. 3. By observing Fig. 3, the IMF obtained by EP-CEEMDAN algorithm. it can be found that the IMF components are arranged from high frequency to low frequency, and the phenomena of EMD mode confusion and endpoint divergence are effectively suppressed. Then, INHT is performed for each IMF component to calculate the actual IMF energy and its proportion in the total signal energy. The energy of each IMF component and the ratio of each IMF to the total energy are listed in Table 1.
In order to highlight that EP-CEEMDAN-INHT algorithm can effectively improve the extraction accuracy of blasting seismic wave time-frequency parameters, and realize the suppression of EMD endpoint effect and modal confusion, EMD is performed on the signal in Fig. 2, and the results shown in Fig. 4. By observing Fig. 4, it can be found that IMF1 and IMF2 are the noises that cannot be removed in the detection. IMF3, IMF4, IMF5 and IMF6 are the key frequency bands. The high frequency mode of IMF3 is confused seriously. IMF4 is relatively stable, but there is modal confusion around 0.35 s. The left end of IMF5 diverges. IMF6 tends to low frequency in 0.3-0.55 s. The right end of IMF7 diverges. Mode splitting occurs in IMF8 and IMF9.
Compared with Figs. 3 and 4, the following conclusions can be drawn: the results of EP-CEEMDAN algorithm can clearly show the signal frequency information contained in the measured blasting vibration data, clearly distinguish the high frequency, medium frequency and low frequency, the modal confusion caused by noise signal and the endpoint effect of the algorithm itself are well suppressed. The unprocessed IMF component in Fig. 4 not only loses some details of real blasting seismic waves, but also reduces the accuracy of time-frequency extraction. The instantaneous frequency obtained by Hilbert transform of IMF component obtained by EMD in Fig. 4 may not have practical physical significance. Therefore, it is necessary to suppress the endpoint effect and modal confusion of EMD. Through the endpoint effect and modal confusion suppression, the decomposition accuracy of IMF will be improved, and the IMF with clearer physical meaning will be obtained.
By observing Table 1, it can be found that the IMF component with the largest blasting energy is IMF2. The envelope of IMF2 is extracted and the envelope line as shown in Fig. 5.
Observing Fig. 5, we can find 5 obvious peaks at 0.0471 s, 0.0734 s, 0.1463 s, 0.2037s and 0.2886 s respectively, which indicates that the blasting is composed of 5 blasting seismic waves. Observing the holes layout of the blasting network of the diversion tunnel in Fig. 1, we can find that the blasting is also divided into 5 sections, which reflects the accuracy of the actual millisecond time calculated by EP-CEEMDAN-INHT algorithm.
Taking the first peak time as the initiation time of the first detonator, the actual millisecond time of each detonator initiation can be calculated as 26.30 ms, 72.90 ms, 57.40 ms and 84.90 ms respectively, which is compared with the theoretical delay time interval of  Table 2.
By observing Table 2, it can be found that the actual millisecond time calculated based on EP-CEEMDAN-INHT algorithm is within the theoretical millisecond time interval in the detonator specification table provided by the manufacturer, which indicates that the performance of this batch of detonators used in this blasting is reliable and can ensure the smooth progress of millisecond blasting of diversion tunnel. Further observation of Fig. 5 shows that there is little  difference in the amplitude of seismic wave in 5 sections of this blasting. According to MATLAB programming, the measured blasting vibration signal in Fig. 5 is divided into five sub signals. Assuming that the waveform amplitude and frequency of each sub signal are roughly the same, the signal shown in Fig. 6 can be used to replace each sub signal (Ling 2004). The reasonable time interval of millisecond is determined by the method of disturbance reduction, which can greatly reduce the damage of blasting    Fig. 7, it can be found that different millisecond time has great influence on blasting vibration intensity. When the millisecond time is less than 2.86 ms, the five sub signals are sent at one time, and the blasting vibration effect reaches the maximum. At this time, the amplitude is the result of the superposition of the five sub signals. When the millisecond time is 54.51-59.75 ms, the amplitude of millisecond blasting is the smallest and the damping effect is the best. When the millisecond time interval is between 2.86 and 120.83 ms, the peak value of vibration velocity increases or weakens in varying degrees after the superposition of five sub signals, which is the result of interference between the sub signals. When the micro difference time is more than 120.83 ms, it can be found that there is not much difference between the peak value of the superimposed signal and the peak value of the sub signal. It shows that the superposition signal is equivalent to the result of each component signal acting alone, which explains why the single hole continuous initiation of the electronic detonator has a good damping effect (Fu et al. 2012).
To sum up, the most reasonable millisecond time of non-electric millisecond detonator initiation in this project is 54.51-59.75 ms and the millisecond blasting with this millisecond time has the smallest peak vibration velocity and the best damping effect. The actual millisecond time obtained by EP-CEEMDAN-INHT algorithm corresponds to the designed blasting section of the initiation network one by one, which reflects that the actual millisecond time obtained by EP-CEEMDAN-INHT algorithm is scientific and authentic. Comparing the calculated actual millisecond time with the theoretical delay time of detonator provided by the manufacturer, we can judge whether the detonator is in normal service during construction, which has important practical significance for blasting safety control.

Conclusions
(1) The decomposition results of EP-CEEMDAN algorithm can clearly display the signal frequency information contained in the measured blasting vibration data, clearly distinguish the high frequency, medium frequency and low frequency, and suppress the modal confusion caused by noise signal and the endpoint effect of the algorithm itself.
(2) The instantaneous frequency obtained by INHT is obtained by the frequency modulation component of IMF, which is no longer limited by Bedrosian theorem, and the obtained instantaneous frequency is more authentic.
(3) The actual millisecond time of non-electric millisecond detonator based on EP-CEEM-DAN-INHT algorithm is within the theoretical millisecond time in the detonator specification table provided by the manufacturer, which shows that the detonator used in this blasting is reliable and in normal service condition in actual blasting construction, and can ensure the smooth progress of millisecond blasting of diversion tunnel. (4) The most reasonable millisecond time of nonelectric millisecond detonator initiation in the diversion tunnel is 54.51-59.75 ms, and the blasting construction with this millisecond time has the smallest signal peak vibration velocity and the best damping effect.