First arrival detection algorithm based on high-resolution time-frequency analysis of seismic signal energy statistics

Abstract

The application discloses a kind of based on high-resolution time-frequency analysis seismic signal energy statistics first break detection algorithm, according to the frequency domain distribution difference of seismic event signal and noise, using narrow-band frequency division technique will seismic time domain waveform be transformed to time-frequency domain, the ratio of local energy mean and global energy mean is calculated in each frequency band and is binarized, when the result is greater than the threshold value set, consider this frequency band the time belongs to the part of seismic event signal, if under a certain threshold value, the binarization result of some time is same, only need to continue to increase threshold value and utilize the method to carry out detection again, until finally only one time point detects the number of frequency band most, the maximum value jump point is considered as initial time.The experimental results show that the detection method disclosed by the application can improve the detection accuracy while solving the high false alarm problem, and has good applicability and robustness under low signal-to-noise ratio signal conditions.

Description

Technical Field

[0001] This invention relates to the field of seismic signal detection technology, specifically to a seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis. Background Technology

[0002] Seismic signal detection plays a crucial role in geological exploration. The principle of seismic wave detection is to identify seismic waves from background noise by utilizing the differences in amplitude, frequency, polarization, and statistical characteristics between seismic signals and background noise. Commonly used seismic signal detection methods can be categorized into three types: those performed in the time domain, those performed in the frequency domain, and those performed in the time-frequency domain. With the rapid development of deep learning, more and more researchers are combining seismic detection with deep learning, achieving good results.

[0003] Among these, the methods performed in the time domain mainly include the Short-Time-Long-Time Energy Ratio (STA / LTA), the Autoregressive Akaike Information Criterion, the Fractal Dimension Method, and the Correlation Method; the methods processed in the frequency domain mainly include the Fourier Transform Method and the Instantaneous Frequency Method; the time-frequency analysis methods in the frequency domain mainly include the Short-Time Fourier Transform Method, the Wavelet Transform Method, and the S-Transform Method; and the deep learning methods mainly utilize convolutional neural networks and U-shaped neural networks to detect the first arrival of seismic waves and identify seismic phases.

[0004] However, the above methods each have their own advantages and limitations: The STA / LTA method uses the dynamic change of the ratio of the average energy of the short-term window to the average energy of the long-term window over time to reflect the change in signal amplitude (or energy) characteristics, and uses this to identify the first arrival time of seismic waves. It is computationally simple, but prone to missed detections for signals with low signal-to-noise ratios (SNR), and has poor robustness to strong interference noise. The AIC (Akaike Information Criterion) criterion first determines the order of the AR model, and then determines the phase based on the change in the AR representation. Similar to STA / LTA, the results of the AIC method also depend on the SNR and the detection interval. STFT is developed based on Fourier transform, and can consider both signal frequency and time. However, the implementation of STFT requires the introduction of a window function, and once the window function is selected, the time and frequency resolution of the entire time-frequency analysis is fixed. Under low SNR conditions, it is difficult to distinguish between signal and noise frequencies. Therefore, the choice of window function has a significant impact on the analysis effect of STFT. In summary, the most prominent problem with the above methods is that they result in a high false alarm rate when processing signals with low signal-to-noise ratios, leading to low reliability of detection results and greatly affecting the accuracy and real-time detection of the seismic phase. Summary of the Invention

[0005] To address the aforementioned problems, this invention provides a seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis.

[0006] The core technical idea of ​​this invention is as follows: This invention uses a time-frequency domain signal energy statistical method to detect seismic signals. First, a high-resolution time-frequency analysis method is used to determine the dominant frequency band of the seismic event signal. Then, fine narrowband frequency division is performed within the dominant frequency band, and wideband frequency division is performed outside the dominant frequency band. Next, the ratio of the average signal energy to the average noise energy within a specific time window is calculated in each frequency band. The times when the ratio is greater than a given threshold are counted, and the results are binarized. Finally, the binarized results in each frequency band are added together to obtain the final result as the detection basis. Since the noise interference signal is distributed outside the dominant frequency band or only partially distributed within the dominant frequency band, while the event signal is widely distributed within the dominant frequency band, this result will enhance the average energy ratio of the seismic signal and suppress the average energy ratio of the noise. Therefore, the binarized cumulative result at the first arrival point of the seismic signal will be much greater than at other times, so this can be used as the basis for judging the seismic signal. If the binarized results are the same at certain times under a certain threshold, the threshold can be increased again and the method can be used for detection until only one time point detects the maximum number of frequency bands.

[0007] The technical solution adopted in this invention is as follows:

[0008] The seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis is characterized by the following steps:

[0009] Step 1: Obtain the raw seismic signal;

[0010] Step 2: Perform continuous wavelet transform on the original seismic signal to determine the dominant frequency band of the seismic signal;

[0011] Step 3: Perform narrowband frequency division on the signal within the dominant frequency band and wideband frequency division on the signal in the non-dominant frequency band to obtain the corresponding narrowband and wideband respectively;

[0012] Step 4: Calculate the mean energy ratio in both narrow and wide frequency bands until the entire seismic data is calculated to obtain the time-frequency band map of the seismic signal;

[0013] Step 5: Obtain the STA / LTA curve from the time-frequency band diagram, perform binarization on the curve based on the given threshold, and sum the binarization results calculated in each frequency band;

[0014] Step 6: Determine the maximum value of the summation result. If there is no maximum value, increase the original threshold by 10% to obtain a new threshold and proceed to step 5; otherwise, proceed to step 7.

[0015] Step 7: Use the obtained maximum value as the initial arrival point and perform earthquake detection through the initial arrival point.

[0016] Furthermore, the specific steps of step 2 include:

[0017] Step 21: For the obtained seismic signal f(t)∈L 2 (R), the wavelet coefficients of the seismic signal are obtained by performing continuous wavelet transform through equation (1);

[0018]

[0019] Where a, b ∈ R, a > 0, a is the scale parameter; b is the translation parameter; W f (a,b) are the wavelet coefficients of the signal f(t); It is a wavelet cluster ψ a,b The conjugate function of (t);

[0020] Step 22: Reconstruct the seismic signal f(t) using equation (2) based on the wavelet coefficients of the seismic signal:

[0021]

[0022] Step 23: Determine the dominant frequency band of the seismic signal based on the time-frequency diagram of the reconstructed seismic signal.

[0023] Furthermore, step 3 includes the following specific steps:

[0024] Step 31: For signals in the dominant frequency band range of N1 to N2 Hz, determine the appropriate number of frequency band intervals M according to the preset frequency division accuracy;

[0025] Step 32: Use a bandpass filter with a bandwidth of A set of narrowband filters decomposes the seismic signal into M consecutive narrowband sub-signals within N1 to N2, thus obtaining the corresponding narrowband.

[0026] Step 33: Perform broadband frequency division on the signals within the frequency ranges of 0 to N1 and N2 to N Hz to obtain the corresponding broadband bands.

[0027] Furthermore, step 4 includes the following specific steps:

[0028] Step 41: Traverse the seismic data and calculate the energy mean ratio in both broadband and narrowband frequencies:

[0029]

[0030] Where STA and LTA are the average energy values ​​of the short-term and long-term windows, respectively; N STA N LTA These are the lengths of the short-term window and the long-term window, respectively; CF is the characteristic function, and CF(i) = Y(i). 2 , where Y(i) is the signal amplitude;

[0031] Step 42: Use the energy mean ratio as the vertical axis and the broadband or narrowband size as the horizontal axis to obtain the time-frequency diagram of the seismic data.

[0032] Furthermore, step 5 includes the following specific steps:

[0033] Step 51: Obtain the STA / LTA curve of the seismic data based on the time-frequency plot of the seismic data;

[0034] Step 52: For a given threshold THR, binarize the STA / LTA curve. Values ​​greater than the threshold THR are identified as seismic event signals and set to 1; conversely, values ​​less than the threshold THR are identified as noise and set to 0.

[0035]

[0036] Step 53: Convert the calculated binarized result X THR (t) is used to sum the results.

[0037] The beneficial effects of this invention are:

[0038] This invention addresses the shortcomings of traditional seismic signal detection methods by proposing a statistical first-arrival detection algorithm based on high-resolution time-frequency analysis. While considering signal energy, it focuses on the dominant frequency band of the seismic event signal. After accurately determining the dominant frequency band through high-resolution time-frequency analysis, waveform bandpass frequency division is performed, effectively suppressing environmental background noise. Furthermore, considering the difference in frequency distribution between noise and event signals, especially when the signal is submerged in noise, this method can still detect seismic signals from high background noise, demonstrating its strong first-arrival detection capability and good applicability and robustness. It improves event detection accuracy while solving the high false alarm problem, and also addresses the stringent threshold selection challenge of the STA / LTA method, enhancing the accuracy and reliability of the detection results. This is more conducive to subsequent research and analysis on event identification, source location, focal mechanism, and geological exploration. Attached Figure Description

[0039] Figure 1 This is a continuous wavelet transform diagram of the seismic signal.

[0040] Figure 2 This is a flowchart of the present invention.

[0041] Figure 3 This is the waveform of the seismic signal received by the Zhoushan station in Zhejiang.

[0042] Figure 4(a) and (b) show the first arrivals of seismic waves picked up by the STA / LTA method when Thr = 2 and Thr = 3.5, respectively.

[0043] Figure 5 This is a continuous wavelet transform diagram of the seismic signal.

[0044] Figure 6 This is the result after frequency division processing of the seismic signal.

[0045] Figure 7 The results are the binarized values ​​for each frequency band.

[0046] Figure 8 The final result is obtained by using this method.

[0047] Figure 9 When THR=2, the average energy ratio of long and short time is used for each frequency band.

[0048] Figure 10 When THR=2, the basis for determining the first arrival of an earthquake is used.

[0049] Figure 11 The waveform of the signal after adding 10dB Gaussian white noise.

[0050] Figure 12 To obtain the initial arrival detection results using this method.

[0051] Figure 13 The waveforms are those of the signal after adding 300dB, 50dB, and 60dB Gaussian white noise.

[0052] Figure 14 The results are from experiments using this algorithm after adding different noise signals.

[0053] Figure 15 The results are from experiments using this algorithm on a signal with 70dB noise added. Detailed Implementation

[0054] To enable those skilled in the art to better understand the technical solutions of the present invention, the technical solutions of the present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0055] 1. High-resolution time-frequency analysis of seismic signals

[0056] Continuous wavelet transform overcomes the limitation of short-time Fourier transform (SFT) where the window size does not change with frequency. It has excellent characterization capabilities for both slowly changing low-frequency and rapidly changing high-frequency signals, meeting the resolution requirements of different signals with varying characteristics. Its core idea is to improve upon the fixed resolution of SFT by employing adaptive, scalable, and shiftable wavelet basis functions, widening the window shape at low frequencies and narrowing it at high frequencies. This invention utilizes continuous wavelet transform for time-frequency analysis of seismic signals to determine the dominant frequency bands, laying the foundation for subsequent signal frequency division processing.

[0057] For a given signal f(t)∈L 2 (R), whose continuous wavelet transform is defined as:

[0058]

[0059] Where a, b ∈ R, a > 0, a is the scale parameter; b is the translation parameter; W f (a,b) are the wavelet coefficients of the signal f(t); It is a wavelet cluster ψ a,b The conjugate function of (t). The inner product operation reflects the local features of the signal and the similarity of the wavelet function; the larger the inner product, the higher the similarity.

[0060] Using the wavelet coefficients W shown in equation (1) f Given (a, b), the original signal f(t) can be reconstructed. The reconstructed signal is:

[0061]

[0062] Among them, C Ψ Assume a wavelet of Fourier transform form satisfies the admissibility condition and For Ψ a,b Fourier transform of (t).

[0063] According to equation (1), we can obtain the following: Figure 1 The continuous wavelet transform time-frequency plot of a certain seismic signal shown is from... Figure 1 The P-wave and S-wave phases of the seismic signal are clearly visible, and their energy is mainly concentrated in the 2-10Hz frequency band. The energy of the signal below 2Hz and above 10Hz is very weak. Therefore, this frequency band is considered to be the dominant frequency band of the seismic signal, and the frequency band outside this band is considered to be noise and other interference signals.

[0064] Frequency division processing is also an important time-frequency analysis method with wide applications in seismic signal denoising. Since seismic signals are non-stationary, frequency division analysis is a powerful tool for studying the structural characteristics of non-stationary signals. Based on the laws of time-frequency analysis, seismic data is rationally divided into frequency bands to understand the distribution patterns of signals and noise within different frequency bands. Targeted processing is then applied to the data in each frequency band. This frequency division noise suppression method can improve signal fidelity. The amplitude and phase of seismic waves are functions of time and frequency. If the boundaries of each narrow frequency band are properly processed, broadband seismic data can be reconstructed, thus achieving frequency division processing. Therefore, frequency division processing can effectively separate seismic signals from background noise. Hence, this invention uses bandpass filtering to perform frequency division processing on the signal.

[0065] 2. Seismic signal energy statistical detection algorithm based on high-resolution time-frequency analysis

[0066] This algorithm is as follows: Figure 2 As shown, this method is based on the feasibility analysis of determining the dominant frequency bands using high-resolution time-frequency analysis of seismic signals and performing frequency division processing. Waveform data for the corresponding frequency bands is extracted using bandpass wavelets. Then, for the waveform data within each frequency band, the ratio of the average signal energy within a specific time window to the average noise energy within the corresponding frequency band is calculated. The times when this ratio exceeds a given threshold are counted, and the results are binarized. Finally, the binarized results for each frequency band are summed to obtain the final result, which serves as the basis for determining the seismic signal.

[0067] For seismic monitoring stations, most of what they record is background noise; the seismic event signal accounts for a very small proportion of the entire time series. Assume that the continuously input single-channel waveform data to be detected x(t) consists of a finite-length infrasound signal s(t) and random noise n(t), that is:

[0068] x(t)=s(t)+n(t)(3)

[0069] Time-frequency analysis determines the frequency range of the data within the channel to be 0–N Hz. This frequency range covers both the background noise and seismic event signal frequencies. The dominant frequency band is identified as N1–N2 Hz. For signals within the 0–N1 and N2–N Hz ranges, broadband coarse frequency division is performed to generate corresponding widebands. Assuming the number of frequency intervals for each division is n1 and n2, for signals within the N1–N2 Hz range, a suitable number of frequency intervals M is determined based on the preset division accuracy. Bandpass filtering is then employed, with a bandwidth of [missing information]. A set of narrowband filters decomposes the seismic event signal into M consecutive narrowband sub-signals within the dominant frequency band;

[0070] Based on numerous experiments, it has been found that when the number of frequency bands in the dominant frequency band is nearly twice or more than the number of other frequency bands (i.e., M≥2*(n1+n2)), the detection accuracy and precision will both achieve good results. At this time, the infrasound signal s(t) is:

[0071] s(t)=s1(t)+s2(t)+s3(t)+…+s M (t)(4)

[0072] As can be seen from equation (4), the original signal is obtained by adding the frequency bands of the frequency-divided signal, which also proves that it is feasible to perform frequency division processing on the seismic signal.

[0073] Next, the energy mean ratio R(t) is calculated for the corresponding wideband and each narrowband sub-signal. The calculation result is used as the ordinate, and the bandwidth as the abscissa to obtain the time-frequency plot (or simply time-frequency plot) for the wideband and narrowband. The formula for calculating R(t) is:

[0074]

[0075] Where STA and LTA are the average energy values ​​of the short-term window and the long-term window, respectively, and N STA N LTA These are the lengths of the short-time window and the long-time window, respectively. CF is a characteristic function that reflects changes in signal amplitude and frequency. In this invention, CF(i) = Y(i) is used. 2 , where Y(i) is the signal amplitude.

[0076] When a seismic event signal has not yet arrived, the ratio approaches 1. When a seismic event arrives, the average energy within the short time window increases, and the ratio also increases accordingly, thus determining the initial arrival of the earthquake. When the long and short time windows traverse the entire seismic data, the STA / LTA curve is obtained. Given an appropriate threshold THR, the curve is binarized: values ​​greater than the threshold THR are considered seismic events and set to 1; values ​​less than the threshold THR are considered noise and set to zero.

[0077]

[0078] Finally, the binarized results X of each frequency band are... THR Summing the results gives the final outcome.

[0079] X THR This reflects information about the frequency and bandwidth of seismic event signals that may be contained in the seismic data. When a seismic event first occurs, X... THRThe value will be significantly greater than the value during the noise period, and this will be used to ultimately determine the first arrival of the earthquake and conduct earthquake detection. Since the earthquake signal is interpreted at the arrival time, only the starting position of the maximum value point is considered, and the maximum value point can be uniquely determined. The maximum value is the first arrival point. If there is no maximum value, the threshold THR is increased. Each time the threshold is increased by 10% based on the previous threshold, so no false alarms will occur, and the interpretation results will be reliable.

[0080] Example

[0081] To further verify the effectiveness of the method proposed in this invention, relevant simulation experiments were conducted and the results were analyzed.

[0082] The STA / LTA method uses the ratio of STA (Short-Time Average) to LTA (Long-Time Average) to reflect changes in signal amplitude, frequency, and other characteristics. When a seismic signal arrives, the STA / LTA value will experience a sudden change. When the ratio exceeds a certain set threshold R, a seismic event is determined to have occurred. The STA / LTA method is fast, simple to apply, and has good real-time performance, making it suitable for detecting large amounts of data. Therefore, it is the main method for real-time seismic wave detection.

[0083] First, an experiment was conducted using a high signal-to-noise ratio (SNR) real seismic signal received from the Zhoushan seismic station in Zhejiang Province. This method primarily utilizes the frequency and relative amplitude information of the seismic signal, thus eliminating the need for data removal to counteract instrument response, reducing data preprocessing steps, and making it more suitable for real-time detection. The sampling frequency was 100Hz, and the manually read signal was 168.30s. Its waveform is shown below. Figure 3 As shown.

[0084] The STA / LTA method was used to pick the first arrival of earthquakes. Following the principle proposed by Liu et al. (2014), the short-term window length is preferably set to 2-3 times the main signal period, and the long-term window is generally 5-10 times the short-term window. The short-term window was set to 3.5 s, and the long-term window to 15 s. Thresholds were set to THR = 2 and THR = 3.5, respectively. The picking results are as follows: Figure 4 As shown.

[0085] from Figure 4 (b) It can be seen that the STA / LTA method can detect the first arrival only when the signal-to-noise ratio is high and the threshold value is appropriate. At this time, the arrival time of the first arrival is 168.58s. However, this method is more sensitive to noise. If the threshold is too small, such as when THR=2 in 4(a), three false alarms occur. The first arrival has already arrived before the maximum value of the STA / LTA curve. The picking result is slightly delayed compared to the actual arrival time. Therefore, the picking accuracy is not high enough.

[0086] Next, this method will be used for first arrival detection. First, continuous wavelet transform will be used to examine the frequency range of the original signal and the dominant frequency band of the seismic event signal. The results are as follows: Figure 5 As shown.

[0087] from Figure 5 The dominant frequency band of the seismic event signal in the identified signal record is 3–8 Hz, i.e., N1 = 3, N2 = 8. Taking n1 = 6, n2 = 2, and M = 12, the signal is divided into 20 narrow frequency bands. The required frequency division accuracy is approximated here. The frequency division result is as follows: Figure 6 As shown.

[0088] Then, the long-short time window energy ratio method is applied in each frequency band. At this point, the threshold selection can be more flexible because after frequency division, noise and event signals are essentially separated. Even if the threshold is small, after binarization and summing the results, the goal of amplifying event information while suppressing noise can still be achieved. Here, the threshold THR is set to 3.5. The results of using the long-short time window energy ratio method and binarization in each frequency band are as follows... Figure 7 As shown.

[0089] Finally, the binarized results of each frequency band are summed to obtain the final criteria for first arrival detection. The point where the value jumps to the maximum value is determined as the first arrival of the earthquake. The results are as follows: Figure 8 As shown.

[0090] Depend on Figure 8 It can be seen that the accumulated value at the initial arrival point is significantly greater than the accumulated values ​​of the other segments, and the reading result is 168.27s. Next, keeping the other parameters unchanged, when the threshold is set to THR=2, the result is compared with the result obtained using the STA / LTA method when THR=2, as shown below. Figure 9 As shown.

[0091] The final result after summing is as follows Figure 10 As shown: by Figure 10 It can be seen that when the threshold is reduced, the cumulative value of the sequence segments other than the first arrival is still significantly different from that at the first arrival. The reading result is 168.43s, which is more favorable for earthquake detection than the traditional STA / LTA method. This proves that the algorithm can overcome the shortcomings of the original STA / LTA method in that the threshold selection is relatively strict. At the same time, it can be applied to the condition of first arrival detection for low signal-to-noise ratio signals, and the algorithm has strong applicability.

[0092] To further verify the performance of this method in picking first arrivals of actual earthquakes, the traditional STA / LTA method and this method were applied to 60 natural earthquake records. 1024 points containing the events were extracted from these records, with the starting point as time zero and a sampling frequency of 40Hz. The picking results of the two algorithms were compared with those obtained manually. The parameter settings for the traditional STA / LTA method continued to follow the principles proposed by Liu et al. (2014), with a short time window of 1.5s and a long time window of 7.5s. The picking results of the two methods are shown in Table 1.

[0093] Table 1 Comparison of initial arrival picking results from two methods with manual results.

[0094]

[0095]

[0096] Table 1 shows that, for the same event, using two different methods for first arrival picking, under given parameters, the average deviation of the traditional STA / LTA method for picking the first arrival time of seismic waves is 0.1079 s, while the standard deviation is 0.1012 s. The average deviation and standard deviation of our proposed method are 0.0383 s and 0.0565 s, respectively. Both indicators are superior to the traditional STA / LTA method, indicating that the detection accuracy of our proposed method is significantly better than that of the STA / LTA method.

[0097] For complex monitoring environments and massive amounts of monitoring data, the ability to accurately detect event signals with significant background noise interference or those submerged in noise is one of the important criteria for evaluating the applicability of an algorithm. The following section examines the applicability and robustness of the algorithm by adding noise interference of varying magnitudes to the original signal.

[0098] First, 10dB of Gaussian white noise is added to the original signal. The waveform of the signal after adding noise is as follows. Figure 11 As shown.

[0099] The results of first-arrival detection using the algorithm presented in this paper are as follows: Figure 12 As shown, from Figure 12 It can be seen that adding 10dB of Gaussian white noise has almost no impact on the detection results.

[0100] Comparative experiments were conducted by adding Gaussian white noise of 30dB, 50dB, and 60dB respectively. Figure 13 The algorithm was used to conduct experiments on three signals from top to bottom, and the results are as follows. Figure 14 As shown. From Figure 14 It can be seen that the initial arrival can still be accurately identified even after adding 60dB of noise, while the signal is as follows when 70dB of Gaussian white noise is added: Figure 15As shown, it can be seen that the signal is almost submerged in noise. When using this method to detect the first arrival, although the value of the binarized accumulation result at the first arrival is not much different from that at other times, it can still clearly determine the first arrival time, thus enabling accurate detection of earthquake events. This also proves the applicability of this algorithm and its robustness under strong interference noise.

[0101] As shown in Table 2, two methods were used to detect signals after adding different amounts of Gaussian white noise. After adding 10dB of Gaussian white noise, both methods could accurately pick up the first arrival. However, when the noise increased to 30dB, the detection accuracy of the STA / LTA method was only 80%, and the false alarm rate also increased significantly. In contrast, the present method still maintained 100% accuracy. When the noise is even greater, the advantage of the present algorithm will be more obvious.

[0102] Table 2 Comparison of detection results between the two methods

[0103]

[0104] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis, characterized in that, Includes the following steps: Step 1: Obtain the raw seismic signal; Step 2: Perform continuous wavelet transform on the original seismic signal to determine the dominant frequency band of the seismic signal; Step 3: Perform narrowband frequency division on the signal within the dominant frequency band and wideband frequency division on the signal in the non-dominant frequency band to obtain the corresponding narrowband and wideband respectively; Step 4: Calculate the mean energy ratio in both narrow and wide frequency bands until the entire seismic data is calculated to obtain the time-frequency band map of the seismic signal; Step 5: Obtain the STA / LTA curve from the time-frequency band diagram, perform binarization on the curve based on the given threshold, and sum the binarization results calculated in each frequency band; Step 6: Determine the maximum value of the summation result. If there is no maximum value, increase the original threshold by 10% to obtain a new threshold and proceed to step 5; otherwise, proceed to step 7. Step 7: Use the obtained maximum value as the initial arrival point and perform earthquake detection through the initial arrival point.

2. The seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis as described in claim 1, characterized in that, Step 2 includes the following specific steps: Step 21: For the obtained seismic signal f(t)∈L 2 (R), the wavelet coefficients of the seismic signal are obtained by performing continuous wavelet transform through equation (1); Where a, b ∈ R, a > 0, a is the scale parameter; b is the translation parameter; W f (a,b) are the wavelet coefficients of the signal f(t); It is a wavelet cluster ψ a,b The conjugate function of (t); Step 22: Reconstruct the seismic signal f(t) using equation (2) based on the wavelet coefficients of the seismic signal: Step 23: Determine the dominant frequency band of the seismic signal based on the time-frequency diagram of the reconstructed seismic signal.

3. The seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis as described in claim 2, characterized in that, Step 3 includes the following specific steps: Step 31: For signals in the dominant frequency band range of N1 to N2 Hz, determine the appropriate number of frequency band intervals M according to the preset frequency division accuracy; Step 32: Use a bandpass filter with a bandwidth of A set of narrowband filters decomposes the seismic signal into M consecutive narrowband sub-signals within N1 to N2, thus obtaining the corresponding narrowband. Step 33: Perform broadband frequency division on the signals within the frequency ranges of 0 to N1 and N2 to N Hz to obtain the corresponding broadband bands.

4. The seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis as described in claim 2, characterized in that, Step 4 includes the following specific steps: Step 41: Traverse the seismic data and calculate the energy mean ratio in both broadband and narrowband frequencies: Where STA and LTA are the average energy values ​​of the short-term and long-term windows, respectively; N STA N LTA These are the lengths of the short-term window and the long-term window, respectively; CF is the characteristic function, and CF(i) = Y(i). 2 , where Y(i) is the signal amplitude; Step 42: Use the energy mean ratio as the vertical axis and the broadband or narrowband size as the horizontal axis to obtain the time-frequency diagram of the seismic data.

5. The seismic signal energy statistical first arrival detection algorithm based on high-resolution time-frequency analysis as described in claim 4, characterized in that, Step 5 includes the following specific steps: Step 51: Obtain the STA / LTA curve of the seismic data based on the time-frequency plot of the seismic data; Step 52: For a given threshold THR, binarize the STA / LTA curve. Values ​​greater than the threshold THR are identified as seismic event signals and set to 1; conversely, values ​​less than the threshold THR are identified as noise and set to 0. Step 53: Convert the calculated binarized result X THR (t) is used to sum the results.

Images