Mine micro-seismic source positioning method and device and micro-seismic event inversion method

Abstract

The invention discloses an underground micro-seismic source positioning method and device and a micro-seismic event inversion method, and belongs to the technical field of mine micro-seismic monitoring. The method comprises the following steps: acquiring micro-seismic signal data under a mine in real time, and performing stream reading; preprocessing the acquired micro-seismic signal data; based on the first arrival signal-to-noise ratio, the ending signal-to-noise ratio and the dominant frequency features, performing preliminary screening on the microseismic waveform window in combination with a feature screening algorithm; picking up the first arrival time and the end time; and micro-seismic source positioning is carried out through combination of a least square method and a Geiger algorithm. Based on the feature screening algorithm and the abnormal window elimination algorithm, the method can effectively extract the micro-seismic event window, eliminates the influence of explosion or other types of interference signals, and facilitates the guarantee of the positioning precision of the micro-seismic source.

Description

Technical Field

[0001] The present invention belongs to the technical field of mine microseismic monitoring, and more specifically, relates to a method and device for locating a microseismic source in a mine and a microseismic event inversion method. Background Art

[0002] Microseismic monitoring in mines is a key technology for geological disaster early warning and mine safety assurance. Microseismic events in mines often indicate changes in geological structure, such as rock fractures and fault activity. These changes can trigger mine disasters such as rockbursts and landslides. Therefore, real-time monitoring and precise location of microseismic events in mines, as well as the effective identification of potential cracks and unstable areas within mines, are crucial for preventing and controlling mine disasters.

[0003] After searching, the Chinese patent application No. 2024108028297 discloses a method for locating the source of microseismic signals induced by overburden fracture based on key layer identification. The application uses a fixed bandpass filter to filter the original signal, then performs wavelet decomposition, and uses the long-short time window comparison method to extract the arrival time of the initial value of the microseismic signal received by each sensor. Although this method is simple and easy to implement, the short time window may be too sensitive to high-frequency noise, especially the fixed range bandpass filter may lead to misjudgment. In addition, the window length of the long-short time window comparison method needs to be set in advance. If the signal characteristics change, the fixed window length may reduce the detection accuracy, especially in the case of large signal interference.

[0004] For example, the Chinese patent application No. 2023105451406 discloses a method, device, computer equipment and storage medium for intelligent positioning of microseismic sources. The application adopts a particle swarm algorithm to solve the preliminary source position, and then uses the preliminary source position as the initial value of the least squares method to solve the precise source position. Although the particle swarm algorithm can avoid falling into the local optimal solution, it is still possible to stay in the local optimal solution when the source location is complex or the data noise is large, thereby affecting the accuracy of source positioning. In addition, if the solution of the particle swarm algorithm is affected by noise or other uncertain factors, the initial value may be inaccurate, which in turn affects the effect of the subsequent least squares method. Summary of the Invention

[0005] In order to solve the technical problem of relatively low accuracy in locating microseismic sources in mines in the prior art, the present invention provides a method and device for locating microseismic sources in mines, thereby realizing real-time monitoring and high-precision positioning of microseismic events.

[0006] The present invention further provides a method for inverting microseismic events in mines, which performs inversion based on the precise positioning results of microseismic sources, thereby providing strong technical support for mine safety monitoring and geological disaster early warning.

[0007] In order to achieve the above object, the technical solution provided by the present invention is:

[0008] A first aspect of the present invention provides a method for locating a microseismic source in a mine, comprising:

[0009] Real-time collection of microseismic signal data in mines, and the use of distributed stream processing and message queues to achieve data streaming and real-time processing;

[0010] Preprocessing the collected microseismic signal data, wherein the preprocessing includes dual filtering processing in the time domain and frequency domain;

[0011] Based on the initial arrival signal-to-noise ratio, final arrival signal-to-noise ratio and main frequency characteristics, and combined with the feature screening algorithm, the microseismic waveform window is preliminarily screened;

[0012] The envelope of the microseismic signal is extracted using Hilbert transform, and the first arrival time and end time of the microseismic signal are picked up by combining the threshold detection method.

[0013] The abnormal window elimination algorithm is used to further eliminate abnormal data;

[0014] The microseismic source is located by combining the least squares method and the Geiger algorithm. First, the least squares method is used to provide the initial iteration point, and then the Geiger algorithm is used to iteratively calculate the source position.

[0015] According to any of the technical solutions described in the first aspect of the present invention, geophones are used to collect microseismic signal data in a mine in real time, and the microseismic signal data collected by each geophone in real time are classified according to topics and published to corresponding topic partitions to achieve effective division and organization of data streams.

[0016] According to any technical solution described in the first aspect of the present invention, the time domain and frequency domain dual filtering processing specifically includes:

[0017] Perform DC removal processing on the data to eliminate the DC component of the signal;

[0018] Perform fast Fourier transform on the DC-removed signal to convert the time domain signal into a frequency domain signal;

[0019] In frequency domain filtering, median filtering is first used to remove shock peaks in the spectrum, smooth the spectrum curve, and eliminate abnormal peaks caused by noise or equipment errors, thereby effectively filtering out high-frequency noise and enhancing signal stability. Then, a peak fitting algorithm is used to correct low-amplitude dips in the spectrum to optimize the signal's spectral shape and improve the accuracy of feature extraction. Finally, the signal's main frequency, low cutoff frequency, and high cutoff frequency are obtained through analysis.

[0020] In time domain filtering, the signal is bidirectionally filtered through an FIR filter based on the signal main frequency, low cutoff frequency, and high cutoff frequency obtained by frequency domain analysis to remove low-frequency and high-frequency noise in the signal and obtain a filtered signal.

[0021] Specifically, the frequency range and order of the FIR filter are set according to the signal main frequency, high and low cutoff frequencies and other characteristics extracted in the frequency domain. By processing the signal through bidirectional filtering (filtfilt), phase delay can be eliminated and useless low-frequency and high-frequency noise can be removed.

[0022] According to any of the technical solutions described in the first aspect of the present invention, when preprocessing collected microseismic signal data, a sliding window mechanism is used to read real-time data to segment the multi-channel microseismic signal data. The windowed data is then subjected to dual filtering in both the time and frequency domains. The size and step size of the segmented window are set as variable parameters, enabling the reading and processing of real-time data using a dynamic windowing approach.

[0023] According to any technical solution described in the first aspect of the present invention, the preliminary screening of the microseismic waveform window based on the initial arrival signal-to-noise ratio, the final arrival signal-to-noise ratio and the main frequency characteristics, combined with a feature screening algorithm, specifically includes:

[0024] After signal preprocessing, the first arrival signal-to-noise ratio of each window is calculated, and the window with the largest first arrival signal-to-noise ratio in each channel is selected;

[0025] Taking the window with the largest initial arrival signal-to-noise ratio as the benchmark, select several windows before and after the window, totaling n candidate windows, to expand the window range; more preferably, select the window before and two windows after the window with the largest signal-to-noise ratio, for a total of four candidate windows;

[0026] In the above candidate windows, the windows are screened according to the dominant frequency characteristics;

[0027] Among the windows that meet the main frequency condition, the window with the largest ending signal-to-noise ratio is selected as the final signal window.

[0028] According to any technical solution described in the first aspect of the present invention, the calculation formulas for the initial arrival signal-to-noise ratio and the final arrival signal-to-noise ratio are respectively as follows:

[0029]

[0030] Where:

[0031] SNR initial represents the first arrival signal-to-noise ratio;

[0032] max_peak is the index of the maximum peak in the array of all peaks;

[0033] p i is the value of the i-th peak;

[0034] SNR end Indicates the end signal-to-noise ratio;

[0035] N is the total number of all peaks.

[0036] According to any technical solution described in the first aspect of the present invention, the use of an abnormal window elimination algorithm to further eliminate abnormal data specifically includes:

[0037] The first arrival time and the end time of the microseismic signal are converted into a global index to map the local index to the corresponding position of the global time series;

[0038] Based on the preliminarily screened multi-channel microseismic event windows, the median M of the first arrival time index and its median absolute difference MAD are calculated; according to the preset threshold T, the data are filtered using MAD as the criterion to obtain a filtered data set;

[0039] Extracting and combining data values ​​from the filtered data set to generate all possible combinations of n data elements, where n is a positive integer of at least 4;

[0040] Outlier detection is performed on each data set based on the first arrival time index of the window. The initial outliers are first marked using the isolation forest algorithm. Then, a greedy method based on sorting is used to check whether the difference between every two data points in the group is less than the preset range L. If so, the combination is returned as a valid combination; otherwise, other combinations are checked.

[0041] If at least one valid combination is found, the corresponding original data point in the data set is retained; otherwise, data processing is terminated and a prompt message indicating that there is no valid data is output;

[0042] Convert the first arrival time index in the valid combination into a timestamp and store the timestamp and microseismic signal geophone coordinates.

[0043] A second aspect of the present invention provides a microseismic source locating device in a mine, comprising:

[0044] The data acquisition and streaming reading module is used to collect microseismic signal data in the mine in real time, and realizes data streaming and real-time processing through distributed stream processing and message queues;

[0045] A microseismic signal data preprocessing module is used to preprocess the collected microseismic signal data, wherein the preprocessing includes dual filtering processing in the time domain and frequency domain;

[0046] The microseismic window selection module is used to preliminarily screen the microseismic waveform window based on the initial arrival signal-to-noise ratio, the final arrival signal-to-noise ratio, and the main frequency characteristics, combined with the feature screening algorithm;

[0047] First arrival time and end time picking module,

[0048] An abnormal window rejection module is used to extract the envelope of the microseismic signal using Hilbert transform and, in combination with the threshold detection method, to pick up the first arrival time and the end time of the microseismic signal; and

[0049] The microseismic source location module is used to locate the microseismic source by combining the least squares method and the Geiger algorithm.

[0050] A third aspect of the present invention provides a microseismic event inversion method, comprising:

[0051] Using any of the methods described in the first aspect of the present invention to locate the microseismic source in a mine;

[0052] Based on the coordinates of microseismic events, the DBSCAN algorithm is used to identify areas with higher density to eliminate noise points;

[0053] The Gaussian mixture model is used to perform cluster analysis on microseismic events based on Gaussian distribution;

[0054] For each identified microseismic event cluster, RANSAC is applied to identify fracture surfaces and faults.

[0055] According to the technical solution described in the third aspect of the present invention, the DBSCAN algorithm is used to identify areas with higher density to eliminate noise points. The specific judgment criteria are: within the neighborhood of radius eps, if the number of neighboring points of a point is ≥ min_samples, it is judged as a core point and participates in clustering; otherwise, it is a boundary point or noise point and is eliminated.

[0056] Furthermore, in terms of parameter selection, we first estimate eps and min_samples through spatial distribution statistics, then use K-distance graphs to assist in judgment. We then conduct multiple rounds of clustering experiments on real data to verify the clustering effect. As a preferred solution, the parameter settings are as follows: eps = 800, min_samples = 50.

[0057] According to the technical solution described in the third aspect of the present invention, the Gaussian mixture model is used to perform cluster analysis on microseismic events based on Gaussian distribution, specifically including: setting the initial distribution number and parameters of the Gaussian mixture model according to the geometric characteristics of the underground mine tunnels and coal seams, prioritizing clustering the roof and floor areas of the tunnels, and training through the expectation maximization algorithm; the trained GMM divides the data into multiple clusters, combines the distribution patterns of faults and cracks in the working face under the mine, dynamically adjusts the number and parameters of the Gaussian distribution, and generates clusters that meet the geological conditions.

[0058] Compared with the prior art, the technical solution provided by the present invention has the following beneficial effects:

[0059] (1) The present invention performs a preliminary screening of microseismic waveform windows based on a feature screening algorithm, and then further eliminates abnormal data using an abnormal window elimination algorithm, thereby automatically identifying and eliminating misidentified windows caused by interference sources such as blasting, accurately extracting microseismic event windows, and thus facilitating the subsequent precise positioning of the microseismic source.

[0060] (2) The present invention performs dual filtering in the time domain and frequency domain on the collected microseismic signal data, thereby effectively improving the quality of the microseismic signal, removing sharp noise and abnormal peaks in the spectrum, improving the smoothness of the spectrum, and improving the accuracy of feature extraction to avoid phase distortion.

[0061] (3) When performing dual filtering in the time and frequency domains on the collected microseismic signal data, the present invention adopts a sliding window mechanism to read real-time data to segment the multi-channel microseismic signal data, and optimizes the calculation formulas for the initial arrival signal-to-noise ratio and the final arrival signal-to-noise ratio, which is conducive to further improving the quality of the microseismic signal.

[0062] (4) The present invention combines the least squares method with the Geiger algorithm to locate the microseismic source. The least squares method provides a preliminary source location, while the Geiger algorithm significantly improves the positioning accuracy through multiple iterative corrections, which helps to ensure that the positioning error is minimized.

[0063] (5) Based on the positioning results of microseismic sources, the present invention innovatively adopts the combination of DBSCAN spatial clustering algorithm and Gaussian mixture model (GMM), which can effectively remove noise points, accurately identify the spatial distribution of microseismic events, and automatically eliminate interference signals, thereby ensuring the accuracy of the inversion results.

[0064] (6) In addition, the present invention significantly reduces manual intervention through automated processing procedures, improves the intelligence level of the method, and makes microseismic monitoring, positioning, and risk identification more efficient and accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 Schematic diagram of the process of microseismic source location and microseismic event inversion method according to an embodiment of the present invention;

[0066] Figure 2 Schematic diagram of the flow of the feature screening algorithm according to an embodiment of the present invention;

[0067] Figure 3 Schematic diagram of the flow of an abnormal window elimination algorithm according to an embodiment of the present invention;

[0068] Figure 4 A schematic diagram of a microseismic positioning process according to an embodiment of the present invention;

[0069] Figure 5 This is a visualization diagram of a small unit of data received by a geophone in a mine in an embodiment of the present invention;

[0070] Figure 6 This is a visualization diagram of the information of the microseismic source location obtained in an embodiment of the present invention. DETAILED DESCRIPTION

[0071] For a further understanding of the present invention, the present invention will now be described in detail with reference to the accompanying drawings and examples. However, it should be understood that these descriptions are exemplary only and are not intended to limit the scope of the present disclosure. In the following detailed description, for ease of explanation, many specific details are set forth to provide a comprehensive understanding of the embodiments of the present disclosure. However, it is apparent that one or more embodiments may also be implemented without these specific details.

[0072] Meanwhile, all terms (including technical and scientific terms) used herein have the meanings commonly understood by those skilled in the art unless otherwise defined. It should be noted that the terms used herein should be interpreted as having a meaning consistent with the context of this specification and should not be interpreted in an idealized or overly rigid manner.

[0073] In addition, the terms used herein are only for describing specific embodiments and are not intended to limit the present disclosure. For example, terms such as "comprises" and "comprising" indicate the presence of the features, steps and / or operations, but do not exclude the presence or addition of one or more other features, steps or operations.

[0074] Combine Figure 1 As shown, the method for locating a microseismic source in a mine according to an embodiment of the present invention includes:

[0075] 1. Real-time data collection and streaming reading

[0076] In a mine, geophones are used to collect microseismic data in real time, generating a large amount of high-frequency time series data. To achieve efficient data transmission and real-time processing, the present invention uses a distributed message queue system.

[0077] Specifically, each channel's geophone, acting as a data producer, categorizes collected microseismic signal data according to topics and stores it in the corresponding partitions. Each geophone publishes real-time microseismic signal data to its corresponding topic partition, effectively partitioning and organizing the data stream. Data consumers subscribe to the corresponding topics, enabling streaming access to real-time data.

[0078] The data input of this method is determined by the sampling rate of the geophone. For example, in the embodiment of the present invention, the sampling rate of the geophone is 2000 Hz. Each time, 16 channels of microseismic signal data are read. Each channel contains 8000 sampling points (corresponding to 4 seconds of data). These data are processed as a small unit in the form of a matrix of 16 columns and 8000 rows. Figure 5 shown.

[0079] 2. Microseismic Signal Data Preprocessing

[0080] In this embodiment of the present invention, a microseismic signal preprocessing module acts as a consumer, acquires data streams from different geophones through partitioning, and performs real-time preprocessing on the received data. This data preprocessing includes: reading real-time data using a sliding window mechanism to segment multi-channel microseismic signal data; extracting key features of the window using dual filtering in the time and frequency domains; and filtering and denoising the window data using an adaptive filter.

[0081] As a further preferred solution of an embodiment of the present invention, when a sliding window mechanism is used to read real-time data, the sliding window size is set to 1000 (500 ms) and the window step is set to 200 (100 ms) according to the characteristics of microseismic events to dynamically divide the data into windows.

[0082] As a further preferred solution of the embodiment of the present invention, when performing dual filtering in the time domain and frequency domain on the window data, the window data is first subjected to DC removal processing to eliminate the DC component of the signal. The formula is:

[0083]

[0084] Where x is the original signal, is the average value of the signal.

[0085] Next, a fast Fourier transform (FFT) is performed on the DC-removed signal to convert the time domain signal into a frequency domain signal to obtain its spectrum characteristics.

[0086] Furthermore, in the frequency domain filtering, median filtering and interpolation smoothing are performed to remove sharp noise and abnormal peaks in the spectrum and improve the smoothness of the spectrum; a peak fitting algorithm is used to correct low-amplitude depressions in the spectrum to optimize the spectral shape of the signal and improve the accuracy of feature extraction; then the signal main frequency, low cutoff frequency and high cutoff frequency are obtained through analysis.

[0087] Furthermore, in time domain filtering, a bandpass finite impulse response (FIR) filter is designed based on the signal main frequency, low cutoff frequency and high cutoff frequency obtained from frequency domain analysis, and the signal is bidirectionally filtered (filtfilt) to obtain a filtered signal, thereby dynamically adjusting the high and low cutoff frequencies to avoid phase distortion, adapt to the complex noise environment in the mine and retain the target signal.

[0088] Furthermore, the design of the FIR filter follows the following parameters:

[0089] Sampling frequency:

[0090] f sampling =2000Hz

[0091] Nyquist frequency:

[0092]

[0093] Since the normalized frequency is used in the filter design process, the actual frequency f low_cut 、f hign_cut Convert to normalized frequency f low 、f high , specifically:

[0094] Low Cutoff Frequency:

[0095]

[0096] High cutoff frequency:

[0097]

[0098] in:

[0099] f low_cut and f hign_cut The low cutoff frequency and high cutoff frequency are automatically extracted from the spectrum;

[0100] f low and f high are the normalized low and high cutoff frequencies used for filter design.

[0101] 3. Selection of microseismic waveform window

[0102] Based on features such as the initial arrival signal-to-noise ratio (SNR), final arrival SNR, and dominant frequency, combined with a "feature screening algorithm," windows containing microseismic events are identified for each channel. Specifically, after signal preprocessing, the initial arrival SNR and final arrival SNR of each window are calculated, and the dominant frequency of each window is determined. Multiple features from multiple channels and windows over a period of time are temporarily stored, and then the "feature screening algorithm" is used to select microseismic waveform windows.

[0103] As a preferred implementation, the calculation formula for the initial arrival signal-to-noise ratio is as follows:

[0104]

[0105] Where:

[0106] SNR initial represents the first arrival signal-to-noise ratio;

[0107] max_peak is the index of the maximum peak in the array of all peaks;

[0108] p i is the value of the i-th peak;

[0109] As a preferred implementation, the calculation formula for the end signal-to-noise ratio is as follows:

[0110]

[0111] Where:

[0112] SNR end Indicates the end signal-to-noise ratio;

[0113] N is the total number of all peaks;

[0114] max_peak is the index of the maximum peak in the array of all peaks;

[0115] p i is the value of the i-th peak.

[0116] Reference Figure 2 As shown in the figure, a feature screening algorithm is used for window selection, including:

[0117] Filter the window with the largest initial arrival signal-to-noise ratio: In each channel, select the window with the largest initial arrival signal-to-noise ratio;

[0118] Expanding the window range: Selecting several windows before and after the current window, for a total of n candidate windows. Specifically, in the embodiment of the present invention, one window is selected before the selected window, and two windows are selected after the selected window, for a total of four windows including the selected window.

[0119] Main frequency screening: In the candidate windows, the windows are screened based on the main frequency characteristics. In the embodiment of the present invention, the windows with main frequencies between 20 Hz and 150 Hz are screened in the candidate windows.

[0120] Filter the window with the largest ending signal-to-noise ratio: Among the windows that meet the main frequency conditions, select the window with the largest ending signal-to-noise ratio as the final signal window.

[0121] 4. Picking up the first arrival time and the end time

[0122] The method of distinguishing significant peaks based on envelope and threshold is to accurately determine the first arrival time and end time of microseismic signals by extracting the envelope of microseismic signals and combining peak detection and threshold setting.

[0123] Specifically, in order to accurately determine the first arrival and end time of a microseismic event, the embodiment of the present invention uses Hilbert transform to extract the envelope of the signal to enhance the significant features of the signal;

[0124] The Hilbert transform formula is:

[0125]

[0126] Where:

[0127] x(t) represents the original signal;

[0128] H(x(t)) represents the Hilbert transform of the original signal;

[0129] z(t) represents the analytical signal;

[0130] E(t) represents the envelope of the signal, which reflects the change of the instantaneous amplitude of the signal over time.

[0131] After obtaining the envelope curve, use the find_peaks function to detect significant peaks in the envelope. Then, set an appropriate threshold for the envelope signal (preferably 3% to 8% of the maximum envelope value to ensure that only strong signal peaks are extracted, and more preferably 5% of the maximum envelope value). Use a peak detection algorithm to find the first peak exceeding the threshold as the first arrival time (the arrival time of the first significant peak, indicating the time when the microseismic signal begins). Similarly, search forward from the end of the signal to find the last peak exceeding the threshold as the end time.

[0132] After determining the local first arrival and end points, we convert them to global indices. This conversion is based on the sliding window step size (step_size) and window size (window). By calculating the relative position of the signal within the window, we map the local index to the corresponding position in the global time series.

[0133] Finally, by using the global index, the corresponding timestamps can be extracted from the original time series to obtain the specific positions of the signal's initial arrival and end points on the time axis, thus preparing for the next step of the "abnormal window removal algorithm."

[0134] 5. Abnormal window removal

[0135] After picking the first arrival and end times, we initially filter out windows containing microseismic events from multiple channels within the same time range. However, due to the presence of blasting or other interfering signals, these windows may be mistakenly identified as microseismic event windows. Therefore, this embodiment of the present invention further employs an abnormal window rejection algorithm to eliminate anomalous data.

[0136] Specifically, due to the high velocity of microseismic signals, the time difference between signals received by detectors in different channels is extremely small. Based on this characteristic, the anomalous window elimination algorithm removes anomalous data by using the multi-channel first arrival time index difference. Furthermore, the algorithm incorporates a combined verification mechanism to filter the microseismic event windows initially identified by the "feature screening algorithm," effectively eliminating anomalous data (e.g., windows of interfering signals from blasting, mining, etc.).

[0137] like Figure 3 As shown, the abnormal window elimination algorithm specifically includes the following steps:

[0138] In the first step, based on the preliminarily screened multi-channel microseismic event windows, the median M of the first arrival time index and its median absolute difference (MAD) are calculated; according to the preset threshold T, the data are filtered using MAD as the criterion to obtain a filtered data set.

[0139] Among them, MAD is a robust statistic, which is the median of the absolute value of the difference between each data point and the median M. It is suitable for anomaly detection that is not affected by extreme values. Furthermore, the specific operations of this step are as follows:

[0140] Calculate the median of event_windows['global_start_index'];

[0141] Calculate the median absolute deviation (MAD) of global_start_index, which is the median of the absolute value of the difference between each data point and the median;

[0142] Based on the set threshold (threshold = 1.7*mad), points that deviate from the median by more than the threshold are removed. This approach can effectively eliminate extreme outliers.

[0143] Step 2: extracting and combining data values ​​from the filtered data set to generate all possible combinations of n data elements, where n is a positive integer of at least 4;

[0144] Step 3: Perform outlier detection on each data set based on the first arrival time index of the window. Specifically, first use the isolation forest to mark preliminary outliers, and then use a sorting-based greedy method to check whether the difference between each two data points in the group is less than a preset range L (for example, 500, exceeding the range indicates that they are not the same microseismic event). If so, return the combination as a valid combination; otherwise, continue to check other combinations;

[0145] Step 4: If at least one valid combination is found, the corresponding original data point in the data set is retained; otherwise, data processing is terminated and a prompt message indicating that there is no valid data is output.

[0146] The fifth step converts the first arrival time index in the valid combination into a timestamp and stores the timestamp and microseismic signal geophone coordinates.

[0147] 6. Microseismic Source Location

[0148] Based on the filtered valid window, the first arrival time and sensor spatial coordinates of each channel are obtained. As a preferred embodiment, the embodiment of the present invention adopts a source location method that combines the least squares method and the Geiger algorithm. First, the least squares method is used to provide the initial iteration point (x0, y0, z0, t0), and then the Geiger algorithm is used to iteratively calculate the source position.

[0149] Specific, combined Figure 4 As shown in Figure 1, during the earthquake source location process, the least squares method is first used to estimate the initial location of the earthquake source. Specifically, by using the distance relationship between the earthquake source and each detector as a constraint, the distance error between each detector and the earthquake source position is minimized, and the preliminary earthquake source location is quickly obtained. The propagation time relationship between the earthquake source and the detector is expressed as follows:

[0150]

[0151] In the formula:

[0152] (x s ,y s ,z s ) is the hypothetical earthquake source location in the mine;

[0153] (x i ,y i ,z i ) is the actual position of the multi-channel detector;

[0154] t iis the propagation time of the microseismic signal received by each geophone;

[0155] t0 is the initial propagation time of the earthquake source;

[0156] v represents the wave velocity;

[0157] The expression of the objective function is as follows:

[0158]

[0159] Among them, (x s ,y s ,z s ) is the source location, t0 is the initial propagation time of the source, N is the number of detectors, and v is the wave velocity.

[0160] After obtaining the preliminary source position (x0, y0, z0, t0) and propagation time t0 through the least squares method, the Geiger algorithm is used to iteratively optimize the source position. By calculating the time error of each sensor, the source position and propagation time are adjusted using a linearization method. The specific calculation is achieved through the following linear equation system:

[0161] AΔp=b

[0162] In the formula:

[0163] A is the Jacobian matrix, which contains the partial derivatives of n sensors;

[0164]

[0165] Δp is the correction vector of the earthquake source position and travel time;

[0166] b is the time error of each sensor;

[0167] By Gaussian elimination or other linear equation solving methods.

[0168] The sum of squared time errors across all sensors is calculated to determine whether the error threshold is met. If the error is less than the threshold, the positioning is considered converged, and the final earthquake source location is the desired one.

[0169] Furthermore, an iterative optimization is performed based on the Geiger algorithm. The specific process is as follows: based on the preliminary source position (x0, y0, z0, t0) and propagation time t0 obtained by the least squares method, the time error of each sensor is calculated using the first-order Taylor expansion. Assuming that the current test point is (x s ,y s ,z s ,t0), calculate the predicted arrival time to the i-th sensor:

[0170]

[0171] Time error Δt i for:

[0172]

[0173] The corrections of the earthquake source position and propagation time are calculated by Taylor expansion. Assume that the correction vectors of the earthquake source position and propagation time are (Δx s ,Δy s ,Δz s ,Δt0), can be calculated by the following linearized equation:

[0174]

[0175] Substituting the above formula for all sensors i yields a system of equations. Using linear equation solving, we determine the correction vector Δp, which is then added to the current source location and propagation time to update the test point.

[0176]

[0177] Calculate whether the updated earthquake source position and propagation time meet the set error threshold. That is, if the error If it is less than a preset threshold, it is considered that the positioning has converged and the iteration is stopped. s ,y s ,z s ,t0) is the desired earthquake source position.

[0178] Figure 6 The figure shows an information visualization diagram of the microseismic source location. Based on the positioning results, the coordinates are projected onto the corresponding stratum. The energy range of the microseismic event is distinguished by small balls of different colors. The greater the energy, the larger the radius of the small ball.

[0179] An embodiment of the present invention further provides a device for locating a microseismic source in a mine, comprising:

[0180] The data acquisition and streaming reading module is used to collect microseismic signal data in the mine in real time, and realizes data streaming and real-time processing through distributed stream processing and message queues;

[0181] A microseismic signal data preprocessing module is used to preprocess the collected microseismic signal data, wherein the preprocessing includes dual filtering processing in the time domain and frequency domain;

[0182] The microseismic window selection module is used to preliminarily screen the microseismic waveform window based on the initial arrival signal-to-noise ratio, the final arrival signal-to-noise ratio, and the main frequency characteristics, combined with the feature screening algorithm;

[0183] First arrival time and end time picking module,

[0184] An abnormal window rejection module is used to extract the envelope of the microseismic signal using Hilbert transform and, in combination with the threshold detection method, to pick up the first arrival time and the end time of the microseismic signal; and

[0185] The microseismic source location module is used to locate the microseismic source by combining the least squares method and the Geiger algorithm.

[0186] A third aspect of the present invention provides a microseismic event inversion method, comprising:

[0187] Using any of the methods described in the first aspect of the present invention to locate the microseismic source in a mine;

[0188] Based on the coordinates of microseismic events, the DBSCAN algorithm is used to identify areas with higher density to eliminate noise points;

[0189] The Gaussian mixture model is used to perform cluster analysis on microseismic events based on Gaussian distribution;

[0190] For each identified microseismic event cluster, RANSAC is applied to identify fracture surfaces and faults.

[0191] The embodiment of the present invention further provides a microseismic event inversion method, comprising the following steps:

[0192] Step 1: spatial clustering and noise point removal.

[0193] First, the calculated microseismic event data (the coordinates of the microseismic events) are processed to remove noise and irrelevant data;

[0194] DBSCAN (density-based spatial clustering algorithm) is then used to identify clusters of microseismic events. DBSCAN identifies clusters of microseismic events by searching for high-density areas in space. It automatically groups closely spaced microseismic events into a single cluster without requiring parameter adjustment. Events that lack sufficient neighborhood density with events in any cluster are considered noise points and are removed.

[0195] Step 2: After removing the noise points, the Gaussian mixture model (GMM) is further used to perform cluster analysis on microseismic events based on Gaussian distribution.

[0196] The algorithm sets the initial number and parameters of a Gaussian mixture model (GMM) based on the geometric characteristics of underground mine tunnels and coal seams. It prioritizes clustering the roof and floor regions of the tunnels and trains the model using an expectation-maximization (EM) algorithm. The trained GMM divides the data into multiple clusters. The algorithm then dynamically adjusts the number and parameters of Gaussian distributions based on the distribution patterns of faults and cracks in the working face of the mine to generate clusters that meet the geological conditions.

[0197] The algorithm combines the known three-dimensional geometric model of the mine tunnel and coal seam, no longer randomly initializes the GMM, and preferentially selects the microseismic coordinates of the tunnel roof, floor and faults as the initial Gaussian distribution center. The number of Gaussian distributions is preferably set to 5 initially; each Gaussian center point μ k Using the tunnel roof and floor as a guide, the initial covariance matrix was set to a diagonal matrix with a standard deviation of 3 meters. The expectation-maximization (EM) algorithm was used to iteratively train all microseismic event coordinates. The trained GMM partitioned the data into multiple clusters. The splitting or merging of Gaussian distributions was determined based on the characteristics of dense areas near the roof, floor, and faults, as well as the covariance shape. The number and parameters of Gaussian distributions were dynamically adjusted to generate three clusters that met the geological conditions.

[0198] Step 3: For each identified microseismic event cluster, RANSAC (random consensus algorithm) is applied to identify fracture surfaces and faults.

[0199] The specific process of this step is as follows: a subset of data points within the cluster is randomly selected, candidate planes are fitted through linear regression, and the fitting error is calculated; the algorithm iteratively optimizes the plane's normal vector and intercept parameters to maximize the number of inliers (i.e., data points with an error below a set threshold); a plane is fitted to the data of each cluster to extract clustered areas that meet the geometric characteristics; finally, the fitted plane parameters are used for three-dimensional visualization to show the spatial distribution of fracture surfaces and faults.

[0200] In summary, the present invention combines time domain and frequency domain automatic filtering technology for signal preprocessing to extract the key features of the signal window; then, through the innovative "feature screening algorithm" and "abnormal window elimination algorithm", the microseismic event window is extracted; in terms of source location, the least squares method is combined with the Geiger algorithm to improve the source location accuracy and computational efficiency; spatial clustering analysis is performed through DBSCAN and GMM algorithms to eliminate noise points and identify the microseismic event clustering area; the RANSAC algorithm is combined to perform plane fitting on the data points in the clustered area to extract the fracture surface and fault structure information, thereby realizing the inversion of microseismic events.

Claims

1. A method for locating microseismic sources in a mine, characterized in that: include: Real-time collection of microseismic signal data in mines, and the use of distributed stream processing and message queues to achieve data streaming and real-time processing; Preprocessing the collected microseismic signal data, wherein the preprocessing includes dual filtering processing in the time domain and frequency domain; Based on the initial arrival signal-to-noise ratio, final arrival signal-to-noise ratio and main frequency characteristics, and combined with the feature screening algorithm, the microseismic waveform window is preliminarily screened; The envelope of the microseismic signal is extracted using Hilbert transform, and the first arrival time and end time of the microseismic signal are picked up by combining the threshold detection method. The abnormal window elimination algorithm is used to further eliminate abnormal data; The microseismic source is located by combining the least squares method and the Geiger algorithm. First, the least squares method is used to provide the initial iteration point, and then the Geiger algorithm is used to iteratively calculate the source position.

2. The method for locating microseismic sources in mines according to claim 1, characterized in that: Geophones are used to collect microseismic signal data in the mine in real time. The microseismic signal data collected by each geophone in real time are classified according to themes and published to the corresponding theme partitions to achieve effective division and organization of data streams.

3. The method for locating microseismic sources in mines according to claim 1, characterized in that: The time domain and frequency domain dual filtering processing specifically includes: Perform DC removal processing on the data to eliminate the DC component of the signal; Perform fast Fourier transform on the DC-removed signal to convert the time domain signal into a frequency domain signal; In the frequency domain filtering, median filtering and interpolation smoothing are performed, and the peak fitting algorithm is used to correct the low-amplitude depressions in the spectrum; then the signal main frequency, low cutoff frequency and high cutoff frequency are obtained through analysis; In time domain filtering, the signal is bidirectionally filtered through an FIR filter according to the signal main frequency, low cutoff frequency and high cutoff frequency obtained by frequency domain analysis to obtain a filtered signal.

4. The method for locating microseismic sources in mines according to claim 3, characterized in that: When preprocessing the collected microseismic signal data, a sliding window mechanism is used to read real-time data to segment the multi-channel microseismic signal data, and then the window data is subjected to dual filtering in the time domain and frequency domain.

5. The method for locating the microseismic source in a mine according to any one of claims 1 to 4, characterized in that: The microseismic waveform window is initially screened based on the initial arrival signal-to-noise ratio, the final arrival signal-to-noise ratio, and the main frequency characteristics in combination with the feature screening algorithm, specifically including: After signal preprocessing, the first arrival signal-to-noise ratio of each window is calculated, and the window with the largest first arrival signal-to-noise ratio in each channel is selected; Taking the window with the largest initial arrival signal-to-noise ratio as the benchmark, select several windows before and after the window, a total of n candidate windows, to expand the window range; In the above candidate windows, the windows are screened according to the dominant frequency characteristics; Among the windows that meet the main frequency condition, the window with the largest ending signal-to-noise ratio is selected as the final signal window.

6. The method for locating microseismic sources in mines according to claim 5, characterized in that: The calculation formulas for the initial signal-to-noise ratio and the final signal-to-noise ratio are as follows: Where: SNR initial represents the first arrival signal-to-noise ratio; max_peak is the index of the maximum peak in the array of all peaks; p i is the value of the i-th peak; SNR end Indicates the end signal-to-noise ratio; N is the total number of all peaks.

7. The method for locating the microseismic source in a mine according to any one of claims 1 to 4, characterized in that: The abnormal window elimination algorithm is used to further eliminate abnormal data, specifically including: The first arrival time and the end time of the microseismic signal are converted into a global index to map the local index to the corresponding position of the global time series; Based on the preliminarily screened multi-channel microseismic event windows, the median M of the first arrival time index and its median absolute difference MAD are calculated; according to the preset threshold T, the data are filtered using MAD as the criterion to obtain a filtered data set; Extracting and combining data values ​​from the filtered data set to generate all possible combinations of n data elements, where n is a positive integer of at least 4; Outlier detection is performed on each data set based on the first arrival time index of the window. The initial outliers are first marked using the isolation forest algorithm. Then, a greedy method based on sorting is used to check whether the difference between every two data points in the group is less than the preset range L. If so, the combination is returned as a valid combination; otherwise, other combinations are checked. If at least one valid combination is found, the corresponding original data point in the data set is retained; otherwise, data processing is terminated and a prompt message indicating that there is no valid data is output; Convert the first arrival time index in the valid combination into a timestamp and store the timestamp and microseismic signal geophone coordinates.

8. A microseismic source locating device in a mine, characterized in that: include: The data acquisition and streaming reading module is used to collect microseismic signal data in the mine in real time, and realizes data streaming and real-time processing through distributed stream processing and message queues; A microseismic signal data preprocessing module is used to preprocess the collected microseismic signal data, wherein the preprocessing includes dual filtering processing in the time domain and frequency domain; The microseismic window selection module is used to preliminarily screen the microseismic waveform window based on the initial arrival signal-to-noise ratio, the final arrival signal-to-noise ratio, and the main frequency characteristics, combined with the feature screening algorithm; First arrival time and end time picking module, The abnormal window rejection module is used to extract the envelope of the microseismic signal using Hilbert transform and, combined with the threshold detection method, to pick up the first arrival time and end time of the microseismic signal; as well as The microseismic source location module is used to locate the microseismic source by combining the least squares method and the Geiger algorithm.

9. A microseismic event inversion method, characterized in that: include: Using the method according to any one of claims 1 to 8 to locate the microseismic source in a mine; Based on the coordinates of microseismic events, the DBSCAN algorithm is used to identify areas with higher density to eliminate noise points; The Gaussian mixture model is used to perform cluster analysis on microseismic events based on Gaussian distribution; For each identified microseismic event cluster, RANSAC is applied to identify fracture surfaces and faults.

10. The microseismic event inversion method according to claim 9, characterized in that: The Gaussian mixture model is used to perform cluster analysis on microseismic events based on Gaussian distribution. Specifically, the following steps are performed: according to the geometric characteristics of underground mine tunnels and coal seams, the initial distribution quantity and parameters of the Gaussian mixture model are set, the roof and floor areas of the tunnels are preferentially clustered, and training is performed using an expectation maximization algorithm; after training, the GMM divides the data into multiple clusters, and dynamically adjusts the number and parameters of the Gaussian distribution based on the distribution patterns of faults and cracks in the working face of the mine to generate clusters that meet geological conditions.

Images