Double-stage seismic source positioning method, device and equipment for coal mine micro-seismic event

Abstract

The invention discloses a double-stage seismic source positioning method, device and equipment for a coal mine micro-seismic event, and relates to the field of seismic source positioning, and the method comprises the steps: determining a first seismic source parameter estimation value of the coal mine micro-seismic event through a moth flame algorithm according to the arrival time data and seismic wave propagation speed of each vibration monitoring device, the first seismic source parameter estimation value comprises a first seismic source estimation value and a first seismic moment estimation value; the first seismic source parameter estimation value serves as the initial state of a Markov chain Monte Carlo method, seismic source parameter samples conforming to posterior distribution are generated through the Markov chain Monte Carlo method, and each seismic source parameter sample comprises a second seismic source estimation value and a second seismic moment estimation value; the joint probability density of each seismic source parameter sample is calculated, the seismic source parameter sample with the maximum joint probability density is selected, the final seismic source parameter estimation value of the coal mine microseismic event is obtained, and positioning efficiency is improved while positioning precision is improved.

Description

Technical Field

[0001] The present application relates to the field of earthquake source location technology, and in particular to a dual-stage earthquake source location method, device and equipment for coal mine microseismic events. Background Art

[0002] Coal mine microseismic positioning is an important tool for preventing dynamic disasters. By quickly and accurately locating the location of the vibration event, that is, the earthquake source, it helps to improve the accuracy and efficiency of mine safety monitoring.

[0003] However, in terms of data monitoring, although a large number of seismic stations theoretically provide more data support for the source location of microseismic events, due to the cumulative time error of each seismic station, the increase in the number of seismic stations does not linearly improve the positioning accuracy. Instead, it will greatly increase the time consumed in the inversion process, resulting in large errors in the source location accuracy of coal mine microseismic events due to the actual distribution of seismic stations. In terms of positioning algorithms, traditional source location methods such as those based on grid search and linear inversion methods have large computational complexity and are prone to falling into local optimal solutions; while source location methods based on genetic algorithms and particle swarm optimization algorithms can globally search for optimal solutions, but their convergence speed is insufficient in high-dimensional parameter spaces (such as three-dimensional spatial coordinates + time). In summary, existing source location methods for coal mine microseismic events are difficult to simultaneously meet the positioning requirements of high efficiency and high precision. Summary of the Invention

[0004] The purpose of this application is to provide a dual-stage source location method, device and equipment for coal mine microseismic events, so as to improve the location efficiency and location accuracy of the source location method for coal mine microseismic events.

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

[0006] In a first aspect, the present application provides a dual-stage source location method for microseismic events in coal mines, comprising:

[0007] Determine, based on the arrival time data and seismic wave propagation velocity observed by each vibration monitoring device, a first source parameter estimate of the coal mine microseismic event using a moth flame algorithm; wherein the first source parameter estimate includes a first source estimate and a first earthquake occurrence time estimate;

[0008] Using the first earthquake source parameter estimate as the initial state of a Markov chain Monte Carlo method, generating earthquake source parameter samples that conform to a posterior distribution through the Markov chain Monte Carlo method, wherein each of the earthquake source parameter samples includes a second earthquake source estimate and a second earthquake occurrence time estimate;

[0009] The joint probability density of each source parameter sample is calculated, and the source parameter sample with the largest joint probability density is selected to obtain the final source parameter estimation value of the coal mine microseismic event.

[0010] Optionally, determining the estimated value of the first source parameter of the coal mine microseismic event by a moth flame algorithm based on the arrival time data and seismic wave propagation velocity observed by each vibration monitoring device specifically includes:

[0011] Perform population initialization operations:

[0012] Initialize the number of moths in the population, the total number of iterations, and randomly initialize the initial position of each moth according to the preset source monitoring range in, is the jth moth m j The initial shock moment, It's a moth j The initial three-dimensional coordinates of , 1≤j≤M, M is the number of moths;

[0013] Perform fitness value evaluation operation:

[0014] Based on the current position of each moth Calculate the current fitness value of each moth through the objective function in, is the number of moths m after the tth iteration j Location is the number of moths m after the tth iteration j The time of earthquake, 0≤t≤A, A is the total number of iterations, is the number of moths m after the tth iteration j The three-dimensional coordinates of

[0015] To perform a flame selection operation:

[0016] Sort the fitness values ​​of the current positions of all moths in ascending order, and select the top N in ascending order f The current position of the moth corresponding to the (t) fitness value is taken as the flame, N f The calculation formula for (t) is:

[0017]

[0018] Among them, g() is the rounding function, and t is the current number of iterations;

[0019] Iteratively executing the steps of position updating, fitness value evaluation, and flame selection in sequence until the deviation of the flame value obtained by multiple consecutive iterations meets the preset requirement or the current iteration number t reaches the total iteration number, and the obtained flame is the first source parameter estimation value;

[0020] The location update operation includes:

[0021] According to the current position of each moth and the corresponding flame, the Euclidean distance between each moth and the corresponding flame is calculated, and the position of each moth is updated according to the Euclidean distance to obtain the updated position of each moth.

[0022] Optionally, the objective function is:

[0023]

[0024]

[0025] Among them, t pred,i is the currently predicted arrival time of the i-th vibration monitoring device, tobs,i is the observation time of the i-th vibration monitoring device, N is the total number of vibration monitoring devices, is the number of moths m after the tth iteration j The moment of earthquake, s i is the three-dimensional coordinate of the i-th vibration monitoring device, v is the velocity model constructed according to the three-dimensional coordinate Determine the speed at which seismic waves propagate.

[0026] Optionally, the calculating the Euclidean distance between each moth and the corresponding flame according to the current position of each moth and the corresponding flame, and updating the position of each moth according to the Euclidean distance to obtain the updated position of each moth specifically includes:

[0027] According to the following formula (4), the Euclidean distance between each moth and the corresponding flame is calculated, and the position of each moth is updated according to the Euclidean distance to obtain the updated position of each moth:

[0028]

[0029] in, is the moth m after the t+1th iteration j Location f k is the number of moths m after the tth iteration j The corresponding k-th flame after ascending sorting, is the number of moths m after the tth iteration j With flames f kThe Euclidean distance is b, the spiral shape control parameter is b, the random parameter is θ∈[a,1], and the attenuation factor is a=-1-t / A.

[0030] Optionally, the Markov Chain Monte Carlo method includes an MH algorithm.

[0031] Optionally, the step of using the first source parameter estimate as an initial state of a Markov chain Monte Carlo method and generating source parameter samples that conform to a posterior distribution by the Markov chain Monte Carlo method specifically includes:

[0032] Based on the current state, use the proposal distribution to generate candidate points:

[0033]

[0034] Among them, x' represents the candidate point, x t′ represents the state of the t′th iteration, t′≥0, x0 is the first source parameter estimate, σp represents the spatial coordinate standard deviation, and diag represents the diagonal covariance matrix;

[0035] For each candidate point x', calculate the current acceptance probability:

[0036]

[0037] Among them, α t represents the acceptance probability of the t′th iteration; L(t obs |x′) represents a given candidate point x', the observed data t obs,i Likelihood value; L(t obs |x t′ ) means given the current state x t′ , observation data t obs,i Likelihood value; p(x′) represents the prior probability of candidate point x', p(x t′ ) represents the current state x t′ The prior probability of q(x t′ |x′) represents the proposal from candidate point x' to the current state x t′ The proposed distribution probability, q(x′|x t′ ) indicates that from the current state x t′ Proposal distribution probability of candidate point x'; where L(t obs |x′) and L(t obs |x t′ )The likelihood function L(t obs )for:

[0038]

[0039] Where σ is the noise standard deviation, t0 is the candidate point x' or the current state x t′ The moment of earthquake, [x, y, z] is the candidate point x' or the current state x t′ The three-dimensional coordinates of [x, y, z] are obtained by the velocity model, and v′ is the seismic wave propagation velocity determined by the three-dimensional coordinates [x, y, z].

[0040] If α t′ ≥Uniform(0,1), let the state x of the t′+1th iteration be t′+1 = x', that is, accept the candidate point x'; otherwise, let x t′+1 =x t′ , that is, keep the current state x t′ ; Among them, Uniform(0,1) is a uniform random number in the interval [0,1];

[0041] The above steps are iterated until the target number of source parameter samples that conform to the posterior distribution is obtained.

[0042] Optionally, the proposed distribution adopts an isotropic Gaussian distribution;

[0043] The method of using the first source parameter estimate as the initial state of the Markov chain Monte Carlo method and generating source parameter samples that conform to the posterior distribution by the Markov chain Monte Carlo method further includes:

[0044] Every Q iterations, the proposal distribution is updated according to the following formula:

[0045] ∑ t′ =Cov(X t′-150:t′ )+10 -6 I; (8)

[0046] q(x t′ |x′)=N(x t′ ,∑ t′ ); (9)

[0047] Among them, X t′-150:t′ is the historical sample matrix, ∑ t′ is a Gaussian distribution N(x t ,∑ t )’s covariance matrix, I represents the identity matrix;

[0048] The method of using the first source parameter estimate as the initial state of the Markov chain Monte Carlo method and generating source parameter samples that conform to the posterior distribution by the Markov chain Monte Carlo method further includes:

[0049] Every Q iterations, the isotropic proposal scale of the isotropic Gaussian distribution is adjusted according to the following formula:

[0050] logσ p,t′+1 =logσ p,t′ +γ t′ (α t′ -P); (10)

[0051] Among them, σ p,t′+1 is the isotropic proposal scale for the t′+1th iteration, σ p,t′ is the isotropic proposal scale for the t′th iteration, γ t′ is the decay learning rate, and P is the optimal acceptance rate.

[0052] Optionally, the calculating the joint probability density of each source parameter sample specifically includes:

[0053] According to the following formulas (12) to (15), the joint probability density of each source parameter sample is calculated:

[0054]

[0055] in, is the joint probability density of a given source parameter sample c, c∈{c i}, c i is the i-th source parameter sample, 1≤i≤n, n is the total number of source parameter samples, H is the bandwidth matrix, d is the parameter space dimension, is the covariance matrix, is the sample mean.

[0056] In a second aspect, the present application provides a dual-stage source location device for coal mine microseismic events, comprising:

[0057] Vibration monitoring device, used to monitor and record the arrival data of microseismic events in coal mines;

[0058] a first optimization module configured to determine, based on the arrival time data and seismic wave propagation velocity observed by each vibration monitoring device, a first source parameter estimate of the vibration event using a moth flame algorithm; wherein the first source parameter estimate includes a first source estimate and a first earthquake onset time estimate;

[0059] A second optimization module is configured to use the first earthquake source parameter estimate as the initial state of a Markov chain Monte Carlo method, and to generate earthquake source parameter samples that conform to a posterior distribution through the Markov chain Monte Carlo method, wherein each of the earthquake source parameter samples includes a second earthquake source estimate and a second earthquake onset time estimate;

[0060] The kernel density estimation module is used to calculate the joint probability density of each source parameter sample and select the source parameter sample with the largest joint probability density to obtain the final source parameter estimation value of the coal mine microseismic event.

[0061] In a third aspect, the present application provides a computer device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the dual-stage source location method for coal mine microseismic events described in any one of the above.

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

[0063] This application provides a dual-stage source location method, device and equipment for coal mine microseismic events. Due to the advantages of high search efficiency and low number of iterations, the moth flame algorithm Therefore, by determining the first source parameter estimate (including the first source estimate and the first earthquake occurrence time estimate) of the coal mine microseismic event based on the arrival data of each vibration monitoring device and the constructed velocity model through the moth flame algorithm, a global search for the first source estimate and the first earthquake occurrence time estimate is realized in the preset search space, with high positioning efficiency and a small number of iterations, and the ability to quickly obtain the initial estimate of the source position and the earthquake occurrence time; by using the first source parameter estimate as the initial state of the Markov chain Monte Carlo method, the initial positioning result of the moth flame algorithm is innovatively used as the high-quality initial value of the Markov chain Monte Carlo method, which reduces the iterations of the Markov chain Monte Carlo method, and generating source parameter samples (including the second source estimate and the second earthquake occurrence time estimate) that conform to the posterior distribution through the Markov chain Monte Carlo method, calculating the joint probability density of each source parameter sample, and selecting the source parameter sample with the largest joint probability density to obtain the final source parameter estimate of the coal mine microseismic event, thereby realizing the dual-stage source positioning of the coal mine microseismic event. Since the Markov chain Monte Carlo method has the ability of local fine optimization, the moth flame algorithm is combined with the Markov chain Monte Carlo method to simultaneously locate the source of coal mine microseismic events from both the global and local levels. Compared with the traditional source location method of coal mine microseismic events that uses a single algorithm, the positioning accuracy is effectively improved; in summary, the present application improves the positioning efficiency while improving the positioning accuracy, and solves the problem that the existing source location method of coal mine microseismic events is difficult to simultaneously meet the positioning requirements of high efficiency and high precision. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0065] Figure 1A schematic flow chart of a dual-stage source location method for a microseismic event in a coal mine provided in one embodiment of the present application;

[0066] Figure 2 A schematic diagram of a vibration monitoring device and monitoring range provided in one embodiment of the present application;

[0067] Figure 3 A schematic diagram of a velocity model provided in one embodiment of the present application;

[0068] Figure 4 A schematic diagram of the positioning effect of a dual-stage source positioning method for a microseismic event in a coal mine provided by another embodiment of the present application;

[0069] Figure 5 A schematic diagram of the positioning effect of a dual-stage source positioning method for a microseismic event in a coal mine provided by another embodiment of the present application;

[0070] Figure 6 A schematic diagram of the positioning effect of a dual-stage source positioning method for a microseismic event in a coal mine provided by another embodiment of the present application;

[0071] Figure 7 A schematic diagram of the positioning effect of a dual-stage source positioning method for a microseismic event in a coal mine provided by another embodiment of the present application;

[0072] Figure 8 A schematic diagram of the positioning effect of a dual-stage source positioning method for a microseismic event in a coal mine provided by another embodiment of the present application;

[0073] Figure 9 A schematic diagram of the distribution and velocity model of the Dongtan Coal Mine vibration monitoring device provided in one embodiment of the present application;

[0074] Figure 10 A schematic diagram of the positioning effect of a dual-stage source location method for coal mine microseismic events provided by another embodiment of the present application on the Dongtan coal mine blasting event;

[0075] Figure 11 A schematic diagram of the functional modules of a dual-stage source locating device for microseismic events in coal mines provided in one embodiment of the present application;

[0076] Figure 12 A schematic diagram of the structure of a computer device provided in one embodiment of the present application. DETAILED DESCRIPTION

[0077] The following will provide a clear and complete description of the technical solutions in the embodiments of this application, in conjunction with the accompanying drawings. It should be understood that the described embodiments are merely a portion of the embodiments of this application, and not all of them. All other embodiments derived by persons of ordinary skill in the art based on the embodiments of this application without inventive effort are intended to fall within the scope of protection of this application.

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

[0079] In an exemplary embodiment, Figure 1 As shown, a dual-stage earthquake source location method for a microseismic event in a coal mine is provided. The method is executed by a computer device, specifically a computer device such as a terminal or a server, and includes the following steps 101 to 103.

[0080] Step 101 : Determine the first source parameter estimate of the coal mine microseismic event using the moth flame algorithm based on the arrival time data of each vibration monitoring device and the constructed velocity model; wherein the first source parameter estimate includes the first source estimate and the first earthquake occurrence time estimate.

[0081] In the embodiments of the present application, the velocity model is a uniform layered medium, and is not specifically limited thereto. An appropriate velocity model can be selected based on actual needs. The arrival time data of the vibration monitoring device refers to the time when the seismic wave reaches each vibration monitoring device. There are no specific limitations on the vibration monitoring device, and it can be selected based on actual needs. For example, a seismic station can be used as the vibration monitoring device. Step 101 achieves the first stage of source location, i.e., coarse location, of the coal mine microseismic event.

[0082] Step 102: Use the first earthquake source parameter estimate as the initial state of the Markov chain Monte Carlo method, and generate earthquake source parameter samples that conform to the posterior distribution through the Markov chain Monte Carlo method, wherein each earthquake source parameter sample includes a second earthquake source estimate and a second earthquake occurrence time estimate.

[0083] Step 103 , calculating the joint probability density of each source parameter sample, and selecting the source parameter sample with the largest joint probability density to obtain the final source parameter estimation value of the coal mine microseismic event.

[0084] In the embodiment of the present application, the second stage of source positioning of the coal mine microseismic event, namely, fine positioning, is achieved through steps 102 to 103.

[0085] In implementing the above steps 101 to 103, since the moth flame algorithm has the advantages of high search efficiency and few iterations, the first source parameter estimation value (including the first source estimation value and the first earthquake occurrence time estimation value) of the coal mine microseismic event is determined by the moth flame algorithm based on the arrival time data of each vibration monitoring device and the constructed velocity model, and the first source estimation value and the first earthquake occurrence time estimation value are globally searched in the preset search space, with high positioning efficiency and few iterations. It can quickly obtain the initial estimated values ​​of the earthquake source position and the time of earthquake occurrence; by taking the first earthquake source parameter estimate as the initial state of the Markov chain Monte Carlo method, the initial positioning result of the moth flame algorithm is innovatively used as the high-quality initial value of the Markov chain Monte Carlo method, so that the iterations of the Markov chain Monte Carlo method are reduced, and the earthquake source parameter samples (including the second earthquake source estimate and the second earthquake occurrence time estimate) that conform to the posterior distribution are generated by the Markov chain Monte Carlo method, the joint probability density of each earthquake source parameter sample is calculated, and the earthquake source parameter sample with the largest joint probability density is selected to obtain the final earthquake source parameter estimate of the coal mine microseismic event, thereby realizing the two-stage earthquake source positioning of the coal mine microseismic event. Since the Markov chain Monte Carlo method has the ability of local fine optimization, the moth flame algorithm is combined with the Markov chain Monte Carlo method, The source location of coal mine microseismic events is carried out at both the global and local levels, which effectively improves the positioning accuracy compared with the traditional source location method of coal mine microseismic events that uses a single algorithm; in summary, the present application improves the positioning efficiency while improving the positioning accuracy, solves the problem that the existing source location method of coal mine microseismic events is difficult to meet the requirements of high-efficiency and high-precision positioning at the same time, and helps to more accurately predict potential dangerous areas in coal mine shafts.

[0086] In another exemplary embodiment of the present application, the above step 101 includes the following steps 201 to 204. Among them:

[0087] Step 201: Perform population initialization operation:

[0088] Initialize the number of moths in the population, the total number of iterations, and randomly initialize the initial position of each moth according to the preset source monitoring range in, is the jth moth m j The initial shock moment, It's a moth j The initial three-dimensional coordinates, 1≤j≤M, M is the number of moths, For moths j The initial X-axis coordinate of For moths j The initial Y-axis coordinate of For moths j The initial Z-axis coordinate of .

[0089] In the embodiment of the present application, when randomly initializing the initial position of each moth, it is ensured that the initial position of the moth covers the possible location of the real earthquake source as much as possible. The preset monitoring range can be set according to actual needs. For example, the earthquake source monitoring range is set to X∈[0,5]km, Y∈[0,5]km, T∈[0,5]s, where X is the value range in the X-axis direction, Y is the value range in the Y-axis direction, Z is the value range in the Z-axis direction, and T is the value range at the time of earthquake, which further improves the accuracy and efficiency of earthquake source positioning.

[0090] Step 202: Execute the fitness value evaluation operation:

[0091] Based on the current position of each moth Calculate the current fitness value of each moth through the objective function in, is the number of moths m after the tth iteration j Location is the number of moths m after the tth iteration j The time of earthquake, 0≤t≤A, A is the total number of iterations, is the number of moths m after the tth iteration j The three-dimensional coordinates of is the number of moths m after the tth iteration j The X-axis coordinate, is the number of moths m after the tth iteration j The Y-axis coordinate, is the number of moths m after the tth iteration j The Z-axis coordinate of .

[0092] In the embodiment of the present application, when the current number of iterations is zero, the current position of each moth is its initial position. When the current number of iterations is greater than zero, the current position is the position updated in each iteration.

[0093] Step 203, perform flame selection operation:

[0094] Sort the fitness values ​​of the current positions of all moths in ascending order, and select the top N in ascending order f The current position of the moth corresponding to the (t) fitness value is taken as the flame (elite solution set), N f The calculation formula for (t) is:

[0095]

[0096] Among them, g() is the rounding function, and t is the current number of iterations.

[0097] In the embodiments of this application, the rounding function, number of moths, and total number of iterations are not specifically limited and can be selected based on actual needs. For example, the rounding function can be set to a rounding function, a rounding-up function, or a rounding-down function. Preferably, the rounding function is a rounding function, the number of moths is 20, and the total number of iterations is 300, further improving the accuracy and efficiency of earthquake source location.

[0098] Flame selection aims to optimize search direction by constructing an elite solution set through a dynamic decay mechanism. This mechanism retains more flames in the early stages of iteration, promoting global exploration through a wide range of elite solutions. In the later stages of iteration, the number of flames is reduced to focus on high-quality solutions and accelerate convergence. This ensures that the moth swarm consistently spirals around the elite solution set, achieving a dynamic balance between exploration and development.

[0099] Step 204, the position update operation, fitness value evaluation operation and flame selection operation are executed in a loop iterative manner until the difference of the flame values ​​obtained by multiple consecutive loop iterations meets the preset requirements or the current iteration number t reaches the total iteration number A, and the obtained flame is the first source parameter estimation value.

[0100] In the embodiment of the present application, the deviation of the flame value obtained through multiple consecutive cyclic iterations satisfies the preset requirement, which means that among the flame values ​​obtained through multiple consecutive cyclic iterations, the absolute value of the difference between the flame value obtained through the first cyclic iteration and the flame value obtained through the last cyclic iteration satisfies the preset requirement, or that among the flame values ​​obtained through multiple consecutive cyclic iterations, the absolute value of the difference between the flame values ​​obtained through any two adjacent cyclic iterations satisfies the preset requirement.

[0101] In this embodiment of the present application, the location update operation includes:

[0102] Step 301 : Calculate the Euclidean distance between each moth and the corresponding flame according to the current position of each moth and the corresponding flame, and update the position of each moth according to the Euclidean distance between each moth and the corresponding flame to obtain the updated position of each moth.

[0103] In the embodiment of the present application, the N selected after the ascending sorting is cyclically assigned to all moths. f (t) flames. If the current number of iterations t is small enough to make N f The value of (t) is M, and the M flames are assigned to the M moths one by one. As the number of iterations t increases, N f When the value of (t) is less than M, follow N f (t) flames are sequenced, and N flames are used in a loop from front to back. f (t) flames are assigned to each of the M moths. For example, if the value of M is 20, N fIf the value of (t) is 15, then first assign the 15 flames one-to-one to any 15 moths, and then assign the first 5 of the 15 flames one-to-one to the remaining 5 moths.

[0104] In another exemplary embodiment of the present application, the objective function is:

[0105]

[0106] Among them, t pred,i is the currently predicted arrival time of the i-th vibration monitoring device, t obs,i is the observed time of the i-th vibration monitoring device. During the test phase, noise (Gaussian noise) can be artificially added to the observed time of the i-th vibration monitoring device to simulate the real observed time; N is the total number of vibration monitoring devices, is the number of moths m after the tth iteration j The moment of earthquake, s i is the three-dimensional coordinate of the i-th vibration monitoring device, v is the velocity model constructed according to the three-dimensional coordinate Determine the speed at which seismic waves propagate.

[0107] In the embodiment of the present application, the number and layout of the vibration monitoring devices can be set according to actual needs. Preferably, 8 vibration monitoring devices are used. Figure 2 The layout is shown in Table 1. The three-dimensional coordinates of the eight vibration monitoring devices are distributed in Figure 2 The monitoring range (cube vertex range) is within the monitoring range to ensure effective monitoring of the monitoring range. The velocity model can be set according to actual needs. For example, using Figure 3 The velocity model of uniform layered medium shown in the figure is bounded by a depth of 2.5 km. If z j Not more than 2.5 km, v = 3.2 km / s, if z j Greater than 2.5km, v=3.9km / s.

[0108] Table 1 Example of three-dimensional coordinate data of vibration monitoring device

[0109]

[0110] In another exemplary embodiment of the present application, the above step 301 includes:

[0111] Step 401: Calculate the Euclidean distance between each moth and the corresponding flame according to the following formula (4), and update the position of each moth according to the Euclidean distance to obtain the updated position of each moth:

[0112]

[0113] in, is the moth m after the t+1th iteration j Location f k is the number of moths m after the tth iteration j The corresponding k-th flame after ascending sorting, is the number of moths m after the tth iteration j With flames f k The Euclidean distance is b, the spiral shape control parameter is b, the random parameter is θ∈[a,1], and the attenuation factor is a=-1-t / A.

[0114] In the embodiment of the present application, b is not specifically limited and can be set according to actual needs. When the source monitoring range is X∈[0,5]km, Y∈[0,5]km, Y∈[0,5]s, it is preferred to set b=1 to further improve the source positioning accuracy and efficiency. This update rule enables the moth to perform a local surround search around the flame (when Small) or a large range (when When the velocity is large), combined with the arrival time-distance linear relationship implied by the velocity model (arrival time is proportional to propagation distance), unreasonable source locations (such as areas where the distance far exceeds the wave velocity × time) can be quickly eliminated, the search for potential solution areas can be concentrated, and invalid calculations can be reduced. In the MFO global search phase, flames (candidate solutions) will be distributed along the direction of distance-time matching, guiding the moth population to move to high-probability areas and avoiding blind random search.

[0115] In another exemplary embodiment of the present application, the above step 301 further includes:

[0116] Step 402: Based on the preset earthquake source monitoring range, the updated position in Truncate.

[0117] In the embodiment of the present application, To truncate, As the final truncated in Substitute into the fitness value evaluation operation of the next step to calculate the updated position The fitness value of Wherein, D is the upper limit of the monitoring range of the preset X axis, E is the upper limit of the monitoring range of the preset Y axis, and F is the upper limit of the monitoring range of the preset Z axis.

[0118] Figure 4This is a schematic diagram of the process of the flame guiding the moth to shrink from the edge of the cube to the center in an embodiment of the present application. The green trajectory is the trajectory of the flame guiding the moth to shrink from the edge of the cube to the center. It can be seen from the figure that the speed / efficiency of locating the earthquake source through MFO (moth flame algorithm) is very high. When the number of iterations reaches 250, the first earthquake source estimate and the first earthquake time estimate located by MFO have reached a stable state, and the calculation time is only 0.2092 seconds.

[0119] In another exemplary embodiment of the present application, the aforementioned Markov Chain Monte Carlo method includes the MH algorithm.

[0120] In the examples of this application, the performance of three classic Markov chain Monte Carlo methods, namely the Hamiltonian Monte Carlo (HMC) algorithm, the No-U-Turn Sampler (NUTS) algorithm, and the Metropolis-Hastings (MH) algorithm, is compared. Through a comprehensive evaluation of the three dimensions of computational efficiency, positioning accuracy, and algorithm stability, the aim is to select the positioning method that is most suitable for engineering test scenarios. Figures 5-7 As shown, compared Figure 5 、 Figure 6 and Figure 7The Hamiltonian Monte Carlo (HMC) method achieved a location error of 0.5206 km and an earthquake onset error of 0.0464 seconds in 2000 iterations, with a computational time of 1.26 seconds. Its location accuracy and computational efficiency were lower than those of the other two methods. The No-U-Turn Sampler (NUTS) method demonstrated higher location accuracy (0.3517 km) and lower earthquake onset error (0.0274 seconds) with the same number of iterations, but its computational time increased significantly to 6.46 seconds, reflecting the additional computational burden introduced by the dynamic adjustment of the integration path. In contrast, the Metropolis-Hastings (MH) method completed the computation in just 0.03 seconds after 2000 iterations, achieving a location error of 0.4644 km, slightly lower than NUTS but significantly better than HMC. Its earthquake onset error also achieved the best of the three methods, at 0.0242 seconds. This result demonstrates that the MH method achieves a better balance between computational efficiency and accuracy. Although the NUTS algorithm has a smaller positioning error, its 6.46-second time consumption is insufficient for real-time processing. The HMC algorithm also offers no advantages in computational efficiency or accuracy. The MH algorithm, which requires no gradient calculations or path integrals and converges rapidly through random walk sampling alone, is more practical in scenarios requiring high real-time performance, such as microseismic monitoring. This is especially true when earthquake-onset accuracy is crucial for early warning systems. The MH algorithm's 0.0242-second error further highlights its engineering applicability. Therefore, for situations with limited computing resources or high-frequency data processing requirements, the MH algorithm offers orders of magnitude faster speeds at the cost of acceptable positioning error, making it a more practical choice.

[0121] In another exemplary embodiment of the present application, the above step 102 includes the following steps 501 to 503. Among them:

[0122] Step 501: Generate candidate points using the proposed distribution based on the current state:

[0123]

[0124] Among them, x' represents the candidate point (new state), x t′ represents the state of the t′th iteration, t′≥0, x0 is the first source parameter estimate, σp represents the standard deviation of the spatial coordinate, preferably setting σp = 0.1 km, diag represents the diagonal covariance matrix, and the time parameter is unconstrained. In this stage, the solution space neighborhood is quickly explored by fixing the step size.

[0125] In the present embodiment, an isotropic Gaussian distribution is preferred. In the first stage, a fixed-step perturbation (0.1 km standard deviation) is applied to the spatial coordinates using an isotropic Gaussian distribution and the time parameter constraints are released. This allows for independent, uniform multi-dimensional search within the earthquake source neighborhood, ensuring both the stability of spatial exploration and rapid convergence of earthquake occurrence times, thereby efficiently covering high-likelihood regions in the solution space.

[0126] Step 502: For each candidate point x', calculate the current acceptance probability:

[0127]

[0128] Among them, α t represents the acceptance probability of the t′th iteration; L(t obs |x′) represents a given candidate point x', the observed data t obs,i Likelihood value; L(t obs |x t′ ) means given the current state x t′ , observation data t obs,i Likelihood value; p(x′) represents the prior probability of candidate point x', p(x t′ ) represents the current state x t′ The prior probability of q(x t′ |x′) represents the proposal from candidate point x' to the current state x t′ The proposed distribution probability, q(x′|x t′ ) indicates that from the current state x t′ Proposal distribution probability of the proposed candidate point x'.

[0129] If α t′ ≥Uniform(0,1), let the state x of the t′+1th iteration be t′+1 = x', that is, accept the candidate point x'; otherwise, let x t′+1 =x t′ , that is, keep the current state x t′ ; Among them, Uniform(0,1) is a uniform random number in the interval [0,1].

[0130] Step 503, iterate the above steps 501 to 502 in a loop until a target number of source parameter samples that conform to the posterior distribution are obtained.

[0131] In the embodiment of the present application, steps 501 to 503 above ensure that the Markov chain accesses the parameter space (a multidimensional space consisting of all possible parameter values, for example, in earthquake source location, the parameter space is a four-dimensional space (x, y, z, t0), each dimension corresponds to a parameter (coordinate or time), and each spatial point represents a possible solution) with a probability proportional to the posterior density, and eventually converges to a stationary distribution. Finally, a posteriori inference is achieved through sample statistics (such as mean and quantiles), and its mathematical essence is the design of transition probabilities that meet detailed balance conditions.

[0132] Preferably, the number of iterations is set to 2000, the initial step size (σp) is set to 0.1, the adjustment interval is set to 50, the theoretical optimal acceptance rate is set to 0.234, Huber_delta is set to 0.15 (1.5σ principle, σ=0.1), the covariance estimation sample size is set to 150, and the number of iterations for starting covariance adaptation is set to 200.

[0133] In another exemplary embodiment of the present application, in the above step 102, in order to enhance the observation noise The robustness of the Markov chain Monte Carlo method is the likelihood function L(t obs )for:

[0134]

[0135] Where σ is the noise standard deviation, preferably set to σ = 0.1 s; t0 is the candidate point x' or the current state x t′ The moment of earthquake, [x, y, z] is the candidate point x' or the current state x t′ The three-dimensional coordinates of [x, y, z] are given by the velocity model, and v′ is the seismic wave propagation velocity determined by the three-dimensional coordinates [x, y, z].

[0136] In the embodiment of the present application, the likelihood function is used in step 502 to calculate L(t obs |x′) and L(t obs |x t′ ) to evaluate the matching degree of the source parameters. The traditional likelihood function uses the least squares function, which is sensitive to outliers and may cause positioning deviations due to noise data of individual vibration monitoring devices (such as outliers in Gaussian noise). The likelihood function of this application enhances robustness through the design of dual-mode penalty and dynamic threshold. The dual-mode penalty refers to the case where the residual |t obs,i -t pred,i When |≤1.5σ, square loss is used to ensure the accuracy of normal data. obs,i -t pred,iWhen |>1.5σ, it switches to linear penalty to suppress the excessive influence of abnormal data of vibration monitoring device (corresponding to Figure 7 The dynamic threshold dynamically sets the switching point based on the noise standard deviation σ (1.5σ = 0.15 seconds), preventing the fixed threshold from being inadequately adaptable to changes in noise levels. High-precision fitting is maintained within the normal data range, while outliers are weighted less, improving positioning stability. By improving the likelihood function, posterior sampling is more focused on high-probability areas, improving convergence efficiency.

[0137] In another exemplary embodiment of the present application, the above step 102 further includes:

[0138] Every Q iterations (preferably 50), the proposal distribution is updated according to the following formula:

[0139] ∑ t′ =Cov(X t′-150:t′ )+10 -6 I; (8)

[0140] q(x t′ |x′)=N(x t′ ,∑ t′ ); (9)

[0141] Among them, X t′-150:t′ is the historical sample matrix, which refers to the matrix composed of all accepted candidate points in the last 150 iterations; t′ is a Gaussian distribution N(x t ,∑ t ) reflects the correlation between parameters; I represents the unit matrix, 10 -6 I is the t′ Regularization term to prevent singular matrix.

[0142] In another exemplary embodiment of the present application, in order to maintain the theoretical optimal acceptance rate (e.g., 23.4%), it is proposed that the distribution adopts an isotropic Gaussian distribution. The above step 102 further includes:

[0143] Every Q iterations, the isotropic proposal scale of the isotropic Gaussian distribution is adjusted according to the following formula:

[0144] logσ p,t′+1 =logσ p,t′ +γ t′ (α t′ -P); (10)

[0145] Among them, σ p,t′+1 is the isotropic proposal scale for the t′+1th iteration, σ p,t′ is the isotropic proposal scale for the t′th iteration, γ t′Let \(\alpha\) be the decay learning rate and \(P\) be the optimal acceptance rate. \(P\) can be set to \(0.234\) to ensure accuracy. This mechanism automatically balances exploration (large step size) and exploitation (small step size), improving the sampling efficiency.

[0146] In another exemplary embodiment of the present application, to ensure that all samples satisfy physical feasibility and avoid invalid calculations, step 102 further includes:

[0147] Before step 502, impose hard constraints on the candidate points:

[0148] Project the spatial coordinates onto the search region:

[0149] x'←max(0,min(x',5)); (11)

[0150] There is no boundary limit for the earthquake origin time. However, due to the setting of the search range for the earthquake origin time, if the earthquake origin time \(t\) of the candidate point a is not within the preset search range \([a, b]\), reflect \(t\) a into the search range \([a, b]\):

[0151] If \(t\) a <a, then reflect it to \(t\) a =a+(a - t a );

[0152] If \(t\) a >b, then reflect it to \(t\) a =b-(t a -b).

[0153] In another exemplary embodiment of the present application, to ensure that all samples satisfy physical feasibility and avoid invalid calculations, in step 202, when generating the source parameter samples that conform to the posterior distribution by the Markov chain Monte Carlo method, it further includes:

[0154] After the warm-up period, retain the Markov chain samples that approximately follow the posterior distribution \(p(x|t\) obs )

[0155] In the embodiment of the present application, going through the warm-up period means discarding the unstable samples in the initial stage of the Markov chain to ensure that the subsequent samples used for analysis truly reflect the target posterior distribution.

[0156] In another exemplary embodiment of the present application, in step 203, calculating the joint probability density of each source parameter sample specifically includes:

[0157] Calculate the joint probability density of each source parameter sample according to the following formulas (12) - (15):

[0158]

[0159]

[0160] in, is the joint probability density of a given source parameter sample c, c∈{c i}, c i is the i-th source parameter sample, 1≤i≤n, n is the total number of source parameter samples, H is the bandwidth matrix, H- 1 is the inverse matrix of the bandwidth matrix, (cc i ) T (cc i ), d is the dimension of parameter space, is the covariance matrix, is the sample mean, for The transpose of .

[0161] Tested with synthetic data, the source monitoring range X, Y, Z is 5km, 8 vibration monitoring devices are set, and the velocity model is used. Figure 3 The layered velocity model shown in the figure is generated theoretically and then 10% random Gaussian noise is added. The positioning results using MFO alone are as follows: Figure 4 As shown, the absolute positioning error is 0.4752km, the earthquake occurrence time error is 0.0399s, and the calculation time is 0.2092s. The positioning results using MH alone are as follows: Figure 7 As shown in the figure, the positioning error is 0.4644 km, the earthquake occurrence time error is 0.0242 s, and the calculation time is 0.02 s; the MFO-MH earthquake source positioning results are as follows Figure 8 As shown in the figure, the positioning error is 0.2518 km, the earthquake occurrence time error is 0.0294 s, and the calculation time is 0.09 s. One point worth explaining is that the initial value of the MFO is not fixed, so the MFO trajectory obtained each time is different.

[0162] A comparison shows that the MFO-MH earthquake source location method improves positioning accuracy while optimizing computational efficiency by integrating the global search of the moth flame algorithm with the local optimization of the adaptive Bayesian (MH) algorithm. The moth flame algorithm, with its flame-guided group parallel search characteristics, quickly approaches the solution space region where the true location of the earthquake source is located, avoiding the initial value sensitivity problem of traditional optimization methods in complex multidimensional space. The subsequent improved Metropolis-Hastings algorithm is based on dynamic covariance adjustment and robust likelihood modeling, and achieves efficient posterior sampling through intelligent step-size control. Its adaptive proposal mechanism effectively captures the correlation between parameters and avoids the inefficient exploration of random walks. The two-stage algorithm forms a relay computing paradigm, an organic combination of global rapid convergence and local fine sampling, which significantly reduces invalid iterations while suppressing observation noise interference, ultimately achieving the coordinated optimization of positioning accuracy and computational efficiency.

[0163] The MFO-MH source location method of the embodiment of the present application was used to conduct actual data field tests on 8 blasting events in Dongtan Coal Mine. The information of the vibration monitoring device is shown in Table 2. The distribution of the vibration monitoring device and the velocity model are shown in Table 2. Figure 9 As shown, if z j Not greater than 0.00, longitudinal wave v = 2.2 km / s, transverse wave v = 1.2717 km / s; if z j Not greater than 0.65, longitudinal wave v = 3.2 km / s, transverse wave v = 1.8497 km / s; if z j Not greater than 1.00, longitudinal wave v = 3.4 km / s, transverse wave v = 1.9653 km / s; if z j Not greater than 4.00, longitudinal wave v = 6.1 km / s, transverse wave v = 3.5300 km / s; if z j Not greater than 5.00, longitudinal wave v = 6.2km / s, transverse wave v = 3.5800km / s.

[0164] Table 2 Dongtan Coal Mine Vibration Monitoring Device Information

[0165]

[0166]

[0167] MFO parameter settings: number of iterations: 50, initial population size (candidate solutions): 20, preset source monitoring range: X∈[0,3]km, Y∈[0,3]km, T∈[0,86400]s. Second stage MH parameter settings: number of iterations: 500, initial step size: 0.2, noise perturbation: 0.1, adaptive adjustment interval: 50, theoretical optimal acceptance rate: 0.234, Huber_delta: 0.15 (1.5σ principle), covariance estimation sample size: 100, adaptive start iteration number: 200. Positioning result example: Figure 10 The 8 complete positioning results are shown in Table 3.

[0168] Table 3 Dongtan Coal Mine blasting event location results

[0169] Event ID Horizontal error / km Depth error / km Absolute error / km Seismic time error / s Calculation time / s 1 0.0210 0.0340 0.0400 0.0162 0.06 2 0.0150 0.0410 0.0437 0.0122 0.06 3 0.0190 0.0330 0.0381 0.0174 0.07 4 0.0180 0.0470 0.0503 0.0201 0.07 5 0.0160 0.0450 0.0478 0.0198 0.06 6 0.0110 0.0210 0.0237 0.0143 0.06 7 0.0370 0.0590 0.0696 0.0158 0.06 8 0.0130 0.0340 0.0364 0.0152 0.06

[0170] It can be seen from the actual data positioning results that the dual-stage source positioning method (MFO-MH source positioning method) of coal mine microseismic events in the embodiment of the present application can accurately locate coal mine microseismic events within 0.1s, and the positioning error complies with the national standard GB / T25217.4-2019 "Impact ground pressure measurement, monitoring and prevention methods Part 4: Microseismic monitoring method", showing excellent performance.

[0171] Based on the same inventive concept, embodiments of the present application also provide a dual-stage source locating device for coal mine microseismic events, which is used to implement the dual-stage source locating method for coal mine microseismic events described above. The solution provided by this device is similar to the solution described in the aforementioned method. Therefore, the specific limitations of one or more embodiments of the dual-stage source locating device for coal mine microseismic events provided below can be found in the aforementioned limitations of the dual-stage source locating method for coal mine microseismic events, and will not be further elaborated here.

[0172] In an exemplary embodiment, Figure 11 As shown, a dual-stage source location device 60 for a microseismic event in a coal mine is provided, comprising:

[0173] Vibration monitoring device 601, used to monitor and record the arrival data of microseismic events in coal mines;

[0174] A first optimization module 602 is configured to determine, based on the arrival time data and seismic wave propagation velocity observed by each vibration monitoring device, a first source parameter estimate of the vibration event using a moth flame algorithm; wherein the first source parameter estimate includes a first source estimate and a first onset time estimate;

[0175] A second optimization module 603 is configured to use the first source parameter estimate as the initial state of a Markov chain Monte Carlo method to generate source parameter samples that conform to a posterior distribution, wherein each source parameter sample includes a second source estimate and a second earthquake onset time estimate;

[0176] The kernel density estimation module 604 is used to calculate the joint probability density of each source parameter sample and select the source parameter sample with the largest joint probability density to obtain the final source parameter estimation value of the coal mine microseismic event.

[0177] In another exemplary embodiment of the present application, the first optimization module 602 is further configured to:

[0178] Perform population initialization operations:

[0179] Initialize the number of moths in the population, the total number of iterations, and randomly initialize the initial position of each moth according to the preset monitoring range in, is the jth moth m j The initial shock moment, It's a moth j The initial three-dimensional coordinates, 1≤j≤M, M is the number of moths, For moths j The initial X-axis coordinate of For moths j The initial Y-axis coordinate of For moths j The initial Z-axis coordinate of

[0180] Perform fitness value evaluation operation:

[0181] Based on the current position of each moth Calculate the current fitness value of each moth through the objective function in, is the number of moths m after the tth iteration j Location is the number of moths m after the tth iteration j The time of earthquake, 0≤t≤A, A is the total number of iterations, is the number of moths m after the tth iteration j The three-dimensional coordinates of is the number of moths m after the tth iteration j The X-axis coordinate, is the number of moths m after the tth iteration j The Y-axis coordinate, is the number of moths m after the tth iteration j The Z-axis coordinate of

[0182] To perform a flame selection operation:

[0183] Sort the fitness values ​​of the current positions of all moths in ascending order, and select the top N in ascending order f The current position of the moth corresponding to the (t) fitness value is taken as the flame (elite solution set), N f The calculation formula for (t) is:

[0184]

[0185] Among them, g() is the rounding function, and t is the current number of iterations;

[0186] The position update operation, fitness value evaluation operation and flame selection operation steps are performed in sequence in a loop iteration until the difference of the flame values ​​obtained by multiple consecutive loop iterations meets the preset requirements or the current iteration number t reaches the total iteration number A. The flame obtained is the first source parameter estimation value.

[0187] In the embodiment of the present application, the above objective function is:

[0188]

[0189] Among them, t pred,i is the predicted arrival time of the ith vibration monitoring device, t obs,i is the superimposed Gaussian noise ∈~N(0,0.1 2 ) is observed, N is the total number of vibration monitoring devices, is the number of moths m after the tth iteration j The moment of earthquake, s i is the three-dimensional coordinate of the i-th vibration monitoring device, v is the velocity model constructed according to the three-dimensional coordinate Determine the speed at which seismic waves propagate.

[0190] Location update operations include:

[0191] According to the current position of each moth and the corresponding flame, the Euclidean distance between each moth and the corresponding flame is calculated, and the position of each moth is updated according to the Euclidean distance between each moth and the corresponding flame to obtain the updated position of each moth.

[0192] In another exemplary embodiment of the present application, the first optimization module 602 is further configured to:

[0193] According to the following formula (4), the Euclidean distance between each moth and the corresponding flame is calculated, and the position of each moth is updated according to the Euclidean distance to obtain the updated position of each moth:

[0194]

[0195] in, is the moth m after the t+1th iteration j Location f k is the number of moths m after the tth iteration j The corresponding k-th flame after ascending sorting, is the number of moths m after the tth iteration j With flames k The Euclidean distance is b, the spiral shape control parameter is b, the random parameter is θ∈[a,1], and the attenuation factor is a=-1-t / A.

[0196] In another exemplary embodiment of the present application, the first optimization module 602 is further configured to:

[0197] According to the preset source monitoring range, the updated position in Truncate.

[0198] In the embodiment of the present application, To truncate, As the final truncated in Substitute into the fitness value evaluation operation of the next step to calculate the updated position The fitness value of Wherein, D is the upper limit of the monitoring range of the preset X axis, E is the upper limit of the monitoring range of the preset Y axis, and F is the upper limit of the monitoring range of the preset Z axis.

[0199] In another exemplary embodiment of the present application, the second optimization module 603 is further configured to:

[0200] It is used to use the first earthquake source parameter estimate as the initial state of the MH algorithm, and generate earthquake source parameter samples that conform to the posterior distribution through the MH algorithm, wherein each earthquake source parameter sample includes a second earthquake source estimate and a second earthquake occurrence time estimate.

[0201] In another exemplary embodiment of the present application, the second optimization module 603 is further configured to:

[0202] Based on the current state, use the proposal distribution to generate candidate points:

[0203]

[0204] Among them, x' represents the candidate point (new state), x t′represents the state of the t′th iteration, t′≥0, t≥0, x0 is the first source parameter estimate, σp represents the standard deviation of the spatial coordinate, preferably setting σp=0.1 km, diag represents the diagonal covariance matrix, and the time parameter is unconstrained;

[0205]

[0206] Among them, α t represents the acceptance probability of the t′th iteration; L(t obs |x′) represents a given candidate point x', the observed data t obs,i Likelihood value; L(t obs |x t′ ) means given the current state x t′ , observation data t obs,i Likelihood value; p(x′) represents the prior probability of candidate point x', p(x t′ ) represents the current state x t′ The prior probability of q(x t′ |x′) represents the proposal from candidate point x' to the current state x t′ The proposed distribution probability, q(x′|x t′ ) indicates that from the current state x t′ Proposal distribution probability of the candidate point x';

[0207] If α t′ ≥Uniform(0,1), let the state x of the t′+1th iteration be t′+1 = x', that is, accept the candidate point x'; otherwise, let x t′+1 =x t′ , that is, keep the current state x t′ ; Among them, Uniform(0,1) is a uniform random number in the interval [0,1];

[0208] The above steps are iterated until the target number of source parameter samples that conform to the posterior distribution is obtained.

[0209] In another exemplary embodiment of the present application, the second optimization module 603 is further configured to:

[0210] The likelihood function L(t obs )for:

[0211]

[0212] Where σ is the noise standard deviation, preferably set to σ = 0.1 s; t0 is the candidate point x' or the current state x t′ The moment of earthquake, [x, y, z] is the candidate point x' or the current state xt′ The three-dimensional coordinates of [x, y, z] are given by the velocity model, and v′ is the seismic wave propagation velocity determined by the three-dimensional coordinates [x, y, z].

[0213] In another exemplary embodiment of the present application, the second optimization module 603 is further configured to:

[0214] Every Q iterations (preferably 50), the proposal distribution is updated according to the following formula:

[0215] ∑ t′ =Cov(X t′-150:t′ )+10 -6 I; (8)

[0216] q(x t′ |x′)=N(x t′ ,∑ t′ ); (9)

[0217] Among them, X t′-150:t′ is the historical sample matrix, which refers to the matrix composed of all accepted candidate points in the last 150 iterations; t′ is a Gaussian distribution N(x t ,∑ t ) reflects the correlation between parameters; I represents the unit matrix, 10 -6 I is the t′ Regularization term to prevent singular matrix.

[0218] In another exemplary embodiment of the present application, the proposed distribution adopts an isotropic Gaussian distribution.

[0219] The second optimization module 603 is further configured to:

[0220] Every Q iterations, the isotropic proposal scale of the isotropic Gaussian distribution is adjusted according to the following formula:

[0221] logσ p,t′+1 =logσ p,t′ +γ t′ (α t′ -P); (10)

[0222] Among them, σ p,t′+1 is the isotropic proposal scale for the t′+1th iteration, σ p,t′ is the isotropic proposal scale for the t′th iteration, γ t′ To decay the learning rate, P is the optimal acceptance rate, which can be set to 0.234 to ensure accuracy. This mechanism automatically balances exploration (large step size) and development (small step size), improving sampling efficiency.

[0223] In another exemplary embodiment of the present application, the above-mentioned second optimization module 603 is further configured to:

[0224] Before calculating the current acceptance probability for each candidate point x', impose a hard constraint on the candidate point:

[0225] Project the spatial coordinates onto the search region:

[0226] x' ← max(0, min(x', 5)); (11)

[0227] There is no boundary limit for the origin time, but since there is a search range for the origin time, therefore, if the origin time t of the candidate point a is not within the preset search range [a, b], reflect t a into the search range [a, b]:

[0228] If t a < a, then reflect it to t a = a + (a - t a );

[0229] If t a > b, then reflect it to t a = b - (t a - b).

[0230] In another exemplary embodiment of the present application, to ensure that all samples meet physical feasibility and avoid invalid calculations, the above-mentioned second optimization module 603 is further configured to:

[0231] After the warm-up period, retain the Markov chain samples approximately following the posterior distribution p(x|t obs )

[0232] In another exemplary embodiment of the present application, the above-mentioned kernel density estimation module 604 is further configured to:

[0233] Calculate the joint probability density of each source parameter sample according to the following formulas (12) to (15):

[0234]

[0235] Where is the joint probability density of the given source parameter sample c, c ∈ {c i}, c​​​​​​​​​), d is the dimension of parameter space, is the covariance matrix, is the sample mean, for The transpose of .

[0236] In an exemplary embodiment, a computer device is provided. The computer device may be a server or a terminal. The internal structure diagram thereof may be as follows: Figure 12 As shown. The computer device includes a processor, a memory, an input / output interface (Input / Output, abbreviated as I / O) and a communication interface. The processor, memory and input / output interface are connected through a system bus, and the communication interface is connected to the system bus through the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program and a database. The internal memory provides an environment for the operation of the operating system and computer program in the non-volatile storage medium. The database of the computer device is used to store dual-stage source location data of coal mine microseismic events. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal through a network connection. When the computer program is executed by the processor, a dual-stage source location method for coal mine microseismic events is implemented.

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

[0238] In an exemplary embodiment, a computer device is further provided, including a memory and a processor. The memory stores a computer program, and the processor implements the steps in the above method embodiments when executing the computer program.

[0239] In an exemplary embodiment, a computer-readable storage medium is provided, storing a computer program. When the computer program is executed by a processor, the steps in the above-mentioned method embodiments are implemented.

[0240] In an exemplary embodiment, a computer program product is provided, including a computer program. When the computer program is executed by a processor, the steps in the above method embodiments are implemented.

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

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

[0243] The databases involved in the various embodiments provided herein may include at least one of a relational database and a non-relational database. Non-relational databases may include, but are not limited to, distributed databases based on blockchains. The processors involved in the various embodiments provided herein may include, but are not limited to, general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic units, data processing logic units based on quantum computing, and the like.

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

[0245] This document uses specific examples to illustrate the principles and implementation methods of this application. The above examples are intended only to facilitate understanding of the method and core concepts of this application. Furthermore, those skilled in the art will appreciate that variations in the specific implementation methods and scope of application are possible based on the concepts of this application. In summary, this specification should not be construed as limiting this application.

Claims

1. A dual-stage source location method for microseismic events in coal mines, characterized by: The dual-stage source location method for coal mine microseismic events includes: Determine, based on the arrival time data and seismic wave propagation velocity observed by each vibration monitoring device, a first source parameter estimate of the coal mine microseismic event using a moth flame algorithm; wherein the first source parameter estimate includes a first source estimate and a first earthquake occurrence time estimate; Using the first earthquake source parameter estimate as the initial state of a Markov chain Monte Carlo method, generating earthquake source parameter samples that conform to a posterior distribution through the Markov chain Monte Carlo method, wherein each of the earthquake source parameter samples includes a second earthquake source estimate and a second earthquake occurrence time estimate; The joint probability density of each source parameter sample is calculated, and the source parameter sample with the largest joint probability density is selected to obtain the final source parameter estimation value of the coal mine microseismic event.

2. The dual-stage source location method for coal mine microseismic events according to claim 1, characterized in that: The method of determining the estimated value of the first source parameter of the coal mine microseismic event by using the moth flame algorithm based on the arrival time data and seismic wave propagation velocity observed by each vibration monitoring device specifically includes: Perform population initialization operations: Initialize the number of moths in the population, the total number of iterations, and randomly initialize the initial position of each moth according to the preset source monitoring range in, is the jth moth m j The initial shock moment, It's a moth j The initial three-dimensional coordinates of , 1≤j≤M, M is the number of moths; Perform fitness value evaluation operation: Based on the current position of each moth Calculate the current fitness value of each moth through the objective function in, is the number of moths m after the tth iteration j Location is the number of moths m after the tth iteration j The time of earthquake, 0≤t≤A, A is the total number of iterations, is the number of moths m after the tth iteration j The three-dimensional coordinates of To perform a flame selection operation: Sort the fitness values ​​of the current positions of all moths in ascending order, and select the top N in ascending order f The current position of the moth corresponding to the (t) fitness value is taken as the flame, N f The calculation formula for (t) is: Among them, g() is the rounding function, and t is the current number of iterations; Iteratively executing the steps of position updating, fitness value evaluation, and flame selection in sequence until the deviation of the flame value obtained by multiple consecutive iterations meets the preset requirement or the current iteration number t reaches the total iteration number, and the obtained flame is the first source parameter estimation value; The location update operation includes: According to the current position of each moth and the corresponding flame, the Euclidean distance between each moth and the corresponding flame is calculated, and the position of each moth is updated according to the Euclidean distance to obtain the updated position of each moth.

3. The dual-stage source location method for coal mine microseismic events according to claim 2, characterized in that: The objective function is: Among them, t pred,i is the currently predicted arrival time of the i-th vibration monitoring device, t obs,i is the observation time of the i-th vibration monitoring device, N is the total number of vibration monitoring devices, is the number of moths m after the tth iteration j The moment of earthquake, s i is the three-dimensional coordinate of the i-th vibration monitoring device, v is the velocity model constructed according to the three-dimensional coordinate Determine the speed at which seismic waves propagate.

4. The dual-stage source location method for coal mine microseismic events according to claim 2, characterized in that: The method further comprises calculating the Euclidean distance between each moth and the corresponding flame according to the current position of each moth and the corresponding flame, and updating the position of each moth according to the Euclidean distance to obtain the updated position of each moth, specifically comprising: According to the following formula (4), the Euclidean distance between each moth and the corresponding flame is calculated, and the position of each moth is updated according to the Euclidean distance to obtain the updated position of each moth: in, is the moth m after the t+1th iteration j Location f k is the number of moths m after the tth iteration j The corresponding k-th flame after ascending sorting, is the number of moths m after the tth iteration j With flames k The Euclidean distance is b, the spiral shape control parameter is b, the random parameter is θ∈[a,1], and the attenuation factor is a=-1-t / A.

5. The dual-stage source location method for coal mine microseismic events according to claim 1, characterized in that: The Markov Chain Monte Carlo method includes the MH algorithm.

6. The dual-stage source location method for coal mine microseismic events according to claim 3, characterized in that: The method of using the first source parameter estimate as the initial state of the Markov chain Monte Carlo method and generating source parameter samples that conform to the posterior distribution by the Markov chain Monte Carlo method specifically includes: Based on the current state, use the proposal distribution to generate candidate points: Among them, x' represents the candidate point, x t′ represents the state of the t′th iteration, t′≥0, x0 is the estimated value of the first source parameter, σ p Represents the standard deviation of spatial coordinates, diag represents the diagonal covariance matrix; For each candidate point x', calculate the current acceptance probability: Among them, α t represents the acceptance probability of the t′th iteration; L(t obs |x′) represents a given candidate point x', the observed data t obs,i Likelihood value; L(t obs |x t′ ) means given the current state x t′ , observation data t obs,i Likelihood value; p(x′) represents the prior probability of candidate point x', p(x t′ ) represents the current state x t′ The prior probability of q(x t′ |x′) represents the proposal from candidate point x' to the current state x t′ The proposed distribution probability, q(x′|x t′ ) indicates that from the current state x t′ Proposal distribution probability of candidate point x'; where L(t obs |x′) and L(t obs |x t′ )The likelihood function L(t obs )for: Where σ is the standard deviation of noise, t0 is the candidate point x' or the current state x t′ The moment of earthquake, [x, y, z] is the candidate point x' or the current state x t′ The three-dimensional coordinates of [x, y, z] are obtained by the velocity model, and v′ is the seismic wave propagation velocity determined by the three-dimensional coordinates [x, y, z]. If α t′ ≥Uniform(0,1), let the state x of the t′+1th iteration be t′+1 = x', that is, accept the candidate point x'; otherwise, let x t′+1 =x t′ , that is, keep the current state x t′ ; Among them, Uniform(0,1) is a uniform random number in the interval [0,1]; The above steps are iterated until the target number of source parameter samples that conform to the posterior distribution is obtained.

7. The dual-stage source location method for coal mine microseismic events according to claim 6, characterized in that: The proposed distribution adopts an isotropic Gaussian distribution; The method of using the first source parameter estimate as the initial state of the Markov chain Monte Carlo method and generating source parameter samples that conform to the posterior distribution by the Markov chain Monte Carlo method further includes: Every Q iterations, the proposal distribution is updated according to the following formula: ∑ t′ =Cov(X t′-150:t′ )+10 -6 I;(8) q(x t′ |x′)=N(x t′ ,∑ t′ );(9) Among them, X t′-150:t′ is the historical sample matrix, ∑ t′ is a Gaussian distribution N(x t ,∑ t )’s covariance matrix, I represents the identity matrix; The method of using the first source parameter estimate as the initial state of the Markov chain Monte Carlo method and generating source parameter samples that conform to the posterior distribution by the Markov chain Monte Carlo method further includes: Every Q iterations, the isotropic proposal scale of the isotropic Gaussian distribution is adjusted according to the following formula: log p,t′+1 =logσ p,t′ +g t′ (a t′ -P);(10) Among them, σ p,t′+1 is the isotropic proposal scale for the t′+1th iteration, σ p,t′ is the isotropic proposal scale for the t′th iteration, γ t′ is the decay learning rate, and P is the optimal acceptance rate.

8. The dual-stage source location method for coal mine microseismic events according to claim 1, characterized in that: The calculation of the joint probability density of each source parameter sample specifically includes: According to the following formulas (12) to (15), the joint probability density of each source parameter sample is calculated: in, is the joint probability density of a given source parameter sample c, c∈{c i }, c i is the i-th source parameter sample, 1≤i≤n, n is the total number of source parameter samples, H is the bandwidth matrix, d is the parameter space dimension, is the covariance matrix, is the sample mean.

9. A dual-stage source location device for microseismic events in coal mines, characterized in that: The dual-stage source location device for coal mine microseismic events includes: Vibration monitoring device, used to monitor and record the arrival data of microseismic events in coal mines; a first optimization module configured to determine, based on the arrival time data and seismic wave propagation velocity observed by each vibration monitoring device, a first source parameter estimate of the vibration event using a moth flame algorithm; wherein the first source parameter estimate includes a first source estimate and a first earthquake onset time estimate; a second optimization module, configured to use the first earthquake source parameter estimate as an initial state of a Markov chain Monte Carlo method, and generate earthquake source parameter samples that conform to a posterior distribution through the Markov chain Monte Carlo method, wherein each of the earthquake source parameter samples includes a second earthquake source estimate and a second earthquake onset time estimate; The kernel density estimation module is used to calculate the joint probability density of each source parameter sample and select the source parameter sample with the largest joint probability density to obtain the final source parameter estimation value of the coal mine microseismic event.

10. A computer device comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the dual-stage source location method for coal mine microseismic events described in any one of claims 1-6.

Images