Three-dimensional in-situ stress microseismic inversion method and system

Abstract

The present application belongs to the field of ground stress monitoring, and particularly relates to a three-dimensional ground stress microseismic inversion method and system. The present application only inverts the longitudinal wave velocity at the profile and plane of the specified position in the three-dimensional space to obtain the microseismic event position; and the inversion distributes the nodes according to the density of the rays. If the ray density of the profile or plane is lower than the set density threshold, the nodes are not distributed on the profile or plane. Therefore, the inversion calculation in the area with smaller ray density and unnecessary area is not needed, the calculation amount is saved as much as possible under the premise of not causing great influence on the inversion result, the calculation efficiency is improved, the inversion result is found in time through the inversion result, the regional stress anomaly is found in time, and the inversion method of the present application performs the grid inversion in the set number of rotation directions, and the inversion results are combined into the same model. Therefore, the error caused by the incomplete node distribution is reduced, the inversion calculation efficiency is improved, and the accuracy of the inversion is ensured.

Description

Technical Field

[0001] This invention belongs to the field of geostress monitoring, specifically relating to a three-dimensional geostress microseismic inversion method and system. Background Technology

[0002] With the rapid development of the domestic economy, the exploitation and mining of shallow resources are reaching saturation, leading to increasingly deeper excavations. Problems encountered during deep mining or tunnel excavation, such as severe roadway deformation, gas outbursts, rock bursts, and a surge in rock bursts and rock pressures, are closely related to geostress. Geostress is an unstable field formed by the combined effects of gravity and tectonic stress, varying with time and space. The geostress field changes in real time due to the superposition of pre-existing stress at the mining face and tectonic stress, causing localized stress concentrations. This complex geostress field can severely impact engineering construction, especially in areas with high geostress concentrations in the coal seam roof, which are frequent sites of rock bursts.

[0003] There are many methods for measuring geostress, which can be broadly divided into absolute geostress measurement and relative geostress measurement. Absolute geostress measurement mainly measures the magnitude and direction of geostress, such as hydraulic fracturing, acoustic emission, borehole collapse, core stress relief, and strain recovery methods. Relative geostress measurement, on the other hand, monitors the dynamic changes of geostress over time, such as rigid hollow inclusion stress gauges, borehole component stress gauges, fiber optic grating-based methods, and microseismic monitoring methods.

[0004] Among these methods, microseismic monitoring offers the advantages of real-time, dynamic, and continuous monitoring of microseismic activity in coal mines, recording data rich in information about the changes in the stress field of coal and rock. Currently, methods for analyzing geostress using tomographic imaging can be divided into active tomography and passive tomography. Active tomography using active sources has seen widespread application; for example, in 2012, Wang Wenshu, Mao Debing, et al., in their paper "Evaluation Model of Rockburst Hazard Based on Seismic CT Technology," initially established an evaluation model for rockburst hazard based on seismic CT technology, with seismic wave velocity anomaly coefficient and wave velocity gradient coefficient as the main factors. Also in 2012, Salnikov et al., in "Transmitted-wave seismictomography for Kuznetsk coal basin: Technology and results," proposed a method for synthesizing projected wave seismic records and applied it to transmitted wave tomography in coal mines, achieving good results. However, passive tomography based on microseismic activity has not yet received sufficient attention and development. In recent years, with the rapid development of microseismic monitoring technology and the demand for real-time monitoring, passive tomography has gradually gained attention. In 2016, Dou Linming, Cai Wu, and others, in their paper "Verification of Shock Hazard Assessment by Seismic Wave Velocity Tomography," combined real-time microseismic monitoring with tomography to invert the P-wave velocity of the working face. Comparing this with traditional monitoring methods, they found that strong seismic events often occur in high-velocity regions. In 2018, Guo Laigong, Dai Guanglong, and others used microseismic travel-time imaging technology to monitor stress anomalies in coal mining faces. Compared to active tomography, microseismic-based passive tomography offers a wider monitoring range, higher safety, lower cost, and larger data volume, meeting the needs of long-term continuous monitoring. This represents a future trend in monitoring geological hazards such as shock hazards both domestically and internationally.

[0005] Furthermore, theoretical research and numerous experiments have demonstrated that the propagation velocity of seismic waves is related to the rock's elastic constants, lithology, density, porosity, and geostress. Simultaneously, influenced by topography and geological structures, the stress dispersion gradually increases with depth. In recent years, many domestic scholars have conducted numerous experiments on this issue. In 2012, Gong Siyuan, Dou Linming, and others, in their work "Experimental Study on the Relationship between P-wave Velocity and Stress in Impact-Distant Coal and Rock," experimentally studied the coupling relationship between stress and P-wave velocity in coal and rock samples under uniaxial conditions. The results showed a power function relationship between stress and P-wave velocity, and that the sensitivity of rock samples to stress via P-wave velocity was greater than that of coal samples. In 2019, Xue Yarong, Song Dazhao, and others, in their work "Correlation between P-wave Velocity and Stress in Outburst-Prone Coal and Rock," conducted uniaxial loading failure experiments on outburst-prone coal samples, demonstrating a positive correlation between P-wave velocity and stress in outburst-prone coal samples. In 2020, Deng Zhigang, Wang Hongwei, and others, in their study "Research on the Correlation between P-wave Velocity and Stress in Coal and Rock Masses with Different Impact Tendencies," selected coal and rock masses with different impact tendencies to test the changes in P-wave velocity at various stress stages under uniaxial and cyclic loading conditions. They found that regardless of whether uniaxial or cyclic loading was performed, wave velocity and stress increased and decreased simultaneously, conforming to a power function relationship. Therefore, the stress state of the coal and rock mass can be indirectly determined based on the magnitude of seismic wave velocity, allowing for the analysis of stress anomalies in the study area.

[0006] However, passive tomography methods for inverting P-wave velocities from working surfaces based on microseismic data typically require processing large amounts of data to achieve accurate P-wave velocity inversion. This often necessitates long processing times and high computing power. In contrast, the detection of stress anomalies in areas prone to rockbursts, such as coal mines, usually requires rapid seismic wave velocity detection and analysis. Inverting P-wave velocities from microseismic data may result in low computational speed and inversion efficiency, affecting the detection and analysis of regional stress anomalies. Summary of the Invention

[0007] The purpose of this invention is to provide a three-dimensional geostress microseismic inversion method and system to solve the problem of low computation speed and inversion efficiency caused by the long processing time of existing methods for inverting P-wave velocity based on microseismic data.

[0008] To achieve the above objectives, the present invention provides a three-dimensional geostress microseismic inversion method, comprising the following steps:

[0009] 1) From the collected microseismic data of the target area, select the data segment corresponding to the maximum microseismic event at a set time interval. The collected microseismic data of the target area includes the coordinates of the acquisition station and the first arrival time of the P-wave.

[0010] A one-dimensional velocity model was established based on the coordinates of the acquisition stations and the first arrival time of the P-waves, and the propagation time of microseismic events at different depths to each acquisition station was calculated.

[0011] 2) Based on the time it takes for microseismic events at different depths to propagate to each acquisition station, establish an objective function for the location probability of microseismic events, and search for the extreme value of the objective function to initially locate microseismic events in the one-dimensional velocity model established in step 1).

[0012] 3) Set a smooth stratigraphic interface to divide the target area into three-dimensional geological models and establish a three-dimensional velocity grid model; set several initial rays at the grid starting point of the three-dimensional velocity grid model, and use the ray solving algorithm to obtain the refraction and reflection paths of the rays corresponding to any two points in the three-dimensional velocity grid model after passing through the stratigraphic interface.

[0013] 4) The longitudinal wave velocity of the three-dimensional velocity grid model is inverted at a specific location in three-dimensional space on a cross section or plane to obtain the location of the microseismic event; the inversion is based on the density of rays to distribute grid nodes. If the ray density of a cross section or plane at a specific location in three-dimensional space is lower than a set density threshold, then no nodes are distributed on that cross section or plane.

[0014] If the microseismic event location obtained by the current inversion does not meet the iteration stopping condition, the three-dimensional velocity grid model is corrected accordingly, and the next inversion is performed according to the corrected three-dimensional velocity grid model until the microseismic event location obtained by the inversion meets the iteration stopping condition. Then, the three-dimensional velocity grid model corresponding to this iteration is taken as the final inversion result.

[0015] The advantages of the above technical solution are: it eliminates the need for inversion calculations in areas with low ray density and unnecessary areas, and can save as much computational load as possible without significantly affecting the inversion results, thereby improving computational efficiency and enabling timely detection of regional stress anomalies through the inversion results.

[0016] Furthermore, the location of the microseismic event is obtained by inverting the P-wave velocity of the three-dimensional velocity grid model at a specified three-dimensional location on a cross-section and plane.

[0017] The three-dimensional mesh nodes are rotated horizontally with the center of the three-dimensional velocity mesh model as the origin. A set number of rotation directions are selected. For each rotation direction, the longitudinal wave velocity of the three-dimensional velocity mesh model is inverted to obtain the inversion results for the corresponding rotation direction. The data corresponding to the inversion results of each rotation direction are added together and averaged to merge the inversion results into the same model.

[0018] The beneficial effects of the above technical solution are: it can reduce the error caused by the incomplete distribution of nodes, that is, reduce the impact of the above inversion method of distributing nodes according to the density of rays on the inversion effect, thereby improving the inversion calculation efficiency while ensuring the accuracy of the inversion.

[0019] Furthermore, the ray initiation parameters of the initial ray described in step 3) are estimated based on a one-dimensional velocity model.

[0020] Furthermore, in step 2), the method for preliminary localization of microseismic events is as follows:

[0021]

[0022] Where G is the objective function for the location probability of microseismic events, and N is the total number of recorded microseismic events. It is the travel time of the grid nodes. and These are the lower and upper limits of the set residual value range, respectively; Represents distance weight;

[0023] The grid search method is used to search for the extrema of the function G, and each extremum corresponds to a microseismic event.

[0024] Furthermore, the longer the ray's path, the smaller the corresponding distance weight value.

[0025] The beneficial effect of the above technical solution is that it can balance the influence of large residuals caused by long-path rays.

[0026] Furthermore, the method for obtaining the refraction and reflection paths of rays corresponding to any two points in the three-dimensional velocity mesh model after passing through the formation interface using the ray-finding algorithm is as follows:

[0027] The formula for calculating the propagation path of a ray in a 3D velocity mesh model corresponding to the ray tracing algorithm is as follows:

[0028]

[0029] Where P x P y P z , , z are the slowness in the x, y, z directions of the three-dimensional velocity mesh model, representing the ray direction; V(x, y, z) is the ray propagation speed; dS is the integration step size;

[0030] The refraction and reflection paths of rays after passing through the geological interface are as follows:

[0031]

[0032] In the formula Let P be the initial slowness of the i-th ray. i N represents the slowness of the i-th ray after refraction or reflection; i It is the normal vector at the intersection of the stratigraphic interface and the i-th ray; and The velocity of the ray on both sides of the geological interface is denoted as .

[0033] Furthermore, if the ray does not encounter a formation interface, the integration step size is adjusted based on the velocity change characteristics around the location of the microseismic event; if the ray encounters a formation interface, the integration step size is increased when the velocity change on both sides of the ray passes through a formation interface that exceeds a set change threshold.

[0034] The beneficial effects of the above technical solution are: when the ray passes through a stratigraphic interface with a significant velocity change, it can quickly approach the next interface by increasing the integration step size, thereby accelerating the solution speed of the refraction and reflection path of the ray after passing through the stratigraphic interface, and thus improving the efficiency of the overall inversion process.

[0035] Furthermore, the iteration stopping condition is: the residual between the microseismic event location obtained by the current inversion and the microseismic event location obtained by the previous inversion is less than a set residual threshold. If the residual is less than the set residual threshold, the iteration stopping condition is determined to be met.

[0036] The present invention also provides a three-dimensional geostress microseismic inversion system, including a processor for executing program instructions to implement the three-dimensional geostress microseismic inversion method as described above.

[0037] This three-dimensional geostress microseismic inversion system can achieve the same beneficial effects as the three-dimensional geostress microseismic inversion method described above. Attached Figure Description

[0038] Figure 1 This is a flowchart of the three-dimensional geostress microseismic inversion method in an embodiment of the present invention;

[0039] Figure 2 This is a schematic diagram of a three-dimensional geological model established in an embodiment of the three-dimensional geostress microseismic inversion method of the present invention;

[0040] Figure 3 This is a schematic diagram of the physical parameters of the strata in the three-dimensional geological model set in the embodiment of the three-dimensional geostress microseismic inversion method of the present invention.

[0041] Figure 4a This is a schematic diagram of the distribution structure of microseismic events uniformly distributed in the space of the coal and rock mass in an embodiment of the three-dimensional geostress microseismic inversion method of the present invention;

[0042] Figure 4b This is a schematic diagram of the inversion results corresponding to the distribution structure of microseismic events uniformly distributed in the space of the coal and rock mass in an embodiment of the three-dimensional geostress microseismic inversion method of the present invention;

[0043] Figure 5aThis is a schematic diagram of the distribution structure of microseismic events concentrated in the space above the coal seam in an embodiment of the three-dimensional geostress microseismic inversion method of the present invention;

[0044] Figure 5b This is a schematic diagram of the inversion results corresponding to the distribution structure of microseismic events concentrated in the space above the coal seam in an embodiment of the three-dimensional geostress microseismic inversion method of the present invention;

[0045] Figure 6a This is a schematic diagram of the distribution structure of microseismic events concentrated in the space below the coal seam in an embodiment of the three-dimensional geostress microseismic inversion method of the present invention;

[0046] Figure 6b This is a schematic diagram of the inversion results corresponding to the distribution structure of microseismic events concentrated in the space below the coal seam in an embodiment of the three-dimensional geostress microseismic inversion method of the present invention. Detailed Implementation

[0047] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0048] Example of a three-dimensional geostress microseismic inversion method

[0049] This embodiment presents a technical solution for a three-dimensional geostress microseismic inversion method, referring to... Figure 1 It includes the following steps:

[0050] 1) From the collected microseismic data of the target area, select the data segment corresponding to the maximum microseismic event at a set time interval (the signals collected by the microseismic instrument are all time-domain signals). The collected microseismic data includes the coordinates of the acquisition station and the first arrival time of the P-wave. The coordinates of the acquisition station correspond to the location where the acquisition instrument (used to collect microseismic data) is placed, and can also be referred to as the station coordinates. In this embodiment, the maximum microseismic events are selected as evenly distributed as possible.

[0051] A one-dimensional velocity model is established based on the coordinates of the acquisition stations and the first arrival time of the P-waves, and the propagation time of microseismic events at different depths to each acquisition station is calculated. The propagation time of microseismic events at different depths to each acquisition station is calculated based on the analytical formulas of the seismic wave propagation ray theory and Fermat's principle, and a travel time table is established.

[0052] 2) Based on the propagation time of microseismic events at different depths to each acquisition station, an objective function for the location probability of microseismic events is established. The extreme values ​​of this objective function are searched to preliminarily locate microseismic events within the one-dimensional velocity model established in step 1), i.e., to estimate the location of microseismic events in the underground space. In this embodiment, the established one-dimensional velocity model is: V=K X represents the model where velocity varies with depth, where V is velocity, K is a coefficient, and X is depth. This is a commonly used, simple one-dimensional model. Since the signals acquired by microseismic instruments are all time-domain signals, while the inversion results belong to the spatial domain, velocity parameters are involved.

[0053] The method for preliminary localization of microseismic events is as follows:

[0054]

[0055] Where G is the objective function for the location probability of microseismic events, and N is the total number of recorded microseismic events. It is the travel time of the grid node (i.e., the time it takes for a ray to travel through the grid node). and These are the lower and upper limits of the set residual value range, respectively. Since the objective function G adopts the conventional iterative calculation method, given an initial value, a calculation result is obtained. The initial value is adjusted and the calculation is continued, and another result is obtained. The residual obtained is the data difference of the objective function obtained from the two iterations. This represents the distance weight; since rays with longer travel distances usually have larger residuals, their weights are smaller than those of rays with shorter travel distances. In other words, the longer the ray travels, the smaller the corresponding distance weight. After dividing the underground three-dimensional space into a grid, underground microseismic events need to pass through many nodes in the grid to propagate to the surface acquisition station. The total time taken is the sum of the travel times of the rays passing through these nodes, i.e., the time it takes for microseismic events at different depths to propagate to each acquisition station.

[0056] According to Fermat's principle, which states that the path of a seismic wave between two points should satisfy the principle of minimum time, microseismic time localization is to determine the grid node position corresponding to the time when the seismic wave takes the least time to propagate. The coordinates of the grid node position are the corresponding positions of the microseismic events. At this time, the objective function should have a minimum value. Therefore, the grid search method is used to search for the extreme values ​​of the function G, and each extreme value corresponds to a microseismic event.

[0057] 3) Set up a three-dimensional geological model corresponding to the target area with a smooth interface segmentation, and use a third-order polynomial to perform linear interpolation of the interlayer grid nodes to establish a three-dimensional velocity grid model; in this embodiment, the grid node spacing in the horizontal and vertical directions is set to be variable to improve the calculation accuracy.

[0058] Several initial rays are set at the starting point of the three-dimensional velocity grid model. The refraction and reflection paths of rays corresponding to any two points in the three-dimensional velocity grid model after passing through the formation interface are obtained using a ray tracing algorithm. Specifically, in this embodiment, a fast calculation method for ray tracing between any two points in three-dimensional space is constructed to obtain the refraction and reflection paths of rays corresponding to any two points in the three-dimensional velocity grid model after passing through the formation interface. That is, a series of initial rays are set at the starting point of the grid, and the following single-point ray tracing algorithm is used:

[0059]

[0060] Where P x P y P z Here, denoted by , V(x, y, z) represents the slowness in the x, y, and z directions of the three-dimensional velocity grid model, signifying the ray direction; V(x, y, z) represents the ray propagation velocity; and dS represents the integration step size. In this embodiment, if the ray does not encounter a geological interface, the integration step size is adjusted based on the velocity variation characteristics around the microseismic event location. If the ray encounters a geological interface, the integration step size is increased when the velocity variation on both sides of the ray (i.e., the velocity difference between the ray and the geological interface) exceeds a set threshold. This formula is used to calculate the seismic wave propagation path, which corresponds to the ray propagation path in the three-dimensional velocity grid model.

[0061] As seismic waves propagate, they encounter geological interfaces, where reflected and refracted waves are generated. Therefore, the refraction and reflection paths of the ray after passing through the geological interface are as follows:

[0062]

[0063] In the formula Let P be the initial slowness of the i-th ray. i N represents the slowness of the i-th ray after refraction or reflection; i It is the normal vector at the intersection of the stratigraphic interface and the i-th ray; and The velocity of the ray on both sides of the geological interface is denoted as .

[0064] The ray initiation parameters of the initial ray are estimated based on a one-dimensional velocity model. In this embodiment, the ray initiation parameters mainly refer to the initial propagation direction of the ray, including the ray tilt angle and azimuth angle.

[0065] 4) To improve computational efficiency, the P-wave velocity of the 3D velocity grid model is inverted only at specific locations in 3D space, on cross-sections and planes, to obtain the location of microseismic events. The inversion distributes grid nodes according to the ray density (here, the grid node is each data point in the corresponding grid of the 3D velocity grid model, and the inversion calculation is to calculate the value corresponding to the location of each grid node). If the ray density of a cross-section or plane at a specific location in 3D space is lower than the set density threshold, nodes are not distributed on that cross-section or plane. Thus, it is not necessary to perform inversion calculations in areas with low ray density and unnecessary areas. This can save as much computational load as possible without significantly affecting the inversion results, thereby improving computational efficiency.

[0066] If the microseismic event location obtained by the current inversion does not meet the iteration stopping condition, the three-dimensional velocity grid model is corrected accordingly, and the next inversion is performed according to the corrected three-dimensional velocity grid model until the microseismic event location obtained by the inversion meets the iteration stopping condition. At this point, the three-dimensional velocity grid model corresponding to this iteration is taken as the final inversion result. In this embodiment, the iteration stopping condition is: the residual between the microseismic event location obtained by the current inversion and the microseismic event location obtained by the previous inversion is less than a set residual threshold. If the residual is less than the set residual threshold, it is determined that the iteration stopping condition is met.

[0067] In this embodiment, the velocity inversion iteration for retrieving the P-wave velocity from the three-dimensional velocity grid model is calculated using the LSQR algorithm. To improve the stability of the inversion calculation, two additional matrices are set for control. One is a diagonal matrix with only one element per row and zero data; increasing its value reduces the amplitude of velocity anomalies. The other is a smooth structure matrix with two equal non-zero elements of opposite signs per row for regularization; increasing its value makes the velocity model smoother.

[0068] Furthermore, the location of the microseismic event is obtained by inverting the P-wave velocity of the three-dimensional velocity grid model at a specified three-dimensional location on a cross-section and plane.

[0069] The 3D mesh nodes are rotated horizontally (counter-clockwise in this embodiment) with the center of the 3D velocity mesh model as the origin. A set number of rotation directions (i.e., rotation angles) are selected, and a rotation is performed once for each rotation direction. The P-wave velocity of the 3D velocity mesh model is inverted for each rotation, obtaining the inversion results for each rotation direction. In a preferred embodiment, the number is set to 5, that is, the P-wave velocity of the 3D velocity mesh model is inverted in five selected rotation directions (0°, 20°, 40°, 60°, 80°), obtaining the inversion results for each rotation direction. The data corresponding to the inversion results for each rotation direction are added together and averaged to merge the inversion results into a single model. This reduces the error caused by incomplete node distribution, that is, it reduces the impact of the above-mentioned inversion method of distributing nodes according to ray density on the inversion effect, thereby improving the inversion calculation efficiency while ensuring the accuracy of the inversion.

[0070] In this embodiment, a three-dimensional geological model with an X-axis of 3000m, a Y-axis of 3000m, and a Z-axis of 1000m is established as the three-dimensional geological model corresponding to the target area. The coal seam is 10m thick and located between 600-610m. A cylindrical high-speed anomaly is inserted at 400-600m in the geological model to represent a high-stress concentration area. Figure 2 And set the physical property parameters of the strata in the three-dimensional geological model, referring to Figure 3 Taking this geological model as an example, the inversion effect of the three-dimensional geostress microseismic inversion method in this embodiment is demonstrated:

[0071] Due to the shielding effect of coal seams on elastic waves, the propagation patterns of seismic waves underground are more complex. Microseismic events caused by factors such as rock fracturing at different locations within the underground rock mass exhibit different propagation patterns of seismic waves within the coal and rock strata. To accurately detect the location of high-stress zones in the coal seam roof, based on a three-dimensional geological model, the spatial locations of microseismic events are divided into areas uniformly distributed within the coal and rock mass (e.g., Figure 4a and Figure 4b ), concentrated above the coal seam (such as Figure 5a and Figure 5b ) and concentrated below the coal seam ( Figure 6a and Figure 6b The tomographic inversion effect is illustrated using three case diagrams. In this embodiment, the inversion result is divided into three parts: the horizontal section, the XZ vertical section, and the YZ vertical section. Figure 4a This represents the distribution structure of microseismic events uniformly distributed within the coal and rock mass. Figure 4b This is the inversion result corresponding to the distribution structure; Figure 5a This represents the distribution structure of microseismic events concentrated in the space above the coal seam. Figure 5bThis is the inversion result corresponding to the distribution structure; Figure 6a This represents the distribution structure of microseismic events concentrated in the space beneath the coal seam. Figure 6b This is the inversion result corresponding to this distribution structure.

[0072] The efficiency of this three-dimensional geostress microseismic inversion method in inversion calculations is significantly better than other existing inversion methods, as detailed below:

[0073] The inversion calculations were performed on a Windows 10 platform using an i7-12700K CPU, as shown in Table 1.

[0074] Table 1

[0075]

[0076] As can be seen from Table 1, the three-dimensional geostress microseismic inversion method (i.e., the fast three-dimensional inversion method) of this embodiment significantly reduces the inversion time (i.e., the corresponding three-dimensional grid calculation time) compared with the conventional three-dimensional inversion method. In particular, the improvement in calculation efficiency is more obvious for larger three-dimensional velocity grid models. Therefore, when facing a large target area, the three-dimensional geostress microseismic inversion method of this embodiment can effectively improve the inversion calculation efficiency.

[0077] Example of a three-dimensional geostress microseismic inversion system

[0078] This embodiment provides a three-dimensional geostress microseismic inversion system, which includes a processor. The processor is used to execute program instructions to implement the three-dimensional geostress microseismic inversion method as described in the above embodiment.

[0079] Since the specific operation mode, principle and example of the three-dimensional geostress microseismic inversion system in this embodiment have been described in detail in the above-mentioned three-dimensional geostress microseismic inversion method embodiment, they will not be repeated here.

[0080] This invention has the following characteristics:

[0081] 1) The P-wave velocity of the three-dimensional velocity grid model is inverted only at specific locations in three-dimensional space, on cross-sections and planes, to obtain the location of microseismic events. The inversion is based on the distribution of nodes according to the ray density. If the ray density of a cross-section or plane at a specific location in three-dimensional space is lower than the set density threshold, nodes are not distributed on that cross-section or plane. Therefore, it is not necessary to perform inversion calculations in areas with low ray density or unnecessary areas. This can save as much computational effort as possible without significantly affecting the inversion results, thereby improving computational efficiency and enabling timely detection of regional stress anomalies through the inversion results.

[0082] 2) When inverting the longitudinal wave velocity of a three-dimensional velocity grid model on a cross-section and plane at a specified location in three-dimensional space, grid inversion is performed separately in a set number of rotation directions, and the inversion results are merged into the same model. This can reduce the error caused by incomplete node distribution, that is, reduce the influence of the above-mentioned inversion method of distributing nodes according to ray density on the inversion effect, thereby improving the inversion calculation efficiency while ensuring the accuracy of the inversion.

[0083] It should be understood that the above-described specific embodiments of the present invention are merely illustrative or explanatory of the principles of the present invention, and do not constitute a limitation thereof.

Claims

1. A three-dimensional microseismic inversion method for geostress, characterized in that, Includes the following steps: 1) From the collected microseismic data of the target area, select the data segment corresponding to the maximum microseismic event at a set time interval. The collected microseismic data of the target area includes the coordinates of the acquisition station and the first arrival time of the P-wave. A one-dimensional velocity model was established based on the coordinates of the acquisition stations and the first arrival time of the P-waves, and the propagation time of microseismic events at different depths to each acquisition station was calculated. 2) Based on the time it takes for microseismic events at different depths to propagate to each acquisition station, establish an objective function for the location probability of microseismic events, and search for the extreme value of the objective function to initially locate microseismic events in the one-dimensional velocity model established in step 1). 3) Set a smooth stratigraphic interface to divide the target area into three-dimensional geological models and establish a three-dimensional velocity grid model; set several initial rays at the grid starting point of the three-dimensional velocity grid model, and use the ray solving algorithm to obtain the refraction and reflection paths of the rays corresponding to any two points in the three-dimensional velocity grid model after passing through the stratigraphic interface. 4) The longitudinal wave velocity of the three-dimensional velocity grid model is inverted at a specific location in three-dimensional space on a cross section or plane to obtain the location of the microseismic event; the inversion is based on the density of rays to distribute grid nodes. If the ray density of a cross section or plane at a specific location in three-dimensional space is lower than a set density threshold, then no nodes are distributed on that cross section or plane. If the microseismic event location obtained by the current inversion does not meet the iteration stopping condition, the three-dimensional velocity grid model is corrected accordingly, and the next inversion is performed according to the corrected three-dimensional velocity grid model until the microseismic event location obtained by the inversion meets the iteration stopping condition. Then, the three-dimensional velocity grid model corresponding to this iteration is taken as the final inversion result.

2. The three-dimensional geostress microseismic inversion method according to claim 1, characterized in that, The location of a microseismic event is obtained by inverting the P-wave velocity of a three-dimensional velocity grid model from a cross-section and plane at a specific location in three-dimensional space. The three-dimensional mesh nodes are rotated horizontally with the center of the three-dimensional velocity mesh model as the origin. A set number of rotation directions are selected. For each rotation direction, the longitudinal wave velocity of the three-dimensional velocity mesh model is inverted to obtain the inversion results for the corresponding rotation direction. The data corresponding to the inversion results of each rotation direction are added together and averaged to merge the inversion results into the same model.

3. The three-dimensional geostress microseismic inversion method according to claim 1 or 2, characterized in that, The ray initiation parameters of the initial ray mentioned in step 3) are estimated based on a one-dimensional velocity model.

4. The three-dimensional geostress microseismic inversion method according to claim 1 or 2, characterized in that, In step 2), the method for preliminary localization of microseismic events is as follows: Where G is the objective function for the location probability of microseismic events, and N is the total number of recorded microseismic events. It is the travel time of the grid nodes. and These are the lower and upper limits of the set residual value range, respectively; Represents distance weight; The grid search method is used to search for the extrema of the function G, and each extremum corresponds to a microseismic event.

5. The three-dimensional geostress microseismic inversion method according to claim 4, characterized in that, The longer the ray's path, the smaller the corresponding distance weight.

6. The three-dimensional geostress microseismic inversion method according to claim 1 or 2, characterized in that, The method for obtaining the refraction and reflection paths of rays corresponding to any two points in a 3D velocity mesh model after passing through the formation interface using the ray tracing algorithm is as follows: The formula for calculating the propagation path of a ray in a 3D velocity mesh model corresponding to the ray tracing algorithm is as follows: Where P x P y P z , , and z represent the slowness in the x, y, and z directions of the three-dimensional velocity mesh model, respectively, and represent the ray direction; V(x,y,z) is the ray propagation speed; dS is the integration step size; The refraction and reflection paths of rays after passing through the geological interface are as follows: In the formula Let P be the initial slowness of the i-th ray. i N represents the slowness of the i-th ray after refraction or reflection; i It is the normal vector at the intersection of the stratigraphic interface and the i-th ray; and The velocity of the ray on both sides of the geological interface is denoted as .

7. The three-dimensional geostress microseismic inversion method according to claim 6, characterized in that, If the ray does not encounter a formation interface, the integration step size is adjusted according to the velocity change characteristics around the location of the microseismic event; if the ray encounters a formation interface, the integration step size is increased when the velocity change on both sides of the ray passes through a formation interface that exceeds a set change threshold.

8. The three-dimensional geostress microseismic inversion method according to claim 1 or 2, characterized in that, The iteration stopping condition is: the residual between the microseismic event location obtained by the current inversion and the microseismic event location obtained by the previous inversion is less than a set residual threshold. If the residual is less than the set residual threshold, the iteration stopping condition is determined to be met.

9. A three-dimensional geostress microseismic inversion system, characterized in that, Includes a processor for executing program instructions to implement the three-dimensional geostress microseismic inversion method as described in any one of claims 1-8.

Images