A water-injection-induced earthquake monitoring and fault activation detection system

Abstract

This invention discloses a water-injection-induced earthquake monitoring and fault activation detection system, belonging to the field of earthquake monitoring and fault activation detection. The method includes: an earthquake observation module for obtaining earthquake monitoring data; a seismic phase data extraction module for identifying microseismic events based on the earthquake monitoring data to obtain seismic phase data; an initial earthquake location module for obtaining an initial earthquake catalog; a fine earthquake location module for performing secondary location to obtain a finely located earthquake catalog; a fault and structure modeling module for integrating the finely located earthquake catalog into a three-dimensional geological structure model for spatial location and visualization of seismic activity; and a fault activation and seismic hazard analysis module for analyzing the temporal variation trend of seismic activity, b-value changes, fault activation, and seismic hazard based on the three-dimensional geological structure model. This invention solves the problem of existing methods lacking monitoring capabilities for fault activation and weakening in shale gas development areas, as well as water-injection-induced earthquakes.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of earthquake monitoring and fault activation detection, and particularly relates to a water injection induced earthquake monitoring and fault activation detection system. BACKGROUND

[0002] The influence of industrial activities on earthquakes, especially the possibility of inducing destructive earthquakes and its mechanism, has become one of the hotspots of international seismology. There are relevant cases that have confirmed that water injection during shale gas exploitation caused induced earthquakes. How to combine earthquake observation data, carry out experiments and physical-based numerical simulation research to understand the physical process of fault activation and weakening and the mechanism of water injection induced earthquakes, and then establish a risk assessment and prevention system for earthquake disasters caused by injection and production, is a common challenge to the future development of international seismology.

[0003] In shale gas and other industrial development areas, due to the need for injection and production for exploitation, during the injection and production process, stress and fluid action often cause nearby seismic activity, especially in the original fault area. Due to the activation of the fault, the seismic activity is enhanced, which may cause a large earthquake, causing the exploitation process to be suspended or stopped, and causing losses. The scale of the activated fault controls the size of the earthquake magnitude, and the size of the earthquake magnitude determines the strength of the earthquake disaster. Induced earthquakes usually occur within 2km of the hydraulic fracturing operation. Based on the existing geological exploration results of the shale gas demonstration area, a 3D geological model is constructed; by designing a seismic station observation, monitoring microseismic activity, carrying out fine structure imaging of buried faults, obtaining the spatiotemporal distribution characteristics of seismic activity, and comprehensively analyzing fault activation. Therefore, it is necessary to carry out induced earthquake monitoring and fault activation detection, and to research a complete set of technical systems to provide technical support for induced earthquake disasters. SUMMARY

[0004] In view of the above shortcomings in the prior art, the water injection induced earthquake monitoring and fault activation detection system provided by the present application solves the problem that the prior art lacks monitoring of fault activation and weakening in shale gas development areas and water injection induced earthquakes.

[0005] In order to achieve the above-mentioned application purposes, the technical scheme adopted by the present application is as follows: a water injection induced earthquake monitoring and fault activation detection system, comprising a seismic observation module, a seismic phase data extraction module, a seismic preliminary positioning module, a seismic precise positioning module, a fault and structure modeling module, and a fault activation and seismic risk analysis module;

[0006] The seismic observation module is used to obtain seismic monitoring data based on the seismic monitoring station arranged in the shale gas development area;

[0007] The seismic phase data extraction module is used to identify microseisms according to the seismic monitoring data to obtain final seismic phase data.

[0008] The earthquake preliminary positioning module is configured to perform earthquake event correlation and preliminary positioning according to the final seismic phase data, and obtain a preliminary positioned earthquake catalog;

[0009] The earthquake precise positioning module is configured to perform secondary positioning according to the preliminary positioned earthquake catalog and corresponding seismic phase data, and obtain a precisely positioned earthquake catalog;

[0010] The fault and structure modeling module is configured to obtain a three-dimensional geological structure model of a shale gas development area, integrate the precisely positioned earthquake catalog into the three-dimensional geological structure model, and perform spatial positioning and visualization of seismic activity;

[0011] The fault activation and earthquake risk analysis module is configured to analyze a time variation trend of seismic activity, a b value variation, fault activation and earthquake risk according to the three-dimensional geological structure model integrated with the earthquake catalog, in combination with fault and seismic data.

[0012] The present application has the beneficial effect that a complete process system is established for shale gas industrial mining area earthquake monitoring, microseismic detection and identification, microseismic positioning, structure modeling and fault activation, a seismic monitoring system is designed, design indicators are given, a seismic detection and seismic positioning scheme is given, the precision and reliability are improved through double methods and grouping positioning, a scheme for discriminating fault activation is given through modeling, and a basis is provided for shale gas development area earthquake disaster prevention.

[0013] Further, the seismic monitoring array covers at least a thousand meters in the well site area and the area near the well site area in the shale gas development area, and extends to cover c kilometers in the well site fault area, and the distance between seismic stations is at least d kilometers.

[0014] The beneficial effect of the above scheme is that the design scheme achieves monitoring of -1 level microseismic, covers the range of earthquake risk areas, and reduces the economic cost of microseismic observation.

[0015] Further, the final seismic phase data is obtained by converting the seismic monitoring data according to a seismic data standard format in units of days to obtain continuous seismic identification input data.

[0016] According to the continuous seismic identification input data, two methods are used for seismic phase identification respectively.

[0017] The first method is to use the PhaseNet seismic phase selection algorithm to perform seismic phase identification according to the continuous seismic identification input data, and obtain first P wave seismic phase identification results and first S wave seismic phase identification results.

[0018] The second method: Based on the continuous earthquake identification input data, the fuzzy K-means clustering algorithm is used to perform seismic phase analysis to obtain the identification results of the second P-wave and the second S-wave.

[0019]

[0020] in, The waveform to be measured is classified into phase type. The probability score; The waveform to be measured; Seismic phase type Cluster center; For multiplication calculation; Seismic phase type Feature weights; Seismic phase type ; Seismic phase type Feature weights; This represents the total number of seismic phase types. The threshold for the accuracy of identifying the corresponding target seismic phase; For column normal form;

[0021] Based on the arrival times of the first motion phases from the identification results of the first P-wave phase, the first S-wave phase, the second P-wave phase, and the second S-wave phase, earthquake phase data within the region where the signal-to-noise ratio of both methods is greater than the signal-to-noise ratio threshold and the travel time difference between the P-wave and S-wave phases is not greater than the travel time difference threshold is retained as the final phase data. The final phase data includes the phase type and the arrival time of the first motion phase.

[0022] The benefits of the previous approach are: improved reliability of seismic phase identification and significantly reduced false identification.

[0023] Furthermore, the obtained initial earthquake catalog specifically includes:

[0024] A regional crustal velocity structure model was initialized using well logging, seismic exploration, and near-seismic tomography.

[0025] Based on the initialization results of the regional crustal velocity structure model, the P-wave first motion phase travel time tables and S-wave first motion phase travel time tables were calculated using the Taup seismic wave travel time calculation software; the grid spacing of the Taup seismic wave travel time calculation software was set to 1 / 10-1 / 5 of the station spacing.

[0026] Based on the final seismic phase data, the P-wave first motion phase travel time table, and the S-wave first motion phase travel time table, the REAL rapid seismic association and location method is used to perform seismic association and preliminary seismic location. This obtains P-wave phase information, S-wave phase information, seismic occurrence event, and seismic location for the same seismic event. Seismic events with a root mean square error greater than the error threshold in the travel time fitting are removed to obtain a preliminary seismic catalog. The preliminary seismic catalog includes several seismic events and the corresponding seismic occurrence time, seismic location, P-wave phase arrival time of each seismic station, and S-wave phase arrival time of each seismic station for each seismic event.

[0027] The beneficial effects of the previous step are: by reasonably setting the velocity model, the reliability of seismic correlation data is improved, providing high-precision input data for the next step.

[0028] Furthermore, the obtained precisely located earthquake catalog specifically includes:

[0029] The initial earthquake catalog was used as input for fine-tuning earthquake localization. For each earthquake event, two methods were used for fine-tuning:

[0030] The first method involves meshing the regional crustal velocity structure model to obtain several model elements, and then using the ray tracing method to calculate the theoretical P-wave and S-wave travel times of seismic waves from each seismic station to each model element node.

[0031] Using the ray tracing method, based on the arrival times of the P-wave phases, the arrival times of the S-wave phases, the theoretical travel times of the P-waves, and the theoretical travel times of the S-waves of each seismic station in the initially located seismic catalog, several source trajectories with arrival time constraints and arrival time difference constraints are constructed.

[0032] The set of intersection points of the earthquake source trajectories is calculated based on the trajectories of each earthquake source;

[0033] Based on the set of intersection points of the source trajectories, determine the distribution area V of the set of intersection points of the source trajectories used for localization;

[0034] Based on the geometric centroid of the set of intersection points of source trajectories within the distribution area V, the observed time residual of each intersection point of source trajectory in the set of intersection points is detected, and the difference is made with the average observed time residual to obtain the observed time difference value of each intersection point of source trajectory. The intersection points of source trajectories with observed time difference values ​​greater than the difference threshold are marked as anomalous intersection points of source trajectories.

[0035] Remove anomalous source trajectory intersections from the set of source trajectory intersections to obtain the set of remaining source trajectory intersections. Calculate the geometric centroid of the set of remaining source trajectory intersections and use it as the source location. Based on the earthquake occurrence locations in the initially located earthquake catalog, obtain the source location error.

[0036] The second method:

[0037] Input the initial positioning earthquake catalog into the regional crustal velocity structure model, and use the hypoDD earthquake event precise positioning method to perform earthquake precise positioning, and based on the earthquake occurrence position in the initial positioning earthquake catalog, obtain the positioning error of the earthquake source position;

[0038] Reserve the earthquakes whose difference of the positioning error of the earthquake source position of the two methods is not greater than the positioning error threshold, to obtain a set of precisely positioned earthquakes;

[0039] Divide the set of precisely positioned earthquakes into several groups based on the depth of shale gas, and then use the hypoDD earthquake event precise positioning method to perform precise positioning to determine the earthquake occurrence position of each precisely positioned earthquake, and update the initial positioning earthquake catalog based on the earthquake occurrence position of each precisely positioned earthquake, to obtain a precisely positioned earthquake catalog.

[0040] The beneficial effect of the above scheme is that the reliability of earthquake positioning is significantly improved, especially the depth error is smaller, and the earthquake positioning result with reliable quality and high stability is selected, which provides reliable support for the next step of fault activation discrimination.

[0041] Further, the spatial positioning and visualization of the earthquake activity includes:

[0042] Data acquisition operation:

[0043] Acquire structural map data, well location and well layer position data, two-dimensional fault line data and topographic data;

[0044] Data preprocessing operation:

[0045] Convert the structural map data into XYZ format which can be spatially positioned and elevation mapped, and filter the elevation values;

[0046] Map the depth and layer position of the well location and well layer position data to ensure that the well information corresponds to the stratigraphic surface and structural profile in the structural map data;

[0047] Numerical three-dimensional coordinate conversion is performed on the two-dimensional fault line data according to the coordinates, depth and dip to form fault XYZ control point data;

[0048] Perform data normalization processing on the precisely positioned earthquake catalog, and form XYZ point data of the earthquake position according to the latitude, longitude and depth;

[0049] Three-dimensional reconstruction operation:

[0050] The fault XYZ control point data is taken as input of a three-dimensional modeling platform, and through combination of fault space expansion, structural surface regularization processing and well data correction technology, three-dimensional reconstruction of the underground structure is realized, specifically: the XYZ control point data of the fault is interpolated through a triangular subdivision interpolation algorithm of adjacent grid points to construct a fine fault surface; the structural profile in the structural map data is smoothed through a regular gridding technology; the depth and horizon of the stratum are corrected through well location and well horizon data to improve the spatial accuracy of the structural profile, so that a complete regional three-dimensional geological structure model is formed.

[0051] Seismic data integration operation:

[0052] The XYZ point data of the seismic position is integrated into the regional three-dimensional geological structure model in the form of point data to realize spatial positioning and visualization of the seismic activity.

[0053] The beneficial effect of the above scheme is that the three-dimensional visualization of the relationship between the earthquake and the structural features is better, and an intuitive basis is provided for judging the fault activation in combination with the comprehensive judgment of the structural data.

[0054] Further, the analysis of the time variation trend of the seismic activity, the b value variation, the fault activation and the seismic risk is specifically:

[0055] Analysis of the time variation trend of the seismic activity: based on the spatial overlap of the fault distribution characteristics in the regional three-dimensional geological structure model and the seismic catalog data, the spatio-temporal position relationship of the earthquake is analyzed according to the spatial distance calculation and projection, the distance between the earthquake position and the fault in different time periods is calculated, the number of earthquake events in the fault neighborhood is determined, and the time variation trend of the seismic activity is obtained;

[0056] Analysis of the b value variation: based on the spatial overlap of the fault distribution characteristics in the regional three-dimensional geological structure model and the seismic catalog data, the b value between the magnitude and the seismic frequency in different time periods is calculated by using the seismological method, and the b value variation is obtained;

[0057] Analysis of the fault activation: the number of earthquakes monitored in the year before the shale gas development is taken as the spatial background seismic activity, and the number of earthquakes in the fault neighborhood in the monitoring period is counted, and when the number of earthquakes in the fault neighborhood in the monitoring period is greater than twice the spatial background seismic activity, the fault is judged to be activated;

[0058] Analysis of the seismic risk: the fault scale is obtained, and the b value between the magnitude and the seismic frequency after the industrial activity is calculated, and when the b value between the magnitude and the seismic frequency after the industrial activity is greater than the b value threshold and the fault scale is greater than the scale threshold, the seismic risk control is performed.

[0059] The beneficial effect of the previous step is that by combining changes in seismic activity with the spatial relationship between earthquakes and tectonic structures, we can more directly and definitively determine fault activation.

[0060] Furthermore, the distance between the earthquake location and the fault is... ;in, This is the distance from the earthquake location to a point at the same depth on the fault plane; This indicates the fault strike direction.

[0061] The beneficial effect of the previous step is that it can provide a more quantitative basis for judging the relationship between earthquakes and faults.

[0062] Furthermore, the calculation of the b-value includes: calculating the cumulative number of seismic events of different magnitudes from a set minimum magnitude to a set maximum magnitude at intervals of 0.1, obtaining the total number of earthquakes, and then calculating the logarithm of the total number of earthquakes.

[0063]

[0064] in, This represents the total number of earthquakes. The minimum magnitude is set. For magnitude index; For the first The range of magnitude [ , The number of earthquakes within ) This is an interval index within the magnitude range; The maximum set magnitude;

[0065] Constantly changing and calculate the corresponding , obtain As the independent variable, with Let b be the fitted line of the dependent variable, and let b be the slope of the fitted line.

[0066] The beneficial effect of the previous step is that it provides a simple and easy-to-use method for assessing seismic hazard. Attached Figure Description

[0067] Figure 1 This is a flowchart of the system operation of the present invention.

[0068] Figure 2 This is a system structure diagram of the present invention.

[0069] Figure 3 This is a schematic diagram of the earthquake monitoring array design in an embodiment of the present invention.

[0070] Figure 4 This is a diagram showing the calculation results of the base noise of the monitoring stations in the southern Sichuan-Zhaotong area in an embodiment of the present invention.

[0071] Figure 5 A seismic correlation diagram for actual observation data in the Chuannan-Zhaotong region in an embodiment of the present application.

[0072] Figure 6 A residual comparison diagram for example microseismic precise positioning in the Chuannan-Zhaotong region in an embodiment of the present application.

[0073] Figure 7 A comparison diagram for depths of precise positioning by hypoDD before and after microseismic grouping in the Chuannan-Zhaotong region in an embodiment of the present application.

[0074] Figure 8 A depth error distribution diagram for precise positioning of microseismic in the Chuannan-Zhaotong region in an embodiment of the present application.

[0075] Figure 9 A structural modeling flowchart in an embodiment of the present application. DETAILED DESCRIPTION

[0076] The specific embodiments of the present application are described below to facilitate the understanding of the present application for those skilled in the art, but it should be clear that the present application is not limited to the scope of the specific embodiments, and for those skilled in the art, it is obvious that various changes are within the spirit and scope of the present application defined and determined by the appended claims, and all the inventions utilizing the concept of the present application are within the scope of protection.

[0077] As shown in Figure 1 and Figure 2 In an embodiment of the present application, a water injection induced seismic monitoring and fault activation detection system includes a seismic observation module, a seismic phase data extraction module, a seismic preliminary positioning module, a seismic precise positioning module, a fault and structure modeling module, and a fault activation and seismic risk analysis module.

[0078] The seismic observation module is configured to obtain seismic monitoring data based on a seismic monitoring array arranged in a shale gas development region.

[0079] The seismic phase data extraction module is configured to perform microseismic identification according to the seismic monitoring data to obtain final seismic phase data.

[0080] The seismic preliminary positioning module is configured to perform seismic event correlation and preliminary positioning according to the final seismic phase data to obtain a preliminary positioned seismic catalog.

[0081] The seismic precise positioning module is configured to perform secondary positioning according to the preliminary positioned seismic catalog and corresponding seismic phase data to obtain a precisely positioned seismic catalog.

[0082] The fault and structure modeling module is configured to obtain a three-dimensional geological structure model of a shale gas development area, integrate the precisely positioned seismic catalog into the three-dimensional geological structure model, and perform spatial positioning and visualization of seismic activities.

[0083] The fault activation and seismic risk analysis module is configured to analyze time variation trends of seismic activities, b value variation, fault activation, and seismic risk based on the three-dimensional geological structure model after integration of the seismic catalog and in combination with fault and seismic data.

[0084] The system includes a technical system of an entire process from observation to final fault activation identification, including shale gas development area microseismic monitoring seismic array layout design, seismic observation scheme, seismic data processing and microseismic detection and identification, seismic positioning, underground structure and fault modeling, and fault identification and fault activation detection, and is particularly directed to microseismic identification, seismic positioning and activation detection technology. The technical system provides technical support for induced seismic monitoring and detection and fault activation in a shale gas development area, and has important significance for earthquake disaster prevention and control and prediction.

[0085] As shown in Figure 3 , the seismic monitoring array covers at least a 5-kilometer range in the well location area and the area near the well location area in the shale gas development area, and extends to cover a 2-kilometer range in the fault area of the well location. The distance between seismic stations is at least 2 kilometers.

[0086] In this embodiment, the well location, well location direction and range of the shale gas development area are obtained, the fault distribution and scale information of the well location area are obtained through geological investigation or detection, the seismic monitoring array is designed according to the shale gas area and the fault distribution, the array covers the well location area and the area within at least 5 km near the well location area, and extends to cover a 2-kilometer range in the fault area. In order to ensure the monitoring ability of monitoring M1 earthquakes, the distance between seismic stations is preferably 2 km. The noise level of the area is investigated (as shown in Figure 4 ), the local noise level is obtained by analyzing the recorded seismic noise data for 3-24 hours using a short-period seismometer, and it is determined which type of seismic station according to the seismic station construction standard (GB / T 1953.1-2004). The seismic station is constructed by using the seismic station on the bedrock or the shallow well method, the noise level is controlled, the seismic station base noise level reaches the III level environmental noise level or above, so that the microseismic can be clearly observed, and the seismic instrument is selected for station erection. The short-period seismometer or the wider band seismometer with a low frequency end of 0.5 Hz-1 Hz and a high frequency end of 20 Hz or above is used for seismic observation.

[0087] The obtained final seismic phase data are specifically obtained by converting the seismic monitoring data according to the seismic data standard format in units of days to obtain continuous seismic identification input data.

[0088] According to the continuous earthquake identification input data, two methods are used for phase identification respectively:

[0089] The first method: according to the continuous earthquake identification input data, the PhaseNet seismic phase selection algorithm is used for phase identification, and the first P wave phase identification result and the first S wave phase identification result are obtained;

[0090] The second method: according to the continuous earthquake identification input data, the fuzzy K-means clustering algorithm is used for phase analysis, and the second P wave phase identification result and the second S wave phase identification result are obtained:

[0091]

[0092] Wherein, is the probability score of the to-be-measured waveform being classified as a phase type; is the to-be-measured waveform; is the cluster center of the phase type; is the multiplication calculation; is the feature weight of the phase type; is the phase type; is the feature weight of the phase type; is the total amount of phase types; is the identification accuracy threshold of the corresponding target phase; is the column norm; is the feature weight of the phase type; is the total amount of phase types; is the identification accuracy threshold of the corresponding target phase; is the column norm;

[0093] According to the first P wave phase identification result, the first S wave phase identification result, the second P wave phase identification result and the second S wave phase identification result, the initial motion phase arrival time, the seismic phase data in the region where the signal-to-noise ratio of the two methods is greater than the signal-to-noise ratio threshold and the P wave phase and S wave phase moveout difference is not greater than the moveout difference threshold is retained as the final phase data; the final phase data includes the phase type and the initial motion phase arrival time.

[0094] In this embodiment, the seismic observation data processing converts the seismic observation data according to the seismic data standard format miniSEED or SAC format according to the day, as the continuous earthquake identification input data; the PhaseNet method AI is used for phase identification, and the P and S wave phases are identified; the algorithm developed based on clustering analysis is used for phase analysis, which includes the steps of phase labeling and intercepting, parameter extraction, parameter screening and waveform clustering, and the P and S phases are confirmed again; according to the P and S phase arrival time, the seismic phase data in the region where the signal-to-noise ratio of the two methods is high is retained, the P and S wave phase moveout difference is not greater than 5 seconds, and the phase data file is output. ​

[0095] The obtained initial positioning earthquake catalog, in particular:

[0096] Initialize the regional crustal velocity structure model by well logging, seismic exploration and near-shock tomography;

[0097] According to the initialization result of the regional crustal velocity structure model, calculate the P-wave initial motion phase travel time table and the S-wave initial motion phase travel time table using the Taup seismic wave travel time calculation software; The grid spacing of the Taup seismic wave travel time calculation software is set to 1 / 10-1 / 5 of the station spacing;

[0098] According to the final seismic phase data, the P-wave initial motion phase travel time table and the S-wave initial motion phase travel time table, the REAL fast seismic correlation and positioning method is used for seismic correlation and preliminary positioning of earthquakes to obtain P-wave phase information, S-wave phase information, earthquake occurrence event and earthquake occurrence position of the same earthquake event, remove the earthquake event with a travel time fitting root mean square error greater than the error threshold, and obtain the initial positioning earthquake catalog; The initial positioning earthquake catalog includes several earthquake events and the corresponding earthquake occurrence time, earthquake occurrence position, P-wave phase arrival time of each seismic station and S-wave phase arrival time of each seismic station.

[0099] As shown in Figure 5 , in this embodiment, the input seismic phase file is in the format of P and S seismic phase travel time, which is the relative time to millisecond precision relative to the reference time (0 point); Determine the regional velocity structure model, perform noise and near-shock tomography, and obtain the velocity structure to determine the regional velocity structure model mainly based on the P-wave velocity of exploration well data and supplemented by noise and near-shock tomography; According to the regional structure velocity model, use Taup to calculate the travel time table, and set the grid to 1 / 10-1 / 5 of the station spacing; Perform seismic correlation, use the REAL method for seismic correlation, and obtain seismic correlation information; Use Hypoinverse for preliminary positioning of earthquakes, obtain an earthquake catalog, and remove earthquake events with large errors.

[0100] The obtained precise positioning earthquake catalog, in particular:

[0101] The initial positioning earthquake catalog is used as the input for precise positioning of earthquakes, and two methods are used for precise positioning of each earthquake event:

[0102] The first method: grid subdivision is performed on the regional crustal velocity structure model to obtain a plurality of model units, and ray tracing method is used to calculate the P-wave theoretical travel time and S-wave theoretical travel time of seismic waves from each seismic station to each model unit node;

[0103] Using the ray tracing method, according to the P wave phase arrival time of each seismic station in the preliminary positioning earthquake catalog, the S wave phase arrival time of each seismic station, the P wave theoretical travel time and the S wave theoretical travel time, a plurality of source trajectory is constructed by time constraint and time difference constraint;

[0104] According to each source trajectory, the source trajectory intersection point set is calculated;

[0105] According to the source trajectory intersection point set, the distribution area V of the source trajectory intersection point set for positioning is determined;

[0106] Based on the geometric center of the source trajectory intersection point set in the distribution area V, the observation arrival time residual of each source trajectory intersection point in the source trajectory intersection point set is detected, and the observation arrival time residual is subtracted from the average value of the observation arrival time residual, to obtain the observation arrival time difference value of each source trajectory intersection point. The source trajectory intersection point with an observation arrival time difference value greater than the difference value threshold is marked as an abnormal source trajectory intersection point;

[0107] The abnormal source trajectory intersection point is removed from the source trajectory intersection point set, to obtain the remaining source trajectory intersection point set, and the geometric center of the remaining source trajectory intersection point set is calculated. The geometric center of the remaining source trajectory intersection point set is taken as the source position, and the source position positioning error is obtained based on the earthquake occurrence position in the preliminary positioning earthquake catalog;

[0108] The second method:

[0109] The preliminary positioning earthquake catalog is input into the regional crustal velocity structure model, and the seismic precise positioning is carried out by using the hypoDD seismic event precise positioning method, and the source position positioning error is obtained based on the earthquake occurrence position in the preliminary positioning earthquake catalog;

[0110] The seismic events with a difference value of the source position positioning error of the two methods not greater than the positioning error threshold are retained, to obtain a set of precisely positioned seismic events;

[0111] The set of precisely positioned seismic events is divided into several groups based on the shale gas depth, and the precise positioning is carried out by using the hypoDD seismic event precise positioning method, to determine the earthquake occurrence position of each precisely positioned seismic event, and the preliminary positioning earthquake catalog is updated based on the earthquake occurrence position of each precisely positioned seismic event, to obtain a precisely positioned earthquake catalog.

[0112] As Figure 6-8As shown, in this embodiment, according to the preliminary positioning directory and the seismic phase data, as the input of fine positioning; using the ray tracing technology-based seismic positioning graphic method for seismic fine positioning, the velocity model of the study area is meshed to obtain a plurality of model units, and the ray tracing method is used to calculate the theoretical travel time of P wave or / and S wave from each seismic station to each model unit node in the study area; using the ray tracing method, a plurality of fine source trajectory constraints are constructed according to all the observed times and the theoretical travel times; the intersection point set of the source trajectory is calculated according to the plurality of source trajectories; the distribution area V of the intersection point set of the source trajectory for positioning in the study area is determined; based on the geometric center of the intersection point set of the source trajectory in the distribution area V, the abnormal source trajectory intersection point is detected from the intersection point set of the source trajectory; the abnormal source trajectory intersection point is removed from the intersection point set of the source trajectory, and the geometric center of the remaining source trajectory intersection point in the intersection point set of the source trajectory in the distribution area V is calculated, and the calculated geometric center of the remaining source trajectory intersection point in the intersection point set of the source trajectory in the distribution area V is taken as the source position, and the seismic fine positioning error is given; using the hypoDD method for seismic fine positioning, constructing the regional structure, inputting Hypoinverse for seismic preliminary positioning results and seismic phase data, and carrying out partition fine positioning on the earthquake, and giving the positioning error; comparing the positioning results and errors of the two methods, separating out the data with higher consistency, and classifying the data with larger positioning error difference; according to the regional velocity structure, the region and depth are grouped, and divided into three categories, the induced earthquake prone horizon within 3 km thickness and the upper and lower three layers, and the partition is carried out, the secondary fine positioning is carried out using hypoDD, and the positioning accuracy is improved.

[0113] The spatial positioning and visualization of the seismic activity include:

[0114] Data acquisition operation:

[0115] Obtain the structure map data, well location and well horizon data, two-dimensional fault line data and topographic data;

[0116] Data preprocessing operation:

[0117] Convert the structure map data into XYZ format which can be spatially positioned and elevation mapped, and filter the elevation values;

[0118] Map the depth and horizon of the well location and well horizon data to ensure that the well information corresponds to the stratigraphic surface and the structure profile in the structure map data;

[0119] Numerical three-dimensional coordinate conversion is performed on the two-dimensional fault line data according to the coordinates, depth and dip to form fault XYZ control point data;

[0120] Data normalization is performed on the fine positioning earthquake catalog, and XYZ point data of earthquake location is formed according to latitude, longitude and depth;

[0121] Three-dimensional reconstruction operation:

[0122] The fault XYZ control point data is taken as the input of the three-dimensional modeling platform, and through the combination of fault space expansion, structural surface regularization processing and well data correction technology, the three-dimensional reconstruction of the underground structure is realized, specifically: the XYZ control point data of the fault is interpolated by the triangular subdivision interpolation algorithm of the adjacent grid points to construct a fine fault surface; the structural profile in the structural map data is smoothed by using the regular gridding technology; the stratum depth and horizon are corrected through the well location and well horizon data to improve the spatial accuracy of the structural profile, so that a complete regional three-dimensional geological structure model is formed;

[0123] Seismic data integration operation:

[0124] The XYZ point data of the earthquake location is integrated into the regional three-dimensional geological structure model in the form of point data to realize the spatial positioning and visualization of the earthquake activity.

[0125] As shown in Figure 9 , in the embodiment, the fault and structure data preparation includes:

[0126] Structural map data: geological reflection layer map and contour map reflecting the structural characteristics of faults, folds and the like in the region;

[0127] Well location and well horizon data: providing well coordinates, well depth, horizon distribution and the like information for structural surface correction and stratum correlation;

[0128] Geological model data: covering fault distribution, fault surface characteristics and the like information to provide a basis for fault feature description;

[0129] Topographic data: regional DEM data for reflecting the surface relief and structural trend;

[0130] Earthquake catalog data: covering fine positioning earthquake information to provide data support for analyzing the relationship between faults and earthquake activity.

[0131] The structural map data is converted into XYZ format which can be spatially positioned and elevation mapped, and the elevation values are reasonably filtered; the well data is depth and horizon mapped to ensure the accuracy of the corresponding relationship between the well information and the structural surface; the two-dimensional fault line data is numerically processed to form standard parameters describing the spatial trend and extension length of the fault; the earthquake catalog data is normalized to form a unified point data set.

[0132] The XYZ control point data of the converted fault is taken as input of a three-dimensional modeling platform, and through combination of fault space expansion, structural surface regularization processing and well data correction, three-dimensional reconstruction of the underground structure is realized. The two-dimensional fault line data is used to construct a fault surface through a near-grid triangular partition interpolation algorithm; the structural surface is constructed according to geological data of each layer, and the model is smoothed through a regular gridding technique; the depth and horizon of the stratum are corrected through well data, the spatial precision of the structural surface is improved, and thus a complete regional three-dimensional geological structure model is formed.

[0133] The XYZ point data of the seismic precise positioning catalog is input, and the specific method comprises: integrating the preprocessed seismic catalog data into the three-dimensional model in the form of point data to realize spatial positioning and visualization of the seismic activity.

[0134] The analysis of the time variation trend of the seismic activity, the b value variation, the fault activation and the seismic risk comprises:

[0135] The time variation trend of the seismic activity is analyzed: based on the spatial overlap of the fault distribution characteristics in the regional three-dimensional geological structure model and the seismic catalog data, the spatio-temporal position relationship of the earthquake is analyzed according to the spatial distance calculation and projection, the distance between the earthquake position and the fault in different time periods is calculated, the number of earthquake events in the fault neighborhood is determined, and the time variation trend of the seismic activity is obtained;

[0136] The b value variation is analyzed: based on the spatial overlap of the fault distribution characteristics in the regional three-dimensional geological structure model and the seismic catalog data, the b value between the magnitude and the seismic frequency in different time periods is calculated by using the seismological method, and the b value variation is obtained;

[0137] The fault activation is analyzed: the number of earthquakes monitored in the year before the shale gas development is taken as the spatial background seismic activity, and the number of earthquakes in the fault neighborhood in the monitoring period is counted, when the number of earthquakes in the fault neighborhood in the monitoring period is greater than twice the spatial background seismic activity, the fault is judged to be activated;

[0138] The seismic risk is analyzed: the fault scale is obtained, and the b value between the magnitude and the seismic frequency after the industrial activity is calculated, when the b value between the magnitude and the seismic frequency after the industrial activity is greater than the b value threshold and the fault scale is greater than the scale threshold, the seismic risk control is performed.

[0139] The distance between the earthquake position and the fault is ; wherein, is the distance from the earthquake position to the point of the same depth of the fault surface; is the strike of the fault.

[0140] The calculation of the b value comprises: calculating the cumulative number of earthquake events of different magnitudes from a set minimum magnitude to a set maximum magnitude at intervals of 0.1, obtaining a total number of earthquakes, and calculating the logarithm of the total number of earthquakes:

[0141]

[0142] wherein, is the total number of earthquakes; is the set minimum magnitude; is the magnitude index; is the number of earthquakes in the th magnitude range , is the interval index in the magnitude range; is the set maximum magnitude; is constantly changing

[0143] , and the corresponding is calculated, to obtain a fitting straight line with as the independent variable and as the dependent variable, and the slope of the fitting straight line is taken as the b value. In this embodiment, based on the spatial overlap of the model fault distribution characteristics and the earthquake catalog data, the spatial and temporal position relationship of the earthquake is analyzed by using a seismological method, the distance between the earthquake position and the fault in a specified time period is calculated, the fault scale and the earthquake positioning accuracy are considered, and the earthquake occurring on the fault within a range of 200 meters from the fault is limited;

[0144]

[0145] The seismic activity of the fault region in the two years before the industrial activity is analyzed to obtain the seismic activity time and spatial background seismic activity characteristics, and the relationship b value between the earthquake magnitude and the earthquake frequency is calculated. The b value calculation: the cumulative number of different magnitudes is calculated from ML -3 magnitude earthquakes to 5 magnitude earthquakes at intervals of 0.1, and then the logarithm of the number of earthquakes is calculated. A straight line with a slope is fitted to the change of the logarithm value of the total number of different magnitudes to the maximum earthquake with the earthquake magnitude, and the slope when the fitting deviation is the smallest is the calculated b value.

[0146] By comparing the background seismic activity, when the seismic activity near the fault is multiplied relative to the background activity, it can be considered as a clear sign of the activation of the fault, and it is judged that the fault is activated.

[0147] ​Seismic hazard analysis is conducted based on the b-value after industrial activity and fault scale. Generally, the b-value for natural earthquakes is less than 1.0. According to the GR law and the statistical results of b-values ​​from southern Sichuan to Zhaotong, if b > 1.2, seismic activity is considered relatively strong. If the b-value fitted line intersects the magnitude axis, the magnitude of the micro-earthquake is the predicted maximum possible earthquake magnitude. If the b-value is small, b < 0.7, the intersection with the horizontal axis may indicate a larger magnitude, potentially leading to a significant increase in seismic hazard. If the fault length is less than a few kilometers, the probability of a large earthquake is considered relatively low, providing support for seismic hazard assessment.

Claims

1. A water-injection-induced earthquake monitoring and fault activation detection system, characterized in that, The system comprises a seismic observation module, a seismic phase data extraction module, a seismic preliminary positioning module, a seismic precise positioning module, a fault and structure modeling module, and a fault activation and seismic risk analysis module. The seismic observation module is configured to obtain seismic monitoring data based on a seismic monitoring array arranged in a shale gas development area. The seismic phase data extraction module is configured to perform microseismic identification based on the seismic monitoring data to obtain final seismic phase data. The final seismic phase data is obtained by converting the seismic monitoring data into continuous seismic identification input data according to a seismic data standard format and in units of days. Two methods are used to perform seismic phase identification based on the continuous seismic identification input data. The first method is to use a PhaseNet seismic phase identification algorithm to perform seismic phase identification based on the continuous seismic identification input data to obtain first P-wave seismic phase identification results and first S-wave seismic phase identification results. wherein, a probability score for classifying the waveform under test as a phase type; a waveform under test; a cluster center for a phase type a multiplication calculation; a feature weight for a phase type a cluster center for a phase type a feature weight for a phase type a total amount of phase types; a recognition accuracy threshold value for a corresponding target phase; a column vector;​​​​​ The second method is to use a fuzzy K-means clustering algorithm to perform seismic phase analysis based on the continuous seismic identification input data to obtain second P-wave seismic phase identification results and second S-wave seismic phase identification results. The final seismic phase data is obtained by retaining seismic phase data in a region where the signal-to-noise ratio of both methods is greater than a signal-to-noise ratio threshold and the travel time difference of P-wave and S-wave is not greater than a travel time difference threshold based on the initial motion phase arrival time of the first P-wave seismic phase identification results, the first S-wave seismic phase identification results, the second P-wave seismic phase identification results, and the second S-wave seismic phase identification results. The seismic preliminary positioning module is configured to perform seismic event correlation and preliminary positioning based on the final seismic phase data to obtain a preliminary positioned seismic catalog. The P-wave initial motion phase travel time table and the S-wave initial motion phase travel time table are calculated using Taup seismic wave travel time calculation software based on the initialization results of the regional crustal velocity structure model. The REAL fast seismic correlation and positioning method is used to perform seismic correlation and seismic preliminary positioning based on the final seismic phase data, the P-wave initial motion phase travel time table, and the S-wave initial motion phase travel time table. The seismic precise positioning module is configured to use the preliminary positioned seismic catalog as an input for seismic precise positioning. The first method is to use a ray tracing method to construct a number of source trajectories based on the P-wave phase arrival time of each seismic station in the preliminary positioned seismic catalog, the S-wave phase arrival time of each seismic station, the P-wave theoretical travel time, and the S-wave theoretical travel time. The second method is to input the preliminary positioned seismic catalog into a regional crustal velocity structure model and use the hypoDD seismic event precise positioning method to perform seismic precise positioning. The source position positioning error is obtained based on the location of the seismic event in the preliminary positioned seismic catalog. The seismic event set is obtained by retaining seismic events whose difference in source position positioning error is not greater than a positioning error threshold. The set of precisely positioned seismic events is divided into several groups based on the depth of shale gas, and then precisely positioned by using a hypoDD seismic event precise positioning method to determine the seismic occurrence positions of the precisely positioned seismic events, and update the initial positioning seismic catalog based on the seismic occurrence positions of the precisely positioned seismic events to obtain a precisely positioned seismic catalog; The fault and structure modeling module is configured to obtain a three-dimensional geological structure model of a shale gas development area, integrate the precisely positioned seismic catalog into the three-dimensional geological structure model, and perform spatial positioning and visualization of seismic activity; The fault activation and seismic risk analysis module is configured to analyze the time variation trend of seismic activity, b value variation, fault activation and seismic risk based on the three-dimensional geological structure model integrated with the seismic catalog and in combination with fault and seismic data.

2. The injection-induced seismic monitoring and fault activation detection system of claim 1, wherein, The initial positioning seismic catalog is obtained by: initializing a regional crustal velocity structure model based on well logging and seismic exploration; calculating P-wave initial motion phase travel time tables and S-wave initial motion phase travel time tables by using Taup seismic wave travel time calculation software according to the initialization results of the regional crustal velocity structure model; setting the grid spacing of the Taup seismic wave travel time calculation software to 1 / 10-1 / 5 of the station spacing; performing seismic correlation and preliminary positioning by using a REAL fast seismic correlation and positioning method according to the final seismic phase data, the P-wave initial motion phase travel time tables and the S-wave initial motion phase travel time tables, obtaining P-wave seismic phase information, S-wave seismic phase information, seismic occurrence events and seismic occurrence positions of the same seismic event, removing seismic events with a root mean square error of travel time fitting greater than an error threshold, and obtaining an initial positioning seismic catalog; the initial positioning seismic catalog includes several seismic events and corresponding seismic occurrence times, seismic occurrence positions, P-wave seismic phase arrival times of each seismic station and S-wave seismic phase arrival times of each seismic station.

3. The injection-induced seismic monitoring and fault activation detection system of claim 1, wherein, The first method of precise positioning is specifically: performing grid division on the regional crustal velocity structure model to obtain several model units, and calculating P-wave theoretical travel time and S-wave theoretical travel time of seismic waves from each seismic station to each model unit node by using a ray tracing method; constructing several seismic source track sets with time constraints and time difference constraints by using the ray tracing method according to the P-wave seismic phase arrival times of each seismic station, the S-wave seismic phase arrival times of each seismic station, the P-wave theoretical travel time and the S-wave theoretical travel time in the initial positioning seismic catalog; calculating a seismic source track intersection point set according to each seismic source track; determining a distribution area V of the seismic source track intersection point set for positioning according to the seismic source track intersection point set; detecting the observation travel time residuals of each seismic source track intersection point in the seismic source track intersection point set based on the geometric center of the seismic source track intersection point set in the distribution area V, and performing subtraction with the average value of the observation travel time residuals to obtain the observation travel time difference values of each seismic source track intersection point, and marking the seismic source track intersection points with observation travel time difference values greater than a difference threshold as abnormal seismic source track intersection points; The abnormal hypocenter intersection points are removed from the hypocenter intersection point set to obtain a remaining hypocenter intersection point set, and a geometric center of the remaining hypocenter intersection point set is calculated, and the geometric center of the remaining hypocenter intersection point set is taken as the hypocenter position, and a hypocenter position positioning error is obtained based on the earthquake occurrence positions in the preliminary positioning earthquake catalog.

4. The injection-induced seismic monitoring and fault activation detection system of claim 1, wherein, The spatial positioning and visualization of the seismic activity comprises: a data acquisition operation: acquiring structural map data, well location and well layer position data, two-dimensional fault line data and terrain data; a data preprocessing operation: converting the structural map data into XYZ format for spatial positioning and elevation mapping, and filtering the elevation values; mapping the depth and layer position of the well location and well layer position data to ensure that the well information corresponds to the stratigraphic surface and the structural profile in the structural map data; numerically converting the two-dimensional fault line data into three-dimensional coordinates according to the coordinates, depth and dip to form fault XYZ control point data; performing data normalization processing on the precisely positioned earthquake catalog to form earthquake position XYZ point data according to the latitude, longitude and depth; a three-dimensional reconstruction operation: taking the fault XYZ control point data as the input of a three-dimensional modeling platform, and realizing the three-dimensional reconstruction of the underground structure by combining fault space expansion, structural surface regularization processing and well data correction technology, specifically: interpolating the fault XYZ control point data by a nearby grid point triangulation interpolation algorithm to construct a fine fault surface; smoothing the structural profile in the structural map data by using a regular gridding technique; correcting the stratigraphic depth and layer position by using the well location and well layer position data to improve the spatial accuracy of the structural profile, thereby forming a complete regional three-dimensional geological structure model; a seismic data integration operation: integrating the earthquake position XYZ point data into the regional three-dimensional geological structure model in the form of point data to realize the spatial positioning and visualization of the seismic activity.

5. The injection-induced seismic monitoring and fault activation detection system of claim 1, wherein, The analysis of the time variation trend of the seismic activity, the b value variation, the fault activation and the seismic risk comprises: analyzing the time variation trend of the seismic activity: based on the spatial overlap of the fault distribution characteristics in the regional three-dimensional geological structure model and the earthquake catalog data, using seismological methods, calculating the spatial and temporal position relationship of the earthquakes according to the spatial distance, calculating the distance between the earthquake positions in different time periods and the faults, determining the number of earthquake events in the fault neighborhood, and obtaining the time variation trend of the seismic activity; analyzing the b value variation: based on the spatial overlap of the fault distribution characteristics in the regional three-dimensional geological structure model and the earthquake catalog data, using seismological methods, calculating the b value between the magnitude and the earthquake frequency in different time periods, and obtaining the b value variation; analyzing the fault activation: counting the number of earthquakes monitored in the year before shale gas development as the spatial background seismic activity, and counting the number of earthquakes in the fault neighborhood in the monitoring period, and when the number of earthquakes in the fault neighborhood in the monitoring period is greater than twice the spatial background seismic activity, the fault is judged to be activated. Analyzing earthquake risk: obtaining fault scale, and calculating the relationship between magnitude and earthquake frequency b value after industrial activities, when the relationship between magnitude and earthquake frequency b value after industrial activities is greater than the b value threshold and the fault scale is greater than the scale threshold, the earthquake risk control is carried out.

6. The injection-induced seismic monitoring and fault activation detection system of claim 5, wherein, The distance of the earthquake location from the fault is ; wherein, is the distance of the earthquake location to a point on the fault plane at the same depth; is the strike of the fault.

7. The waterflood-induced seismic monitoring and fault activation detection system of claim 5, wherein, The calculation of the b value comprises: calculating the cumulative number of earthquake events of different magnitudes from the set minimum magnitude to the set maximum magnitude at an interval of 0.1, obtaining the total number of earthquakes, and calculating the logarithm of the total number of earthquakes: wherein, is the total number of earthquakes; is the set minimum magnitude; is the magnitude index; is the number of earthquakes in the th magnitude range , is the interval index in the magnitude range; is the set maximum magnitude;​ constantly changing and calculates the corresponding , obtains a fitting straight line with as an independent variable and as a dependent variable, and takes the slope of the fitting straight line as the b value.

Images