A method for automatic identification of seismic event first arrival based on multi-module filter integration

Abstract

This invention discloses an automatic identification method for the initial arrival of vibration events based on multi-module filter integration. The main steps are as follows: First, the continuously recorded raw signal is preprocessed and then segmented. For each segmented signal, the effective frequency distribution of vibration energy is determined. Then, operations are sequentially performed, including embedding an anti-aliasing low-pass filter, downsampling the signal, embedding an adaptive detrending term filter, embedding a median filter, and embedding a borderless delay filter. Next, the signal envelope energy is calculated based on the borderless delay filter. Finally, a discriminant filter is embedded, and the signal is normalized to obtain the final identification feature function. This invention establishes an automatic picking algorithm integrated with multiple filters. While reducing noise, it constructs a feature function, requiring only the most basic signal parameters and effectively reducing noise interference. This improves the accuracy of automatically picking the occurrence of disaster events.

Description

Technical Field

[0001] This invention relates to the field of disaster early warning and prediction technology, and in particular to an automatic identification method for the initial arrival of vibration events based on multi-module filter integration. Background Technology

[0002] Currently, disaster early warning and prediction increasingly rely on the observation of dynamic characteristics. Natural disasters such as earthquakes, volcanoes, landslides, rockfalls, and debris flows, as well as man-made disasters such as mine rockfalls, pipeline ruptures, and high-pressure gas leaks, are all accompanied by the vibrational propagation of energy. Therefore, by recording the temporal characteristics of these vibrations with instruments, and based on mathematical physics methods, the eruptive characteristics of the disaster source can be analyzed. The most fundamental analysis is determining the moment of the event, i.e., the initial arrival. Theoretically, an array of observation instruments (broad / narrow frequency seismographs, long-period seismographs, accelerometers, acoustic sensors, etc.) can obtain characteristic information such as the initial arrival time of the event, the energy scale of the event, and the azimuth of the energy source over time. However, in actual monitoring, the number of detection instruments is often limited. Therefore, how to obtain the initial arrival time of an event with a small number of monitoring instruments has always been a research hotspot in various fields, and also presents many technical challenges. Finding the source energy carried by a disaster event from complex and random vibration records is of great significance for disaster early warning and the formation of multi-scale disasters.

[0003] The main purpose of methods using dynamic vibration signals to reveal disaster events at different scales is to accurately monitor the spatiotemporal distribution characteristics of disaster events, based on detecting the initial arrival of the disaster event, i.e., the moment the disaster event occurs. Currently, the acquisition of the initial arrival moment is mainly based on the analysis of dynamic time-series signals. Typically, in the early stages of a disaster event, the time-series information recorded by the sensors is essentially background noise. At the moment of the disaster event, there is usually a relatively strong pulse signal, which then becomes the source propagating in the medium until it is received by sensors deployed at different distances. As time continues, the signal strength gradually attenuates until the vibration amplitude decreases to the scale of background noise. Initial arrival acquisition is generally based on the envelope strength, absolute amplitude, or power spectrum characteristics of the signal in the time or frequency domain. Currently popular algorithms are the Short-Time Average / Long-Time Average Ratio (STA / LTA) algorithm and the Power Spectrum Determination (PSD) method.

[0004] While the short-time average / long-time average ratio algorithm can pick up the initial arrival of an event, it relies too heavily on pre-set parameters, although some studies have modified the algorithm's sensitivity by changing the feature function. Furthermore, the length of the short-time and long-time windows affects the threshold; a low threshold may lead to many false arrivals (false positives), while a high threshold may cause weak events to be missed (false negatives). To improve upon the short-time average / long-time average ratio algorithm's over-reliance on set parameters, the power spectrum thresholding method has emerged. Background noise typically has a relatively uniform spectral content across frequency bands, while the spectral content of disaster events is very high in certain specific frequency bands. Therefore, the core of the power spectrum thresholding algorithm is to find moments with abnormally high frequency content and pick up the initial arrival of a disaster event.

[0005] While power spectrum thresholding algorithms require fewer initial parameter settings, they are highly sensitive to signal-to-noise ratio (SNR). For signals with high SNR, such as strong earthquake energy, both methods demonstrate good recognition capabilities. However, for vibrations with low SNR, or where the signal energy is obscured by background noise, such as micro-seismic events, volcanic tremors, and pipe bursts, these algorithms suffer from significant misidentification. Therefore, developing a more stable algorithm that incorporates both noise reduction and feature extraction modules is of significant practical importance for initial vibration signal acquisition and, consequently, for early warning of disasters. Summary of the Invention

[0006] To address the problems of numerous dependent parameters and low accuracy in recognizing low signal-to-noise ratio signals in the two existing vibration event first-arrival picking algorithms, this invention provides an automatic vibration event first-arrival recognition method based on multi-module filter integration.

[0007] The core of the method of this invention consists of two parts. The first part is mainly to reduce the noise of the signal, that is, to suppress noise and peak distortion. The second part is mainly to construct a feature function based on the envelope function and the discriminant function to complete the initial arrival picking.

[0008] The automatic identification method for the first arrival of vibration events based on multi-module filter integration provided by this invention has the following specific steps:

[0009] S1. The continuously recorded raw signal is preprocessed and then segmented. The effective frequency distribution of vibration energy is determined segment by segment for each segment of the segmented signal. Specifically:

[0010] S11. Preprocess the continuously recorded raw signals, including baseline leveling, mean removal, and trend removal steps;

[0011] S12. Cut the preprocessed original signal by overlapping and cutting it.

[0012] S13. For each segment of the cutting signal, determine the effective impact time of vibration energy and set the duration of the impact:

[0013] d min : The shortest duration of the disaster; d max : The longest duration of the disaster; N d The number of time windows that can be decomposed during the duration;

[0014] S14. For each segment of the cutting signal, determine the effective frequency distribution of vibration energy, select the low-frequency band, and set the frequency band range:

[0015] f low : Low-frequency cutoff frequency; f high : High-frequency cutoff frequency; N f : The number of frequency-decomposable bandwidth windows.

[0016] S2. An embedded anti-aliasing low-pass filter is shown as follows:

[0017]

[0018] In the formula: n is the filter order, ω c ω is the cutoff angular frequency; p This refers to the passband edge frequency.

[0019] S3. Signal Downsampling: Based on the low-pass filtering in the previous step, downsampling is performed, introducing a downsampling scaling factor R, expressed as:

[0020]

[0021] fsp = 2f high

[0022] In the formula, fs is the original sampling frequency of the signal; fsp is the sampling frequency after downsampling; f high This is a pre-set high-frequency cutoff frequency.

[0023] S4, Embedded adaptive detrending term filter; this filter is expressed as:

[0024] B = FIR(M,ω) c )

[0025] M = round(d max ×fsp)

[0026]

[0027] trend = filtfit(B, 1, x r )

[0028] x rd =xr -trned

[0029] In the formula, B is the filter, M is the cutoff angular frequency order, and ω c It is the cutoff angular frequency, d max The maximum duration of the disaster as preset; f low The preset low-frequency cutoff frequency; trend means de-trending, indicating that this step is de-trending filtering; x r The signal is downsampled; x rd It is the signal after trend term filtering.

[0030] S5. Embedded median filtering, median filtering is expressed as:

[0031] y = medfilt(x rd ,round(fsp×d min / 10))

[0032] In the formula, d min The minimum duration of a disaster as preset.

[0033] S6, Embedded edge-free delay filter; this filter is expressed as:

[0034] y′=filter(H,1,y)

[0035] In the formula, y' is the signal after edge-delay filtering; H is the filter parameter; and y is the signal after median filtering.

[0036] The filter parameter H is set as follows:

[0037] H = acos(2π × fc × th)

[0038]

[0039] a = 1: length(th): 0

[0040]

[0041] k = 1, 2, ..., N f

[0042] In the formula, f c For [f low f high The center frequency of [], th is the time series, a is the amplitude transformation coefficient representing the range from 0 to 1, the number of th, k is the number of frequency points, and N f The number of frequency band windows that can be decomposed from a predefined frequency.

[0043] S7. Calculate the signal envelope energy: Based on edge-free delay filtering, the signal is divided into sub-signals of different frequency bands according to different center frequencies, and then the envelope energy E of each signal segment is calculated. k The calculation formula is as follows:

[0044]

[0045] k = 1, 2, ..., N f

[0046] In the formula, H is the Hilbert transform of the signal; y' is the signal after edge-delay filtering.

[0047] S8. First, embed the discriminant filter, then establish the feature function;

[0048] The expression for the discrimination filter is as follows:

[0049] D(m) = filter(T,1,E) k )

[0050] m = 1, 2, ..., N d

[0051] k = 1, 2, ..., N f

[0052] In the formula, D is the signal after discrimination filtering; T is an intermediate variable containing the penalty factor β for non-pulse signals; E k It is the signal envelope energy, where m represents the sequence of natural numbers, indicating that the values ​​are taken from the first to the Nth. d 1, N d It is the number of time windows that can be decomposed over time;

[0053] The intermediate variable T is represented as:

[0054] T=β×((1:l)-(l+1))

[0055]

[0056] L = round(D(m) × fsp)

[0057] D = [d min N d ,d max ]

[0058]

[0059] In the formula, β is the penalty factor for non-pulse signals, and L represents the data length of D.

[0060] The signal is normalized to obtain the final recognition feature function, which is expressed as follows:

[0061] CF = 10 × (log 10 (1+D / V) / log 10 (2))

[0062] In the formula, CF is the initial arrival recognition feature function; V is the maximum amplitude of the signal; and D is the signal after discrimination filtering.

[0063] Compared with the prior art, the advantages of the present invention are:

[0064] (1) The method of the present invention can pick up stable first arrival times within a wide range of parameters such as vibration time window, signal frequency, and penalty factor. Different settings only affect the picking resolution of the event, but the picking results for important events are basically consistent. It effectively avoids the dependence of the feature function on the threshold parameter.

[0065] (2) The present invention incorporates multiple filters, which avoids the influence of background noise, high-frequency instrument noise and low-frequency environmental noise, effectively improving the quality of the recorded signal and improving the accuracy of the calculated feature parameters.

[0066] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Attached Figure Description

[0067] Figure 1 The flowchart of the automatic vibration event initial arrival design calculation method based on the multi-module filter integration method of the present invention.

[0068] Figure 2 An example diagram of automatic vibration event picking based on a multi-module filter ensemble algorithm is shown. In the diagram, a is the original vibration signal time series; b is the signal time series processed by anti-aliasing low-pass filtering; c is the signal time series processed by downsampling; d is the signal time series processed by adaptive detrending term; e is the signal time series processed by median filtering; f is the signal time series processed by edge-delay-free filtering; g is the energy envelope curve of the multi-band signal; h is the superposition of energy envelope curves; and i is the normalized characteristic function.

[0069] Figure 3 Characteristic functions and first arrival time extraction for different filter frequency bands (d min =1s,d max =2s, N d =1, N f =1, β=1).

[0070] Figure 4 Characteristic functions and first arrival times for different vibration durations (f) low =1Hz, fhigh =2Hz, N d =1, N f =1, β=1).

[0071] Figure 5 Characteristic functions and first arrival time extraction for different number of frequency bands (Nf) (f low =1Hz, f high =3Hz, d min =1s,d max =3s, N d =1, β=1).

[0072] Figure 6 Characteristic functions and first arrival times for different penalty factors β (f low =1Hz, f high =3Hz, d min =1s,d max =3s, N d =1, N f =1).

[0073] Figure 7 Characteristic functions and first arrival times of signals (SNR=1,2,3,5) under different noise intensities.

[0074] Figure 8 Vibration timing record and initial arrival time of an underground pipeline rupture event.

[0075] Figure 9 Sequence of vibrations and first arrival time of earthquake events. Detailed Implementation

[0076] The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0077] like Figure 1-9 As shown, the automatic identification method for the initial arrival of vibration events based on multi-module filter integration provided by the present invention includes the following specific steps:

[0078] S1. Set the basic parameters of the timing signals: as follows:

[0079] S11. Preprocess the continuously recorded raw signals, including steps such as baseline leveling, mean removal, and trend removal. Figure 2 As shown in a.

[0080] S12. The preprocessed original signal is cut. To ensure the continuity of the cut signal and minimize energy leakage, an overlapping method is used for cutting, with an overlap of more than 50%.

[0081] S13. For segment-by-segment vibration energy impact time determination of the cutting signal, since the energy triggered by a disaster is usually pulse energy, the duration of the impact time is set here:

[0082] d min : The shortest duration of the disaster; d max : The longest duration of the disaster; N d The number of time windows that can be decomposed during the duration;

[0083] S14. For each segment of the cutting signal, the effective frequency distribution of vibration energy is determined. Typically, the initial arrival of disaster-induced vibration energy is a compressional body wave phase (P-wave). Body waves do not exhibit dispersion effects, meaning the energy propagation speed is the same across all frequency bands. Therefore, a low-frequency band can be selected to ensure a simple and clear signal. Frequency band range is set here:

[0084] f low : Low-frequency cutoff frequency; f high : High-frequency cutoff frequency; N f : The number of frequency-decomposable bandwidth windows.

[0085] The following steps will be explained using a segmented signal as an example.

[0086] S2. Embedding an anti-aliasing FIR filter: To eliminate the influence of high-frequency noise on the vibration signal generated by the disaster energy, this step designs an embedded low-frequency filter (such as...). Figure 2 As shown in b), to ensure the signal's frequency response curve is as flat as possible in the passband frequency domain, while gradually approaching zero in the stopband frequency domain, a Butterworth low-pass filter can be used to ensure significant information within the signal's observation frequency band. The filter is represented as:

[0087]

[0088] In the formula: n is the filter order, ω c ω is the cutoff angular frequency, which is the frequency at which the amplitude drops to -3dB; p This refers to the passband edge frequency. The conversion between angular frequency and line frequency is: ω = 2πf. For example, if the signal sampling frequency (fs) is 100Hz, the maximum frequency of the signal is 50Hz. To retain energy below 10Hz, i.e., filter energy above 10Hz, ω... c =2*π*10=62.8rad / s.

[0089] S3. Signal Downsampling: To simplify the computational workload of signal analysis and reduce signal complexity as much as possible, and to retain important initial arrival information of the signal, downsampling is performed based on the low-pass filtering in the previous step (e.g., ...). Figure 2 As shown in c), the downsampling scaling factor R is introduced and expressed as:

[0090]

[0091] fsp = 2f high

[0092] In the formula, fs is the original sampling frequency of the signal; fsp is the sampling frequency after downsampling; f high This is a pre-set high-frequency cutoff frequency.

[0093] S4. Embedded adaptive de-trending filter;

[0094] To eliminate the impact of low-frequency noise (weather noise or slow instrument drift) on signal monitoring, this step involves designing an embedded detrending term filter. This is because low-frequency noise primarily causes a continuous deviation of the signal from the baseline (e.g., Figure 2 (As shown in d). This filter consists of a linear phase filter and a zero-phase-shift filter function. The filter is represented as:

[0095] B = FIR(M,ω) c )

[0096] M = round(d max ×fsp)

[0097]

[0098] trend = filtfit(B, 1, x r )

[0099] x rd =x r -trned

[0100] In the formula, B is the filter, M is the cutoff angular frequency order, and ω c It is the cutoff angular frequency, d max The maximum duration of the disaster as preset; f low The preset low-frequency cutoff frequency; x r The signal is downsampled; x rd It is the signal after trend term filtering.

[0101] S5, Median filter;

[0102] To further eliminate spike noise in the signal caused by monitoring instrument malfunctions or data transmission, this step involves embedding a median filter to further smooth the signal (e.g., ...). Figure 2 As shown in e), median filtering is expressed as:

[0103] y = medfilt(x rd ,round(fsp×d min / 10))

[0104] In the formula, d min The minimum duration of a disaster as preset.

[0105] S6. Embedded No-edge delay filter;

[0106] Disasters, acting as vibration sources, release energy, often in the form of pulse signals. To enhance the prominence of these pulse signals, this step embeds a borderless delay filter to further reduce background noise interference with the signal (e.g., Figure 2 (as shown in f); the filter is represented as:

[0107] y′=filter(H,1,y)

[0108] In the formula, y' is the signal after edge-delay filtering; H is the filter parameter; and y is the signal after median filtering.

[0109] The filter parameter H is set as follows:

[0110] H = acos(2π × fc × th)

[0111]

[0112] a = 1:th:0

[0113]

[0114] k = 1, 2, ..., N f

[0115] In the formula, f c For [f low f high The center frequency of N f This is the number of frequency band windows that can be predefined for frequency decomposition. For example, if the analysis frequency band of a signal is defined as [1Hz, 3Hz], then N... f If the value is 3, then f c The result is [1.33Hz, 2Hz, 2.67Hz].

[0116] S7. Calculate the signal envelope energy:

[0117] Based on edge-free delay filtering, the signal is divided into sub-signals of different frequency bands according to different center frequencies, and then the envelope energy E of each signal segment is calculated. k (like Figure 2 (as shown in g), the calculation formula is as follows:

[0118]

[0119] k = 1, 2, ..., N f

[0120] In the formula, H is the Hilbert transform of the signal; y' is the signal after edge-delay filtering.

[0121] S8. First, embed the discriminant filter, then establish the feature function;

[0122] (1) In order to ensure the identification accuracy of the initial vibration energy caused by the disaster event, this step embeds a discrimination filter, which enhances signal attenuation by rapidly increasing and gradually decreasing speed, while further reducing the influence of impulse noise. In addition, a penalty factor β for non-impulse signals is introduced to improve the filter's ability to suppress background noise. The expression of the discrimination filter is as follows:

[0123] D(m) = filter(T,1,E) k )

[0124] m = 1, 2, ..., N d

[0125] k = 1, 2, ..., N f

[0126] In the formula, D is the signal after discrimination filtering; T is an intermediate variable containing the penalty factor β for non-pulse signals; E k It is the signal envelope energy, N d It is the number of time windows that can be decomposed over time.

[0127] The intermediate variable T is represented as:

[0128] T = β × ((1:L) - (L+1))

[0129]

[0130] L = round(D(m) × fsp)

[0131] D = [d min N d ,d max ]

[0132]

[0133] In the formula, β is the penalty factor for non-pulse signals.

[0134] (2) To reduce the loss of initial arrival amplitude intensity in the filtered signal, the signal is normalized to obtain the final recognition feature function (e.g., Figure 2 As shown in i), the expression for the characteristic function is as follows:

[0135] CF = 10 × (log 10 (1+D / V) / log 10 (2))

[0136] In the formula, CF is the initial arrival recognition feature function; V is the maximum amplitude of the signal; and D is the signal after discrimination filtering.

[0137] The automatic vibration event first-arrival design calculation process based on the multi-module filter ensemble algorithm of this invention is shown below. Figure 1 As shown.

[0138] To verify the performance of the recognition algorithm proposed in this invention in terms of both parameter dependence and noise suppression, the following experiments were conducted:

[0139] To investigate the dependence of the identification feature function on parameter settings, orthogonal verification experiments were conducted. This can be achieved by controlling the signal duration characteristics (d... min d max N d Parameters), signal frequency characteristics (f low f high N f The vibration first-arrival time was compared among three categories of parameters: parameters ), the penalty factor β for non-pulse signals, and d. Regarding the duration characteristics of the signal, the first-arrival time of vibration energy signals triggered by disaster events is generally within 10 seconds; therefore, d is set... min =1s,d max =3~10s, N d =1~5. Disaster event signals do not have dispersion characteristics; the energy propagation speed is the same across all frequency bands. Therefore, a more stable and clearer low-frequency band can be selected. l0w =1Hz, f high =3~10Hz, N f =1 to 5. The penalty factor β for non-pulse signals can be 1 to 3. The parameter settings are shown in Table 1.

[0140] Table 1. Characteristic parameters of orthogonal verification experiments

[0141]

[0142]

[0143] Figure 3This study demonstrates how changing the signal bandwidth reveals the characteristics of the characteristic function while keeping two parameter variables—signal duration characteristics and the penalty factor β for non-pulse signals—fixed. Overall, as the bandwidth increases, the peak values ​​of the identified characteristic function gradually increase, indicating a growing number of picked-up vibration events. However, with varying bandwidth, major events carrying more vibrational energy can be identified. The bandwidth (1-10Hz) selected in this experiment encompasses the energy frequency bands released by most natural disasters (volcanoes, earthquakes, landslides) and engineering disasters (pipeline ruptures, high-pressure gas leaks).

[0144] Figure 4 This demonstrates how, with fixed parameters for the signal filtering frequency band and the penalty factor β for non-pulse signals, the characteristics of the characteristic function can be identified by varying the signal duration. As the duration changes, the peak value of the characteristic function gradually becomes delayed, especially for durations longer than 7 seconds. Characteristic functions with shorter durations maintain a relatively consistent shape. For vibrations released by a disaster source, the energy release is a high-frequency pulse, resulting in a short vibration duration. For durations of 7 seconds or longer, the corresponding low-frequency (<1Hz) vibration source is not the type of disaster characterized by explosive seismic sources discussed in this invention. Therefore, the choice of duration does not significantly affect the acquisition of the initial arrival time of the vibration event. Figure 5 This displays the number of different frequency bands (N) of the signal. f Regarding the effect of frequency bands on the characteristic function, for a given frequency band width, different numbers of frequency bands have almost no impact on the shape of the characteristic function, especially at the initial arrival time of the main event. However, fewer frequency bands can reduce the algorithm's computation time. The same situation occurs with different durations (N). d The situation is as follows. Figure 6 This demonstrates the effect of the penalty factor β on the characteristic function of non-pulse signals while keeping the signal duration and frequency characteristics unchanged. For larger β values, only the initial arrival of strong energy events can be detected, while as the β value gradually decreases, more weak energy events are identified. The selection of pulse signals also becomes more stringent as the penalty factor β increases.

[0145] To assess the stability of the identification feature function under varying noise levels, verification tests were conducted using the same signal but with different noise levels. This experiment mixed the original signal with Gaussian random white noise, setting signal-to-noise ratios (SNR) of 5, 3, 2, and 1 respectively. The signal timing and feature function are shown below. Figure 7As shown, when the signal-to-noise ratio (SNR) is 5, the noise intensity has little impact on the signal, and the vibration characteristics of the signal are clear. Based on the characteristic function, the arrival times of the four main events can be clearly picked out, corresponding to the signal timing. As the SNR decreases to 3, the noise content in the signal increases, but the shape of the characteristic function remains basically unchanged, and the arrival times of the four main events are still obvious. When the SNR further decreases to 2, the characteristic curve only changes slightly. When the SNR is 1, that is, the amplitude of the signal is equal to the amplitude of the noise. Such signals are usually of low quality, and the vibration timing characteristics change significantly. The signal is almost annihilated by the noise, and the shape of the characteristic curve changes significantly compared to before. However, important events can still be picked out, especially the arrival time of the first event, which remains consistent with the high SNR case. This result shows that the noise reduction module in the algorithm of this invention effectively restores the signal shape.

[0146] To further apply this algorithm in extracting the first arrival of vibration events, it was used to pick up the first arrival of vibration signals induced by different seismic sources, and compared with the traditional short-time average / long-time average ratio (STA / LTA) algorithm. Figure 8 As shown, this is a microseismic signal generated when an underground water pipe ruptures. The energy arrival time is extracted to be 5.03 s using the multi-filter module integrated algorithm proposed in this invention. Simultaneously, the arrival time is obtained as 5.16 s using the commonly used third-order threshold algorithm of the traditional Short-Time Average / Long-Time Average Ratio (STA / LTA) algorithm, while the arrival time is obtained as 5.27 s using the fourth-order threshold algorithm. This signal has a good signal-to-noise ratio, thus the noise interference is relatively small, and the calculation results of the two algorithms are basically consistent.

[0147] like Figure 9 As shown, this is a near-field earthquake of magnitude 3.1 in Iburiku, Hokkaido, Japan. Because the seismograph was buried in the loose overburden near a nearby sewage treatment plant, the earthquake time series record suffered significant interference, with the earthquake event almost completely obscured by factory noise, resulting in a low signal-to-noise ratio. The energy arrival time was extracted to be 3.17s using the multi-filter module integrated algorithm proposed in this invention. However, the results calculated by the traditional Short-Time Average / Long-Time Average Ratio (STA / LTA) algorithm using different parameter thresholds showed significant differences: the third-order threshold algorithm yielded an arrival time of 1.29s, while the fourth-order threshold algorithm yielded an arrival time of 6.53s. The distance between the epicenter and the seismic station is about 3 km. The bedrock in the large area is mainly strongly weathered tuff. The first arrival P-wave velocity is 800-1000 m / s, and the theoretical first arrival time is 3-3.75 s. Although this theoretical value is calculated by a one-dimensional velocity model, it is close to the calculation result of the multi-filter module integrated algorithm proposed in this invention. Therefore, the algorithm of this invention is more reasonable.

[0148] In summary, this invention establishes an automatic picking algorithm integrating multiple filters. While reducing noise, it constructs a feature function, requiring only the most basic signal parameters and effectively reducing noise interference. This allows for the extraction of the initial arrival time of disaster events, further improving the accuracy of automatically picking up the occurrence of disaster events.

[0149] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. An automatic identification method for the initial arrival of vibration events based on multi-module filter integration, characterized in that, The steps are as follows: S1. The continuously recorded raw signal is preprocessed and then segmented. The effective frequency distribution of vibration energy is determined for each segment of the segmented signal. S2, Embedded anti-aliasing low-pass filter; S3. Signal Downsampling: Based on the low-pass filtering in the previous step, downsampling is performed, introducing a downsampling scaling factor R, expressed as: fsp=2f high In the formula, fs is the original sampling frequency of the signal; fsp is the sampling frequency after downsampling; f high The preset high-frequency cutoff frequency; S4, Embedded adaptive detrending term filter; this filter is expressed as: B=FIR(M,ω) c ) M=round(d max ×fsp) trend=filtfit(B,1,x r ) x rd =x r -trned In the formula, B is the filter, M is the cutoff angular frequency order, and ω c It is the cutoff angular frequency, d max The maximum duration of the disaster as preset; f low The preset low-frequency cutoff frequency; trend is the detrending filter, x r The signal is downsampled; x rd It is the signal after trend term filtering; S5. Embedded median filtering, median filtering is expressed as: y=medfilt(x rd ,round(fsp×d min / 10)) In the formula, d min The minimum duration of the disaster as preset; S6, Embedded edge-free delay filter; this filter is expressed as: y′=filter(H,1,y) In the formula, y' is the signal after edge-delay filtering; H is the filter parameter; and y is the signal after median filtering. S7. Calculate the signal envelope energy; S8. First, embed the discriminant filter, then establish the feature function. The expression of the feature function is as follows: CF=10×(log 10 (1+D / V) / log 10 (2)) In the formula, CF is the initial arrival recognition feature function; V is the maximum amplitude of the signal; and D is the signal after discrimination filtering. To assess the stability of the identification feature function under varying noise levels, verification tests were conducted using the same signal but with different noise levels. The original signal was mixed with Gaussian random white noise, and signal-to-noise ratios (SNRs) of 5, 3, 2, and 1 were set. At an SNR of 5, the noise intensity had minimal impact on the signal, and the vibration characteristics were clear. Based on the feature function, the arrival times of the four main events could be captured, corresponding to the signal timing. As the SNR decreased to 3, the noise content in the signal increased, but the shape of the feature function remained largely unchanged, and the arrival times of the four main events were still evident. When the SNR further decreased to 2, the feature curve showed only a slight change. When the SNR was 1, meaning the signal amplitude was equal to the noise amplitude, the shape of the feature curve changed significantly, but the arrival time of the first event could still be captured.

2. The automatic identification method for the initial arrival of vibration events based on multi-module filter integration as described in claim 1, characterized in that, In step S6, the filter parameter H is set as follows: H = acos(2π × fc × th) a = 1: length(th): 0 k=1,2,……,N f In the formula, f c For [f low f high The center frequency of [], th is the time series, a is the amplitude transformation coefficient, representing the number of th values ​​from 0 to 1; k is the number of frequency points, N f The number of frequency band windows that can be decomposed from a predefined frequency.

3. The automatic identification method for the initial arrival of vibration events based on multi-module filter integration as described in claim 1, characterized in that, In step S8, the expression for the discrimination filter is as follows: D(m)=filter(T,1,E k ) m=1,2,……,N d k=1,2,……,N f In the formula, D is the signal after discrimination filtering; T is an intermediate variable containing the penalty factor β for non-pulse signals; E k It is the signal envelope energy, where m is a sequence of natural numbers, representing values ​​from the first to the Nth. d 1, N d It is the number of time windows that can be decomposed over time; The intermediate variable T is represented as: T = β × ((1:L) - (L+1)) L = round(D(m) × fsp) D=[d min ,N d ,d max ] In the formula, β is the penalty factor for non-pulse signals, and L represents the data length of D.

4. The automatic identification method for the initial arrival of vibration events based on multi-module filter integration as described in claim 1, characterized in that, The specific steps of S1 are as follows: S11. Preprocess the continuously recorded raw signals, including baseline leveling, mean removal, and trend removal steps; S12. Cut the preprocessed original signal by overlapping and cutting it. S13. For each segment of the cutting signal, determine the effective impact time of vibration energy and set the duration of the impact: d min : The shortest duration of the disaster; d max : The longest duration of the disaster; N d The number of time windows that can be decomposed during the duration; S14. For each segment of the cutting signal, determine the effective frequency distribution of vibration energy, select the low-frequency band, and set the frequency band range: f low : Low-frequency cutoff frequency; f high : High-frequency cutoff frequency; N f : The number of frequency-decomposable bandwidth windows.

5. The automatic identification method for the initial arrival of vibration events based on multi-module filter integration as described in claim 4, characterized in that, In step S2, the anti-aliasing low-pass filter is represented as follows: In the formula: n is the filter order, ω c ω is the cutoff angular frequency; p This refers to the passband edge frequency.

6. The automatic identification method for the initial arrival of vibration events based on multi-module filter integration as described in claim 5, characterized in that, Specifically, step S7 involves: based on edge-free delay filtering, dividing the signal into sub-signals of different frequency bands according to different center frequencies, and then calculating the envelope energy E of each signal segment. k The calculation formula is as follows: k=1,2,……,N f In the formula, H is the Hilbert transform of the signal; y' is the signal after edge-delay filtering.

Images