P-wave first arrival time automatic picking method for multi-channel seismic signals based on first arrival time window

Abstract

The present application relates to a kind of based on first arrival variable time window's multiple channel seismic signal P-wave arrival time automatic picking method.Seismic P-wave first arrival picking is the key technology of unfolding seismic exploration, earthquake warning and other fields work.Facing the production requirement of the data quantity of acquisition is more and more big, need to carry out processing to multiple channel seismic signal in actual detection process.Manual picking method commonly used is time-consuming and not accurate enough.STA / LTA and other methods need to be input and adjusted respectively to each channel parameter, cannot realize full-automatic picking multiple first arrival time, therefore a kind of based on first arrival variable time window's multiple channel seismic signal P-wave arrival time automatic picking method is proposed.Aiming at the relationship between parameter and waveform, first arrival time window is optimized to select.Compared with existing natural seismic P-wave picking method, the present application not only has the help to improve positioning accuracy, and algorithm real-time performance is strong, accuracy is high, can realize multiple P-wave first arrival simultaneous picking in the case of zero input parameter.

Description

Technical Field

[0001] This invention belongs to the field of earthquake detection, and specifically refers to a method for automatically picking up the arrival time of P-waves of multi-channel seismic signals based on a first-arrival variable time window. Background Technology

[0002] Faced with ever-increasing production demands for data acquisition, the actual detection process requires processing multiple microseismic signals. For example, in natural earthquakes, due to differences in the properties of the subsurface medium, there can be significant errors in calculating the arrival time difference between channels. Furthermore, the large number and uneven distribution of seismic stations, along with varying sampling rates, complicates first-arrival acquisition. Manually acquiring large amounts of seismic data would be extremely wasteful of manpower. Methods such as the long-short window ratio (STA / LTA), envelope function method, and Akaike Information Criterion (AIC) are common P-wave arrival time acquisition methods that play a crucial role in practical applications. They not only quickly and effectively acquire P-wave first arrivals but also offer lower errors and higher accuracy compared to manual acquisition. However, for first-arrival acquisition of multiple seismic signals, these methods require separate parameter input and adjustment for each channel, resulting in low computational efficiency and failing to meet the practical needs of multi-channel seismic signal acquisition. Therefore, to achieve rapid first-arrival acquisition of multiple seismic signals, it is necessary to dynamically change the parameters of each channel based on the actual seismic wave conditions, achieving automatic, simple, and highly accurate multi-channel P-wave first-arrival acquisition.

[0003] CN103995290 proposes a high-precision automatic P-wave phase first arrival acquisition method for microseismic signals. This method utilizes the Hilbert transform to solve for the envelope of the microseismic signal, sets a threshold for the envelope to coarsely acquire the P-wave first arrival, and then selects a time window and uses the AIC function to accurately calculate the P-wave arrival time. However, the method suffers from subjective factors in parameter selection, which significantly impacts the final acquisition results. Therefore, parameter selection should be quantified as much as possible. Furthermore, the Hilbert transform for signal envelope acquisition exhibits endpoint effects, severely affecting the coarse acquisition results and consequently the accuracy of first arrival acquisition in multi-channel automatic acquisition. Additionally, this method can only effectively acquire one seismic signal with a single parameter input, and cannot simultaneously acquire the P-wave arrival times of multiple seismic signals.

[0004] In summary, there is a need for a first-arrival acquisition method that can overcome the end-point effect while ensuring P-wave acquisition accuracy and allows for parameter self-selection. To this end, this invention proposes an automatic first-arrival acquisition method for multi-channel seismic signal P-waves based on a variable first-arrival time window. This method divides acquisition into pre-acquisition and precise acquisition. Pre-acquisition selects different time windows of varying sizes near the first arrival, tailored to the characteristics of different channels, and performs precise acquisition of signals within the time window range. In multi-channel acquisition, this method adjusts the time window size according to the changing characteristics of different channels, achieving both high acquisition accuracy and simultaneous multi-channel acquisition, thus solving the problem of low acquisition efficiency in existing technologies for multi-channel microseismic P-wave acquisition. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of the prior art by providing an automatic acquisition method for the arrival time of P-waves in multi-channel seismic signals based on a first-arrival variable time window.

[0006] The key technical points of this invention are as follows: First, the single-channel seismic signal is transformed by Hilbert to obtain an envelope signal. Since the first peak point of the envelope signal is closest to the first arrival of the earthquake, a peak detection threshold is used to find the first peak point. During this process, the endpoint effect can affect the selection of the peak point, so it is necessary to determine whether the endpoint effect affects the picking method and to develop a method to smooth the endpoint effect. Then, based on the characteristics of different channels, the size of the picking window is designed using a right boundary detection threshold of the picking range to ensure that the first arrival time can be accurately picked within the window. The arrival time of the first arrival P-wave is picked using the AIC method to improve the accuracy of precise picking. Finally, the P-wave picking algorithm traverses each channel of microseismic signal to achieve multi-channel P-wave picking with varying first arrival time windows. Experiments show that this method can achieve multi-channel P-wave picking with varying first arrival time windows. Compared with existing P-wave picking methods, this method not only has high computational accuracy and strong real-time performance, but also can achieve simultaneous picking of multiple P-wave first arrivals with zero input parameters.

[0007] The objective of this invention is achieved through the following technical solution:

[0008] An automatic acquisition method for P-wave arrival times of multi-channel seismic signals based on first-arrival variable time windows includes the following steps:

[0009] a. Read the vertical component seismic signal X(t) received by multiple seismic detector arrays, with a total number of traces of K. Then X(t) can be represented by a set as X(t) = {x1(t), x2(t), ..., x...} K (t)}, where column vector x k (t) represents the k-th signal, k = 1, 2, ..., K, t = 1, 2, ..., N, where N is the number of sampling points for each seismic signal, and the sampling rate is Fs. Additionally, x... k The number of background noise sampling points collected before the earthquake in (t) is denoted as b.k ;

[0010] b. Let k = 1;

[0011] c. Regarding x k (t) Perform Hilbert transform:

[0012]

[0013]

[0014]

[0015] Among them, y (k)max With y (k)min y k The maximum and minimum values ​​of (t);

[0016] d. Define the peak detection threshold

[0017]

[0018] e. Column vector C has a size of K×1, and all elements have a value of 0.05. Defined...

[0019] y k "(t)=y k ′(t)-C (5)

[0020] If y k "(1) < 0, then it is not necessary to adjust y." k ″(t) is processed; otherwise, y is searched. k The first value less than 0 in (t) is denoted as the sampling point number corresponding to that value. k m k ∈[1, N) and

[0021]

[0022] f, regarding y k Binarize ″(t) and define

[0023]

[0024] D k The position where four or more consecutive "1"s appear in (t) is denoted as x. k Find the valid signal position of (t), and locate the start and end sampling point numbers s of the first valid signal. k and e k s k e k ∈(1, N) and

[0025] g, search y k "(t) in [s k e k The maximum value within the range is denoted as P. k The corresponding sampling point number is r k ;

[0026] h. Define the right boundary detection threshold of the picking range.

[0027] Q k =0.5(P) k +0.05)-0.05 (8)

[0028]

[0029] Looking for y k The value in (t) that satisfies all the conditions of equation (9) is the sampling point number corresponding to that value. k ;

[0030] i. Let Δn k =r k -u k ;

[0031] j, let z k =3, "·" indicates multiplication, calculated in [u k -z k ·Δn k u k -(z k -1)·Δn k Within the range of -1], y k "(t) is greater than 0.1Q" k Number of sampling points v k ,

[0032]

[0033] If v k If equation (10) is not satisfied, then z k =z k +1, recalculate v k And determine whether equation (10) is satisfied, until equation (10) is satisfied, if v k If equation (10) is satisfied, proceed to step k;

[0034] k, the range of the initial arrival and arrival times [u] k -z k ·Δn k u k Define a one-dimensional vector.

[0035] W k=[x k (u k -z k ×Δn k ), x k (u k -z k ×Δn k +1), ..., x k (u k )] T (11)

[0036] l. Using AIC functions for precise picking

[0037] AIC k (h)=h·lg(Var(W k [1, h]))+(f k -h-1)·lg(Var(W k [h+1, f k (12)

[0038] Among them, f k For W k length, f k =z k ·Δn k +1, h = 1, 2, ..., f k Var(W k [1, h]) refers to W k (1) To W k The variance of h sampling points between (h);

[0039] m, Search AIC k The sampling point number a corresponding to the minimum value of the finite rational number in (h) k Then the sampling point number of the first arrival of the longitudinal wave is (u k -z k ·Δn k +z k -1);

[0040] n. Definition

[0041]

[0042] Here, Time k It is x k (t) When the longitudinal wave first arrives;

[0043] o. If k < K, let k = k + 1 and repeat steps c ~ o. If k = K, then the multi-channel seismic P-wave pickup is completed.

[0044] Beneficial effects:

[0045] Through experimentation, this invention discloses an automatic P-wave arrival time acquisition method based on a variable first-arrival time window for multi-channel seismic signals. This method can acquire the first-arrival phase and arrival time of P-waves from multiple seismic signals. Existing high-precision methods for acquiring first-arrival P-waves are mainly geared towards single-channel microseismic first-arrival acquisition. While acquiring first-arrival phases from single-channel seismic signals offers relatively high accuracy, it suffers from low computational efficiency. Complex algorithms can acquire first-arrival phases from multiple channels, but their computational complexity is high. This new method, while ensuring P-wave acquisition accuracy, automatically adjusts the acquisition time window and acquires first-arrival phases from multiple seismic P-waves, achieving a balance between acquisition accuracy and simultaneous multi-channel acquisition. This addresses the problem of low acquisition efficiency in existing technologies for seismic P-wave acquisition. Attached Figure Description

[0046] Figure 1 The results of this method, which extracts two data points after picking up the first arrival of 60 signals from an actual natural earthquake set, are compared with the results of STA / LTA picking and the arrival results provided by SAGE (c), and the data near the first arrival are magnified (a and b).

[0047] Figure 2 The picking errors of the method of this invention compared to STA / LTA and SAGE on actual natural earthquake sets are shown in the figure (circles represent the picking error of the method of this invention; inverted triangles represent the picking error of STA / LTA; and crosses represent the picking error calculated based on the arrival time provided by the SAGE official website). Detailed implementation method:

[0048] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0049] An automatic acquisition method for P-wave arrival times of multi-channel seismic signals based on first-arrival variable time windows includes the following steps:

[0050] a. A magnitude 6.5 earthquake signal that occurred in the Hindu Kush Region, Afghanistan, at 16:47:23 on March 21, 2023, was acquired by NSF's Seismological Facility for the Advancement of Geoscience (SAGE). The vertical component seismic signal X(t) received by multiple seismic detector arrays was read, with a total trace count of K = 60. X(t) can be represented by the set X(t) = {x1(t), x2(t), ..., x...}. K (t)}, where column vector x k(t) represents the k-th signal, k = 1, 2, ..., K, t = 1, 2, ..., N, where N is the number of sampling points for each seismic signal, and the sampling rate is Fs. Additionally, x... k The number of background noise sampling points collected before the earthquake in (t) is denoted as b. k ;

[0051] b. Let k = 1;

[0052] c. Regarding x k (t) Perform Hilbert transform:

[0053]

[0054]

[0055]

[0056] Among them, y (k)max With y (k)min y k The maximum and minimum values ​​of (t);

[0057] d. Define the peak detection threshold

[0058]

[0059] e. Column vector C has a size of K×1, and all elements have a value of 0.05. Defined...

[0060] y k "(t)=y k ′(t)-C (5)

[0061] If y k "(1) < 0, then it is not necessary to adjust y." k ″(t) is processed; otherwise, y is searched. k The first value less than 0 in (t) is denoted as the sampling point number corresponding to that value. k m k ∈[1, N) and

[0062]

[0063] f, regarding y k Binarize ″(t) and define

[0064]

[0065] D k The position where four or more consecutive "1"s appear in (t) is denoted as x. kFind the valid signal position of (t), and locate the start and end sampling point numbers s of the first valid signal. k and e k s k e k ∈(1, N) and

[0066] g, search y k "(t) in [s k e k The maximum value within the range is denoted as P. k The corresponding sampling point number is r k ;

[0067] h. Define the right boundary detection threshold of the picking range.

[0068] Q k =0.5(P) k +0.05)-0.05 (8)

[0069]

[0070] Looking for y k The value in (t) that satisfies all the conditions of equation (9) is the sampling point number corresponding to that value. k ;

[0071] i. Let Δn k =r k -u k ;

[0072] j, let z k =3, "·" indicates multiplication, calculated in [u k -z k ·Δn k u k -(z k -1)·Δn k Within the range of -1], y k "(t) is greater than 0.1Q" k Number of sampling points v k ,

[0073]

[0074] If v k If equation (10) is not satisfied, then z k =z k +1, recalculate v k And determine whether equation (10) is satisfied, until equation (10) is satisfied, if v k If equation (10) is satisfied, proceed to step k;

[0075] k, the range of the initial arrival and arrival times [u] k -z k ·Δn k u k Define a one-dimensional vector.

[0076] W k =[x k (u k -z k ×Δn k ), x k (u k -z k ×Δn k +1), ..., x k (u k )] T (11)

[0077] l. Using AIC functions for precise picking

[0078] AIC k (h)=h·lg(Var(W k [1, h]))+(f k -h-1)·lg(Var(W k [h+1, f k (12)

[0079] Among them, f k For W k length, f k =z k ·Δn k +1, h = 1, 2, ..., f k Var(W k [1, h]) refers to W k (1) To W k The variance of h sampling points between (h);

[0080] m, Search AIC k The sampling point number a corresponding to the minimum value of the finite rational number in (h) k Then the sampling point number of the first arrival of the longitudinal wave is (u k -z k ·Δn k +a k -1);

[0081] n. Definition

[0082]

[0083] Here, Time k It is x k (t) When the longitudinal wave first arrives;

[0084] o. If k < K, let k = k + 1 and repeat steps c ~ o. If k = K, then the multi-channel seismic P-wave pickup is completed.

[0085] After acquiring the first arrival of 60 seismic signals, two seismic signals were extracted from the 60 signals, and the corresponding seismic station names were MAKZ and KURK. Figure 1 c represents the acquisition result of extracting two channels after multi-channel first arrival acquisition of 60 channels using this method, and is compared with the acquisition results of STA / LTA and the arrival results provided by SAGE. To observe the acquisition effect, magnification is performed near the first arrival of the two stations. Figure 1 (the part in box c), respectively corresponding to Figure 1 b and Figure 1 a. It can be clearly seen that this invention can not only achieve multi-channel first arrival pickup, but also has high pickup accuracy, and the pickup results are better than the pickup results of STA / LTA and the arrival results provided by SAGE.

[0086] To more intuitively represent the results of first arrival picking for each P-wave, the errors of the picking results obtained by the method of this invention, the picking results obtained by STA / LTA, and the first arrival results provided by the SAGE website, as well as the results obtained manually, were calculated. Figure 2 Table 1 shows the basic statistical parameters of the absolute error of the natural earthquake dataset. By comparison, the picking results of the method proposed in this invention in the natural earthquake dataset are basically close to those of the manual picking results, relatively stable and generally better than the first arrival time and STA / LTA results provided by SAGE, and are more suitable for P-wave picking of multichannel natural earthquakes.

[0087] Table 1. Basic statistical parameters of the absolute error of natural earthquake sets (in seconds).

[0088] method Minimum error Maximum error median average Standard deviation Method of the present invention 0 0.0200 0.0033 0.0047 0.0045 STA / LTA 0 0.0279 0.0065 0.0072 0.0061 SAGE 4.1667e-04 1 0.0373 0.0672 0.1756

Claims

1. A method for automatically picking up the arrival time of P-waves in multi-channel seismic signals based on a first-arrival variable time window, characterized in that, Includes the following steps: a. Read the vertical component seismic signal X(t) received by multiple seismic detector arrays, where t is the sampling point number and the total number of traces is K. Then X(t) can be represented by a set as X(t) = {x1(t), x2(t), ..., x...} K (t)}, where column vector x k (t) represents the k-th signal, k = 1, 2, ..., K, t = 1, 2, ..., N, where N is the total number of sampling points for each seismic signal, and the sampling rate is Fs. Additionally, x... k The number of background noise sampling points collected before the earthquake in (t) is denoted as b. k ; b. Let k = 1; c. Regarding x k (t) Perform Hilbert transform: Where τ is the integration variable, y (k)max With y (k)min y k The maximum and minimum values ​​of (t); d. Define the peak detection threshold e. Column vector C has a size of K×1, and all elements have a value of 0.

05. Defined... y k ″(t)=y k ′(t)-C (5) If y k "(1) < 0, then it is not necessary to adjust y." k ″(t) is processed; otherwise, y is searched. k The first value less than 0 in (t) is denoted as the sampling point number corresponding to that value. k m k ∈[1, N) and f, regarding y k Binarize ″(t) and define D k The position where four or more consecutive "1"s appear in (t) is denoted as x. k Find the valid signal position of (t), and locate the start and end sampling point numbers s of the first valid signal. k and e k s k e k ∈(1, N) and g, search y k "(t) in [s k e k The maximum value within the range is denoted as P. k The corresponding sampling point number is r k ; h. Define the right boundary detection threshold of the picking range. Q k =0.5(P k +0.05)-0.05 (8) Looking for y k The value in (t) that satisfies all the conditions of equation (9) is the sampling point number corresponding to that value. k ; i, order Δn k = r k -u k ; j, let z k =3, "·" indicates multiplication, calculated in [u k -z k ·Δn k u k -(z k -1)·Δn k Within the range of -1], y k "(t) is greater than 0.1Q" k Number of sampling points v k , If v k If equation (10) is not satisfied, then z k =z k +1, recalculate v k And determine whether equation (10) is satisfied, until equation (10) is satisfied, if v k If equation (10) is satisfied, proceed to step k; k, the range of the initial arrival and arrival times [u] k -z k ·Δn k u k Define a one-dimensional vector. W k =[x k (u k -z k ×Δn k ),x k (u k -z k ×Δn k +1),...,x k (u k )] T (11) 1. Use AIC functions for precise picking AIC k (h)=h·lg(Var(W k [1,h]))+(f k -h-1)·lg(Var(W k [h+1,f k ])) (12) Among them, f k For W k length, f k =z k ·Δn k +1, h = 1, 2, ..., f k Var(W k [1, h]) refers to W k (1) To W k The variance of h sampling points between (h); m, Search AIC k The sampling point number a corresponding to the minimum value of the finite rational number in (h) k Then the sampling point number of the first arrival of the longitudinal wave is (u k -z k ·Δn k +a k -1); n. Definition Here, Time k It is x k (t) When the longitudinal wave first arrives; o. If k < K, let k = k + 1 and repeat steps c ~ o. If k = K, then the multi-channel seismic P-wave pickup is completed.

Images