Mine earthquake P wave pickup method, system and equipment based on adaptive characteristic function

Abstract

The invention discloses a mine earthquake P wave pickup method, system and device based on a self-adaptive characteristic function, and relates to the technical field of micro-earthquake signal processing.The method comprises the steps that on the basis of original mine earthquake records, an ICEEMDAN algorithm is adopted, and a mode function component without noise and a mode function component with noise are determined; processing the noise-containing modal function component by adopting a wavelet threshold to obtain a denoised modal function component; obtaining a denoised mine earthquake record based on the noise-free mode function component and the denoised mode function component; obtaining a self-adaptive characteristic function sequence based on the denoised mine earthquake record, and calculating the movement standard deviation of the sequence to obtain a movement standard deviation sequence; and determining an AIC window length based on the mobile standard deviation sequence and a preset pickup threshold value, calculating an AIC value of the adaptive feature function sequence corresponding to the AIC window length, and taking a moment corresponding to the minimum AIC value as the P wave arrival time of the original mine earthquake record. According to the invention, the accuracy and reliability of micro-seismic monitoring are improved.

Description

Technical Field

[0001] The present application relates to the technical field of microseismic signal processing, and in particular to a method, system and device for picking up mine seismic P waves based on an adaptive characteristic function. Background Art

[0002] As a reliable and effective monitoring tool, microseismic monitoring technology has been widely used in many fields such as rock burst, hydraulic fracturing monitoring and mine safety. The accuracy of arrival time recognition of microseismic events directly affects the accuracy of microseismic location results, and high-quality arrival time estimation can significantly improve the efficiency of microseismic monitoring. At present, there are many mature pickers for evaluating the arrival time of high-quality mine earthquake records, including Akaike Information Criterion (AIC) and other improved pickers. In addition, a variety of methods have been developed to suppress noise in noisy mine earthquake records, and a joint picking method has been designed on this basis.

[0003] However, existing technologies still face many challenges and limitations in practical applications. First, digital filters usually need to set specific parameters according to the different main frequencies of microseismic data during use, which not only increases the complexity of operation, but also may lead to instability of filtering effects and suboptimal performance. Secondly, in high-noise environments, existing technologies have limited ability to predict the main frequency of mine earthquake records, resulting in reduced signal processing accuracy. Especially for complex coal mine environments, traditional machine learning methods rely on manually designed features, and these feature extraction processes cannot fully capture key information in the data, which in turn affects the accuracy and practicality of the model. Deep learning models are limited by scarce high-quality training data in coal mine environments, which affects their generalization ability under different conditions and exhibits instability when processing high-noise, non-stationary signals. Although the model has improved in picking accuracy, in some cases, the performance still fails to surpass human experts.

[0004] In summary, existing technologies in the field of microseismic monitoring and phase picking, especially in coal mine environments, face problems such as the accuracy of signal arrival recognition, the complexity of filter parameter setting, the feature design limitations of machine learning methods, the computing resource requirements of deep learning models and their dependence on high-quality data. Therefore, it is urgent to develop a technical solution that can efficiently and stably extract distortion-free mine earthquake records in complex and high-noise coal mine environments, so as to improve the application effect and reliability of microseismic monitoring. Summary of the invention

[0005] The purpose of this application is to provide a method, system and equipment for picking up mine earthquake P waves based on an adaptive characteristic function, thereby improving the accuracy and reliability of microseismic monitoring.

[0006] To achieve the above objectives, this application provides the following solutions:

[0007] In a first aspect, the present application provides a method for picking up mine earthquake P waves based on an adaptive characteristic function, comprising:

[0008] Obtain original mine earthquake records;

[0009] Based on the original mine earthquake records, the ICEEMDAN algorithm is used to determine the modal function components without noise and the modal function components with noise.

[0010] The wavelet threshold is used to perform denoising on the modal function components containing noise, and the denoised modal function components are obtained;

[0011] Based on the modal function components without noise and the modal function components after noise reduction, the denoised mine earthquake records are obtained;

[0012] Based on the denoised mine earthquake records and the adaptive characteristic function, an adaptive characteristic function sequence is obtained; the adaptive characteristic function is constructed based on the amplitude sequence of the denoised mine earthquake records;

[0013] Calculate the moving standard deviation of the adaptive characteristic function sequence to obtain the corresponding moving standard deviation sequence;

[0014] Based on the moving standard deviation sequence and the preset picking threshold, the corresponding AIC window length is determined;

[0015] Within the AIC window length range, the AIC value corresponding to the adaptive characteristic function sequence is calculated, and the time corresponding to the minimum AIC value is used as the P wave arrival time of the original mine earthquake record.

[0016] Optionally, based on the original mine earthquake record, the ICEEMDAN algorithm is used to determine the modal function components without noise and the modal function components with noise, specifically including:

[0017] The ICEEMDAN algorithm is used to decompose the original mine earthquake records to obtain multiple intrinsic mode function components;

[0018] Calculate the sample entropy value of each intrinsic mode function component;

[0019] Determine whether the sample entropy value of each intrinsic mode function component is greater than a preset threshold;

[0020] If so, the corresponding intrinsic mode function component is taken as the noisy mode function component;

[0021] If not, the corresponding intrinsic modal function component is taken as the modal function component without noise.

[0022] Optionally, obtaining a denoised mine earthquake record based on the noise-free modal function component and the denoised modal function component specifically includes:

[0023] Calculate the sample entropy value of each denoised modal function component;

[0024] Determine whether the sample entropy value of each denoised modal function component is greater than a preset threshold;

[0025] If so, the corresponding denoised modal function component is discarded;

[0026] If not, the corresponding denoised modal function components and the noise-free modal function components are reconstructed to obtain the denoised mine earthquake record.

[0027] Optionally, an adaptive feature function is constructed, specifically including:

[0028] Based on the amplitude sequence of the denoised mine earthquake records, the weight parameters corresponding to each sampling point are determined;

[0029] Based on the weight parameters, the relative time series is determined;

[0030] Based on the relative time series and weight parameters, a relative energy coefficient is determined;

[0031] Based on the relative time series and relative energy coefficients, an adaptive characteristic function is constructed.

[0032] Optionally, the expression of the weight parameter W is:

[0033]

[0034] Where W(i) represents the weight parameter of the i-th sampling point; x(j) represents the amplitude value of the j-th sampling point in the amplitude sequence of the mine earthquake record;

[0035] The expression of the relative time series S is:

[0036] S(i)=|x(i) / X max |;

[0037] Where S(i) represents the relative time series of the i-th sampling point; X max It is the maximum absolute amplitude in the amplitude sequence of the mine earthquake record;

[0038] The expression of relative energy coefficient α is:

[0039]

[0040] The expression of adaptive characteristic function ACF is:

[0041] ACF(i)=x α (i)+W·[x(i)-x(i-1)] 2 .

[0042] Optionally, based on the moving standard deviation sequence and a preset picking threshold, the corresponding AIC window length is determined, specifically including:

[0043] Determining the time interval based on the moving standard deviation sequence and a preset picking threshold;

[0044] Based on the time interval, the corresponding AIC window length is determined.

[0045] Optionally, the moving standard deviation is calculated as:

[0046]

[0047] Where M(i) represents the moving standard deviation of the i-th sampling point; N represents the size of the sliding window; ACF(j) represents the corresponding value of the j-th sampling point in the adaptive feature function; u i is the moving average of the i-th sampling point in the sliding window.

[0048] Optionally, the expression of the AIC function is:

[0049] AIC(i)=(i-2)log(var(ACF[1,i]))+(L-2-i)log(var(ACF[1,i]));

[0050] Where AIC(i) represents the AIC value of the i-th sampling point; var(ACF[1,i]) represents the variance of the ACF from the 1st sampling point to the i-th sampling point; and L represents the AIC window length.

[0051] In a second aspect, the present application provides a mine shock P wave picking system based on an adaptive characteristic function, the mine shock P wave picking system based on an adaptive characteristic function is used to implement the mine shock P wave picking method based on an adaptive characteristic function, and the mine shock P wave picking system based on an adaptive characteristic function includes:

[0052] A data acquisition unit, used for acquiring original mine earthquake records;

[0053] A modal function component division unit is used to determine the modal function components without noise and the modal function components with noise based on the original mine earthquake records by using the ICEEMDAN algorithm;

[0054] A denoising unit, used for performing denoising processing on the noisy modal function component by using a wavelet threshold to obtain a denoised modal function component;

[0055] A denoised mining earthquake record determination unit, used to obtain a denoised mining earthquake record based on a modal function component without noise and a modal function component after denoising;

[0056] An adaptive characteristic function sequence determination unit is used to obtain an adaptive characteristic function sequence based on the denoised mine earthquake record and the adaptive characteristic function; the adaptive characteristic function is constructed based on the amplitude sequence of the denoised mine earthquake record;

[0057] A moving standard deviation sequence determination unit is used to calculate the moving standard deviation of the adaptive characteristic function sequence to obtain a corresponding moving standard deviation sequence;

[0058] An AIC window length determination unit, used to determine a corresponding AIC window length based on a moving standard deviation sequence and a preset picking threshold;

[0059] The P wave arrival time determination unit is used to calculate the AIC value corresponding to the adaptive characteristic function sequence within the AIC window length range, and use the time corresponding to the minimum AIC value as the P wave arrival time of the original mine earthquake record.

[0060] In a third aspect, the present application provides a computer device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement any of the above-mentioned methods for picking up mine earthquake P waves based on the adaptive characteristic function.

[0061] According to the specific embodiments provided in this application, this application has the following technical effects:

[0062] The present application discloses a method, system and equipment for picking up P-waves of mine earthquakes based on an adaptive characteristic function. The method is based on original mine earthquake records and adopts the ICEEMDAN-SamEn-WTD method to denoise the original mine earthquake records, thereby significantly improving the accuracy of subsequent P-wave picking. Then, the MSD-ACF-AIC method is adopted to pick up the P-wave arrival time of the denoised mine earthquake records, thereby significantly improving the accuracy of P-wave arrival time picking and laying a good foundation for subsequent microseismic event monitoring. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0064] Figure 1 A schematic diagram of a process flow of a mine earthquake P wave picking method based on an adaptive characteristic function provided in an embodiment of the present application;

[0065] Figure 2A schematic diagram of IMF components obtained by decomposing original mine earthquake records using the ICEEMDAN method provided in an embodiment of the present application;

[0066] Figure 3 A schematic diagram showing the comparison between the original mine earthquake record and the mine earthquake record after denoising provided in one embodiment of the present application;

[0067] Figure 4 A schematic diagram comparing the results of picking the P-wave arrival time of synthetic mine earthquake records by the present application and the M-AIC picking method provided in one embodiment of the present application;

[0068] Figure 5 A schematic diagram showing a comparison between an original mine earthquake record and a mine earthquake record after denoising provided in another embodiment of the present application;

[0069] Figure 6 A schematic diagram of the distribution of P-wave picking deviations provided in another embodiment of the present application;

[0070] Figure 7 A schematic diagram of functional modules of a mine earthquake P-wave picking system based on an adaptive characteristic function provided in one embodiment of the present application;

[0071] Figure 8 A schematic diagram of the structure of a computer device provided in one embodiment of the present application.

[0072] Reference numerals:

[0073] Data acquisition unit-1, modal function component division unit-2, noise reduction unit-3, denoised mine earthquake record determination unit-4, adaptive characteristic function sequence determination unit-5, moving standard deviation sequence determination unit-6, AIC window length determination unit-7, P wave arrival time determination unit-8. DETAILED DESCRIPTION

[0074] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0075] The present application provides a method, system and device for picking up P waves of mine earthquakes based on adaptive characteristic functions to improve the accuracy of picking up P waves in high-noise mine earthquake records, which is mainly divided into two aspects: denoising and P-wave signal extraction. Among them, in terms of denoising, a method combining improved complete adaptive noise ensemble empirical mode decomposition (Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise, ICEEMDAN), sample entropy (SamEn) and wavelet threshold denoising (Wavelet Threshold Denoising, WTD) is used to extract undistorted useful signals from high-noise microseismic data, significantly improve the accuracy of P wave picking, and ensure that key mine earthquake records can still be effectively extracted in a high-noise environment. In terms of P-wave signal extraction, relative time series and relative energy coefficients are introduced to construct an improved adaptive characteristic function (ACF). In order to further improve the picking accuracy, the moving standard deviation (MSD) is creatively introduced to calculate the improved adaptive characteristic function. The window length of the autoregressive Akaike Information Criterion (Akaike Information Criterion, AIC) is determined by the picking threshold, and then the AIC within the window length is calculated, and then the minimum value is determined as the initial P-wave arrival time of the original mine earthquake record, which further improves the accuracy of P-wave picking.

[0076] In order to make the above-mentioned objects, features and advantages of the present application more obvious and easy to understand, the present application is further described in detail below in conjunction with the accompanying drawings and specific implementation methods.

[0077] In an exemplary embodiment, Figure 1 As shown, a method for picking up mining earthquake P waves based on an adaptive characteristic function is provided, comprising the following steps S1 to S6. Among them:

[0078] Step S1, obtaining original mine earthquake records.

[0079] Step S2, based on the original mine earthquake record, using the ICEEMDAN algorithm to determine the modal function components without noise and the modal function components with noise.

[0080] Step S3, using a wavelet threshold to perform denoising on the noisy modal function component to obtain a denoised modal function component.

[0081] Step S4, obtaining a denoised mine earthquake record based on the noise-free modal function component and the denoised modal function component.

[0082] Step S5, obtaining an adaptive characteristic function sequence based on the denoised mine earthquake records and the adaptive characteristic function; the adaptive characteristic function is constructed based on the amplitude sequence of the denoised mine earthquake records.

[0083] Step S6, calculating the moving standard deviation of the adaptive characteristic function sequence to obtain the corresponding moving standard deviation sequence.

[0084] Step S7, determining the corresponding AIC window length based on the moving standard deviation sequence and the preset picking threshold.

[0085] Step S8, within the AIC window length range, calculate the AIC value corresponding to the adaptive characteristic function sequence, and use the time corresponding to the minimum AIC value as the P wave arrival time of the original mine earthquake record.

[0086] Specifically, in microseismic monitoring, determining the arrival time of the P wave is a key step in identifying the location of the earthquake source and determining the parameters of the earthquake source. By accurately identifying the arrival time of the P wave, the spatial location of the microseism, the intensity of the earthquake source and its activity pattern can be further calculated. Using the arrival time of the P wave and the propagation speed of seismic waves, the location of the earthquake source can be calculated by the differential method or other algorithms; the earthquake source position calculated by the arrival time of the P wave, combined with other earthquake source information such as the arrival time and amplitude of the S wave, can be further calculated. Parameters such as the magnitude and stress release of the earthquake source, and then these data are used to analyze the activity characteristics and potential risks of the earthquake source; by continuously monitoring the arrival time of the P wave, real-time monitoring of microseisms can be achieved. Once microseismic activity is identified and its location and intensity are determined, the monitoring system can issue an early warning in time to avoid possible accidents; the collected P wave arrival data can be accumulated to conduct in-depth analysis of the stress state and fault activity of the mining area or engineering area, identify potential risk areas, and provide reference for subsequent operations.

[0087] As an optional implementation, step S2 specifically includes:

[0088] Step S21, using the ICEEMDAN algorithm to decompose the original mine earthquake record to obtain multiple intrinsic mode function (IMF) components, and at the same time, multiple residual components.

[0089] Step S22, calculate the sample entropy value of each intrinsic mode function component. After the IMF components are obtained by ICEEMDAN decomposition, it is necessary to determine which IMF components are of substantial help in reconstructing new mine earthquake records. Specifically, the smaller the sample entropy, the higher the self-similarity of the signal and the lower the noise; on the contrary, the larger the sample entropy, the greater the noise.

[0090] Step S23, determining whether the sample entropy value of each intrinsic mode function component is greater than a preset threshold; if so, executing step S24, if not, executing step S25.

[0091] Step S24: taking the corresponding intrinsic modal function component as the modal function component containing noise.

[0092] Step S25: taking the corresponding intrinsic modal function component as the modal function component without noise.

[0093] Specifically, the original mine earthquake records were decomposed using the ICEEMDAN algorithm. ICEEMDAN is an advanced signal processing method that improves the traditional empirical mode decomposition (EMD) technique and its variants (such as EEMD and CEEMDAN), aiming to more accurately decompose complex nonlinear and non-stationary signals. Compared with the modal aliasing problem in EEMD, ICEEMDAN makes the separation of each intrinsic mode function (IMF) clearer by adaptively adding noise, significantly improving the accuracy and stability of signal processing. Therefore, ICEEMDAN can provide more accurate and robust results when processing complex signals, effectively reduce modal aliasing, and make signal decomposition more accurate.

[0094] As an optional implementation, step S21 specifically includes:

[0095] Step S211, set the noise amplitude parameter ε0, generate a noise sequence n(i), and add Gaussian white noise to the original mine earthquake record to obtain X(i).

[0096] X(i)=X+ε0E1(n(i)),i=1,2,...,G (1)

[0097] Where X is the original mine earthquake record, E1(·) is the function for calculating IMF1, and i is the sampling point number.

[0098] Step S212, obtaining the first residual component R1:

[0099]

[0100] Where A(·) is the local average of the original mine earthquake records generated by EMD, and G is the number of time series.

[0101] Step S213, define k as the kth intrinsic mode function component generated by ICEEMDAN decomposition. When k=1, the first intrinsic mode function component F1 is obtained:

[0102]

[0103] When k=2,3,4,...,K, calculate the kth residual component R k and the kth intrinsic mode function component F k :

[0104]

[0105] F k =R k-1 -R k (5)

[0106] Among them, E k (·) is the function that generates the EMD decomposition model.

[0107] Step S214, repeatedly calculate to obtain all IMF components until the final residual component R n It satisfies monotonicity and cannot be decomposed any further. The original mine earthquake record is decomposed into two parts:

[0108]

[0109] Furthermore, in order to verify the effectiveness of ICEEMDAN decomposition, a mine earthquake record with a sampling rate of 500 and a sample number of 2000 was synthesized using Ricker wavelet, and then 70% random noise was randomly added to it. The ICEEMDAN method was used to decompose it into 9 IMF components. The decomposed IMF components correspond to the different frequency components of the original mine earthquake record in descending order. The decomposed IMF components and residual components are shown in Figure 2 As shown, from Figure 2 It can be seen that ICEEMDAN decomposition can help extract useful signals from complex mine earthquake records.

[0110] As an optional implementation, step S22 specifically includes:

[0111] Step S221, assuming that C f (i) is the IMF sequence obtained by ICEEMDAN decomposition, where i = 1, 2, 3, ..., G, f represents the total number of IMFs, and C f (i) Construct an m-dimensional block sequence Y(i):

[0112] Y(i)=|C f (i),C f (i+1),...,C f (i+m-1)| (7)

[0113] Where i=1,2,...,k, and k=G-m+1.

[0114] Step S222, calculate the distance between the sequence Y(i) and all the sequences Y(j), which is defined as:

[0115] d ij =max|C f (i+l)-C f (j+l)| (8)

[0116] Among them, i≠j, l=0,1,2,...,m-1.

[0117] Step S223, define T = r × SD, where r ranges from 0.1 to 0.25, SD represents the standard deviation of the calculated sequence, and then calculate d ij The number greater than T is denoted by

[0118] Step S224, define Represents the mean of the i-th m-dimensional block sequence, which is used to calculate the average similarity of the sample sequence in the current dimension.

[0119]

[0120] Step S225, all The average value of is defined as φ m , which represents the overall average similarity under the current dimension.

[0121]

[0122] Step S226, let m=m+1, and then repeat steps S221 to S225, that is, during the calculation process, all relevant distances and average values ​​will be calculated again.

[0123] Step S227, define sample entropy SamEn as:

[0124] SamEn=lnφ m -lnφ m+1 (11)

[0125] Among them, φ m and φ m+1 It is the average value calculated in m dimension and m+1 dimension.

[0126] Among them, m and T are two important parameters. The parameter m is set to 1 and T is set to 0.2 times the standard deviation of the IMF to distinguish the different characteristics of the decomposed IMF. Subsequently, the SampEn values ​​of different IMFs are calculated to illustrate their respective complexities. In order to extract useful IMF components, it is necessary to determine a reasonable threshold to eliminate redundant IMF components. The noise characteristics of the signal are judged by randomness, and the threshold is set to 0.2.

[0127] As an optional implementation, step S4 specifically includes:

[0128] Step S41, calculating the sample entropy value of each denoised modal function component.

[0129] Step S42, determining whether the sample entropy value of each denoised modal function component is greater than a preset threshold; if so, executing step S43; if not, executing step S44.

[0130] Step S43: discard the corresponding denoised modal function components.

[0131] Step S44, reconstructing the corresponding denoised modal function components and the noise-free modal function components to obtain the denoised mine earthquake record.

[0132] Specifically, since the original mine earthquake records contain some medium and high frequency components, the original ICEEMDAN algorithm and sample entropy denoising mode will directly discard significant IMF components, which will reduce noise while also causing loss of useful signals. To this end, this application has made improvements to this, performing wavelet threshold denoising on the IMF components that are higher than the preset threshold, and calculating the sample entropy value of the denoised IMF components again. If the sample entropy value is higher than the preset threshold, it will be discarded; otherwise, it will be retained. Figure 2 The sample entropy values ​​of IMF components 1 to 4 are all greater than 0.2. Wavelet threshold denoising is performed on IMF components 1 to 4, and the sample entropy values ​​of the denoised IMF components are less than 0.2. Finally, the noise-free IMF components and the denoised IMF components are reconstructed to obtain effective mine earthquake records, as shown in the schematic diagram. Figure 3 As shown, Figure 3 The blue curve in the middle is the original mine earthquake record, and the red curve is the mine earthquake record after denoising. The above method is named ICEEMDAN-SamEn-WTD, and the mine earthquake record denoised by this method is used as the input of the subsequent P-wave picking method.

[0133] As an optional implementation, in step S5, constructing an adaptive feature function specifically includes:

[0134] Step S511, based on the amplitude sequence of the denoised mine earthquake record, determine the weight parameter corresponding to each sampling point.

[0135] Step S512, determining the relative time series based on the weight parameter.

[0136] Step S513, determining the relative energy coefficient based on the relative time series and the weight parameter.

[0137] Step S514, constructing an adaptive characteristic function based on the relative time series and the relative energy coefficient.

[0138] Specifically, the arrival time of the mine earthquake P wave can be indicated by amplifying and analyzing the changes in amplitude and frequency in the earthquake sequence. This is the main task of the characteristic function. By establishing a suitable and sensitive characteristic function, the resolution between the signal and the noise will be further improved, thereby improving the accuracy of the arrival time. In this application, relative time series and relative energy coefficients are introduced to construct an adaptive characteristic function.

[0139] As an optional implementation, in step S511, the expression of the weight parameter W is:

[0140]

[0141] Among them, W(i) represents the weight parameter of the i-th sampling point, and the parameter W is a weight parameter that changes with the sampling rate and the characteristics of stationary noise; x(j) represents the amplitude value of the j-th sampling point in the amplitude sequence of the mine earthquake record.

[0142] The expression of the relative time series S is:

[0143] S(i)=|x(i) / X max | (13)

[0144] Where S(i) represents the relative time series of the i-th sampling point; X max It is the maximum absolute amplitude in the amplitude sequence of the mine earthquake record and is used to adjust the time series in the mine earthquake record.

[0145] The expression of relative energy coefficient α is:

[0146]

[0147] The expression of adaptive characteristic function ACF is:

[0148] ACF(i)=x α (i)+W·[x(i)-x(i-1)] 2 (15)

[0149] Among them, the adaptive characteristic function ACF is modified from Allen's characteristic function.

[0150] As an optional implementation, step S5 specifically includes:

[0151] The waveform record obtained after denoising the above-mentioned noisy synthetic mine earthquake record is used as the input for P wave arrival time picking. The ACF calculation is performed on it using formula (15) to obtain the adaptive characteristic function sequence, as shown in Figure 4 As shown in (c), Figure 4 (a) is the original mine earthquake record. Figure 4 (b) in the figure is the denoised mine earthquake record. Figure 4 It can be seen from (a) to (c) that the obtained adaptive characteristic function sequence further suppresses the noise in the mine earthquake signal.

[0152] As an optional implementation, step S6 specifically includes:

[0153] Calculate the moving standard deviation of each sampling point in the adaptive characteristic function sequence. In order to further improve the accuracy of picking up the P wave arrival time, the MSD and ACF-AIC methods are combined. Use MSD for ACF, and manually set the sliding window and picking threshold E0. In the application, the sliding window is set to 0.25s, and the picking threshold M0 is set to 0.5. The calculation formula for the moving standard deviation of each sampling point is:

[0154]

[0155] Where M(i) represents the moving standard deviation of the i-th sampling point; N represents the size of the sliding window; ACF(j) represents the corresponding value of the j-th sampling point in the adaptive feature function; u i is the moving average of the i-th sampling point in the sliding window.

[0156] As an optional implementation, step S7 specifically includes:

[0157] Step S71, based on the moving standard deviation sequence and the preset picking threshold, determine the time interval. A basic time interval L0 is determined by M(i and M0, such as Figure 4 As shown in (d), the length between the two black dotted lines is the basic time interval L0, and the basic time interval L0 is determined by the time t1 and t2 corresponding to the two points where the red dotted line intersects the picking threshold M0 and M(i). In addition, the maximum value of M(i) is also included in this interval.

[0158] Step S72: determining the corresponding AIC window length based on the time interval.

[0159] The AIC window length is determined based on t1 determined by MSD and the basic time interval L0. The AIC window length should be large enough to include the prominent waveform. In this application, the AIC window L is controlled by N*L0. In this application, N is set to 2, and the AIC window is Figure 4 The distance L between the two black dashed lines shown in (c).

[0160] The AIC function is corrected by using ACF to obtain the corrected AIC function, and the corrected AIC function is used as the final AIC function, whose expression is:

[0161] AIC(i)=(i-2)log(var(ACF[1,i]))+(L-2-i)log(var(ACF[1,i])) (18)

[0162] Where AIC(i) represents the AIC value of the i-th sampling point; var(ACF[1,i]) represents the variance of the ACF from the 1st sampling point to the i-th sampling point, i = 1, 2, .., L; L represents the AIC window length.

[0163] Further, in step S8, based on the above AIC function, within the AIC window length range, the AIC value corresponding to the adaptive feature function sequence is calculated, such as Figure 4 As shown by the red dotted line in (e), the time corresponding to the minimum AIC value is taken as the P wave arrival time of the original mine earthquake record.

[0164] The above-mentioned method of automatic P wave picking is named MSD-ACF-AIC. Figure 4 The comparison chart shows the results of ICEEMDAN-SamEn-WTD denoising the entire original mine earthquake record and applying the MSD-ACF-AIC method to the reconstructed denoised mine earthquake record and executing the conventional M-AIC picker. Figure 4 In (f), Ori represents the synthetic signal of mining earthquake containing noise, and the blue solid line represents the result of M-AIC pickup for the noisy signal; Figure 4 In (g), Pre represents the denoised mine earthquake signal obtained by executing the conventional M-AIC picker, and the green solid line represents the picking result of the M-AIC picker on the denoised signal. By comparison, it can be found that the denoised mine earthquake record obtained after reconstruction can reliably replace the original mine earthquake signal for P wave first arrival picking. Especially in Figure 4 (f) and Figure 4 As can be seen from (g), the denoising strategy of ICEEMDAN-SamEn-WTD of the present application significantly improves the accuracy of P-wave arrival time picking in a high noise environment. In addition, Figure 4 The MSD-ACF-AIC method shown in (b) is more accurate than the traditional M-AIC method, indicating that the MSD-ACF-AIC method has obvious advantages in seismic phase picking, especially in P-wave first arrival picking. Therefore, the MSD-ACF-AIC method can provide higher accuracy and robustness in complex noise environments, and the ICEEMDAN-SamEn-WTD denoising strategy effectively improves the signal quality, so that the denoised signal can be better used for subsequent seismic signal analysis.

[0165] Furthermore, the mine earthquake P-wave picking method based on the adaptive characteristic function of the present application is applied to the manually picked mine earthquake data set of Dongtan Coal Mine (a total of 343 mine earthquake records) to further verify the technical solution of the present application.

[0166] First, the ICEEMDAN-SamEn-WTD method of this application is used to denoise the manually picked mine earthquake dataset of Dongtan Coal Mine to obtain the denoised mine earthquake records. Figure 5 It shows the mine earthquake records before and after denoising. Figure 5 (a) in the figure is the mine earthquake record before denoising (i.e. the original mine earthquake record). Figure 5 (b) in the figure is the denoised mine earthquake record. It can be seen that compared with the original mine earthquake record, the noise of the denoised mine earthquake record obtained by denoising using the ICEEMDAN-SamEn-WTD method is significantly reduced, and the signal-to-noise ratio is increased from 5.9dB of the original mine earthquake record to 10.8dB, indicating that the denoising algorithm can effectively enhance the signal quality, thereby improving the accuracy of picking up the P wave arrival time.

[0167] Furthermore, the MSD-ACF-AIC method is used to pick the P-wave arrival time of the denoised mine earthquake records, and the autocorrelation coefficient sequence obtained after ACF calculation is as follows: Figure 5 As shown in (c) in the figure, the schematic diagram of the P wave arrival time manually calibrated by professionals is as follows Figure 5 As shown in (d) in the figure, the P wave arrival time automatically calibrated using the MSD-ACF-AIC method is as follows Figure 5 As shown in (e) in the figure, it can be seen that the P wave arrival time (i.e. Figure 5 The red dashed line in (e) is different from the P wave arrival time manually calibrated by professionals (i.e. Figure 5 The error between the green dotted line in (d) is zero, indicating that the accuracy of the automated method in this embodiment has reached a level comparable to that of manual calibration.

[0168] In addition, the test results show that in 94.7% of the data (i.e., 325 records), the error between the automatically picked P wave arrival time and the manually calibrated time is within the range of 0 to 3 sample intervals. Specifically, these errors range from + / -0ms to + / -6ms. Figure 6 shown. Figure 6 The distribution of the P-wave picking deviation is further shown. Most of the errors are concentrated around 0ms, and the distribution of the deviation shows a typical normal distribution. This result shows that this method has extremely high accuracy and stability in automatically picking the P-wave arrival time.

[0169] In general, the experimental results verify that the mine earthquake P-wave picking method based on adaptive characteristic functions in this application not only performs well in noise suppression, but also can accurately estimate the arrival time of noisy mine earthquake records, greatly improving the efficiency and accuracy of mine earthquake monitoring, providing a solid technical foundation for the automated application of mine earthquake data processing, and also demonstrating its broad application prospects in coal mine safety monitoring under complex geological conditions.

[0170] Based on the same inventive concept, the embodiment of the present application also provides a mine shock P wave picking system based on an adaptive characteristic function for implementing the mine shock P wave picking method based on an adaptive characteristic function. The implementation scheme for solving the problem provided by the system is similar to the implementation scheme recorded in the above method, so the specific limitations in one or more embodiments of the mine shock P wave picking system based on an adaptive characteristic function provided below can refer to the limitations of the mine shock P wave picking method based on an adaptive characteristic function in the above text, and will not be repeated here.

[0171] In an exemplary embodiment, Figure 7 As shown, a mine earthquake P wave picking system based on an adaptive characteristic function is provided, comprising:

[0172] The data acquisition unit 1 is used to acquire original mine earthquake records.

[0173] The modal function component division unit 2 is used to determine the modal function components without noise and the modal function components with noise based on the original mine earthquake records by using the ICEEMDAN algorithm.

[0174] The noise reduction unit 3 is used to perform noise reduction processing on the noisy modal function component by using a wavelet threshold to obtain the noise-reduced modal function component.

[0175] The denoised mining earthquake record determination unit 4 is used to obtain the denoised mining earthquake record based on the modal function component without noise and the modal function component after denoising.

[0176] The adaptive characteristic function sequence determination unit 5 is used to obtain an adaptive characteristic function sequence based on the denoised mine earthquake record and the adaptive characteristic function; the adaptive characteristic function is constructed based on the amplitude sequence of the denoised mine earthquake record.

[0177] The moving standard deviation sequence determining unit 6 is used to calculate the moving standard deviation of the adaptive characteristic function sequence to obtain the corresponding moving standard deviation sequence.

[0178] The AIC window length determining unit 7 is used to determine the corresponding AIC window length based on the moving standard deviation sequence and the preset picking threshold.

[0179] The P wave arrival time determination unit 8 is used to calculate the AIC value corresponding to the adaptive characteristic function sequence within the AIC window length range, and take the time corresponding to the minimum AIC value as the P wave arrival time of the original mine earthquake record.

[0180] In an exemplary embodiment, a computer device is provided, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement a mine earthquake P-wave picking method based on an adaptive characteristic function.

[0181] In an exemplary embodiment, a computer device is provided. The computer device may be a server or a terminal. The internal structure diagram thereof may be as follows: Figure 8 As shown. The computer device includes a processor, a memory, an input / output interface (Input / Output, referred to as I / O) and a communication interface. The processor, the memory and the input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal through a network connection. When the computer program is executed by the processor, a method for picking up mine earthquake P waves based on an adaptive characteristic function is implemented.

[0182] Those skilled in the art will understand that Figure 8 The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0183] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, stored data, displayed data, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of relevant data must comply with relevant regulations.

[0184] Those of ordinary skill in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to the memory, database or other medium used in the embodiments provided in the present application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM may be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0185] The database involved in each embodiment provided in this application may include at least one of a relational database and a non-relational database. The non-relational database may include a distributed database based on blockchain, etc., but is not limited thereto. The processor involved in each embodiment provided in this application may be a general-purpose processor, a central processing unit, a graphics processor, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, etc., but is not limited thereto.

[0186] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0187] This article uses specific examples to illustrate the principles and implementation methods of this application. The description of the above embodiments is only used to help understand the method and core ideas of this application. At the same time, for those skilled in the art, according to the ideas of this application, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting this application.

Claims

1. A method for picking up P waves of mine earthquakes based on an adaptive characteristic function, characterized in that: The method for picking up mining earthquake P waves based on the adaptive characteristic function includes: Obtain original mine earthquake records; Based on the original mine earthquake records, the ICEEMDAN algorithm is used to determine the modal function components without noise and the modal function components with noise. The wavelet threshold is used to perform denoising on the modal function components containing noise, and the denoised modal function components are obtained; Based on the modal function components without noise and the modal function components after noise reduction, the denoised mine earthquake records are obtained; Based on the denoised mine earthquake records and the adaptive characteristic function, an adaptive characteristic function sequence is obtained; the adaptive characteristic function is constructed based on the amplitude sequence of the denoised mine earthquake records; Calculate the moving standard deviation of the adaptive characteristic function sequence to obtain the corresponding moving standard deviation sequence; Based on the moving standard deviation sequence and the preset picking threshold, the corresponding AIC window length is determined; Within the AIC window length range, the AIC value corresponding to the adaptive characteristic function sequence is calculated, and the time corresponding to the minimum AIC value is used as the P wave arrival time of the original mine earthquake record.

2. The method for picking up mining earthquake P waves based on adaptive characteristic function according to claim 1, characterized in that: Based on the original mine earthquake records, the ICEEMDAN algorithm is used to determine the modal function components without noise and the modal function components with noise, including: The ICEEMDAN algorithm is used to decompose the original mine earthquake records to obtain multiple intrinsic mode function components; Calculate the sample entropy value of each intrinsic mode function component; Determine whether the sample entropy value of each intrinsic mode function component is greater than a preset threshold; If so, the corresponding intrinsic mode function component is taken as the noisy mode function component; If not, the corresponding intrinsic modal function component is taken as the modal function component without noise.

3. The method for picking up mining earthquake P waves based on adaptive characteristic function according to claim 1, characterized in that: Based on the modal function components without noise and the modal function components after noise reduction, the denoised mine earthquake records are obtained, which specifically include: Calculate the sample entropy value of each denoised modal function component; Determine whether the sample entropy value of each denoised modal function component is greater than a preset threshold; If so, the corresponding denoised modal function component is discarded; If not, the corresponding denoised modal function components and the noise-free modal function components are reconstructed to obtain the denoised mine earthquake record.

4. The method for picking up mining earthquake P waves based on adaptive characteristic function according to claim 1, characterized in that: Construct an adaptive feature function, including: Based on the amplitude sequence of the denoised mine earthquake records, the weight parameters corresponding to each sampling point are determined; Based on the weight parameters, the relative time series is determined; Based on the relative time series and weight parameters, a relative energy coefficient is determined; Based on the relative time series and relative energy coefficients, an adaptive characteristic function is constructed.

5. The method for picking up mining earthquake P waves based on adaptive characteristic function according to claim 4, characterized in that: The expression of the weight parameter W is: Where W(i) represents the weight parameter of the i-th sampling point; x(j) represents the amplitude value of the j-th sampling point in the amplitude sequence of the mine earthquake record; The expression of the relative time series S is: S(i)=|x(i) / X max |; Among them, S(i) represents the relative time series of the i-th sampling point, X max It is the maximum absolute amplitude in the amplitude sequence of the mine earthquake record; The expression of relative energy coefficient α is: The expression of adaptive characteristic function ACF is: ACF(i)=x α (i)+W·[x(i)-x(i-1)] 2 。 6. The method for picking up mining earthquake P waves based on adaptive characteristic function according to claim 1, characterized in that: Based on the moving standard deviation sequence and the preset picking threshold, the corresponding AIC window length is determined, including: Determining the time interval based on the moving standard deviation sequence and a preset picking threshold; Based on the time interval, the corresponding AIC window length is determined.

7. The method for picking up mining earthquake P waves based on adaptive characteristic function according to claim 6, characterized in that: The formula for calculating the moving standard deviation is: Where M(i) represents the moving standard deviation of the i-th sampling point; N represents the size of the sliding window; ACF(j) represents the corresponding value of the j-th sampling point in the adaptive feature function; u i is the moving average of the i-th sampling point in the sliding window.

8. The method for picking up mining earthquake P waves based on adaptive characteristic function according to claim 6, characterized in that: The expression of the AIC function is: AIC(i)=(i-2)log(var(ACF[1,i]))+(L-2-i)log(var(ACF[1,i])); Where AIC(i) represents the AIC value of the i-th sampling point; var(ACF[1, i]) represents the variance of the ACF from the 1st sampling point to the i-th sampling point; and L represents the AIC window length.

9. A mine earthquake P wave picking system based on adaptive characteristic function, characterized in that: The mine shock P wave picking system based on the adaptive characteristic function is used to implement the mine shock P wave picking method based on the adaptive characteristic function according to any one of claims 1 to 8, and the mine shock P wave picking system based on the adaptive characteristic function comprises: A data acquisition unit, used for acquiring original mine earthquake records; A modal function component division unit is used to determine the modal function components without noise and the modal function components with noise based on the original mine earthquake records by using the ICEEMDAN algorithm; A denoising unit, used for performing denoising processing on the noisy modal function component by using a wavelet threshold to obtain a denoised modal function component; A denoised mining earthquake record determination unit, used to obtain a denoised mining earthquake record based on a modal function component without noise and a modal function component after denoising; An adaptive characteristic function sequence determination unit is used to obtain an adaptive characteristic function sequence based on the denoised mine earthquake record and the adaptive characteristic function; the adaptive characteristic function is constructed based on the amplitude sequence of the denoised mine earthquake record; A moving standard deviation sequence determination unit is used to calculate the moving standard deviation of the adaptive characteristic function sequence to obtain a corresponding moving standard deviation sequence; An AIC window length determination unit, used to determine a corresponding AIC window length based on a moving standard deviation sequence and a preset picking threshold; The P wave arrival time determination unit is used to calculate the AIC value corresponding to the adaptive characteristic function sequence within the AIC window length range, and use the time corresponding to the minimum AIC value as the P wave arrival time of the original mine earthquake record.

10. A computer device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the mine earthquake P-wave picking method based on the adaptive characteristic function as described in any one of claims 1 to 8.

Images