Microseismic signal classification method based on genetic algorithm optimization

Abstract

The invention discloses a micro-seismic signal classification method based on genetic algorithm optimization. The method comprises the following steps: data acquisition: acquiring waveform data containing multiple types of micro-seismic event forms; preprocessing: carrying out mean value removal and band-pass filtering preprocessing on the obtained waveform data, and removing baseline drift and high-frequency noise interference; envelope calculation: respectively calculating a time domain envelope and a frequency spectrum envelope for the preprocessed microseismic signal; feature extraction: extracting a plurality of statistics and time frequency parameters from the time domain envelope and the frequency spectrum envelope to form a feature complete set; feature selection and classifier optimization: utilizing a genetic algorithm to perform joint search optimization on structure parameters and feature subsets of the multi-layer perceptron; and training and verification: training the multilayer perceptron by using the optimized feature subset and the network structure, and evaluating the classification accuracy on a verification set. The method can realize high-precision classification, has robustness to noise, and can be widely applied to the fields of microseism monitoring, blasting event monitoring, other vibration source identification and the like.

Description

Technical Field

[0001] The present invention relates to the technical field of microseismic signal processing and identification, and in particular to a microseismic signal classification method based on genetic algorithm optimization. Background Art

[0002] Microseismic monitoring networks record a large number of microseismic signals from various source types daily, including local earthquakes (LE), distant earthquakes (DE), quarry blasts (QB), and noise (NS). Without automated and effective classification of this massive amount of data, research and early warning efforts will face enormous challenges. To this end, researchers around the world have developed and proposed numerous automated microseismic signal classification methods. These methods can be roughly divided into two categories: those that rely on individual feature extraction and those that do not. The latter category, such as convolutional neural networks, performs well in large data sets but are often difficult to interpret and require massive amounts of labeled data. Summary of the Invention

[0003] The purpose of the present invention is to overcome the above technical problems and propose a microseismic signal classification method based on genetic algorithm optimization, which can achieve high-precision automatic classification in the case of strong noise, limited data or diverse signal types.

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

[0005] A microseismic signal classification method based on genetic algorithm optimization includes the following steps:

[0006] S1. Data acquisition: Acquire waveform data containing multiple types of microseismic event morphologies;

[0007] S2. Preprocessing: Perform de-averaging and bandpass filtering on the acquired waveform data to remove baseline drift and high-frequency noise interference;

[0008] S3. Envelope calculation: Calculate the time domain envelope and spectrum envelope of the preprocessed microseismic signal;

[0009] S4. Feature extraction: Extract various statistics and time-frequency parameters from the time domain envelope and spectral envelope to form a complete feature set;

[0010] S5. Feature Selection and Classifier Optimization: A genetic algorithm is used to jointly search and optimize the multilayer perceptron structural parameters and feature subsets. The genetic algorithm selects and combines features from the full set of features. The optimized and selected feature combination is the feature subset.

[0011] S6. Training and Validation: Train the multilayer perceptron using the optimized feature subset and network structure, and evaluate the classification accuracy on the validation set.

[0012] Optionally, in step S3, the calculation formula of the time domain envelope is:

[0013]

[0014] Where z(t) is the real-valued microseismic waveform function, H[z(t)] is the Hilbert transform result of z(t), and E t (t) is the time domain envelope, t is time;

[0015] The calculation formula of the spectrum envelope is:

[0016] E f (f)=Smooth|Z(f)|

[0017] Among them, Z(f) is the result of fast Fourier transform of z(t), Smooth represents the median filter operation, and E f (f) is the spectrum envelope, and f is the frequency.

[0018] Optionally, the features in step S4 include but are not limited to:

[0019] Maximum value of the time domain envelope:

[0020] Time domain envelope duration: T es =t e -t s , where t s Indicates the time point when the time domain envelope exceeds the set threshold for the first time, t e Indicates the time point when the time domain envelope falls below the threshold;

[0021] Spectral envelope centroid: where f n Represents discrete frequency points, E f (f n ) indicates that at frequency f n The spectrum envelope value at ;

[0022] Envelope gradient / energy ratio: Among them E t (t n ) represents the time point t n The time domain envelope value at , E t (t n-1 ) represents the time point t n-1 The time domain envelope value at , N represents the total number of sampling points.

[0023] Optionally, step S5 includes the following steps:

[0024] (a) Chromosome encoding: Encode the feature selection into a binary string, where "1" indicates the selected feature and "0" indicates the discarded feature; it also encodes the network structure parameters of the multilayer perceptron;

[0025] (b) Population initialization: randomly generate the initial population;

[0026] (c) Fitness evaluation: The fitness of each chromosome is evaluated based on the classification accuracy and model complexity;

[0027] (d) Selection operation: selecting excellent individuals;

[0028] (e) Crossover operation: perform single-point or double-point crossover on the selected chromosome;

[0029] (f) Mutation operation: performing gene mutation on chromosomes;

[0030] (g) Repeat steps (c) to (f) until the termination condition is met and output the optimal solution.

[0031] Optionally, in step (c), the fitness function is:

[0032]

[0033] Among them, α and β are weight factors, ||x||1 indicates how many features are selected, n indicates the number of hidden layer neurons or network size, and C v Indicates the number of correctly classified samples, N v Indicates the validation set size.

[0034] Optionally, the multilayer perceptron structure includes an input layer, a hidden layer and an output layer, the input layer is used to receive the selected feature vector; the hidden layer includes neurons of a number optimized by a genetic algorithm; and the output layer includes output neurons of a number equal to the number of microseismic event types.

[0035] Optionally, in step S1, the microseismic event types include local microseismic events, teleseismic events, quarry or mine blasting, and noise interference.

[0036] Optionally, the bandpass filtering frequency range in step S2 is 1-45 Hz.

[0037] In summary, the technical effects and advantages of the present invention are as follows: the feature extraction method of the method of the present invention is simple and effective: a series of features based on the time domain envelope and spectral envelope of microseismic signals are proposed, and the optimal feature set can be automatically selected; the use of genetic algorithms: the multi-layer perceptron structure and feature subsets are jointly searched and optimized through genetic algorithms, which makes it easier to find the global optimal solution and avoid manual repeated parameter adjustment; reducing noise interference: the envelope features are more robust to noise, and can retain the necessary classification information even under the interference of microseismic sub-phases such as P waves, S waves and coda waves. DETAILED DESCRIPTION

[0038] The following is a clear and complete description of the technical solutions in the embodiments of the present invention. Obviously, the embodiments described are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making any creative efforts are within the scope of protection of the present invention.

[0039] In addition, the technical features involved in the different embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0040] This embodiment provides a microseismic signal classification method based on genetic algorithm optimization, comprising the following steps:

[0041] S1. Data Acquisition: Waveform data containing various types of microseismic events are acquired from the microseismic monitoring array. These events include localized microseismic events, teleseismic events, quarry or mine blasting, and noise disturbances.

[0042] S2. Preprocessing: Perform basic data preprocessing such as de-averaging and bandpass filtering to remove significant baseline drift and high-frequency noise. The bandpass filter frequency range is 1-45 Hz.

[0043] S3. Calculation of time domain envelope and spectrum envelope: Calculate the time domain envelope and spectrum envelope of the preprocessed microseismic signal respectively, which serve as the core of subsequent feature extraction.

[0044] Time domain envelope calculation:

[0045] Let the microseismic signal waveform be z(t). First construct the complex analytical signal:

[0046] z + (t)=z(t)+j·H[z(t)] (1)

[0047] in,

[0048] H[z(t)] represents the Hilbert transform result of z(t);

[0049] j represents the imaginary unit.

[0050] The time domain envelope can be obtained from the modulus of the complex analytic signal, as shown in Formula 2:

[0051]

[0052] in,

[0053] z(t) is a function of the real-valued microseismic waveform, and the physical quantity is the vibration velocity that changes with time t.

[0054] H[z(t)] is the Hilbert transform result of z(t), which can be understood as the compensation of the phase information of the microseismic signal.

[0055] E t (t) is the time domain envelope, and its physical meaning is the smooth outer envelope of the signal energy distribution over time.

[0056] Spectral envelope calculation:

[0057] Perform fast Fourier transform (FFT) on z(t) to get the spectrum Z(f). To reduce the influence of noise, Z(f) can be smoothed by median filtering to get the spectrum envelope E. f (f), namely:

[0058] E f (f)=Smooth|Z(f)| (3)

[0059] Where f is the frequency, and the physical quantity is the energy distribution of the signal in different frequency bands. Z(f) is the result of fast Fourier transform of z(t), Smooth represents the median filter operation, and E f (f) is the spectrum envelope. This envelope avoids directly using the sharp peaks in the spectrum graph and instead extracts the overall outer envelope, making it easier to distinguish the main frequency range or energy distribution characteristics of different types of signals.

[0060] S4. Feature extraction: Extract various statistics and time-frequency parameters from the envelope to form a full set of features.

[0061] Envelope feature extraction:

[0062] Based on the time domain envelope E t (t) and the spectrum envelope E f (f) Features that can distinguish different microseismic types can be extracted. Features include but are not limited to the maximum value of the time domain envelope, the duration of the time domain envelope, the centroid of the spectrum envelope, and the envelope gradient / energy ratio.

[0063] Maximum value of the time domain envelope:

[0064]

[0065] This feature measures the energy peak of the signal in the time domain.

[0066] Time domain envelope duration:

[0067] T es =t e -t s (5)

[0068] where t s Indicates the time point when the time domain envelope exceeds the set threshold for the first time, t e It indicates the time point when the time domain envelope falls below the threshold, and is used to measure the duration of effective events in continuous microseismic records.

[0069] Spectral envelope centroid:

[0070]

[0071] Among them, f n represents discrete frequency points,

[0072] E f (f n ) indicates that at frequency f n The spectral envelope value at . This feature is used to reflect the centroid of the overall frequency distribution.

[0073] Envelope gradient / energy ratio:

[0074]

[0075] Among them, E t (t n ) represents the time point t n The time domain envelope value at , E t (t n-1 ) represents the time point t n-1 The time domain envelope value at , N represents the total number of sampling points.

[0076] This type of feature measures the smoothness or steepness of the envelope curve over time or frequency, and can distinguish whether the rising segment of the blasting signal and the natural microseismic signal has a sudden change.

[0077] The above features can form a candidate feature set, which is then automatically selected by the genetic algorithm.

[0078] S5. Feature Selection and Classifier Optimization (Genetic Algorithm): Genetic Algorithm (GA) is used to automatically search for multi-layer perceptron (MLP) structural parameters and feature subsets to obtain the optimal classifier configuration.

[0079] Among them, the feature subset is part of the full feature set. In the S4 feature extraction step, by extracting multiple statistics and time-frequency parameters from the time domain envelope and spectral envelope, a set containing all candidate features is formed, which is the full feature set. In the S5 feature selection and classifier optimization step, the genetic algorithm will screen and combine from the full feature set to select a set of optimal feature combinations for subsequent classification tasks; the feature combination after optimization and screening is the feature subset. The output of the S4 feature extraction step (full feature set) is the input of the S5 feature selection and classifier optimization step. The goal of S4 is to extract as many potentially useful features as possible to form a comprehensive feature pool (full feature set). The goal of S5 is to automatically find the most effective and streamlined feature combination (feature subset) for the classification result from this large and potentially redundant feature pool, while optimizing the structural parameters of the classifier itself. This process can remove invalid or interfering features, reduce model complexity, and improve classification accuracy and robustness.

[0080] Multilayer Perceptron (MLP) classifier design:

[0081] This embodiment uses a multi-hidden layer MLP as a classifier, and the network structure is shown as follows:

[0082] Input layer: receives the selected d-dimensional feature vector, x∈R d . where x represents a single feature vector input to the multilayer perceptron (MLP), and R d Represents a d-dimensional real vector space, where R refers to the set of real numbers and d refers to the dimension of the vector (i.e., the number of features contained in the feature subset).

[0083] Hidden layer: Contains a number of neurons optimized by a genetic algorithm, for example, n neurons, each of which stores a set of trainable weights and biases; the activation function can be ReLU, Sigmoid, or Tanh, etc.

[0084] Output layer: contains output neurons equal to the number of microseismic event types. In this embodiment, it contains 4 output neurons, corresponding to the 4 microseismic event types (LE, DE, QB, NS). The activation function generally uses Sigmoid to output probability distribution.

[0085] Mathematically, the output of neurons in layer l can be described by the following equation:

[0086] a (l) =f(W (l) a (l-1) +b (l) ) (8)

[0087] in,

[0088] a(l) is the output vector of the lth layer,

[0089] W (l) is the weight matrix from the previous layer to the lth layer,

[0090] a (l-1) is the output vector of the l-1th layer,

[0091] b (l) is the bias vector,

[0092] f(·) is the activation function.

[0093] Genetic algorithm (GA) is used for feature selection and network structure optimization:

[0094] This embodiment uses a genetic algorithm to search for feature subsets and network hyperparameters of a multilayer perceptron structure. The specific process includes the following:

[0095] (a) Chromosome encoding, binary representation of feature subsets, encoding feature selection as a binary string, where "1" indicates a selected feature and "0" indicates a discarded feature; it also encodes the network structure parameters of the multilayer perceptron;

[0096] The candidate feature dimension is mapped to a binary bit string, with 0 or 1 indicating whether the feature is selected.

[0097] In addition, the chromosome is added with the encoding of hyperparameters such as the number of hidden layer neurons and learning rate.

[0098] Assume that there are D optional features, and use a binary string containing D gene bits to represent the feature selection result:

[0099] chrom=(x1,x2,...,x D ),x i ∈{0,1} (9)

[0100] Among them, when x i =1 means the i-th feature is selected. When x i = 0, it means that the i-th feature is discarded.

[0101] Encoding of multilayer perceptron hyperparameters:

[0102] If you need to automatically search for the appropriate number of hidden layer neurons n or other hyperparameters (such as learning rate, etc.), you can add some special bits or integer coding segments to the chromosome. For example, use a 5-bit binary number to represent the number of hidden layer neurons n in [1,2 5 ] within the range of .

[0103] In this way, each chromosome in the genetic algorithm uniquely corresponds to: a feature subset and a multilayer perceptron network configuration.

[0104] (b) Population initialization: randomly generate an initial population; randomly generate several chromosomes, each chromosome represents a set of feasible solutions.

[0105] Set the population size to Npop. GA first randomly or semi-randomly generates Npop chromosomes. Each chromosome represents a potential solution (i.e., a set of feature selection + network hyperparameter combination).

[0106] At this time, the crossover probability p can be set c , mutation probability p m Wait for GA super parameters to prepare for subsequent evolution and iteration.

[0107] (c) Fitness evaluation: The fitness of each chromosome is evaluated based on the classification accuracy and model complexity;

[0108] Assemble the feature subset corresponding to the chromosome and the multi-layer perceptron structure, train it on the training set, and then test its classification accuracy on the validation set or cross-validation. Let the number of correctly classified samples be Cv and the validation set size be Nv, then the fitness can be defined as:

[0109]

[0110] Considering the complexity of the model, it is also necessary to add a penalty term to the fitness. For example, a penalty coefficient is added to the number of selected features ||x||1 (that is, the number of 1s in the binary string) or to the network size. The fitness function is formed as follows:

[0111]

[0112] Among them, α, β are trade-off factors,

[0113] ||x||1 indicates how many features are selected.

[0114] n represents the number of hidden layer neurons or network size;

[0115] This method can prevent excessive feature subsets or excessive network size, thus ensuring the simplicity of the model.

[0116] (d) Selection operation, selecting excellent individuals;

[0117] Common methods include roulette wheel selection and tournament selection. The essence is to make chromosomes with high fitness values ​​more likely to enter the next generation.

[0118] Roulette wheel selection: The fitness of all individuals is normalized and regarded as a "probability distribution", and the chance of randomly selecting an individual is proportional to its fitness.

[0119] Tournament selection: Randomly select 2-3 individuals for comparison, and select the one with the highest fitness to enter the next generation, and repeat this process.

[0120] Goal: To retain the chromosomes with excellent performance in the current population and eliminate the weak ones.

[0121] In this embodiment, the tournament selection is selected.

[0122] (e) Crossover operation, which performs single-point or double-point crossover on the selected chromosomes with a certain probability;

[0123] Crossover is the main source of GA to obtain potential excellent offspring. The typical approach is to randomly pair the selected parent chromosomes with a crossover probability p. c implement:

[0124] Single-point crossover: A cut is made at a random position on the chromosome, and the first half is taken from father 1 and the second half is taken from father 2, or vice versa:

[0125] Double-point crossover: Randomly select two cut points on the chromosome to form more recombination methods.

[0126] After executing the crossover, new offspring chromosomes are obtained, and their feature subset encoding and network structure encoding may undergo combinatorial changes and be inherited to the next generation.

[0127] (f) Mutation operation, which performs gene mutation on chromosomes with a small probability;

[0128] In order to maintain diversity and prevent the population from converging prematurely, GA will perform mutation on offspring genes with a probability p after crossover. m Perform mutation:

[0129] If a gene is 0, it is flipped to 1; if it is 1, it is flipped to 0;

[0130] For integer coding segments (such as the number of neurons), operations such as random increase or decrease of 1 can be performed.

[0131] Mutation can generate new search directions and avoid falling into local optimality.

[0132] (g) Termination condition and optimal solution: repeat steps (c) to (f) until the termination condition is met and output the optimal solution;

[0133] Repeat the selection-crossover-mutation iterative process until one or more of the following conditions are met:

[0134] Reach the maximum number of iterations N gen ;

[0135] The population fitness no longer improves over a long period of time;

[0136] The fitness reaches the pre-set target accuracy.

[0137] Finally, the chromosome with the highest fitness is regarded as the optimal solution, and the feature subset and hyperparameter configuration it contains are the selected / designed microseismic signal classifier.

[0138] S6. Training and Validation: The constructed multilayer perceptron (MLP) is trained using the optimized feature subset and network structure, and the classification accuracy is evaluated on the validation set as the fitness of the genetic algorithm (GA); finally, the generalization performance is tested on the test set.

[0139] S7. Online application: Deploy the optimized classifier to the microseismic monitoring network to realize automatic microseismic event identification.

[0140] This embodiment addresses the practical problems of limited data at many monitoring sites, significant noise interference, and complex manual labeling. The artificial neural network (especially the multi-layer perceptron structure) optimized based on the genetic algorithm is combined with targeted feature selection to achieve high-precision, low-complexity microseismic signal classification under limited training set conditions. The method proposed in the present invention is not only suitable for real-time online monitoring, but also takes into account the applicability and interpretability in small-scale databases. The present invention uses the time domain envelope and the spectrum envelope as the basis, extracts a variety of statistics and time-frequency parameters from them to form a full set of features, and uses the genetic algorithm to jointly optimize the feature subset and the multi-layer perceptron structure to obtain the optimal classifier configuration. The present invention can achieve high-precision classification in the case of strong noise, limited data or diverse signal types, is robust to noise, reduces repeated parameter adjustment work, is suitable for real-time or quasi-real-time deployment, and can be widely used in microseismic monitoring, blasting event monitoring, and other vibration source identification fields.

[0141] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the common and similar parts between the various embodiments. For the systems disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the method description.

[0142] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The above examples are only intended to help understand the method and core concept of the present invention. At the same time, those skilled in the art will find that the specific implementation methods and application scopes may vary based on the concept of the present invention. In summary, the contents of this specification should not be construed as limiting the present invention.

Claims

1. A microseismic signal classification method based on genetic algorithm optimization, characterized in that: The following steps are involved: S1. Data acquisition: Acquire waveform data containing multiple types of microseismic event morphologies; S2. Preprocessing: Perform de-averaging and bandpass filtering on the acquired waveform data to remove baseline drift and high-frequency noise interference; S3. Envelope calculation: Calculate the time domain envelope and spectrum envelope of the preprocessed microseismic signal; S4. Feature extraction: Extract various statistics and time-frequency parameters from the time domain envelope and spectral envelope to form a complete feature set; S5. Feature Selection and Classifier Optimization: A genetic algorithm is used to jointly search and optimize the multilayer perceptron structural parameters and feature subsets. The genetic algorithm selects and combines features from the full set of features. The optimized and selected feature combination is the feature subset. S6. Training and Validation: Train the multilayer perceptron using the optimized feature subset and network structure, and evaluate the classification accuracy on the validation set.

2. The microseismic signal classification method based on genetic algorithm optimization according to claim 1 is characterized in that: In step S3, the calculation formula of the time domain envelope is: Where z(t) is the real-valued microseismic waveform function, H[z(t)] is the Hilbert transform result of z(t), and E t (t) is the time domain envelope, t is time; The calculation formula of the spectrum envelope is: E f (f)=Smooth|Z(f)| Among them, Z(f) is the result of fast Fourier transform of z(t), Smooth represents the median filter operation, and E f (f) is the spectrum envelope, and f is the frequency.

3. The microseismic signal classification method based on genetic algorithm optimization according to claim 1 is characterized in that: The features in step S4 include but are not limited to: Maximum value of the time domain envelope: Time domain envelope duration: T es =t e -t s , where t s Indicates the time point when the time domain envelope exceeds the set threshold for the first time, t e Indicates the time point when the time domain envelope falls below the threshold; Spectral envelope centroid: where f n Represents discrete frequency points, E f (f n ) indicates that at frequency f n The spectrum envelope value at ; Envelope gradient / energy ratio: Among them E t (t n ) represents the time point t n The time domain envelope value at , E t (t n-1 ) represents the time point t n-1 The time domain envelope value at , N represents the total number of sampling points.

4. The microseismic signal classification method based on genetic algorithm optimization according to claim 1 is characterized in that: Step S5 includes the following steps: (a) Chromosome encoding: Encode the feature selection into a binary string, where "1" indicates the selected feature and "0" indicates the discarded feature; at the same time, it encodes the network structure parameters of the multilayer perceptron; (b) Population initialization: randomly generate the initial population; (c) Fitness evaluation: The fitness of each chromosome is evaluated based on the classification accuracy and model complexity; (d) Selection operation: selecting excellent individuals; (e) Crossover operation: perform single-point or double-point crossover on the selected chromosome; (f) Mutation operation: performing gene mutation on chromosomes; (g) Repeat steps (c) to (f) until the termination condition is met and output the optimal solution.

5. The microseismic signal classification method based on genetic algorithm optimization according to claim 4 is characterized in that: In step (c), the fitness function is: Among them, α and β are weight factors, ||x||1 indicates how many features are selected, n indicates the number of hidden layer neurons or network size, and C v Indicates the number of correctly classified samples, N v Indicates the validation set size.

6. The microseismic signal classification method based on genetic algorithm optimization according to claim 1 is characterized in that: The multilayer perceptron structure includes an input layer, a hidden layer and an output layer, wherein the input layer is used to receive the selected feature vector; the hidden layer includes neurons of a number optimized by a genetic algorithm; and the output layer includes output neurons of a number equal to the number of microseismic event types.

7. The microseismic signal classification method based on genetic algorithm optimization according to claim 1 is characterized in that: In step S1 , the types of microseismic events include local microseismic events, teleseismic events, quarry or mine blasting, and noise disturbance.

8. The microseismic signal classification method based on genetic algorithm optimization according to claim 1 is characterized in that: The frequency range of the bandpass filter in step S2 is 1-45 Hz.

Images