Micro-seismic time sequence prediction method and system based on energy filling

Abstract

The invention discloses a micro-seismic time sequence prediction method and system based on energy filling, and relates to the technical field of mine safety monitoring and geological disaster early warning. The method comprises the following steps: performing time difference positioning analysis on a non-uniform wave velocity model through acquired microseismic waveform data and microseismic sensor coordinates to obtain three-dimensional coordinates and initial time of a seismic source; correcting the original energy based on the path characteristics, the bedding attenuation law and the spherical diffusion effect to obtain corrected seismic source energy; calculating a correction factor and a confidence coefficient based on the corrected seismic source energy and the original energy, and screening the corrected seismic source energy according to the confidence coefficient to obtain representative energy; based on the initial time, the representative energy and the confidence coefficient, constructing a microseismic event sequence feature tensor; inputting the micro-seismic event sequence feature tensor into a pre-trained micro-seismic energy prediction model to obtain a maximum energy prediction value; according to the method, the problem of energy deviation caused by rock anisotropy is solved, and the accuracy of micro-seismic monitoring and early warning is improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of mine safety monitoring and geological disaster early warning, more specifically, the present application relates to a microseismic time series prediction method and system based on energy filling. BACKGROUND

[0002] In the incubation and occurrence process of coal mine dynamic disasters such as rock burst and coal and gas outburst, the concentration and release of stratum stress are often accompanied, and a large number of microseismic events are generated in this process. Through time series analysis of microseismic monitoring data including spatial coordinates, occurrence time and energy, early warning of high-risk areas can be realized. The existing related technology usually extracts the intensity information of the seismic source from the time series signal by using the energy calculation model based on the elastic wave data received by the sensor, and then establishes a dangerous event discrimination model. Specifically, this kind of technology mainly depends on the energy parameters calculated from the waveform, and then regards the calculated energy value as the representation of the release intensity of the seismic source, so as to complete the identification and early warning of the dangerous event.

[0003] However, the core design logic of the existing technology depends on the energy attenuation law of the elastic wave generated by the microseismic in the uniform or ideal medium model, and does not fully consider the key physical scene of the azimuth angle energy deviation caused by the anisotropy of the rock. In actual rock mass, the sedimentary rock with obvious bedding structure or the rock mass with developed joints, the elastic wave attenuates slowly along the bedding direction, and attenuates fast vertically to the bedding direction. When the direction of the line connecting the seismic source and the sensor is different from the angle of the rock layer, the same seismic source is received by different azimuth sensors, and the calculated energy value may be significantly different. At this time, the fluctuation of the time series data may not be caused by the change of the source intensity, but by the directional attenuation difference of the propagation path relative to the bedding. The existing technology performs energy calculation and comparison under the premise of assuming that the medium is isotropic, and the evaluation criteria cannot distinguish whether the energy difference comes from the change of the source intensity or from the directional attenuation of the propagation path. Therefore, it is difficult to accurately judge whether the low energy calculation value is caused by the weak source intensity or by the fast attenuation of the wave propagation path perpendicular to the bedding.

[0004] The above problems will cause system error of the existing technology when facing the rock mass scene with significant anisotropy, which is specifically manifested as underestimation of the intensity of the dangerous event from a specific propagation direction. This kind of situation will reduce the reliability and early warning accuracy of the microseismic monitoring system, so that part of the real threat cannot be effectively identified, thereby affecting the actual effectiveness of the rock mass stability monitoring system.

[0005] In view of this, the present application proposes a microseismic time series prediction method and system based on energy filling to solve the above problems. SUMMARY

[0006] To overcome the aforementioned shortcomings of the prior art and achieve the above objectives, the present invention provides the following technical solution: a microseismic time series prediction method based on energy filling, comprising:

[0007] A non-uniform wave velocity model was built based on the collected rock geological exploration data. The time difference positioning analysis of the non-uniform wave velocity model was performed by the collected microseismic waveform data and microseismic sensor coordinates to obtain the three-dimensional coordinates of the source and the initial time.

[0008] The raw energy and path characteristics of microseismic waveform data are extracted. Based on the path characteristics, bedding attenuation law and spherical diffusion effect, the raw energy is corrected to obtain the corrected source energy.

[0009] The correction factor is calculated based on the corrected source energy and the original energy. The confidence coefficient is calculated based on the correction factor. The corrected source energy is then selected based on the confidence coefficient to obtain the representative energy.

[0010] Based on initial time, representative energy, and confidence coefficient, a feature tensor for the microseismic event sequence is constructed.

[0011] The maximum energy prediction value is obtained by inputting the feature tensor of the microseismic event sequence into the pre-trained microseismic energy prediction model.

[0012] Furthermore, methods for correcting the original energy include:

[0013] Path characteristics include total propagation distance, actual propagation distance of each rock layer, and path direction angle;

[0014] Based on path characteristics and bedding attenuation laws, the initial attenuation parameters of each rock mass are calculated according to the path direction angle.

[0015] The initial attenuation parameters are weighted and summed using the ratio of the actual propagation distance of each rock mass to the total propagation distance to obtain the comprehensive attenuation parameters.

[0016] Based on the total propagation distance and spherical diffusion effect, the geometric attenuation factor is calculated using the classic spherical diffusion attenuation factor formula;

[0017] Multiplying the geometric attenuation factor by the comprehensive attenuation parameter yields the total attenuation compensation factor;

[0018] Divide the original energy by the total attenuation compensation factor to obtain the corrected source energy.

[0019] Furthermore, the method for calculating the initial attenuation parameters of each rock mass based on path direction angle decomposition includes:

[0020] The projections of the path direction angle onto the preset vertical and parallel bedding directions are used as the initial weight coefficients for the vertical and parallel bedding directions, respectively.

[0021] The initial weighting coefficients for the perpendicular and parallel bedding directions are normalized.

[0022] Based on the normalized initial weight coefficients for the vertical bedding direction and the parallel bedding direction, the actual propagation distance is decomposed into the effective propagation distance in the vertical bedding direction and the effective propagation distance in the parallel bedding direction. For each rock mass layer traversed by the propagation path, the target vertical bedding attenuation parameter corresponding to the effective propagation distance in the vertical bedding direction is extracted from the pre-constructed vertical bedding attenuation parameter table, and the target parallel bedding attenuation parameter corresponding to the effective propagation distance in the parallel bedding direction is extracted from the pre-constructed parallel bedding attenuation parameter table.

[0023] The initial attenuation parameter of the current rock layer is obtained by adding the target vertical bedding attenuation parameter to the target parallel bedding attenuation parameter.

[0024] Furthermore, methods for calculating correction factors based on the corrected source energy and the original energy include:

[0025] Calculate the ratio of the corrected source energy to the original energy to obtain the basic correction coefficient;

[0026] Divide the basic correction factor by the geometric attenuation factor to obtain the correction factor.

[0027] Furthermore, methods for calculating confidence coefficients based on correction factors include:

[0028] The correction factor for each channel is standardized to obtain the standardized correction factor for each channel.

[0029] The standardized correction factor is mapped to the [0,1] interval using a preset mapping function to obtain the confidence coefficient.

[0030] Furthermore, methods for screening and correcting seismic source energy based on confidence coefficients include:

[0031] Eliminate corrected seismic source energies with confidence coefficients lower than the preset confidence threshold;

[0032] The corrected source energies after elimination are sorted in descending order, and the maximum corrected source energy is selected as the representative energy of the corresponding microseismic event.

[0033] Furthermore, the training process of the microseismic energy prediction model employs a weighted loss function to optimize the microseismic energy prediction model;

[0034] Collect a dataset of historical microseismic event sequences, which includes the confidence coefficients and true values ​​of the actual source energy of the historical microseismic events.

[0035] The actual source energy is incremented by 1 and then logarithmically transformed to obtain the potential damage value. The potential damage value is added to the preset base weight to obtain the risk-oriented weight.

[0036] Calculate the mean of the confidence coefficients corresponding to all historical microseismic events to obtain the physical confidence weights;

[0037] The risk-oriented weight is coupled with the physical confidence weight to obtain a comprehensive weight;

[0038] The overall weight is used as the weight of the weighted loss function.

[0039] Furthermore, methods for coupling risk-oriented weights with physical confidence weights include:

[0040] Physical uncertainty is defined as 1 minus the physical confidence weight;

[0041] The actual penalty term is obtained by multiplying the preset adjustment factor by the physical uncertainty;

[0042] Subtract the actual penalty term from 1 to obtain the linear adjustment coefficient;

[0043] The risk-oriented weight is multiplied by the linear adjustment coefficient to obtain the comprehensive weight.

[0044] Furthermore, methods for constructing non-uniform wave velocity models based on collected rock mass geological exploration data include:

[0045] A reference coordinate system is established with the horizontal plane of the monitoring area as the xy plane and the z-axis pointing vertically upwards, thus constructing a three-dimensional geological framework.

[0046] The geological exploration data of the rock mass includes the strike and dip of the bedding planes; the unit vector of the intersection of the strike and the horizontal plane is taken as the strike vector; the unit vector perpendicular to the strike vector is selected as the dip vector, and the z-axis component of the dip vector corresponds to the dip angle of the bedding planes; the bedding normal vector is obtained by the cross product operation of the strike vector and the dip vector.

[0047] The three-dimensional geological framework is divided into multiple three-dimensional units according to a preset grid size. Within each three-dimensional unit, a local orthogonal coordinate system is formed by the strike line vector, dip line vector, and bedding normal vector.

[0048] In a local orthogonal coordinate system, the pre-tested vertical and parallel bedding wave velocities are mapped to the x, y, and z directions of the reference coordinate system through coordinate transformation. The equivalent directional wave velocity components of the corresponding three-dimensional elements in the reference coordinate system are obtained. The equivalent directional wave velocity components of all three-dimensional elements are then combined to form a non-uniform wave velocity model.

[0049] A microseismic time series prediction system based on energy filling includes:

[0050] The microseismic location module builds a non-uniform wave velocity model based on the collected rock geological exploration data. It performs time difference location analysis on the non-uniform wave velocity model by collecting microseismic waveform data and microseismic sensor coordinates to obtain the three-dimensional coordinates of the seismic source and the initial time.

[0051] The attenuation correction module extracts the original energy and path characteristics of the microseismic waveform data, and corrects the original energy based on the path characteristics, bedding attenuation law and spherical diffusion effect to obtain the corrected source energy.

[0052] The energy filling module calculates a correction factor based on the corrected source energy and the original energy, calculates a confidence coefficient based on the correction factor, and filters the corrected source energy according to the confidence coefficient to obtain representative energy;

[0053] The serialization module constructs a feature tensor for the microseismic event sequence based on the initial time, representative energy, and confidence coefficient.

[0054] The microseismic prediction module inputs the feature tensor of the microseismic event sequence into the pre-trained microseismic energy prediction model to obtain the maximum energy prediction value.

[0055] Compared with existing technologies, the technical effects and advantages of the energy-filling-based microseismic time series prediction method and system of the present invention are as follows:

[0056] This invention integrates rock mass geological exploration data, microseismic sensor coordinates, and microseismic waveform data to build a non-uniform wave velocity model that matches actual geological conditions. It combines field tests to obtain vertical and parallel bedding attenuation parameter tables and completes the original energy correction based on propagation path characteristics and spherical diffusion effects. Addressing the shortcomings of existing technologies that cannot distinguish the causes of energy differences and tend to underestimate the intensity of dangerous events, this invention achieves accurate restoration of the true energy of the seismic source, eliminates systematic errors caused by rock mass anisotropy, and improves the accuracy of microseismic energy characterization.

[0057] This invention constructs a microseismic energy prediction model based on a bidirectional long short-term memory neural network and an attention mechanism. It captures temporal evolution patterns in the forward direction and supplements precursor correlations in the reverse direction. At the same time, it associates confidence coefficients through attention weights. It adopts a weighted loss function that couples risk-oriented weights and physical confidence weights, focusing on high-energy catastrophic events to reduce the risk of underreporting. The penalty intensity is dynamically adjusted according to data reliability to avoid fitting noise. This approach balances the generalization ability of the microseismic energy prediction model with the needs of risk prevention and control, thereby improving the accuracy of microseismic maximum energy prediction.

[0058] The present invention is adapted to complex rock mass scenarios with anisotropy, improves the reliability and practical effectiveness of microseismic monitoring systems, reduces the engineering safety risks of dynamic disasters such as rock bursts, coal and gas outbursts, and ensures the continuous and effective operation of rock mass stability monitoring systems. Attached Figure Description

[0059] Figure 1 This is a schematic diagram of a microseismic timing prediction system based on energy filling, according to an embodiment of the present invention.

[0060] Figure 2 This is a flowchart of the energy-filled microseismic time series prediction method according to an embodiment of the present invention;

[0061] Figure 3 This is a flowchart illustrating the initial attenuation parameters for synthesizing each rock layer according to an embodiment of the present invention.

[0062] Figure 4 This is a flowchart illustrating the process of obtaining representative energy according to an embodiment of the present invention. Detailed Implementation

[0063] The technical solutions of the embodiments of the present invention will be described in detail, clearly, and completely below with reference to the accompanying drawings. It should be particularly noted that the specific embodiments described below are only for better illustrating and explaining the technical solutions of the present invention, and are intended to enable those skilled in the art to better understand and implement the present invention, and should not be construed as limiting the scope of protection of the present invention. Without departing from the spirit and substance of the present invention, those skilled in the art can modify, adjust, or make equivalent substitutions based on the content disclosed in the present invention, and these should all be considered within the scope of protection of the present invention.

[0064] Example 1:

[0065] Please see Figure 1 As shown, this embodiment discloses a microseismic time series prediction system based on energy filling, including a microseismic location module, an attenuation correction module, an energy filling module, a serialization module, and a microseismic prediction module. Each module is connected by wired and / or wireless means to realize data transmission.

[0066] The microseismic location module, based on the collected microseismic waveform data, rock geological exploration data, and microseismic sensor coordinates, builds a non-uniform wave velocity model, performs time difference location analysis on the non-uniform wave velocity model, and obtains the three-dimensional coordinates and initial time of the seismic source.

[0067] Microseismic waveform data is acquired synchronously by an array of microseismic sensors deployed in the monitoring area. This data includes the raw elastic wave signals of the microseismic event captured by each sensor at the current sampling time. Due to the significant anisotropic characteristics of the rock mass, the propagation velocity and attenuation characteristics of the elastic waves generated by the microseismic event vary significantly in different spatial directions, resulting in differences in arrival time and energy characteristics of the same seismic source in different microseismic sensor channels. To avoid misinterpreting the propagation differences caused by the anisotropy of the propagation medium as differences in the properties of the seismic source itself, a non-uniform wave velocity model that reflects the spatial anisotropy of the rock mass needs to be introduced during the microseismic localization process.

[0068] Rock mass geological exploration data was obtained through preliminary drilling, geological mapping, and geophysical exploration. This data includes the bedding strike, dip angle, and spacing of the rock mass in the monitoring area. Bedding strike characterizes the direction of bedding extension in the horizontal plane and is a fundamental parameter for determining the relationship between the horizontal propagation path of microwaves and the rock mass structure. Bedding dip angle describes the spatial inclination of bedding planes and is an important basis for calculating the spatial orientation and normal vector of bedding planes. Bedding spacing quantifies the density of bedding distribution in the rock mass and is an important reference parameter for constructing spatially partitioned wave velocity models.

[0069] The coordinates of the microseismic sensors are collected during the installation and commissioning phase. High-precision positioning equipment is used to obtain the precise coordinates of each microseismic sensor in a unified monitoring coordinate system, providing a spatial coordinate basis for subsequent propagation path length calculation and propagation direction analysis.

[0070] First, the microseismic sensor array is started to carry out synchronous acquisition and obtain microseismic waveform data. The arrival time of the first arrival wave is extracted for each microseismic sensor channel. The arrival time of the first arrival wave is the time point at which the elastic wave generated by the microseismic event is first reliably identified by the microseismic sensor.

[0071] Based on the clearly defined bedding strike, dip angle, and spacing in the rock mass geological exploration data, the wave velocities perpendicular to and parallel to bedding planes are determined. The wave velocity perpendicular to bedding planes represents the equivalent propagation speed of elastic waves propagating approximately perpendicular to the bedding planes, while the wave velocity parallel to bedding planes represents the equivalent propagation speed of elastic waves propagating approximately parallel to the bedding planes. Both the wave velocities perpendicular to and parallel to bedding planes are obtained through controlled elastic wave testing and are used to characterize the main propagation velocity features under anisotropic conditions in the rock mass.

[0072] Using the unified monitoring coordinate system of the monitoring area as a reference, the horizontal plane is defined as the xy plane, and the z-axis is vertically upward, to establish a reference coordinate system.

[0073] Define bedding strike and bedding dip. Bedding strike is the azimuth angle of the intersection of the bedding plane and the xy-plane, measured clockwise from true north, ranging from 0 to 360 degrees. Bedding dip is the angle between the bedding plane and the xy-plane, ranging from 0 to 90 degrees. Bedding strike and bedding dip are the only independent parameters characterizing the spatial attitude of bedding. Bedding strike is the spatial extension direction of the intersection of the bedding plane and the horizontal plane of the reference coordinate system; bedding dip is the angle between the bedding plane and the horizontal plane of the reference coordinate system. Together, they determine the spatial orientation of the bedding plane.

[0074] The unit direction vector of the intersection line between the bedding plane and the xy plane is calculated based on the bedding strike, thus obtaining the strike vector. Based on the spatial component decomposition rules of trigonometric functions, the x-component of the strike vector is the cosine of the azimuth angle, the y-component is the sine of the azimuth angle, and the z-component is 0. Geologically, the dip line is a straight line on the bedding plane perpendicular to the strike line and pointing downhill. The bedding dip angle is the angle between the bedding plane and the horizontal plane, directly determining the steepness of the bedding plane's dip. As the downhill direction line on the bedding plane, the z-axis component of the dip line must reflect the degree of bedding plane dip. Therefore, the unit direction vector located within the bedding plane and perpendicular to the strike vector is calculated based on the bedding dip angle, thus obtaining the dip vector. The x-component of the dip vector is the product of the negative sine of the azimuth angle and the cosine of the bedding dip angle; the y-component is the product of the cosine of the azimuth angle and the cosine of the bedding dip angle; and the z-component is the negative sine of the bedding dip angle. Since the strike line vector and the dip line vector are not parallel vectors, the bedding normal vector is obtained by performing a cross product operation between them. The bedding normal vector is used to characterize the spatial normal direction of the bedding plane in the reference coordinate system.

[0075] Based on the bedding normal vector, known spatial points on the bedding plane are selected to construct the bedding plane spatial equation, which is used to describe the distribution of bedding in three-dimensional space.

[0076] Using a reference coordinate system as a unified spatial reference, a three-dimensional geological framework for the monitoring area is constructed, and this framework is divided into multiple three-dimensional units according to a preset grid size. Each three-dimensional unit serves as the basic spatial unit of the non-uniform wave velocity model, representing the main structural orientation characteristics of the rock mass within that area. The grid size is set based on the minimum geological structure thickness that the three-dimensional geological model needs to resolve, the spatial distribution density of geological data points, and the scale of the number of three-dimensional units that the wave velocity simulation calculation can handle. For example, in a monitoring area with well-developed bedding, if the minimum target rock layer thickness is 2 meters, the average spacing of the main geological exploration data points on the horizontal plane is 10 meters, and the maximum number of three-dimensional units that the wave velocity simulation calculation hardware can efficiently process is 1 million, then, after comprehensively considering the minimum rock layer thickness, the average spacing of the geological exploration data points, and the maximum number of three-dimensional units, the global basic grid size for the monitoring area can be preset as follows: a horizontal grid side length of 10 meters and a vertical grid thickness of 0.67 meters. Based on a horizontal grid side length of 10 meters and a vertical grid thickness of 0.67 meters, for a monitoring area 1000 meters long, 500 meters wide, and 200 meters high, the total number of three-dimensional elements obtained is approximately 1.5 million. This exceeds the upper limit of 1 million three-dimensional elements and requires optimization. The vertical grid thickness in non-core areas can be adjusted to 1 meter, or, after evaluation, the global horizontal grid side length can be adjusted to 12 meters. This would reduce the total number of three-dimensional elements to within the 1 million upper limit, satisfying both geological resolution requirements and computational efficiency constraints.

[0077] Based on the bedding plane spatial equation, the bedding strike and bedding dip angle within each three-dimensional unit are determined, so that each three-dimensional unit can reflect the main rock mass structural characteristics of its spatial location.

[0078] Within each three-dimensional unit, based on the wave velocities perpendicular to and parallel to the bedding planes, and combined with the bedding normal vector direction corresponding to that unit, the wave velocities perpendicular to and parallel to the bedding planes are mapped to the x, y, and z directions of the reference coordinate system using coordinate transformation methods. This yields the equivalent directional wave velocity components of the three-dimensional unit in the reference coordinate system. The specific process is based on the determination of bedding geometric features, the construction of a local coordinate system, and orthogonal coordinate transformation. Using a unified reference coordinate system for the monitoring area as a reference, the spatial orientation of the bedding plane in three-dimensional space is calculated according to the bedding strike and dip angle corresponding to the three-dimensional unit. The bedding normal vector is then calculated from the bedding strike and dip angle. The bedding normal vector is normalized to obtain a unit direction vector perpendicular to the bedding plane direction, used to characterize the spatial direction characteristics of elastic wave propagation along a direction approximately perpendicular to the bedding plane. Within the bedding plane, a unit direction vector is calculated based on the bedding strike direction. This unit direction vector is orthogonal to the bedding normal vector and is used to characterize the horizontal extension characteristics of the bedding plane. By performing a cross product operation between the bedding normal vector and the strike direction unit direction vector, the dip direction unit direction vector, located within the bedding plane and orthogonal to the strike direction unit direction vector, is obtained. This vector characterizes the spatial directional characteristics of another orthogonal direction within the bedding plane. The strike direction unit direction vector, dip direction unit direction vector, and bedding normal vector together constitute a local orthogonal coordinate system within the three-dimensional unit, used to describe the spatial directional relationships under bedding-controlled conditions. In this local orthogonal coordinate system, the wave velocity perpendicular to bedding is defined as the propagation velocity along the direction of the bedding normal vector, and the wave velocity parallel to bedding is defined as the propagation velocity along the plane containing the strike and dip direction unit direction vectors. It is assumed that the propagation velocities in all directions within the bedding plane are consistent, used to characterize the main anisotropic propagation characteristics of the rock mass within the three-dimensional unit under bedding-controlled conditions. Propagation characteristic expressions corresponding to the wave velocities perpendicular to and parallel to bedding are established in the local orthogonal coordinate system, making the propagation characteristics of the perpendicular and parallel bedding directions mathematically independent. The components of the unit direction vectors along the strike line, dip line, and bedding normal vector in the reference coordinate system are organized to construct an orthogonal coordinate transformation matrix from the local orthogonal coordinate system to the reference coordinate system. By applying the orthogonal coordinate transformation matrix to the propagation characteristic expression established in the local orthogonal coordinate system, the wave velocities perpendicular to and parallel to the bedding are mapped to the reference coordinate system, resulting in equivalent propagation characteristic descriptions along the x, y, and z directions in the reference coordinate system. After completing the coordinate transformation, the equivalent propagation velocity components corresponding to the x, y, and z directions are extracted from the equivalent propagation characteristic descriptions in the reference coordinate system to obtain the equivalent directional wave velocity components of the three-dimensional element in the reference coordinate system.The equivalent directional wave velocity component is used to approximate the equivalent propagation capability of elastic waves when they propagate along different spatial directions in the reference coordinate system within a three-dimensional element. Under the premise of keeping the anisotropic propagation characteristics of rock mass reflected by the wave velocity perpendicular to bedding and the wave velocity parallel to bedding unchanged, it provides a unified and queryable propagation velocity constraint for the construction of non-uniform wave velocity models and the numerical calculation of propagation time during microseismic location.

[0079] The equivalent directional wave velocity components of all three-dimensional elements are uniformly summarized, and spatial integration of these components is performed using finite difference numerical methods or finite element numerical methods to construct a non-uniform wave velocity model covering the entire monitoring area. The non-uniform wave velocity model stores the equivalent propagation velocity information corresponding to each spatial location in the form of a three-dimensional mesh. Each mesh element in the three-dimensional mesh corresponds to an equivalent directional wave velocity component of a three-dimensional element, providing propagation velocity constraints for calculating elastic wave propagation time during microseismic location.

[0080] By utilizing microseismic sensor coordinates, first-arrival wave arrival time, and a non-uniform wave velocity model, a constraint relationship between the source location and arrival time is established using a time-difference localization method. The three-dimensional coordinates of the source and the initial time are set as unknowns in the localization calculation. For each microseismic sensor involved in the calculation, a corresponding propagation time constraint equation is established to describe the relationship between the source location, the source occurrence time, and the first-arrival wave arrival time of the microseismic sensor.

[0081] In the location calculation process, based on the currently assumed three-dimensional coordinates of the seismic source, a straight path between the seismic source and the microseismic sensor is used as an engineering approximation of the elastic wave propagation path. Along this straight path, the three-dimensional mesh elements traversed by the path are sequentially sampled in the non-uniform wave velocity model. The equivalent directional wave velocity components at the corresponding spatial locations are extracted, and the equivalent propagation velocity along the straight path is calculated based on the direction information of the straight path in the reference coordinate system. By numerically integrating the equivalent propagation velocity at each spatial location along the path, the propagation time of the elastic wave from the seismic source to the microseismic sensor is numerically estimated. The straight path approximation is used to perform engineering modeling of the actual propagation process while ensuring computational efficiency.

[0082] By linearizing the propagation time constraint equations corresponding to all microseismic sensors, an overdetermined set of equations is constructed. The least squares fitting algorithm is used to solve the overdetermined set of equations to obtain the optimal solution for the three-dimensional coordinates and initial time of the seismic source.

[0083] The construction and least squares fitting solution process of the overdetermined system of equations includes:

[0084] Using microseismic sensor coordinates, first-arrival wave arrival time, and a non-uniform wave velocity model, the three-dimensional coordinates and initial time of the seismic source are solved. The three-dimensional coordinates of the source serve as the unknown spatial location of the source within the monitoring area, while the initial time serves as the unknown moment of the microseismic event. Together, the three-dimensional coordinates and initial time constitute the unknown vector in the location calculation.

[0085] For each microseismic sensor, a propagation time constraint is established based on the arrival time of the corresponding first arrival wave. The actual propagation time of the elastic wave from the seismic source to the microseismic sensor is obtained by subtracting the initial time from the arrival time of the first arrival wave of the corresponding microseismic sensor.

[0086] In the reference coordinate system, based on the currently assumed three-dimensional coordinates of the seismic source and the corresponding microseismic sensor coordinates, the straight-line distance from the seismic source to the microseismic sensor is calculated using the Euclidean distance formula. This straight-line distance is used as the approximate propagation distance of the corresponding propagation path of the microseismic sensor.

[0087] By combining the non-uniform wave velocity model, the corresponding equivalent propagation velocity information is extracted on the propagation path from the source to the microseismic sensor. Distance constraint equations are established between propagation distance, propagation time and equivalent propagation velocity to describe the propagation relationship that should be satisfied between propagation time and propagation distance under given three-dimensional coordinates of the source and initial time conditions.

[0088] Since the propagation distance has a nonlinear relationship with the three-dimensional coordinates of the seismic source, the distance constraint equation for each microseismic sensor is also a nonlinear equation. A first-order Taylor expansion is performed on the nonlinear distance constraint equation near the currently assumed three-dimensional coordinates of the seismic source and the initial time. This linearization process yields the linearized distance equation for the corresponding microseismic sensor.

[0089] A linearized distance equation is established for each microseismic sensor involved in the location calculation. When the number of microseismic sensors involved in the location calculation is greater than or equal to four, the number of linearized distance equations exceeds the number of unknowns in the unknown vector, thus forming an overdetermined system of equations.

[0090] All linearized distance equations are uniformly organized into matrix form, resulting in a system of linear equations. In this system, the unknown vectors include the three-dimensional coordinates of the source and the initial time. The coefficient matrix consists of the partial derivatives of each linearized distance equation with respect to the three-dimensional coordinates of the source and the initial time. The constant term vector is determined by the arrival time of the first wave, the equivalent propagation velocity, and the current assumed values ​​for each microseismic sensor.

[0091] The least squares fitting algorithm fits the residuals of all linearized distance equations in the overdetermined system of equations. The residual is defined as the difference between the product of the coefficient matrix and the unknown vector and the constant term vector, and is used to quantify the degree of inconsistency between the propagation time constraints of each microseismic sensor under the current three-dimensional coordinates of the seismic source and the initial time assumptions.

[0092] The objective function of residual sum of squares is constructed by squaring and summing all residuals. A least-squares fitting algorithm is then used to minimize this objective function, yielding the optimal solution for the unknown vector that minimizes the residual sum of squares. This optimal solution is obtained by calculating the inverse of the product of the transpose and the coefficient matrix, and then multiplying it sequentially by the transpose of the coefficient matrix and the constant term vector.

[0093] The optimal solution of the unknown vector obtained by the least squares fitting algorithm corresponds to the estimation results of the three-dimensional coordinates of the seismic source and the initial time. Based on the estimation results of the three-dimensional coordinates of the seismic source and the initial time, the propagation time residual corresponding to each microseismic sensor is recalculated, and all propagation time residuals are statistically analyzed. The residual error is calculated as an evaluation index of positioning accuracy.

[0094] When the residual error is less than the preset residual error threshold, the estimated results of the source's three-dimensional coordinates and initial time are deemed to meet the positioning accuracy requirements, and the source's three-dimensional coordinates and initial time are output as the final positioning result. When the residual error is greater than the preset residual error threshold, the currently estimated source's three-dimensional coordinates and initial time are used as new initial assumptions, the nonlinear distance constraint equation is re-linearized, and the least squares fitting solution process is repeated until the residual error meets the accuracy requirements or reaches the preset iteration termination condition. The residual error threshold is set based on the time acquisition accuracy of the microseismic sensor, the propagation time calculation accuracy of the non-uniform wave velocity model, and the engineering requirements for consistency of microseismic positioning time. For example, in microseismic monitoring engineering, if the maximum acquisition error of the first arrival time of the microseismic sensor is 2 milliseconds, the maximum propagation time error introduced by the non-uniform wave velocity model in the equivalent propagation time calculation process is 3 milliseconds, and the engineering application requires that the time deviation between the propagation time calculation results of each microseismic sensor and the first arrival time does not exceed 5 milliseconds, then after comprehensively considering the above sources of time error, the residual error threshold can be preset to 4 milliseconds.

[0095] The final output of the three-dimensional source coordinates is used to characterize the spatial location of the microseismic event within the monitoring area, and the final output of the initial time is used to characterize the actual occurrence time of the microseismic event on the time axis. The three-dimensional source coordinates and the initial time serve as the basic input parameters for the subsequent attenuation correction module and energy filling module.

[0096] The positioning results are quality controlled by calculating residual error. When the residual error is less than the preset residual error threshold, the current solution result is determined to be the effective three-dimensional coordinates and initial time of the seismic source. When the residual error is greater than the preset threshold, the input microseismic waveform data or non-uniform wave velocity model is verified and the solution is iterated again until the accuracy requirements are met.

[0097] The attenuation correction module extracts the original energy and path characteristics from the microseismic waveform data, and corrects the original energy based on the path characteristics combined with the attenuation law of stratification and the spherical diffusion effect to obtain the corrected source energy.

[0098] Based on the clearly defined bedding strike, dip angle, and other bedding distribution characteristics from the rock mass geological exploration data, 3 to 5 typical test areas are selected within the monitoring area. Typical test areas refer to regions that represent the overall rock mass characteristics of the monitoring area. Each area should focus on covering both perpendicular and parallel bedding directions, with multiple test points at different distances set up under each direction. Test points should avoid geologically anomalous areas such as fracture zones to ensure that the test results can characterize the general attenuation characteristics of the rock mass perpendicular and parallel to the bedding directions within the area. Specifically, perpendicular bedding direction refers to the direction where the angle between the propagation path of the test point and the normal vector of the bedding plane is less than or equal to 10°, and parallel bedding direction refers to the direction where the angle between the propagation path and the normal vector of the bedding plane is greater than or equal to 80°.

[0099] A controllable elastic wave excitation source is deployed at each test point, and the precise location of the excitation source and the standard excitation energy released are recorded. The precise location of the excitation source is determined using a total station, and the coordinates of the excitation source and the microseismic sensor are unified to the same monitoring coordinate system. Microseismic test signals after the elastic wave propagation are collected using microseismic sensors already deployed in the surrounding area. The microseismic test signals are then subjected to square integration to obtain the signal energy value of the corresponding channel. Test data for calculating the vertical and parallel bedding attenuation parameters is obtained, including the standard excitation energy, excitation source coordinates, microseismic sensor coordinates, and microseismic signal energy value.

[0100] Methods for performing square integral operations on microseismic test signals include:

[0101] For any microseismic sensor channel, the microseismic test signal collected under the excitation conditions at the corresponding test point is obtained. The microseismic test signal is a continuous vibration amplitude sequence recorded by the microseismic sensor during the sampling process. The change of vibration amplitude over time reflects the propagation characteristics of elastic waves in this channel.

[0102] Using the triggering time of the controllable elastic wave excitation source as a time reference, and combining it with the arrival time of the first arrival wave in the microseismic test signal, a preset effective signal time range is established. This effective signal time range is used to cover the period when the main energy of the elastic wave is concentrated, in order to avoid interference from environmental background noise and wake wave signals on the energy calculation results. The effective signal time range is set based on the rock mass layer structure parameters, interlayer wave impedance differences, elastic wave multipath propagation effects, background noise levels, and testing experience with similar rock mass layers. For example, in a hammer-driven microseismic test of sandstone-shale interbedded rock mass, if the propagation velocity of the elastic wave in sandstone is 4500 m / s and that in shale is 2500 m / s, as measured by cross-hole testing in the field, and the direct path length between the excitation source and the sensor is 15 meters, the calculated arrival time of the first arrival wave is 3.3 milliseconds, the average amplitude of the background noise is 0.08 volts, and a threshold of invalid signal is set for a signal amplitude less than twice the background noise amplitude, and the statistical value of the duration of the main energy of the elastic wave under the same working conditions in similar sandstone-shale interbedded rock masses is 30 to 50 milliseconds, then after comprehensively considering the characteristics of the layered rock mass and the test parameters, the effective signal time range can be preset to an interval starting from the arrival time of the first arrival wave of 3.3 milliseconds and ending at the arrival time of the first arrival wave plus 40 milliseconds, i.e., 43.3 milliseconds.

[0103] Within the effective signal time range, the microseismic test signal is subjected to DC removal processing. Specifically, the average value of all vibration amplitudes within the effective signal time range is calculated, and the vibration amplitude at each moment is subtracted from the average value to eliminate the influence of possible DC bias and low-frequency drift in the test signal on the energy calculation results.

[0104] After the DC removal process is completed, the micro-vibration test signal within the effective signal time range is squared point by point. The vibration amplitude at each moment is squared to obtain an energy sequence reflecting the instantaneous energy intensity at each moment. The squaring process is used to ensure that the energy calculation results are not affected by the positive or negative change of the vibration direction.

[0105] After squaring, the energy sequence is integrated over the effective signal time range. Specifically, the instantaneous energy at each moment is accumulated according to the signal sampling time interval to obtain the signal energy value received by the corresponding microseismic sensor channel under the corresponding test point conditions.

[0106] Based on pre-collected test data, the vertical bedding attenuation parameters and parallel bedding attenuation parameters are calculated at corresponding distances. For all propagation distances and their corresponding vertical bedding attenuation parameters at the vertical bedding azimuth, the data are sorted in ascending order of propagation distance to form a vertical bedding attenuation parameter table. Similarly, for all propagation distances and their corresponding parallel bedding attenuation parameters at the parallel bedding azimuth, the data are sorted in ascending order of propagation distance to form a parallel bedding attenuation parameter table.

[0107] Methods for calculating the vertical bedding attenuation parameter and the parallel bedding attenuation parameter at the corresponding distances include:

[0108] For test points at different distances under the vertical bedding orientation, test data with an angle of less than or equal to 10° between the propagation path and the normal vector of the bedding plane are selected for each test point. The energy values ​​of the microseismic signal are extracted from three repeated excitations for each test point, and the average value of the microseismic signal energy value is calculated to obtain the average signal energy value.

[0109] Based on the coordinates of the excitation source and the microseismic sensor, the propagation distance is calculated using the Euclidean distance formula. The ratio of the standard excitation energy to the average signal energy is calculated to obtain the energy attenuation ratio. The natural logarithm of the energy attenuation ratio is then performed to obtain the logarithmic energy attenuation at the corresponding propagation distance. The logarithmic energy attenuation is used as the vertical bedding attenuation parameter at the corresponding propagation distance. The vertical bedding attenuation parameter is a dimensionless quantity and is only used to characterize the exponential energy attenuation characteristics caused by the rock mass medium during the propagation of elastic waves in the vertical bedding direction.

[0110] For test points at different distances under the orientation of parallel bedding, test data with an angle greater than or equal to 80° between the propagation path of each test point and the normal vector of the bedding plane are selected. Using the same method as obtaining the vertical bedding attenuation parameter at the corresponding distance, the parallel bedding attenuation parameter at the corresponding distance is calculated.

[0111] For the microseismic waveform data acquired by each microseismic sensor, the raw energy corresponding to a single channel is calculated through square integration. Path characteristics of the propagation path from the seismic source to each microseismic sensor are extracted, including the actual propagation distance and path direction angle. Based on the three-dimensional coordinates of the seismic source and the precise coordinates of each microseismic sensor, the total actual length of the entire propagation path is calculated using the Euclidean distance formula to obtain the total propagation distance. For each rock layer, the coordinates of the start and end points of the propagation path within the corresponding rock layer are located, and the actual length between the start and end points is calculated using the Euclidean distance formula to obtain the actual propagation distance of the corresponding rock layer. The sum of the actual propagation distances of each rock layer equals the total propagation distance. The normal vector of the corresponding bedding plane is calculated based on the spatial equation of the bedding plane of each rock layer. For each rock layer, the angle between the spatial vector of the propagation path within the corresponding rock layer and the normal vector of the bedding plane is calculated. The acute angle corresponding to the angle between the spatial vector and the normal vector of the bedding plane is taken as the path direction angle of the corresponding rock layer. The path direction angle is used to reflect the degree to which the direction of elastic wave propagation is controlled by the anisotropy of the rock mass structure. The cosine value of the vector angle is obtained by calculating the dot product of the spatial vector of the propagation path and the normal vector of the bedding plane, and then converting it into an angle value through inverse cosine operation.

[0112] Because the rock mass bedding structure exhibits mechanical anisotropy, it directly affects the degree of energy loss during elastic wave propagation. The original energy data needs to be specifically corrected for the bedding attenuation law, which can accurately distinguish the energy attenuation differences in different propagation directions and avoid misjudgment of the source energy caused by directional factors.

[0113] Based on path characteristics, initial attenuation parameters are synthesized for each rock layer. The actual propagation distance is decomposed into effective propagation distances perpendicular to bedding planes and parallel to bedding planes based on the path direction angle corresponding to the rock layer. These effective propagation distances perpendicular to bedding planes and parallel to bedding planes describe the relative projection scale of the propagation path in different structural control directions, serving only as the directional weighting basis for attenuation parameter matching. Within the current rock layer, the actual propagation distance of the corresponding rock layer is considered as a complete path integration interval, including weights for perpendicular and parallel bedding planes, used to characterize the attenuation contribution ratio of the corresponding propagation path in different structural control directions. Based on the effective propagation distances perpendicular to and parallel to bedding planes, matching is performed from the vertical and parallel bedding plane attenuation parameter tables to obtain target vertical and parallel bedding plane attenuation parameters. Based on these target vertical and parallel bedding plane attenuation parameters, initial attenuation parameters for each rock layer are synthesized. After completing the initial attenuation parameter synthesis process for all rock layers, the overall attenuation parameter of the entire path can be obtained by combining the initial attenuation parameters of each layer, thereby ensuring that the overall attenuation parameter is accurately matched with the actual propagation state.

[0114] Please see Figure 3 As shown, the method for synthesizing initial attenuation parameters for each rock mass layer includes:

[0115] Based on the path direction angle corresponding to the current rock stratum, the cosine value of the path direction angle is calculated as the initial weighting coefficient for the perpendicular bedding direction, and the sine value of the path direction angle is calculated as the initial weighting coefficient for the parallel bedding direction. The initial weighting coefficients for the perpendicular bedding direction and the initial weighting coefficients for the parallel bedding direction are normalized. Specifically, the initial weighting coefficients for the perpendicular bedding direction and the initial weighting coefficients for the parallel bedding direction are divided by their sum to obtain the weighting coefficients for the perpendicular bedding direction and the weighting coefficients for the parallel bedding direction. The sum of the normalized weighting coefficients for the perpendicular bedding direction and the weighting coefficients for the parallel bedding direction is one, which is used to characterize the relative attenuation contribution ratio of the propagation path in the perpendicular bedding direction and the parallel bedding direction, thereby ensuring that the weighting process of attenuation effect in different directions does not introduce additional energy amplification or reduction.

[0116] The effective propagation distance in the vertical bedding direction is obtained by multiplying the actual propagation distance of the current rock stratum by the initial weighting coefficient in the direction perpendicular to the bedding; the effective propagation distance in the parallel bedding direction is obtained by multiplying the actual propagation distance of the current rock stratum by the initial weighting coefficient in the direction parallel to the bedding; the equivalent propagation distance in the vertical bedding direction is used to characterize the equivalent attenuation scale of the propagation path in the direction controlled by the vertical bedding while keeping the actual propagation path length unchanged; the equivalent propagation distance in the parallel bedding direction is used to characterize the equivalent attenuation scale of the propagation path in the direction controlled by the parallel bedding.

[0117] From the vertical bedding attenuation parameter table and the parallel bedding attenuation parameter table, the propagation distances consistent with the effective propagation distance in the vertical bedding direction are found, and the corresponding vertical bedding attenuation parameters are extracted as target vertical bedding attenuation parameters. Similarly, the propagation distances consistent with the effective propagation distance in the parallel bedding direction are found, and the corresponding parallel bedding attenuation parameters are extracted as target parallel bedding attenuation parameters. The target vertical bedding attenuation parameters and the target parallel bedding attenuation parameters are added together to obtain the initial attenuation parameters of the current rock layer. The initial attenuation parameters are the equivalent representation of the propagation path attenuation effect within the corresponding rock layer after path integration, and are used to characterize the cumulative influence of the corresponding rock layer on the attenuation of elastic wave energy.

[0118] If the effective propagation distance perpendicular to the bedding direction or the effective propagation distance parallel to the bedding direction exceeds the range of the vertical bedding attenuation parameter table or the parallel bedding attenuation parameter table, for the effective propagation distance in the vertical bedding direction or the effective propagation distance in the parallel bedding direction that exceeds the range, select the vertical or parallel bedding attenuation parameter corresponding to the known propagation distance that is closest to the corresponding effective propagation distance from the corresponding vertical bedding attenuation parameter table or the parallel bedding attenuation parameter table.

[0119] Methods for integrating the initial attenuation parameters of each layer include:

[0120] Calculate the ratio of the actual propagation distance to the total propagation distance to obtain the attenuation weight;

[0121] Multiply the initial attenuation parameter by the attenuation weight to obtain the attenuation contribution value;

[0122] The attenuation contribution values ​​corresponding to each rock layer are added together to obtain the comprehensive attenuation parameter.

[0123] Based on comprehensive attenuation parameters, single-channel energy correction is completed, and the true source energy is inferred. Since the propagation of microseismic elastic waves in rock mass can be approximated as spherical diffusion, the corresponding geometric attenuation calculation adopts the classic spherical diffusion attenuation factor formula in the field of microseismic monitoring:

[0124] ;

[0125] Where r is the total propagation distance, the spherical diffusion attenuation factor formula is used to quantify the degree of energy attenuation caused by the diffusion effect.

[0126] Substituting the total propagation distance into the spherical diffusion attenuation factor formula, we obtain the geometric attenuation factor;

[0127] Multiplying the geometric attenuation factor by the comprehensive attenuation parameter yields the total attenuation compensation factor, which simultaneously covers both geometric diffusion attenuation and rock mass anisotropic attenuation. The reason for using the product form is that, under the condition that the amplitude of the signal energy received by the microseismic sensor is affected by both geometric diffusion caused by propagation distance and directional attenuation caused by rock mass anisotropy, using only a uniform attenuation model would confuse the distance effect and the azimuth effect, leading to distortion of the source energy inversion. Therefore, it is necessary to compensate for both geometric attenuation and anisotropic attenuation simultaneously through the total attenuation compensation factor in the product form in order to restore the true source energy.

[0128] Divide the original energy by the total attenuation compensation factor to obtain the corrected source energy of the corresponding channel.

[0129] The energy filling module calculates a correction factor based on the corrected source energy and the original energy, calculates a confidence coefficient based on the correction factor, and selects representative energies based on the confidence coefficient.

[0130] Please see Figure 4 As shown, the ratio of the corrected source energy to the original energy is calculated to obtain the basic correction coefficient, which characterizes the adjustment range of the attenuation correction on the original energy.

[0131] The ratio of the basic correction coefficient to the geometric attenuation coefficient is calculated to obtain the correction factor, which reflects the degree of matching between the corrected source energy and the theoretical energy after deducting the geometric attenuation from the original energy. It is the core quantitative indicator for judging the authenticity of the data.

[0132] Based on the correction factors corresponding to each channel, the confidence coefficient of the single-channel corrected source energy data is calculated. The confidence coefficient is the core quantitative indicator that characterizes the reliability of single-channel data and provides a clear reliability basis for subsequent representative energy selection.

[0133] Methods for calculating the confidence coefficient of single-channel corrected source energy data include:

[0134] Calculate the mean and standard deviation of the correction factors for all channels to obtain the corrected mean and corrected variance;

[0135] The correction factor for each channel is standardized, for example, by Z-score standardization, to obtain the standardized correction factor for each channel, thus eliminating the influence of different dimensions and distribution ranges of the correction factor.

[0136] Since the standardized correction factor is an arbitrary real number, it cannot be directly used as the confidence coefficient. The standardized correction factor must be mapped to the [0,1] interval to obtain the confidence coefficient. For example, the Sigmoid function can be used, with the natural constant as the base and the negative standardized correction factor as the power, to calculate the negative exponent value. One is added to the negative exponent value to obtain the base negative exponent cumulative value. One is then divided by the base negative exponent cumulative value to obtain the confidence coefficient. The closer the confidence coefficient is to one, the higher the consistency between the corresponding channel's correction factor and the overall channel correction factor, and the stronger the data reliability. If the standard deviation of all calculated channel correction factors is zero, it indicates that all channel correction factors are completely consistent, and in this case, the confidence coefficient is uniformly set to one.

[0137] The maximum value among the source energies of all microseismic sensor channels is selected as the representative energy of the corresponding microseismic event. This avoids underestimation of energy caused by local signal attenuation and sensor deviation, accurately pinpoints the maximum energy level of the microseismic event, and provides a risk upper limit reference for disaster risk early warning.

[0138] The corrected source energy corresponding to all microseismic sensor channels output by the attenuation correction module is extracted, and the corrected source energy with a confidence coefficient less than the preset confidence threshold is removed to avoid abnormal data interfering with the screening results. The confidence threshold is preset based on the accuracy of the microseismic sensor and the actual data quality requirements of the project. For example, in the early warning of tunnel surrounding rock stability, a high-precision microseismic sensor with a positioning error of less than or equal to 3cm and an initial arrival wave arrival time acquisition error of less than or equal to 1ms is used. Since the project needs to rely on highly reliable data to avoid false early warning judgments, the confidence threshold is preset to 0.7.

[0139] The corrected source energies after removal are sorted in descending order, and the corrected source energy with the highest value at the top of the list is selected as the representative energy of the microseismic event. If multiple channels have the same maximum corrected source energy, the corrected source energy with the highest confidence coefficient is selected as the representative energy, ensuring that the representative energy has both maximum value characteristics and high reliability. The corrected source energies after single-channel correction are all the true energy of the source, but different channels correspond to different propagation paths. The maximum corrected energy is the true upper limit of the energy released by the source, directly corresponding to the maximum impact force of the microseismic event on the rock mass, and is also the most dangerous signal of rock mass instability. Selecting the corrected source energy with the largest value as the representative energy can avoid the problem of high risk being masked by averaging multiple channel data, ensuring that even if strong energy in a certain direction is corrected and restored due to propagation attenuation, it can still be accurately captured and will not be hidden due to channel differences.

[0140] The serialization module serializes microseismic events based on initial time, representative energy, and confidence coefficients to obtain the microseismic event sequence feature tensor.

[0141] All the microseismic events that occurred in succession were sorted according to their initial time sequence;

[0142] For each sorted microseismic event, the representative energy and the corresponding confidence coefficient are extracted, and the initial time difference between the microseismic event and the previous microseismic event is calculated as the time interval; the representative energy, confidence coefficient and time interval together constitute the microseismic characteristics of the corresponding microseismic event.

[0143] The microseismic features are standardized by normalization to eliminate the dimensional differences between different microseismic features. For example, the maximum and minimum value normalization method is used to obtain standardized feature vectors, ensuring that all standardized feature values ​​are in the range of 0 to 1.

[0144] The standardized feature vectors are arranged sequentially according to their initial time to form a feature vector set. The feature vector set is then transformed into a microseismic event sequence feature tensor through tensor reconstruction.

[0145] Methods for transforming a set of feature vectors into a feature tensor of a microseismic event sequence through tensor reconstruction include:

[0146] A unique association identifier is added to each generated standardized feature vector. This identifier directly corresponds to the initial time of the microseismic event to which it belongs, ensuring that each standardized feature vector is accurately bound to the time attribute of the original microseismic event and avoiding misalignment between features and time in the subsequent sorting process.

[0147] Based on the initial time, all standardized feature vectors bound to time identifiers are sorted in ascending order, starting from the standardized feature vector corresponding to the earliest microseismic event and proceeding sequentially to the standardized feature vector corresponding to the latest microseismic event, forming an ordered set of feature vectors. The set of feature vectors is a linear list structure, and each element in the list is a three-dimensional standardized feature vector of a single microseismic event. The vector contains three standardized feature values ​​representing energy, confidence coefficient, and time interval, respectively.

[0148] Based on the structure of the feature vector set, the three dimensions of the microseismic event sequence feature tensor are defined. The first dimension is the sequence length, which corresponds to the total number of standardized feature vectors in the feature vector set, i.e., the total number of microseismic events participating in the sorting. The second dimension is the feature dimension, which corresponds to the three features included in the standardized feature vector: representative energy, confidence coefficient, and time interval. The third dimension is the batch dimension, with a default value of one. If multiple microseismic event sequences need to be processed simultaneously, the value of the batch dimension can be adjusted according to the actual number of microseismic event sequences to be processed, in order to adapt to the batch calculation requirements of the microseismic energy prediction model.

[0149] Each standardized feature vector in the ordered set of feature vectors is sequentially filled into a tensor frame of pre-defined dimensions to form a well-structured microseismic event sequence feature tensor. The microseismic event sequence feature tensor fully preserves the initial temporal sequence of microseismic events, as well as the three core feature information of each microseismic event: representative energy, confidence coefficient, and time interval. It can be directly used as input data for the microseismic energy prediction model.

[0150] The microseismic prediction module inputs the feature tensor of the microseismic event sequence after mask filling and alignment into the pre-trained microseismic energy prediction model to obtain the maximum energy prediction value.

[0151] Controllable elastic wave excitation sources are pre-distributed evenly at key test points in the monitoring area, and a microseismic sensor array is simultaneously deployed throughout the monitoring area. By precisely controlling the release of elastic waves from the excitation sources, the standard energy released is recorded to obtain the true value of the actual source energy. Simultaneously, the microseismic sensor array collects historical microseismic waveform data during the propagation of the corresponding elastic waves. The historical microseismic event sequence dataset includes historical microseismic event sequence data and corresponding microseismic result data. The historical microseismic event sequence data is the historical microseismic event sequence feature tensor obtained after processing the collected historical microseismic waveform data through a microseismic localization module, an attenuation correction module, an energy filling module, and a serialization module. The microseismic result data is the true value of the actual source energy, which serves as the basis for calculating the loss function during model training.

[0152] Training methods for microseismic energy prediction models include:

[0153] The sequence length of historical microseismic event data varies depending on the frequency and duration of the events. To adapt to batch computing requirements, a mask-filling alignment method is used to process the historical microseismic event sequence data. Pre-defined, meaningless, specific identifiers are padded at the end of shorter sequences, and a corresponding binary mask matrix is ​​generated. In the binary mask matrix, 1 indicates valid microseismic feature data, and 0 indicates the padded portion. This ultimately yields an aligned feature tensor. The mask-filling alignment method not only fully preserves the original unequal length characteristics of each sequence but also achieves accurate differentiation between valid microseismic feature data and the padded portion through the binary mask matrix. This ensures that the microseismic energy prediction model can stably receive input data and avoids interference from the padded portion during training.

[0154] The standardized historical microseismic event sequence dataset is divided into a training set, a validation set, and a test set in a 7:2:1 ratio. The training set is augmented by random rotation of the coordinate system to simulate the propagation characteristics of microseismic events under different anisotropic geological conditions, thereby improving the generalization ability of the microseismic energy prediction model. The validation set is used to adjust the model hyperparameters in real time and monitor overfitting during the training process. The test set strictly uses the latest data split in chronological order to avoid interference from future data on model evaluation and ensure the fairness and authenticity of the evaluation results.

[0155] A microseismic energy prediction model is constructed using a bidirectional long short-term memory neural network based on an attention mechanism. The microseismic energy prediction model includes an input layer, a feature extraction layer, an attention mechanism layer, and an output layer.

[0156] The input layer is used to receive the sequence of aligned feature tensors.

[0157] The feature extraction layer consists of a forward long short-term memory (LSTM) network layer and a backward LSM network layer. The forward LSM network layer traverses the microseismic events in ascending chronological order, from the earliest to the latest microseismic event. This traversal aligns the feature tensor sequence to capture the forward temporal evolution pattern, which refers to the sequential influence and overall trend between features of earlier and later microseismic events. The backward LSM network layer traverses the microseismic events in reverse chronological order. That is, from the latest microseismic event to the earliest microseismic event, the alignment feature tensor sequence is traversed to capture the inverse dependency. The inverse dependency refers to the precursor correlation between the features of the later microseismic event and the features of the earlier microseismic event. The forward long short-term memory network layer and the backward long short-term memory network layer output the forward temporal feature sequence and the inverse temporal feature sequence, respectively. The forward temporal feature sequence and the inverse temporal feature sequence are spliced ​​and fused in the feature dimension to form a comprehensive feature sequence containing bidirectional temporal information, which is then passed to the attention mechanism layer.

[0158] The attention mechanism layer correlates the comprehensive feature sequence with the confidence coefficients embedded in the corresponding aligned feature tensor sequence. The attention weight is obtained by calculating the correlation between the comprehensive feature sequence and the corresponding confidence coefficient. For example, the correlation score between the comprehensive feature vector of each microseismic event and the corresponding confidence coefficient is calculated using dot product operation. The correlation scores of all microseismic events are normalized to obtain the attention weight. The comprehensive feature sequence and the attention weight are multiplied element-wise, and then all weighted feature vectors are summed to obtain the final weighted feature vector. At the same time, the attention mechanism layer combines the binary mask matrix passed from the input layer to assign zero weight to the features of the padding part, ensuring that the padding data does not participate in the feature fusion process.

[0159] The output layer consists of a fully connected layer that receives the weighted feature vector output from the attention mechanism layer. Through linear transformation, the high-dimensional feature vector is mapped to a single scalar output, which is the maximum energy prediction value of microseismic events within a future time interval, such as predicting the maximum energy value of microseismic events in the next hour.

[0160] A weighted loss function is used to optimize the microseismic energy prediction model, and a comprehensive weight is constructed by coupling risk-oriented weights and physical confidence weights.

[0161] To differentiate the penalty for microseismic events of different energy levels, the risk-oriented weights are calculated using the following formula:

[0162] ;

[0163] in, Risk-oriented weighting This represents the true value of the actual seismic source energy. Based on the weights.

[0164] In practical seismology, the Richter scale magnitude, based on the logarithmic relationship of seismic energy, directly correlates with the potential destructive power of seismic events. Taking the logarithm of the actual source energy value aligns with the corresponding physical meaning. Furthermore, the energy of microseismic events typically spans multiple orders of magnitude. Directly using energy might lead to weights dominating extremely high-energy events while being insensitive to low-energy events. Taking the logarithm compresses the scale, making the differences in weights between events of different energy levels smoother. Adding 1 to the logarithm ensures that risk-oriented weights are still defined even when the actual source energy value is zero.

[0165] The base weights are determined based on the distribution characteristics of low-energy events in the historical microseismic event sequence dataset. Low-energy events are defined as microseismic events whose actual source energy is less than 50% of the relevant engineering safety threshold. The average energy of the actual source energy of low-energy events in the historical microseismic event sequence dataset is statistically analyzed, and the base weights are preset to 0.3 to 0.5 times the average energy of the actual source energy of low-energy events. This ensures that the base weights provide sufficient penalty attention for low-energy events during training without interfering with the focus on learning high-energy catastrophic events due to excessive weights. High-energy catastrophic events, defined as microseismic events whose actual source energy reaches or exceeds the relevant engineering safety threshold, have a risk of triggering rock mass instability disasters. When the microseismic energy prediction model has a large prediction deviation for high-energy catastrophic events, the loss value will increase significantly, thus forcing the microseismic energy prediction model to focus on learning the precursor features of high-energy events and avoiding missed detections. The base weights, on the other hand, ensure the training focus on low-energy events and prevent the microseismic energy prediction model from ignoring the evolutionary patterns of low-energy events.

[0166] The physical confidence weight adjusts the penalty intensity based on the overall confidence level of the input aligned feature tensor, which is the average of the confidence coefficients of all microseismic events. If the overall confidence is low, it indicates that the proportion of noisy data in the aligned feature tensor is high. In this case, the penalty coefficient is reduced to prevent the microseismic energy prediction model from forcibly fitting noisy data, which would lead to a decrease in generalization ability.

[0167] A pre-defined adjustment factor is introduced, and multiplied by the physical uncertainty to obtain the actual penalty term. Subtracting the actual penalty term from 1 constructs a linear adjustment coefficient acting on the risk-oriented weight. Multiplying the risk-oriented weight by the linear adjustment coefficient yields the comprehensive weight. The comprehensive weight couples the risk-oriented weight and the physical confidence weight through a product. Simultaneously, it linearly combines the physical confidence weight (reflecting the reliability of the alignment feature tensor sequence) with the scenario-appropriate adjustment factor. The influence of the adjustment factor on the physical confidence weight is used to regulate the comprehensive weight, ensuring that the comprehensive weight retains the risk-oriented weight's differentiated focus on energy levels while incorporating the physical confidence weight's judgment on the reliability of the alignment feature tensor. The adjustment factor is directly preset to a fixed value based on the safety level requirements of the microseismic monitoring scenario; specifically, it is preset to 0.6 for key protection scenarios and 0.4 for routine protection scenarios. Key protection scenarios refer to areas with poor rock mass stability, dense engineering personnel, or severe disaster consequences, which require stronger physical confidence weighting to improve the model's learning accuracy on reliable data; conventional protection scenarios refer to areas with good rock mass stability and low disaster risk, where the influence of physical confidence weighting can be appropriately weakened.

[0168] The method for calculating the deviation between the predicted results and the actual source energy values ​​based on comprehensive weights is as follows: Mean squared error is used. The difference between the maximum predicted energy value and the actual source energy value for each historical microseismic event sequence is calculated to obtain the loss difference. This loss difference is then squared to obtain the loss variance. The loss variance is multiplied by the comprehensive weights corresponding to the historical microseismic event sequence data to obtain the weighted loss value. The average of these weighted loss values ​​is then calculated to obtain the overall loss value, which is the quantified result of the deviation between the predicted results and the actual source energy values. Based on this overall loss value, the parameters of all layers in the microseismic energy prediction model are iteratively updated using a gradient descent algorithm. The parameters of all layers include all learnable parameters contained in the input layer, feature extraction layer, attention mechanism layer, and output layer of the microseismic energy prediction model. This guides the iterative optimization of the microseismic energy prediction model parameters, continuously improving the model's prediction accuracy.

[0169] Configure a suitable optimizer, learning rate decay strategy, and batch size for the microseismic energy prediction model to balance training efficiency and stability. For example, the AdamW optimizer can be configured with weight decay set to... The learning rate is set to... It is paired with cosine annealing attenuation; the batch size is set to 64 to ensure training efficiency and stability.

[0170] The microseismic event sequence feature tensor is input into the trained microseismic energy prediction model to obtain the maximum energy prediction value. The maximum energy prediction value is then directly compared with a preset engineering safety threshold. The engineering safety threshold is set comprehensively based on the rock mass stability requirements of the monitoring scenario, the engineering protection capacity, and historical disaster data. If the maximum energy prediction value is greater than or equal to the engineering safety threshold, a rock mass instability early warning is generated, which includes suggestions for engineering personnel to take reinforcement measures, adjust construction plans, or carry out emergency evacuation. If the maximum energy prediction value is less than the engineering safety threshold, a safety warning is output, and subsequent changes in the microseismic event sequence are continuously monitored.

[0171] Example 2:

[0172] Please see Figure 2 As shown, this embodiment provides a microseismic time series prediction method based on energy filling, including:

[0173] A non-uniform wave velocity model was built based on the collected rock geological exploration data. The time difference positioning analysis of the non-uniform wave velocity model was performed by the collected microseismic waveform data and microseismic sensor coordinates to obtain the three-dimensional coordinates of the source and the initial time.

[0174] The raw energy and path characteristics of microseismic waveform data are extracted. Based on the path characteristics, bedding attenuation law and spherical diffusion effect, the raw energy is corrected to obtain the corrected source energy.

[0175] The correction factor is calculated based on the corrected source energy and the original energy. The confidence coefficient is calculated based on the correction factor. The corrected source energy is then selected based on the confidence coefficient to obtain the representative energy.

[0176] Based on initial time, representative energy, and confidence coefficient, a feature tensor for the microseismic event sequence is constructed.

[0177] The maximum energy prediction value is obtained by inputting the feature tensor of the microseismic event sequence into the pre-trained microseismic energy prediction model.

[0178] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

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

Claims

1. A microseismic time series prediction method based on energy filling, characterized in that, include: A non-uniform wave velocity model was built based on the collected rock geological exploration data. The time difference positioning analysis of the non-uniform wave velocity model was performed by the collected microseismic waveform data and microseismic sensor coordinates to obtain the three-dimensional coordinates of the source and the initial time. The raw energy and path characteristics of microseismic waveform data are extracted. Based on the path characteristics, bedding attenuation law and spherical diffusion effect, the raw energy is corrected to obtain the corrected source energy. The correction factor is calculated based on the corrected source energy and the original energy. The confidence coefficient is calculated based on the correction factor. The corrected source energy is then selected based on the confidence coefficient to obtain the representative energy. Based on initial time, representative energy, and confidence coefficient, a feature tensor for the microseismic event sequence is constructed. The maximum energy prediction value is obtained by inputting the feature tensor of the microseismic event sequence into the pre-trained microseismic energy prediction model.

2. The microseismic time series prediction method based on energy filling according to claim 1, characterized in that, Methods for correcting the original energy include: Path characteristics include total propagation distance, actual propagation distance of each rock layer, and path direction angle; Based on path characteristics and bedding attenuation laws, the initial attenuation parameters of each rock mass are calculated according to the path direction angle. The initial attenuation parameters are weighted and summed using the ratio of the actual propagation distance of each rock mass to the total propagation distance to obtain the comprehensive attenuation parameters. Based on the total propagation distance and spherical diffusion effect, the geometric attenuation factor is calculated using the classic spherical diffusion attenuation factor formula; Multiplying the geometric attenuation factor by the comprehensive attenuation parameter yields the total attenuation compensation factor; Divide the original energy by the total attenuation compensation factor to obtain the corrected source energy.

3. The microseismic time series prediction method based on energy filling according to claim 2, characterized in that, Methods for calculating the initial attenuation parameters of each rock mass based on path direction angle decomposition include: The projections of the path direction angle onto the preset vertical and parallel bedding directions are used as the initial weight coefficients for the vertical and parallel bedding directions, respectively. The initial weighting coefficients for the perpendicular and parallel bedding directions are normalized. Based on the normalized initial weight coefficients for the vertical bedding direction and the parallel bedding direction, the actual propagation distance is decomposed into the effective propagation distance in the vertical bedding direction and the effective propagation distance in the parallel bedding direction. For each rock mass layer traversed by the propagation path, the target vertical bedding attenuation parameter corresponding to the effective propagation distance in the vertical bedding direction is extracted from the pre-constructed vertical bedding attenuation parameter table, and the target parallel bedding attenuation parameter corresponding to the effective propagation distance in the parallel bedding direction is extracted from the pre-constructed parallel bedding attenuation parameter table. The initial attenuation parameter of the current rock layer is obtained by adding the target vertical bedding attenuation parameter to the target parallel bedding attenuation parameter.

4. The microseismic time series prediction method based on energy filling according to claim 2, characterized in that, Methods for calculating correction factors based on corrected source energy and original energy include: Calculate the ratio of the corrected source energy to the original energy to obtain the basic correction coefficient; Divide the basic correction factor by the geometric attenuation factor to obtain the correction factor.

5. The microseismic time series prediction method based on energy filling according to claim 4, characterized in that, Methods for calculating confidence coefficients based on correction factors include: The correction factor for each channel is standardized to obtain the standardized correction factor for each channel. The standardized correction factor is mapped to the [0,1] interval using a preset mapping function to obtain the confidence coefficient.

6. The microseismic time series prediction method based on energy filling according to claim 5, characterized in that, Methods for screening and correcting earthquake source energy based on confidence coefficients include: Eliminate corrected seismic source energies with confidence coefficients lower than the preset confidence threshold; The corrected source energies after elimination are sorted in descending order, and the maximum corrected source energy is selected as the representative energy of the corresponding microseismic event.

7. The microseismic time series prediction method based on energy filling according to claim 1, characterized in that, The training process of the microseismic energy prediction model uses a weighted loss function to optimize the microseismic energy prediction model; Collect a dataset of historical microseismic event sequences, which includes the confidence coefficients and true values ​​of the actual source energy of the historical microseismic events. The actual source energy is incremented by 1 and then logarithmically transformed to obtain the potential damage value. The potential damage value is added to the preset base weight to obtain the risk-oriented weight. Calculate the mean of the confidence coefficients corresponding to all historical microseismic events to obtain the physical confidence weights; The risk-oriented weight is coupled with the physical confidence weight to obtain a comprehensive weight; The overall weight is used as the weight of the weighted loss function.

8. The microseismic time series prediction method based on energy filling according to claim 7, characterized in that, Methods for coupling risk-oriented weights with physical confidence weights include: Physical uncertainty is defined as 1 minus the physical confidence weight; The actual penalty term is obtained by multiplying the preset adjustment factor by the physical uncertainty; Subtract the actual penalty term from 1 to obtain the linear adjustment coefficient; The risk-oriented weight is multiplied by the linear adjustment coefficient to obtain the comprehensive weight.

9. The microseismic time series prediction method based on energy filling according to claim 1, characterized in that, Methods for constructing non-uniform wave velocity models based on collected rock mass geological exploration data include: A reference coordinate system is established with the horizontal plane of the monitoring area as the xy plane and the z-axis pointing vertically upwards, thus constructing a three-dimensional geological framework. The geological exploration data of the rock mass includes the strike and dip of the bedding planes; the unit vector of the intersection of the strike and the horizontal plane is taken as the strike vector; the unit vector perpendicular to the strike vector is selected as the dip vector, and the z-axis component of the dip vector corresponds to the dip angle of the bedding planes; the bedding normal vector is obtained by the cross product operation of the strike vector and the dip vector. The three-dimensional geological framework is divided into multiple three-dimensional units according to a preset grid size. Within each three-dimensional unit, a local orthogonal coordinate system is formed by the strike line vector, dip line vector, and bedding normal vector. In a local orthogonal coordinate system, the pre-tested vertical and parallel bedding wave velocities are mapped to the x, y, and z directions of the reference coordinate system through coordinate transformation. The equivalent directional wave velocity components of the corresponding three-dimensional elements in the reference coordinate system are obtained. The equivalent directional wave velocity components of all three-dimensional elements are then combined to form a non-uniform wave velocity model.

10. A microseismic time series prediction system based on energy filling, used to implement the microseismic time series prediction method based on energy filling as described in any one of claims 1-9, characterized in that, include: The microseismic location module builds a non-uniform wave velocity model based on the collected rock geological exploration data. It performs time difference location analysis on the non-uniform wave velocity model by collecting microseismic waveform data and microseismic sensor coordinates to obtain the three-dimensional coordinates of the seismic source and the initial time. The attenuation correction module extracts the original energy and path characteristics of the microseismic waveform data, and corrects the original energy based on the path characteristics, bedding attenuation law and spherical diffusion effect to obtain the corrected source energy. The energy filling module calculates a correction factor based on the corrected source energy and the original energy, calculates a confidence coefficient based on the correction factor, and filters the corrected source energy according to the confidence coefficient to obtain representative energy; The serialization module constructs a feature tensor for the microseismic event sequence based on the initial time, representative energy, and confidence coefficient. The microseismic prediction module inputs the feature tensor of the microseismic event sequence into the pre-trained microseismic energy prediction model to obtain the maximum energy prediction value.

Images