High-robustness P-wave seismic phase first arrival picking method

Abstract

The invention discloses a high-robustness P-wave seismic phase first arrival picking method, which is used for identifying a P-wave seismic phase by adopting a three-step picking method and comprises the following steps of: 1, performing one-dimensional multi-scale wavelet decomposition on an original vertical acceleration record, and reconstructing the acceleration record; step 2, respectively searching moments greater than a trigger threshold value for the original and reconstructed acceleration records by adopting an STA / LTA method, comparing differences, judging whether a signal is a P-wave seismic phase or an interference signal, and preliminarily picking up a P-wave first arrival moment; and 3, taking a section of record from the first arrival time point to the front and back of the original vertical acceleration record, and accurately picking up the P-wave first arrival time by applying the AIC criterion. In order to eliminate the influence of an interference signal on a first arrival P wave seismic phase, the acceleration record is reconstructed by using a large-scale wavelet coefficient, and the difference between the first arrival time point of the original acceleration record and the first arrival time point of the reconstructed acceleration record is compared, so that the accuracy of the pickup method is improved. The method can effectively eliminate the influence of various interference signals, is simple in calculation, and meets the requirement of earthquake early warning timeliness.

Description

Technical Field

[0001] The invention relates to the technical field of earthquake early warning, and in particular to a method for picking up the first arrival of a P-wave seismic phase with strong robustness. Background Art

[0002] In earthquake early warning technology, the identification of the first arrival time of the P-wave phase is the basis for subsequent earthquake location, magnitude measurement and other work. Therefore, reliably and accurately identifying the first arrival time of the P-wave phase is the first key problem faced by earthquake early warning technology. The identification of the first arrival time of the P-wave phase is mainly based on the different characteristics of the seismic waveform and the interference signal. At present, commonly used methods include polarization analysis, energy analysis, autoregression analysis, wavelet transform analysis, etc. These methods extract different features of the seismic waveform, while the same features of the interference signal are relatively weak, so that the first arrival time of the P-wave phase can be identified.

[0003] However, the single method in the prior art only extracts a certain feature of the seismic waveform signal and is easily affected by interference signals. From the perspective of the timeliness of the recognition of the first arrival moment of the P-wave seismic phase, the method for recognizing the first arrival moment of the seismic phase needs to be computationally efficient. Therefore, there is ample room for research on how to combine different recognition methods, stably analyze the differences in the characteristics of the seismic waveform signal and the interference signal from multiple angles, and then conduct a comprehensive analysis and judgment, and finally reliably, accurately and quickly recognize the first arrival moment of the P-wave seismic phase. Summary of the invention

[0004] Purpose of the invention: To solve the problems mentioned in the background technology, the present invention provides a robust P-wave seismic phase first arrival picking method. The wavelet transform is comprehensively utilized to analyze and directly decompose the p-wave and use it as a limiting condition, and then the characteristics of the seismic waveform signal and the interference signal are distinguished by energy analysis and autoregressive analysis methods, so as to reliably and accurately identify the first arrival time of the P-wave seismic phase quickly and accurately. By using the integration method of the present invention, the result can be obtained by inputting the original data without complex calculations, thereby ensuring its robustness.

[0005] Technical solution:

[0006] The present invention discloses a robust P-wave seismic phase first arrival picking method, the method comprising the following steps:

[0007] S1; Perform one-dimensional multi-scale wavelet decomposition on the original vertical acceleration record to obtain multi-scale wavelet coefficients, and select the highest-scale wavelet coefficient with the least noise influence to reconstruct the acceleration record.

[0008] S2: Calculate the ratio of the short-time window average value of the characteristic function to the long-time window average value of the original vertical acceleration record and the reconstructed acceleration record, and take the moment when the ratio is greater than a preset threshold as the initial estimated P wave first arrival time.

[0009] S3: The estimated P wave first arrival time is limited, and the duration window and critical threshold are added. If the STA / LTA ratio is greater than the critical threshold within the duration window, the next step of picking is performed, otherwise it is judged as an interference signal.

[0010] S4: Compare the difference between the initial arrival time points of the original vertical acceleration record and the reconstructed acceleration record. If the initial arrival time difference between the two is less than the set value, proceed to the next step of picking up, otherwise it is judged as an interference signal.

[0011] S5; Taking the initially estimated P-wave first arrival time as a reference, take a certain time length forward and backward records of the original vertical acceleration records, and use the AIC criterion to find the minimum value as the final P-wave first arrival time.

[0012] Furthermore, the specific steps of S1 are as follows:

[0013] The original vertical acceleration record is subjected to one-dimensional multi-scale wavelet decomposition using the sym4 wavelet function to obtain the wavelet decomposition structure, and the wavelet coefficients at the highest scale are extracted. The acceleration record is reconstructed using the wavelet coefficients at the highest scale to suppress the influence of noise. The wavelet decomposition processing formula is as follows:

[0014]

[0015] Among them, H is a low-pass filter, G is a high-pass filter, j is the number of decomposition layers, c is the scale coefficient, d is the wavelet coefficient, and the original acceleration record is C0.

[0016] The wavelet coefficient reconstruction process is shown as follows:

[0017] C J =H * C J+1 +G * D J+1 ,J=j-1,j-2,...,1

[0018] C J Represents the reconstructed record, H * and G * represents the dual operator of H and G.

[0019] Furthermore, the specific steps of S2 are as follows:

[0020] For the original vertical acceleration record and the reconstructed acceleration record, the characteristic function CFp is obtained, and the processing formula is as follows:

[0021] CF p (k) = x(k) 2 +[x(k)-x(k-1)] 2

[0022] Among them, CFp is the characteristic function, x(k) represents the acceleration record, k represents the current time point of the record, and two records of different lengths are taken forward from the current time point as the boundary to calculate STA(k) / LTA(k) of the characteristic function. The processing formula is as follows:

[0023]

[0024] Wherein, k1 and k2 are a certain time point before the current time i, requiring k2<k1<i, and the time when the STA(i) / LTA(i) ratio is greater than the set threshold is taken as the initial estimated P wave first arrival time.

[0025] Furthermore, the picking of the first arrival time of the P wave in S2 is added to the limitation of S3, and the specific steps are as follows:

[0026] The initially estimated P wave arrival time is judged, and the duration window and critical threshold constraints are added. The initially estimated P wave arrival time is taken as the starting point, and a certain length of records are taken as the duration window. Within the duration window, if the ratio of STA(i) / LTA(i) is greater than the critical threshold, this moment is taken as the P wave arrival time for the next calculation, otherwise this moment is considered to be an interference signal.

[0027] Furthermore, the picking of the first arrival time of the P wave in S2 is added to the limitation of S4, and the specific steps are as follows:

[0028] Compare the P wave first arrival time of the original vertical acceleration record and the reconstructed acceleration record. If the difference between the first arrival time of the two is less than the set value, this time is used as the P wave first arrival time; otherwise, the first arrival time is judged to be an interference signal.

[0029] Furthermore, the specific steps of S5 are:

[0030] The first arrival time of the P wave determined in the above steps is used as a reference. Records of a certain length of time are taken forward and backward from the original acceleration record to calculate the value of the AIC criterion. The processing formula is as follows:

[0031] AIC(k)=klog{var(x[1,k])}+(Nk-1)log{(var(x[k+1,N])}

[0032] Among them, k represents the recorded time point, var represents the variance, and the point corresponding to the minimum AIC value is taken as the final P wave first arrival time.

[0033] Beneficial effects:

[0034] The present invention utilizes the vertical acceleration records observed by the station, adopts a three-step picking method to identify the first arrival time of the P-wave seismic phase, and comprehensively utilizes wavelet transform analysis, energy analysis, and autoregressive analysis methods to extract the characteristics of the seismic waveform signal. By applying three limiting conditions, the influence of interference signals is eliminated. The method is simple and efficient in calculation, and does not require complex calculations. It provides a reference for the identification of the first arrival time of the P-wave seismic phase in the implementation of earthquake early warning technology. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 is a flow chart of the method of the present invention;

[0036] Figure 2 This is a schematic diagram of the original vertical acceleration recording of the first embodiment of the present invention;

[0037] Figure 3 This is the original vertical acceleration record diagram of 0-15s of the first embodiment of the present invention;

[0038] Figure 4 This is the 0-15s reconstructed acceleration record diagram of the first embodiment of the present invention;

[0039] Figure 5 Schematic diagram of the values ​​of STA / LTA and trigger judgment results of the original acceleration record and the reconstructed acceleration record in the first embodiment of the present invention;

[0040] Figure 6 This is a schematic diagram of the AIC value of the first embodiment of the present invention;

[0041] Figure 7 This is a diagram showing the initial arrival and final identification result of the P-wave seismic phase according to the first embodiment of the present invention;

[0042] Figure 8 This is a schematic diagram of the original vertical acceleration recording of the second embodiment of the present invention;

[0043] Fig. 9 This is the original vertical acceleration record diagram of 0-18s of the second embodiment of the present invention;

[0044] Fig.10 This is the 0-18s reconstructed acceleration record diagram of the second embodiment of the present invention;

[0045] Fig.11 Schematic diagram of the values ​​of STA / LTA and trigger judgment results of the original acceleration record and the reconstructed acceleration record in the second embodiment of the present invention;

[0046] Fig.12 This is a schematic diagram of the AIC value of the second embodiment of the present invention;

[0047] Fig.13 This is a diagram showing the first arrival and final identification result of the P-wave seismic phase according to the second embodiment of the present invention; DETAILED DESCRIPTION

[0048] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. The following embodiments are only for illustrating the technical concept and features of the present invention, and their purpose is to enable people familiar with the technology to understand the content of the present invention and implement it accordingly, and they cannot be used to limit the protection scope of the present invention. Any equivalent or modification made according to the spirit of the present invention should be included in the protection scope of the present invention.

[0049] The first embodiment of the present invention, as Figure 1 As shown, the P-wave seismic phase first arrival picking method includes the following steps:

[0050] S1; performing one-dimensional multi-scale wavelet decomposition on the original vertical acceleration record to obtain multi-scale wavelet coefficients, and selecting the highest-scale wavelet coefficient with the least noise influence to reconstruct the acceleration record;

[0051] The wavelet decomposition processing formula is as follows:

[0052]

[0053] Among them, H is a low-pass filter, G is a high-pass filter, j is the number of decomposition layers, c is the scale coefficient, d is the wavelet coefficient, and the original acceleration record is C0.

[0054] The wavelet coefficient reconstruction process is shown as follows:

[0055] C J =H * C J+1 +G * D J+1 ,J=j-1,j-2,...,1

[0056] C J Represents the reconstructed record, H * and G * represents the dual operator of H and G.

[0057] In this embodiment, if Figure 2 As shown in the figure, the sym4 wavelet function is used to perform one-dimensional 5-scale wavelet decomposition on the vertical acceleration record, as shown in the figure. Figure 3 As shown, in order to show the interference situation, the vertical acceleration record only retains the original vertical acceleration record 15s ago. This record contains the P wave first arrival phase and interference situation, and the acceleration amplitude is limited to between -0.1 and 0.1.

[0058] In order to show the effect of high-scale wavelet coefficient reconstruction on suppressing interference, only the first 15 seconds of the reconstructed acceleration record are taken, and the amplitude of the acceleration is limited to between -0.1 and 0.1. The reconstructed acceleration record at scale 5 is as follows: Figure 4 shown.

[0059] S2; calculating the ratio of the short-time window average value of the characteristic function to the long-time window average value for the original vertical acceleration record and the reconstructed acceleration record, and taking the time when the ratio is greater than a preset threshold as the initial estimated P wave first arrival time;

[0060] The specific steps are as follows:

[0061] For the original vertical acceleration record and the reconstructed acceleration record, obtain the characteristic function CF p , the processing formula is as follows:

[0062] CF p (k) = x(k) 2 +[x(k)-x(k-1)] 2

[0063] Among them, CF p is the characteristic function, x(k) represents the acceleration record, and k represents the current time point of the record. Taking the current time point as the boundary, take two records of different lengths forward and calculate STA(k) / LTA(k) of the characteristic function. The processing formula is as follows:

[0064]

[0065] Where k1 and k2 are time points before the current time i, and k2 < k1 < i is required. The time when the STA(i) / LTA(i) ratio is greater than the set threshold is taken as the initial estimated P wave first arrival time.

[0066] S3: limit the initial estimated P wave arrival time, add the duration window and critical threshold judgment, if within the duration window, the STA / LTA ratio is greater than the critical threshold, then proceed to the next step of picking, otherwise it is judged as an interference signal; that is, the picking of the P wave arrival time in S2 is added to the limitation of S3, the specific steps are:

[0067] The estimated P wave first arrival time is determined, and the duration window and critical threshold constraints are added:

[0068] The preset threshold is 5, the critical threshold is 5, and the duration window is 0.5s. If the STA / LTA ratio is greater than the critical threshold 5 within 0.5s after the initial estimation time, the next step of picking is performed. The STA / LTA values ​​of the original acceleration record and the reconstructed acceleration record and the results of the comprehensive judgment for each trigger are as follows: Figure 5 As shown in the figure, the estimated P wave first arrival time is taken as the starting point, and a certain length of records are taken as the duration window. Within the duration window, if the ratio of STA(i) / LTA(i) is greater than the critical threshold, this moment is taken as the P wave first arrival time for the next calculation. Otherwise, this moment is considered to be an interference signal.

[0069] S4: Compare the difference between the initial arrival time points of the original vertical acceleration record and the reconstructed acceleration record. If the initial arrival time difference between the two is less than the set value, the next step of picking is performed, otherwise it is judged as an interference signal; that is, the picking of the initial arrival time of the P wave in S2 is added to the limitation of S4, and the specific steps are:

[0070] Compare the P wave first arrival time of the original vertical acceleration record and the reconstructed acceleration record. If the difference between the first arrival time of the two is less than the set value, this time is used as the P wave first arrival time. Otherwise, the first arrival time is judged to be an interference signal.

[0071] In this embodiment, if the difference between the two is less than 1s, the next step of picking up is performed. Otherwise, it is judged as an interference signal, such as Figure 5 shown.

[0072] S5; Taking the initially estimated P-wave first arrival time as a reference, take a certain time length forward and backward records of the original vertical acceleration records, and use the AIC criterion to find the minimum value as the final P-wave first arrival time.

[0073] The specific steps are:

[0074] Taking the first arrival time of the P wave determined in the above steps as a reference, take a certain length of time forward and backward records of the original acceleration record to calculate the value of the AIC criterion. The processing formula is as follows:

[0075] AIC(k)=klog{var(x[1,k])}+(Nk-1)log{(var(x[k+1,N])}

[0076] Where k represents the recorded time point, var represents the variance, and the point corresponding to the minimum AIC value is taken as the final P wave first arrival time.

[0077] In this embodiment, the estimated P wave first arrival time is taken as a reference, 1.5s is taken forward and 1s is taken backward, and the minimum value is obtained by applying the AIC criterion as the final P wave first arrival time. The AIC value obtained by intercepting a section of the original acceleration record with the initial P wave phase first arrival time as a reference is as follows: Figure 6 As shown in Figure 1, in order to show the characteristics of the P-wave phase, only the original vertical acceleration record before 15 seconds is retained, and the amplitude of the acceleration is limited to between -0.1 and 0.1. The obtained P-wave phase initial and final identification results are as follows: Figure 7 shown.

[0078] The second embodiment of the present invention is different from the first embodiment in that Figure 8 As shown in the figure, the sym4 wavelet function is used to perform one-dimensional 5-scale wavelet decomposition on the vertical acceleration record, as shown in the figure. Fig. 9As shown, the vertical acceleration record only retains the original vertical acceleration record before 18 seconds, which contains the initial arrival phase of the P wave and the interference situation, and limits the acceleration amplitude to between -0.05 and 0.05;

[0079] like Fig.10 As shown, in order to show the effect of high-scale wavelet coefficient reconstruction on suppressing interference, only the first 18s of the reconstructed acceleration record is taken, and the amplitude of the acceleration is limited to between -0.003 and 0.003;

[0080] The STA / LTA values ​​of the original acceleration record and the reconstructed acceleration record of the second embodiment and the result of comprehensive judgment for each trigger are as follows: Fig.11 As shown in the figure, the AIC value obtained by cutting a section of the original acceleration record with the initial arrival time of the P-wave phase as the reference is as follows: Fig.12 shown.

[0081] In order to show the characteristics of the P-wave phase, only the original vertical acceleration record before 18 seconds is retained, and the amplitude of the acceleration is limited to between -0.05 and 0.05. The final identification results of the P-wave phase are as follows: Fig.13 shown.

[0082] The above embodiments are only for illustrating the technical concept and features of the present invention, and their purpose is to enable people familiar with the technology to understand the content of the present invention and implement it accordingly, and they cannot be used to limit the protection scope of the present invention. Any equivalent transformation or modification made according to the spirit of the present invention should be included in the protection scope of the present invention.

Claims

1. A robust P-wave seismic phase first arrival picking method, characterized in that: The method comprises the following steps: S1; performing one-dimensional multi-scale wavelet decomposition on the original vertical acceleration record to obtain multi-scale wavelet coefficients, and selecting the highest-scale wavelet coefficient with the least noise influence to reconstruct the acceleration record; S2; calculating the ratio of the short-time window average value of the characteristic function to the long-time window average value for the original vertical acceleration record and the reconstructed acceleration record, and taking the time when the ratio is greater than a preset threshold as the initial estimated P wave first arrival time; S3: limit the initial estimated P wave arrival time, add the duration window and critical threshold determination, if within the duration window, the STA / LTA ratio is greater than the critical threshold, then proceed to the next step of picking, otherwise it is judged as an interference signal; S4; comparing the difference between the initial arrival time points of the original vertical acceleration record and the reconstructed acceleration record, if the initial arrival time difference between the two is less than the set value, proceed to the next step of picking up, otherwise it is judged as an interference signal; S5; Taking the initially estimated P-wave first arrival time as a reference, take a certain time length forward and backward records of the original vertical acceleration records, and use the AIC criterion to find the minimum value as the final P-wave first arrival time.

2. A robust P-wave seismic phase first arrival picking method according to claim 1, characterized in that: The specific steps of S1 are as follows: The original vertical acceleration record is subjected to one-dimensional multi-scale wavelet decomposition using the sym4 wavelet function to obtain the wavelet decomposition structure, and the wavelet coefficients at the highest scale are extracted. The acceleration record is reconstructed using the wavelet coefficients at the highest scale to suppress the influence of noise. The wavelet decomposition processing formula is as follows: Among them, H is a low-pass filter, G is a high-pass filter, j is the number of decomposition layers, c is the scale coefficient, d is the wavelet coefficient, and the original acceleration record is C0. The wavelet coefficient reconstruction process is shown as follows: C J =H * C J+1 +G * D J+1 ,J=j-1,j-2,...,1 C J Represents the reconstructed record, H * and G * represents the dual operator of H and G.

3. A robust P-wave seismic phase first arrival picking method according to claim 1, characterized in that: The specific steps of S2 are as follows: For the original vertical acceleration record and the reconstructed acceleration record, obtain the characteristic function CF p , the processing formula is as follows: CF p (k)=x(k) 2 +[x(k)-x(k-1)] 2 Among them, CF p is the characteristic function, x(k) represents the acceleration record, k represents the current time point of the record, and two records of different lengths are taken forward from the current time point as the boundary to calculate STA(k) / LTA(k) of the characteristic function. The processing formula is as follows: Wherein, k1 and k2 are a certain time point before the current time i, requiring k2<k1<i, and the time when the STA(i) / LTA(i) ratio is greater than the set threshold is taken as the initial estimated P wave first arrival time.

4. A robust P-wave seismic phase first arrival picking method according to claim 1 or 3, characterized in that: The picking of the first arrival time of the P wave in S2 is added to the limitation of S3. The specific steps are: The initially estimated P wave arrival time is judged, and the duration window and critical threshold constraints are added. The initially estimated P wave arrival time is taken as the starting point, and a certain length of records are taken as the duration window. Within the duration window, if the ratio of STA(i) / LTA(i) is greater than the critical threshold, this moment is taken as the P wave arrival time for the next calculation, otherwise this moment is considered to be an interference signal.

5. A robust P-wave seismic phase first arrival picking method according to claim 1 or 3, characterized in that: The picking of the first arrival time of the P wave in S2 is added to the limitation of S4. The specific steps are: Compare the P wave first arrival time of the original vertical acceleration record and the reconstructed acceleration record. If the difference between the first arrival time of the two is less than the set value, this time is used as the P wave first arrival time; otherwise, the first arrival time is judged to be an interference signal.

6. A robust P-wave seismic phase first arrival picking method according to claim 1, characterized in that: The specific steps of S5 are: The first arrival time of the P wave determined in the above steps is used as a reference. Records of a certain length of time are taken forward and backward from the original acceleration record to calculate the value of the AIC criterion. The processing formula is as follows: AIC(k)=klog{var(x[1,k])}+(Nk-1)log{(var(x[k+1,N])} Among them, k represents the recorded time point, var represents the variance, and the point corresponding to the minimum AIC value is taken as the final P wave first arrival time.

Images