Microseismic signal classification method based on VMD feature extraction and HOA optimized random forest

Abstract

The present application relates to a kind of microseismic signal classification method based on VMD feature extraction and HOA optimization random forest.The present application method is to solve the problems of feature extraction failure in strong noise environment, poor self-adaptive ability of classification model parameters and low recognition accuracy of existing microseismic signal recognition method, by linear detrending and amplitude normalization preprocessing of the obtained original microseismic signal, eliminating direct current drift and non-stationary trend;Then extract the multi-dimensional physical feature vector of fusion time energy duration, frequency domain robust main frequency and VMD-based main seismic phase maximum energy time, to overcome the limitation of single feature;On the basis of the foregoing, the hyperparameters of random forest classifier are adaptively optimized using HOA algorithm, and the optimal classification model is constructed.The present application can effectively distinguish artificial blasting, mine earthquake and natural earthquake signal, and still maintain high classification accuracy and robustness under low signal-to-noise ratio conditions, suitable for real-time microseismic monitoring in complex mine environment.

Description

Technical Field

[0001] This invention belongs to the field of mine safety monitoring and geophysical signal processing technology, specifically involving a microseismic signal classification method based on variational mode decomposition (VMD) physical feature extraction and hiking optimization algorithm (HOA) optimized random forest. Background Technology

[0002] During deep mining in coal and non-coal mines, the stress distribution in the rock mass undergoes complex changes with increasing depth and intensity, easily inducing dynamic disasters such as fault slippage and roof fracture. Microseismic monitoring technology, acting as a "stethoscope" for sensing changes in rock mass stability, analyzes vibration signals to identify events such as mine tremors, artificial blasting, and natural earthquakes caused by mining activities, playing a crucial role in ensuring safe mine production.

[0003] Current microseismic signal classification methods mainly rely on manual experience or traditional machine learning algorithms. However, existing technologies have the following significant drawbacks in practical applications:

[0004] (1) Feature extraction failure under strong noise environment: The background noise such as electromagnetic interference and mechanical vibration in the mine is extremely strong, resulting in low signal-to-noise ratio of micro-vibration signal. Traditional time-frequency analysis methods (such as short-time Fourier transform STFT) are limited by the uncertainty principle and it is difficult to obtain high-precision time-frequency resolution at the same time; and existing methods often ignore DC drift and non-stationary trend in the signal, which causes the extracted frequency features (such as the main frequency) to be often masked by false low-frequency components near 0Hz, which seriously affects the classification accuracy.

[0005] (2) Lack of physical interpretability of features: Although deep learning methods (such as CNN and RNN) have improved the recognition rate to some extent, their "black box" characteristics make the model decision-making process lack geological and physical interpretation. For example, the energy lag characteristics of natural earthquakes (S waves arrive later than P waves) and the instantaneous impact characteristics of artificial blasting are difficult to be explicitly captured and utilized in traditional neural networks.

[0006] (3) Difficulty in optimizing classification model parameters: Random Forest (RF) is a commonly used classifier, and its performance is highly dependent on the setting of hyperparameters such as the number of decision trees and the maximum depth. Existing methods mostly use manual trial and error or grid search to adjust parameters, which is not only computationally inefficient, but also prone to getting trapped in local optima and difficult to adapt to the complex and ever-changing data distribution characteristics of different mines.

[0007] Therefore, there is an urgent need to develop a microseismic signal classification method that can integrate a clear physical mechanism, has strong noise resistance, and can adaptively optimize model parameters. Summary of the Invention

[0008] This invention aims to address the problems of existing microseismic signal identification methods, such as feature extraction failure in high-noise environments, lack of physical interpretability, and poor adaptive capability of classification model parameters. It provides a microseismic signal classification method based on variational mode decomposition (VMD) feature extraction and hiking optimization (HOA) algorithm-optimized random forest. Furthermore, this invention also relates to a microseismic signal classification system and a computer-readable storage medium for implementing the above classification method.

[0009] To achieve the above objectives, the present invention employs the following technical solutions:

[0010] In a first aspect, the present invention provides a microseismic signal classification method based on VMD physical feature extraction and HOA-optimized random forest, the method comprising the following steps:

[0011] Step 1: Deploy a microseismic sensor array in the mine monitoring area to acquire the raw microseismic monitoring signal data to be identified;

[0012] Step 2: Perform physical constraint preprocessing on the original microseismic monitoring signal data. The preprocessing includes linear detrending processing and amplitude normalization processing to eliminate DC drift and non-stationary trends, and obtain a pure microseismic signal.

[0013] Step 3: Perform multi-domain robust physical feature extraction on the pure microseismic signal to construct a multi-dimensional feature vector; the multi-dimensional feature vector includes at least time-domain energy duration features, frequency-domain robust dominant frequency features, and VMD-based maximum energy moment features of the main seismic phase;

[0014] Step 4: Construct a microseismic classifier based on random forest, and use the HOA algorithm to adaptively optimize the hyperparameter combination of the random forest to obtain the optimal microseismic classification model;

[0015] Step 5: Input the multidimensional feature vector extracted in Step 3 into the optimal microseismic classification model trained in Step 4, and output the category label of the microseismic signal. The category label includes artificial blasting, mining tremors, and natural earthquakes.

[0016] Furthermore, in step 2, the linear detrending process involves fitting the linear trend term using the least squares method and performing detrending using the following formula:

[0017]

[0018] Where, x detrend (t) represents the detrended signal data; x raw (t) represents the original microseismic monitoring signal data; a represents the slope of the linear trend term; b represents the intercept of the linear trend term; t represents the time sampling point.

[0019] Furthermore, the energy integration method is used to extract the time-domain energy duration feature in step 3. First, the cumulative energy curve E(i) of the pure microseismic signal is calculated, as follows:

[0020]

[0021] Where E(i) represents the cumulative energy at the i-th sampling point; x(k) represents the amplitude of the pure microseismic signal at the k-th sampling point;

[0022] Then, the time span during which the accumulated energy increases from 5% to 95% of the total energy is calculated, i.e., the time-domain energy duration characteristic, calculated using the following formula:

[0023]

[0024] Among them, T duration Represented as the time-domain energy duration characteristic; t 95% This represents the index of the sampling point when the accumulated energy reaches 95% of the total energy; t 5% The index represents the sampling point index corresponding to when the accumulated energy reaches 5% of the total energy; ∆t represents the sampling interval of the signal. This method can accurately define the duration of the effective signal in a noisy environment.

[0025] Furthermore, step 3, which involves extracting the maximum energy moment feature of the mainshock phase based on VMD, specifically includes:

[0026] A constrained variational model is constructed to decompose the pure microseismic signal into K intrinsic modal components. The model formula is as follows:

[0027]

[0028] Among them, {u k} represents the k-th intrinsic mode component obtained from the decomposition; {ω k} represents the center frequency of each modal component; δ(t) represents the Dirac distribution function; * represents the convolution operation; ∂ t This represents the partial derivative with respect to time; j represents the imaginary unit.

[0029] After solving, the mode with the highest energy is selected as the dominant seismic phase mode u. best The envelope A(t) is extracted using the Hilbert transform:

[0030]

[0031] Where A(t) represents the instantaneous amplitude envelope; u best (t) represents the main seismic phase mode signal; best (t) is represented as u bestThe Hilbert transform result of (t);

[0032] By searching for the global maximum point T of A(t) maxE To capture the differences in physical mechanisms between different earthquake sources (such as the instantaneous impact of artificial blasting and the energy lag of natural earthquakes), that is: to extract the moment of maximum energy:

[0033]

[0034] Where, argmax t This represents the variable operation corresponding to the maximum value, that is, obtaining the time point corresponding to the maximum amplitude of the instantaneous amplitude envelope.

[0035] Furthermore, in step 4, the HOA algorithm is used for parameter optimization based on the Tobler trekking function. The algorithm dynamically adjusts the search speed v according to the terrain slope s between the current position and the target position, and its core position update formula is based on the Tobler trekking function:

[0036]

[0037] Where v represents the hiker's movement speed in the solution space; s represents the terrain slope between the current position and the target position, which is defined as the ratio of the difference in fitness function values ​​to the distance in the solution space, and is calculated using the following formula:

[0038]

[0039] Where f(·) represents the fitness function; X target and X old ϵ represents the hyperparameter combination of the target position and the current position, respectively; ϵ is a small constant to prevent the denominator from being zero.

[0040] The formula for updating a hiker's new location is:

[0041]

[0042] Among them, X new It is represented as the updated hyperparameter combination vector; r is represented as a random factor between [0, 1].

[0043] Through this speed control mechanism, the algorithm can achieve fast convergence and global optimization in the non-convex parameter space of the random forest.

[0044] Furthermore, the input features of the random forest classifier also include frequency domain robust dominant frequency features, the calculation formula of which is:

[0045]

[0046] Among them, Fdom This represents the robust dominant frequency characteristic in the frequency domain; FFT(·) represents the Fast Fourier Transform operation; Ω represents the effective frequency range, and Ω={f|f>f} th}, f th is a preset low-frequency cutoff threshold used to filter out low-frequency DC interference through masking operations; f represents the frequency variable.

[0047] The present invention also relates to a microseismic signal classification system, 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 the steps of the microseismic signal classification method as described above.

[0048] Furthermore, the present invention also relates to a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the aforementioned microseismic signal classification method.

[0049] The beneficial effects of this invention are as follows:

[0050] (1) Strong noise resistance and physical interpretability: This invention eliminates DC drift and background noise interference in the mining environment in principle through a mathematically defined linear detrending and energy integration mechanism. In particular, the introduction of the variational constraint model of VMD separates the effective vibration components from the complex background and explicitly quantifies the physical differences in energy release timing between artificial blasting and natural earthquakes, ensuring the authenticity of the features under low signal-to-noise ratio conditions.

[0051] (2) Adaptive optimization of model parameters: This invention utilizes the unique slope-driven speed update mechanism (Tobler function) of the HOA algorithm to effectively solve the problem of difficulty in tuning hyperparameters of random forests in discrete, non-convex spaces. Compared with traditional grid search or particle swarm optimization algorithms, this method has a faster convergence speed and a stronger ability to escape local optima, significantly improving the generalization performance and recognition accuracy of the classification model under different mining geological conditions. Attached Figure Description

[0052] Figure 1 This is a data processing flowchart of the microseismic signal classification method of the present invention.

[0053] Figure 2 This is a VMD decomposition diagram in the microseismic signal classification method of the present invention.

[0054] Figure 3 This is a violin plot showing the characteristic distribution in the microseismic signal classification method of the present invention. Detailed Implementation

[0055] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0056] like Figure 1 As shown, the present invention provides a microseismic signal classification method based on VMD physical feature extraction and HOA-optimized random forest, which mainly includes five core steps: data acquisition, physical constraint preprocessing, multi-domain robust feature extraction, HOA-RF model construction and optimization, and online classification.

[0057] Step 1: Acquire raw microseismic signals

[0058] To verify the effectiveness of the method of this invention under complex geological conditions, this embodiment first constructs a dataset through high-fidelity forward simulation. Forward simulation is performed using Sofia2D software based on the viscoelastic wave equation. To simulate the heterogeneity of real mine strata, a Von Karman-type stochastic medium model is introduced, with a correlation length set as a. x =a z =50m, Hearst exponent H=0.3, disturbance standard deviation of 10%. Three typical vibration events were simulated:

[0059] Artificial blasting: The seismic source is an expansion source, and the main frequency is set to 100Hz to simulate the high-frequency impact signal of underground blasting operations;

[0060] Mine seismic activity: The seismic source is a shear source with a dominant frequency of 20Hz to simulate the vibration generated by rock strata fracturing;

[0061] Natural earthquake: A dual-couple source was used as the seismic source, with a dominant frequency of 8Hz to simulate low-frequency, long-period signals transmitted from the far field. The data sampling rate was set to 10kHz. To simulate the harsh electromagnetic and mechanical noise environment downhole, colored noise with a signal-to-noise ratio (SNR) of 10dB was superimposed on all raw signals.

[0062] Step 2: Physical Constraint Preprocessing

[0063] like Figure 1 As shown in step 2, the acquired raw signal often contains zero drift and a trend term. Specifically, for each segment of the microseismic signal x... raw (t), fit a straight line y=at+b using the least squares method, and then calculate x(t)=x raw (t)-(at+b).

[0064] Results: After detrending, the time-domain waveform of the signal oscillates around the zero axis again, eliminating the DC component caused by instrument temperature drift or ground tilt. In the frequency domain, spurious high-energy peaks near 0 Hz are completely eliminated, providing an accurate basis for subsequent main frequency extraction. Subsequently, the signal amplitude is normalized to the [-1, 1] interval to eliminate amplitude differences caused by different source distances.

[0065] Step 3: Multi-domain Robust Feature Extraction

[0066] This embodiment extracts four features with clear physical meaning and constructs a feature vector V=[T duration ,F dom ,T maxE ,F maxE ].

[0067] (1) Duration of time-domain energy (T) duration )

[0068] To address the difficulty in defining the "start and end times" of noisy signals, this embodiment employs the energy integration method.

[0069] First, calculate the cumulative energy curve E(i) of the microseismic signal. The calculation formula is:

[0070]

[0071] Where E(i) represents the cumulative energy at the i-th sampling point; x(k) represents the amplitude of the microseismic signal at the k-th sampling point.

[0072] Define the total energy of the signal as E total Find the sampling point index t corresponding to when the accumulated energy reaches 5% of the total energy. 5% And the sampling point index t corresponding to reaching 95%. 95% .

[0073] The final time-domain energy duration characteristic is:

[0074]

[0075] Where Δt represents the sampling interval of the signal.

[0076] Physical significance: Artificial blasting releases energy extremely rapidly, T duration Typically less than 0.2s, while natural earthquakes are affected by the superposition of multiple wave phases, T duration It is usually greater than 0.5s.

[0077] (2) Frequency domain robust main frequency (F dom )

[0078] The detrended signal is subjected to a Fast Fourier Transform (FFT). A frequency mask is used to prevent low-frequency pulsation interference. The frequency corresponding to the point of maximum amplitude is searched within the effective frequency range Ω, and its calculation formula is as follows:

[0079]

[0080] Where Ω = {f|f>f} th}, f th This is the preset low-frequency cutoff threshold.

[0081] Physical significance: Experimental statistics show that the main frequency of artificial blasting is concentrated in the range of 300Hz-600Hz (high frequency region on logarithmic scale), while mine tremors and natural earthquakes are concentrated in the range of 10Hz-50Hz.

[0082] (3) Time-frequency domain maximum energy characteristics (T) based on VMD maxE )

[0083] like Figure 2 As shown, this is the core step of the present invention.

[0084] VMD decomposition process: Construct a variational problem to decompose the pure microseismic signal into K intrinsic modal components. The constrained variational model formula is as follows:

[0085]

[0086] Among them, {u k} represents the k-th intrinsic mode component obtained from the decomposition; {ω k} represents the k-th intrinsic mode component obtained from the decomposition {ω}. k} represents the center frequency of each modal component; δ(t) represents the Dirac distribution function; * represents the convolution operation; ∂ t This represents the partial derivative with respect to time; j represents the imaginary unit.

[0087] In this example, by introducing a penalty factor α and Lagrange multipliers λ, the variational model is transformed into the following augmented Lagrangian function for solution, thereby making the decomposition process more robust to noise. The overall formula is as follows:

[0088]

[0089] Set the decomposition level K=3, penalty factor α=2000, and noise margin to 0. Decompose the signal into 3 intrinsic mode components (IMF1, IMF2, IMF3). Figure 2 As shown in (b), the mode with the highest energy is selected as the dominant phase mode u. best (t) (IMF1 in this embodiment).

[0090] Feature extraction: The instantaneous amplitude envelope A(t) of the analytic signal is obtained by performing a Hilbert transform on the main shock phase mode. The calculation formula is as follows:

[0091]

[0092] Among them, û best (t) is u best The Hilbert transform result of (t).

[0093] The moment of maximum energy extraction:

[0094] Physical significance: This feature effectively captures the energy distribution pattern of the waveform. Explosive signals are instantaneously excited, with the energy peak occurring early; natural earthquake signals, due to the velocity difference between P-waves and S-waves, often have a delayed energy peak. Experiments show that the T-wave velocity of natural earthquakes... maxE Significantly greater than 0.3s.

[0095] Step 4: HOA-RF Model Construction and Optimization

[0096] Construct a random forest classifier and use the hiking optimization algorithm (HOA) to find the optimal combination of hyperparameters.

[0097] (1) Optimization target:

[0098] Number of decision trees: Search range [10, 200].

[0099] Maximum depth of the tree: search range [2, 20].

[0100] Minimum number of split samples for internal nodes: search range [2, 10].

[0101] (2) HOA Implementation Steps

[0102] Initialization: Set the population size N=10 and the maximum number of iterations T=15. Randomly initialize the positions of 10 "hikers" in the solution space.

[0103] Fitness evaluation: The accuracy of 3-fold cross-validation is used as the fitness function f(·).

[0104] Location Update: HOA simulates hikers adjusting their speed based on terrain slope by first calculating the terrain slope s between the current location and the target location.

[0105]

[0106] If the fitness of the target location (global optimum) is much higher than that of the current location (steeper slope), the hiker will move at a faster speed v, following Tobler's rule:

[0107]

[0108] The formula for updating a hiker's new location is:

[0109]

[0110] Here, r represents a random factor between [0, 1]. This allows the algorithm to quickly approach the optimal region in the early stages and perform a fine-grained search in the later stages.

[0111] Termination condition: After reaching the maximum number of iterations, output the globally optimal parameter combination.

[0112] After HOA optimization, the model converged in the 8th generation. The optimal parameter combination is: number of trees = 100, maximum depth = 12, minimum number of splits = 4.

[0113] Step 5: Online classification and output

[0114] The method of this embodiment was applied to a test set (300 samples in total, 100 samples per class), and the output results are as follows. Figure 3 As shown. By Figure 3 As can be seen, after introducing logarithmic coordinates, artificial blasting (Blast) exhibits a clear linearly separable boundary with the other two types in the frequency domain (Blast > 200Hz, Others < 50Hz). In the time domain, natural earthquakes exhibit a significant long tail characteristic. This verifies the effectiveness of the features proposed in this invention.

[0115] Combination Figure 2 , Figure 3 This indicates that, in the microseismic signal classification method based on VMD physical feature extraction and HOA optimized random forest of the present invention, the HOA-RF model exhibits strong stability and superiority in the microseismic signal classification task, and has very satisfactory classification ability.

[0116] In summary, this invention addresses the issues of feature failure and physical interpretability under strong noise by introducing physical constraint preprocessing and VMD feature extraction; it also solves the problem of difficult parameter optimization in classification models by implementing adaptive optimization of the classifier through the HOA algorithm. This method has low computational cost (feature extraction only requires milliseconds) and high accuracy, making it highly suitable for deployment in mine microseismic monitoring terminals with limited computing power.

[0117] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A microseismic signal classification method based on VMD feature extraction and HOA-optimized random forest, characterized in that, The method includes the following steps: Step 1: Deploy a microseismic sensor array in the mine monitoring area to acquire the raw microseismic monitoring signal data to be identified; Step 2: Perform physical constraint preprocessing on the original microseismic monitoring signal data. The preprocessing includes linear detrending processing and amplitude normalization processing to eliminate DC drift and non-stationary trends, and obtain a pure microseismic signal. Step 3: Feature extraction based on VMD: Perform multi-domain robust feature extraction on the pure microseismic signal to construct a multi-dimensional feature vector containing physical constraints; the multi-dimensional feature vector includes time-domain energy duration features, frequency-domain robust dominant frequency features, and the maximum energy moment of the main seismic phase features; Step 4: Construct a microseismic classifier based on random forest, and use the HOA algorithm to adaptively optimize the hyperparameter combination of the random forest to obtain the optimal microseismic classification model; Step 5: Input the multidimensional feature vector extracted in Step 3 into the optimal microseismic classification model trained in Step 4, and output the category label of the microseismic signal. The category label includes artificial blasting, mining tremor and natural earthquake.

2. The microseismic signal classification method based on VMD feature extraction and HOA-optimized random forest according to claim 1, characterized in that, In step 2, the linear detrending process involves fitting the linear trend term using the least squares method. The specific formula is as follows: ； Where, x detrend (t) represents the detrended signal data; x raw (t) represents the original microseismic monitoring signal data; a represents the slope of the linear trend term; b represents the intercept of the linear trend term; t represents the time sampling point; The slope a and intercept b are obtained by fitting the original microseismic monitoring signal using the least squares method, and are used to eliminate the DC component and linear drift in the original microseismic monitoring signal.

3. The microseismic signal classification method based on VMD feature extraction and HOA-optimized random forest according to claim 1, characterized in that, Step 3 extracts the time-domain energy duration feature using the energy integration method, as follows: First, calculate the cumulative energy curve E(i) of the pure microseismic signal: ； Where E(i) represents the cumulative energy at the i-th sampling point; x(k) represents the amplitude of the pure microseismic signal at the k-th sampling point; Then, the time-domain energy duration characteristic is calculated using the following formula: ； Among them, T duration Represented as the time-domain energy duration characteristic; t 95% This represents the index of the sampling point when the accumulated energy reaches 95% of the total energy; t 5% This represents the sampling point index when the accumulated energy reaches 5% of the total energy; ∆t represents the sampling interval of the signal.

4. The microseismic signal classification method based on VMD feature extraction and HOA-optimized random forest according to claim 1, characterized in that, Step 3, the VMD-based feature extraction process, includes: A variational problem is constructed, decomposing the pure microseismic signal into K intrinsic modal components. The model formula is as follows: ； Among them, {u k } represents the k-th intrinsic mode component obtained from the decomposition; {ω k } represents the center frequency of each modal component; δ(t) represents the Dirac distribution function; * represents the convolution operation; ∂ t This represents the partial derivative with respect to time; j represents the imaginary unit.

5. The microseismic signal classification method based on VMD physical feature extraction and HOA optimized random forest according to claim 1, characterized in that, In step 3, during the feature extraction process based on VMD, the method for extracting the maximum energy moment feature of the mainshock phase is as follows: First, calculate the energy of each modal component, and then select the mode with the highest energy from the energy of each modal component as the dominant seismic phase mode; Then, the instantaneous amplitude envelope of the analytic signal is obtained by performing a Hilbert transform on the main phase mode: ； Where A(t) represents the instantaneous amplitude envelope; u best (t) represents the main seismic phase mode signal; û best (t) is represented as u best The Hilbert transform result of (t); The moment of maximum energy extraction: ； Where, argmax t This represents the variable operation corresponding to the maximum value, that is, obtaining the time point corresponding to the maximum amplitude of the instantaneous amplitude envelope.

6. The microseismic signal classification method based on VMD feature extraction and HOA-optimized random forest according to claim 1, characterized in that, In step 4, the HOA algorithm is used to optimize the hyperparameters of the random forest, and its core location update formula is based on the Tobler trekking function: ； Where v represents the hiker's movement speed in the solution space; s represents the terrain slope between the current position and the target position, which is defined as the ratio of the difference in fitness function values ​​to the distance in the solution space, and is calculated using the following formula: ； Where f(·) represents the fitness function; X target and X old ϵ represents the hyperparameter combination of the target position and the current position, respectively; ϵ is a small constant to prevent the denominator from being zero. The formula for updating a hiker's new location is: ； Among them, X new It is represented as the updated hyperparameter combination vector; r is represented as a random factor between [0, 1].

7. The microseismic signal classification method based on VMD feature extraction and HOA-optimized random forest according to claim 1, characterized in that, In step 4, the input features of the random forest microseismic classifier include the frequency domain robust dominant frequency feature, which is calculated using the following formula: ； Among them, F dom This represents the robust dominant frequency characteristic in the frequency domain; FFT(·) represents the Fast Fourier Transform operation; Ω represents the effective frequency range, and Ω={f|f>f} th }, f th is a preset low-frequency cutoff threshold used to filter out low-frequency DC interference through masking operations; f represents the frequency variable.

8. A microseismic signal classification system, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method as described in any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1 to 7.

Images