Water injection induced earthquake monitoring and fault activation detection system

Abstract

The invention discloses a water injection induced earthquake monitoring and fault activation detection system, and belongs to the field of earthquake monitoring and fault activation detection. The earthquake observation module is used for obtaining earthquake monitoring data; the seismic phase data extraction module is used for carrying out micro-seismic identification according to the seismic monitoring data to obtain seismic phase data; the earthquake initial positioning module is used for obtaining an earthquake directory of initial positioning; the earthquake fine positioning module is used for carrying out secondary positioning to obtain an earthquake directory of fine positioning; the fault and structure modeling module is used for integrating the finely positioned seismic catalog into a three-dimensional geologic structure model and carrying out spatial positioning and visualization of seismic activities; and the fault activation and earthquake risk analysis module is used for analyzing the earthquake activity time change trend, b value change, fault activation and earthquake risk according to the three-dimensional geologic structure model. The method solves the problem that an existing method lacks monitoring of fault activation and weakening of the shale gas development area and water injection induced earthquakes.

Description

Technical Field

[0001] The present invention belongs to the field of earthquake monitoring and fault activation detection, and in particular relates to a water injection-induced earthquake monitoring and fault activation detection system. Background Art

[0002] The impact of industrial activities on earthquakes, particularly the potential for inducing destructive earthquakes and their mechanisms, has become a hot topic of intense concern in the international seismological community. Cases have confirmed induced earthquakes caused by water injection during shale gas extraction. Combining seismic observation data with experimental and physics-based numerical simulations to understand the physical processes of fault activation and weakening and the mechanisms of water injection-induced earthquakes, and ultimately to establish a risk assessment and prevention system for water injection-induced earthquake disasters, remains a common challenge for the future development of international seismology.

[0003] In industrial development areas, such as shale gas, production requires injection and production. During the injection and production process, stress and fluid interaction often induce nearby seismic activity, particularly in areas of pre-existing faults. Fault activation increases seismic activity, potentially causing large earthquakes, halting or even stopping production, and causing losses. The size of the activated fault controls the magnitude of the earthquake, which in turn determines the severity of the hazard. Induced earthquakes typically occur within 2 km of hydraulic fracturing operations. Based on existing geological exploration results in shale gas demonstration areas, a 3D geological model was constructed. Seismic array observations were designed to monitor microseismic activity, conduct fine structural imaging of hidden faults, and determine the spatial and temporal distribution of seismicity, allowing for a comprehensive analysis of fault activation. Therefore, it is necessary to develop induced earthquake monitoring and fault activation detection, and to develop a comprehensive technical system to provide technical support for induced earthquake disasters. Summary of the Invention

[0004] In response to the above-mentioned deficiencies in the prior art, the present invention provides a water injection-induced seismic monitoring and fault activation detection system that solves the problem that the existing methods lack monitoring of fault activation and weakening in shale gas development areas, as well as water injection-induced seismicity.

[0005] In order to achieve the above-mentioned purpose, the technical solution adopted by the present invention is as follows: a water injection-induced earthquake monitoring and fault activation detection system, including an earthquake observation module, a seismic phase data extraction module, an earthquake preliminary positioning module, an earthquake precise positioning module, a fault and structure modeling module, and a fault activation and earthquake hazard analysis module; The seismic observation module is used to obtain seismic monitoring data based on the seismic monitoring array deployed in the shale gas development area; The seismic phase data extraction module is used to identify microseisms based on seismic monitoring data to obtain final seismic phase data; The earthquake preliminary location module is used to perform earthquake event correlation and preliminary location based on the final seismic phase data to obtain a preliminary located earthquake catalog; The earthquake precise positioning module is used to perform secondary positioning based on the initially positioned earthquake catalog and the corresponding seismic phase data to obtain a precisely positioned earthquake catalog; The fault and structure modeling module is used to obtain and construct a three-dimensional geological structure model of the shale gas development area, and integrate the precisely located earthquake catalog into the three-dimensional geological structure model to perform spatial positioning and visualization of earthquake activities; The fault activation and earthquake hazard analysis module is used to analyze the temporal trend of earthquake activity, b-value changes, fault activation and earthquake hazard based on the three-dimensional geological structure model integrated from the earthquake catalog, combined with fault and earthquake data.

[0006] The beneficial effects of the present invention are as follows: the present invention establishes a complete process system for earthquake monitoring, microseismic detection and identification, microseismic positioning, structural modeling and fault activation in shale gas industrial mining areas, designs an earthquake monitoring system, provides design indicators, and simultaneously provides earthquake detection and earthquake positioning schemes. The accuracy and reliability are improved through dual methods and group positioning, and a scheme for distinguishing fault activation is provided through modeling, providing a basis for earthquake disaster prevention in shale gas development areas.

[0007] Furthermore, the seismic monitoring array covers the well location area in the shale gas development area and the area within a kilometer near the well location area, and extends to cover the well location fault area for c kilometers. The distance between seismic stations is at least d kilometers.

[0008] The beneficial effects of the previous plan are: the design plan can monitor microseismic events of magnitude -1, cover the earthquake-prone area, and reduce the economic cost of microseismic observation.

[0009] Furthermore, the final seismic phase data is obtained by converting the seismic monitoring data according to a standard seismic data format in units of days to obtain continuous earthquake identification input data; According to the continuous earthquake identification input data, two methods are used to identify the seismic phase: The first method: Based on the continuous earthquake recognition input data, the PhaseNet seismic phase selection algorithm is used to perform seismic phase recognition, and the first P-wave seismic phase recognition results and the first S-wave seismic phase recognition results are obtained; 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 second P-wave seismic phase identification results and the second S-wave seismic phase identification results:

[0010] in, Classify the waveform to be measured into the seismic phase type The probability score of is the waveform to be measured; Phase type The cluster center of To calculate multiplication; Phase type The feature weights of Phase type ; Phase type The feature weights of is the total number of seismic phase types; is the recognition accuracy threshold corresponding to the target seismic phase; is the column normal form; 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 arrival time of the first-motion phase of the second S-wave phase identification result, the earthquake phase data in the area where the signal-to-noise ratio of the two methods is greater than the signal-to-noise ratio threshold and the travel time difference between the P-wave phase and the S-wave phase is not greater than the travel time difference threshold are retained as the final phase data; the final phase data includes the phase type and the arrival time of the first-motion phase.

[0011] The beneficial effects of the previous solution are: improving the reliability of seismic phase identification and significantly reducing misidentification.

[0012] Furthermore, the earthquake catalog obtained by preliminary location is specifically: Initialize regional crustal velocity structure models from well logging, seismic exploration, and near-seismic tomography; Based on the initialization results of the regional crustal velocity structure model, the Taup seismic travel time calculation software was used to calculate the P-wave first-motion phase travel time table and the S-wave first-motion phase travel time table. The grid spacing of the Taup seismic travel time calculation software was set to 1 / 10-1 / 5 of the station spacing. Based on the final seismic phase data, P-wave initial motion phase travel time table and S-wave initial motion phase travel time table, the REAL rapid earthquake correlation and positioning method is used to perform earthquake correlation and preliminary earthquake positioning, and obtain the P-wave phase information, S-wave phase information, earthquake occurrence events and earthquake occurrence locations belonging to the same earthquake event. Earthquake events with travel time fitting root mean square errors greater than the error threshold are removed to obtain a preliminary located earthquake catalog; the preliminary located earthquake catalog includes several earthquake events and the corresponding earthquake occurrence time, earthquake occurrence location, P-wave phase arrival time at each seismic station and S-wave phase arrival time at each seismic station.

[0013] The beneficial effect of the previous solution is: by reasonably setting the velocity model, the reliability of earthquake correlation data is improved, providing high-precision input data for the next step.

[0014] Furthermore, the precisely located earthquake catalog is specifically: The earthquake catalog of the initial location is used as the input for earthquake precise location. Two methods are used to precisely locate each earthquake event: The first method is to mesh the regional crustal velocity structure model to obtain several model units. The theoretical travel time of P waves and S waves from each seismic station to each model unit node is calculated using the ray tracing method. Using the ray tracing method, we construct several earthquake source trajectories with arrival time constraints and arrival time difference constraints based on the P-wave phase arrival time, S-wave phase arrival time, P-wave theoretical travel time, and S-wave theoretical travel time of each seismic station in the initially located earthquake catalog. Calculate the intersection point set of the earthquake source trajectories according to the earthquake source trajectories; Determine a distribution area V of the set of intersection points of the earthquake source trajectories used for positioning according to the set of intersection points of the earthquake source trajectories; Based on the geometric centroid of the set of earthquake source trajectory intersection points in the distribution area V, the observed arrival time residual of each earthquake source trajectory intersection point in the set of earthquake source trajectory intersection points is detected, and the observed arrival time residual is subtracted from the average value of the observed arrival time residual to obtain the observed arrival time difference value of each earthquake source trajectory intersection point. The earthquake source trajectory intersection point with an observed arrival time difference value greater than the difference threshold is marked as an abnormal earthquake source trajectory intersection point; Abnormal source trajectory intersection points are removed from the source trajectory intersection point set to obtain a remaining source trajectory intersection point set, and the geometric center of gravity of the remaining source trajectory intersection point set is calculated. The geometric center of gravity of the remaining source trajectory intersection point set is used as the source position, and the source position positioning error is obtained based on the earthquake occurrence position in the initially located earthquake catalog; Second method: The initially located earthquake catalog is input into the regional crust velocity structure model, and the hypoDD earthquake event precise location method is used to accurately locate the earthquake. Based on the earthquake occurrence location in the initially located earthquake catalog, the source location error is obtained. The earthquake events whose source location errors of the two methods differ by no more than the location error threshold are retained to obtain a set of precisely located earthquake events. The set of precisely located earthquake events is divided into several groups based on the shale gas depth, and then precisely located using the hypoDD earthquake event precise positioning method to determine the earthquake occurrence location of each precisely located earthquake event. The initially located earthquake catalog is updated based on the earthquake occurrence location of each precisely located earthquake event to obtain a precisely located earthquake catalog.

[0015] The beneficial effects of the previous scheme are: significantly improving the reliability of earthquake positioning, especially the smaller depth error, selecting earthquake positioning results with reliable quality and high stability, and providing reliable support for the next step of fault activation judgment.

[0016] Furthermore, the spatial positioning and visualization of seismic activity includes: Data acquisition operation: Obtain structural map data, well location and well layer data, two-dimensional fault line data and topographic data; Data preprocessing operations: Convert the construction map data into XYZ format for spatial positioning and elevation mapping, and filter the elevation values; Perform depth and layer mapping on well locations and well layer data to ensure that the well information corresponds to the stratigraphic layers and structural sections in the structural map data; The two-dimensional fault line data is converted into three-dimensional coordinates according to the coordinates, depth and dip angle to form the fault XYZ control point data; Normalize the data of the precisely located earthquake catalog to form XYZ point data of the earthquake location according to latitude, longitude and depth; 3D reconstruction operation: Using fault XYZ control point data as input for the 3D modeling platform, the 3D reconstruction of underground structures is achieved by combining fault spatial expansion, structural surface regularization, and well data correction technology. Specifically, the fault XYZ control point data is interpolated using a neighboring grid triangulation algorithm to construct a fine fault surface; regular gridding technology is used to smooth the structural sections in the structural map data; and the stratigraphic depth and layer position are corrected using well location and well layer position data to improve the spatial accuracy of the structural section, thereby forming a complete regional 3D geological structure model. Seismic data integration operations: The earthquake location XYZ point data is integrated into the regional three-dimensional geological structure model in the form of point data to achieve spatial positioning and visualization of earthquake activities.

[0017] The beneficial effects of the previous solution are: better three-dimensional visualization of the relationship between earthquakes and structural characteristics, providing an intuitive basis for judging fault activation combined with structural data.

[0018] Furthermore, the analysis of temporal trends in seismic activity, b-value changes, fault activation, and seismic risk is specifically as follows: Analyze the temporal trend of earthquake activity: Based on the spatial overlap of fault distribution characteristics in the regional 3D geological structure model and earthquake catalog data, seismological methods are used to analyze the temporal and spatial positional relationship of earthquakes through spatial distance calculation and projection. The distance between earthquake locations and faults in different time periods is calculated, and the number of earthquake events in the fault neighborhood is determined to obtain the temporal trend of earthquake activity. Analyze b-value changes: Based on the spatial overlap of fault distribution characteristics in the regional 3D geological structure model and earthquake catalog data, seismological methods are used to calculate the b-value relationship between earthquake magnitude and earthquake frequency in different time periods to determine the b-value changes; Analyze fault activation: Count the number of earthquakes detected in the year before shale gas development as the spatial background seismicity, and count the number of earthquakes in the fault neighborhood during the monitoring period. When the number of earthquakes in the fault neighborhood during the monitoring period is greater than twice the spatial background seismicity, it is determined that the fault has been activated. Analyze earthquake hazard: Obtain the fault scale and calculate the b-value of the relationship between the magnitude and earthquake frequency after industrial activities. When the b-value of the relationship between the magnitude and earthquake frequency after industrial activities is greater than the b-value threshold and the fault scale is greater than the scale threshold, earthquake hazard control is carried out.

[0019] The beneficial effect of the previous step is that it can more directly and accurately judge fault activation by combining changes in seismic activity with the spatial relationship between earthquakes and structures.

[0020] Furthermore, the distance between the earthquake location and the fault is ;in, is the distance from the earthquake location to the point at the same depth on the fault plane; The fault direction.

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

[0022] Furthermore, the calculation of the b value includes: 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 to obtain the total number of earthquakes, and calculating the logarithm of the total number of earthquakes:

[0023] in, is the total number of earthquakes; is the minimum magnitude set; is the magnitude index; For the Magnitude range[ , ) number of earthquakes within; is the interval index within the magnitude range; is the maximum magnitude set; Constantly changing , and calculate the corresponding , get As the independent variable, is the fitted straight line of the dependent variable, and the slope of the fitted straight line is taken as the b value.

[0024] The beneficial effect of the previous solution is to provide a simple and easy judgment method for earthquake hazard assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 This is a flow chart of the system operation of the present invention.

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

[0027] Figure 3 Schematic diagram of the earthquake monitoring array design in an embodiment of the present invention.

[0028] 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.

[0029] Figure 5 Schematic diagram of earthquake correlation of actual observation data in the southern Sichuan-Zhaotong area in an embodiment of the present invention.

[0030] Figure 6 This is a schematic diagram of the comparison of residuals of microseismic precise positioning in the southern Sichuan-Zhaotong area in an embodiment of the present invention.

[0031] Figure 7 Schematic diagram of comparison of depth of precise positioning using hypoDD before and after microseismic grouping in the Southern Sichuan-Zhaotong example area in an embodiment of the present invention.

[0032] Figure 8 Schematic diagram of depth error distribution after microseismic precise positioning in the Southern Sichuan-Zhaotong example area in an embodiment of the present invention.

[0033] Figure 9 This is a flow chart of structural modeling in an embodiment of the present invention. DETAILED DESCRIPTION

[0034] The specific embodiments of the present invention are described below to facilitate understanding of the present invention by those skilled in the art. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the appended claims, these changes are obvious, and all inventions and creations utilizing the concepts of the present invention are protected.

[0035] like Figure 1 and Figure 2 As shown, in one embodiment of the present invention, a water injection-induced earthquake monitoring and fault activation detection system includes an earthquake observation module, a seismic phase data extraction module, an earthquake preliminary positioning module, an earthquake precise positioning module, a fault and structure modeling module, and a fault activation and earthquake hazard analysis module; The seismic observation module is used to obtain seismic monitoring data based on the seismic monitoring array deployed in the shale gas development area; The seismic phase data extraction module is used to identify microseisms based on seismic monitoring data to obtain final seismic phase data; The earthquake preliminary location module is used to perform earthquake event correlation and preliminary location based on the final seismic phase data to obtain a preliminary located earthquake catalog; The earthquake precise positioning module is used to perform secondary positioning based on the initially positioned earthquake catalog and the corresponding seismic phase data to obtain a precisely positioned earthquake catalog; The fault and structure modeling module is used to obtain and construct a three-dimensional geological structure model of the shale gas development area, and integrate the precisely located earthquake catalog into the three-dimensional geological structure model to perform spatial positioning and visualization of earthquake activities; The fault activation and earthquake hazard analysis module is used to analyze the temporal trend of earthquake activity, b-value changes, fault activation and earthquake hazard based on the three-dimensional geological structure model integrated from the earthquake catalog, combined with fault and earthquake data.

[0036] This system includes the layout design of microseismic monitoring seismic arrays in shale gas development areas, seismic observation plans, seismic data processing and microseismic detection and identification, earthquake positioning, underground structure and fault modeling, fault identification, fault activation detection, and other technical systems covering the entire process from observation to final fault activation identification. It is particularly aimed at microseismic identification, earthquake positioning and activation detection technologies. The technical system provides technical support for induced earthquake monitoring and detection, and fault activation in shale gas development areas, and is of great significance to the prevention, control and prediction of earthquake disasters.

[0037] like Figure 3 As shown in the figure, the seismic monitoring array covers the well location area in the shale gas development area and the area within a kilometer around the well location area, and extends to c kilometers to the well location fault area. The distance between seismic stations is at least d kilometers.

[0038] In this embodiment, the well location, well strike, and range in the shale gas development area are obtained; the fault distribution and scale information in the well area are obtained through geological survey or exploration; a seismic monitoring array is designed based on the shale gas area and fault distribution, with the array covering the well area and the surrounding area within at least 5 km, and extending coverage to 2 km in the fault area. To ensure the monitoring capability of detecting -1 magnitude earthquakes, the distance between seismic stations is preferably 2 km; a noise level survey is conducted in the area (such as Figure 4(as shown), use short-period seismometers to observe for 3-24 hours, analyze the recorded seismic noise data, determine the local noise level, and determine the type of seismic station according to the Seismic Station Construction Standard (GB / T 1953.1-2004). When constructing a seismic station, control the noise level by installing a seismic station on bedrock or using shallow wells, so that the noise level at the seismic station base reaches Class III environmental noise level or above, so that microearthquakes can be clearly observed. For station installation, select seismographs with a low-frequency end of the operating band within 0.5 Hz to 1 Hz and a high-frequency end of 20 Hz or above, or seismic instruments with a wider band, for earthquake observations.

[0039] The final seismic phase data is obtained by converting the seismic monitoring data into a day-based unit according to a standard seismic data format to obtain continuous earthquake identification input data; According to the continuous earthquake identification input data, two methods are used to identify the seismic phase: The first method: Based on the continuous earthquake recognition input data, the PhaseNet seismic phase selection algorithm is used to perform seismic phase recognition, and the first P-wave seismic phase recognition results and the first S-wave seismic phase recognition results are obtained; 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 second P-wave seismic phase identification results and the second S-wave seismic phase identification results:

[0040] in, Classify the waveform to be measured into the seismic phase type The probability score of is the waveform to be measured; Phase type The cluster center of To calculate multiplication; Phase type The feature weights of Phase type ; Phase type The feature weights of is the total number of seismic phase types; is the recognition accuracy threshold corresponding to the target seismic phase; is the column normal form; 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 arrival time of the first-motion phase of the second S-wave phase identification result, the earthquake phase data in the area where the signal-to-noise ratio of the two methods is greater than the signal-to-noise ratio threshold and the travel time difference between the P-wave phase and the S-wave phase is not greater than the travel time difference threshold are retained as the final phase data; the final phase data includes the phase type and the arrival time of the first-motion phase.

[0041] In this embodiment, earthquake observation data is processed by converting the earthquake observation data into daily format according to the standard earthquake data format miniSEED or SAC format as continuous earthquake identification input data; using the PhaseNet method AI to identify phases and identify P and S wave phases; using an algorithm developed based on cluster analysis to perform phase analysis, the method includes phase labeling and interception, parameter extraction, parameter screening, and waveform clustering steps to secondary confirm the P and S phases; according to the arrival time of the P and S phases, retaining data with a high signal-to-noise ratio of the two methods, retaining earthquake phase data in the area where the travel time difference between the P and S wave phases is not more than 5 seconds, and outputting the phase data file.

[0042] The earthquake catalog obtained by preliminary positioning is specifically: Initialize regional crustal velocity structure models from well logging, seismic exploration, and near-seismic tomography; Based on the initialization results of the regional crustal velocity structure model, the Taup seismic travel time calculation software was used to calculate the P-wave first-motion phase travel time table and the S-wave first-motion phase travel time table. The grid spacing of the Taup seismic travel time calculation software was set to 1 / 10-1 / 5 of the station spacing. Based on the final seismic phase data, P-wave initial motion phase travel time table and S-wave initial motion phase travel time table, the REAL rapid earthquake correlation and positioning method is used to perform earthquake correlation and preliminary earthquake positioning, and obtain the P-wave phase information, S-wave phase information, earthquake occurrence events and earthquake occurrence locations belonging to the same earthquake event. Earthquake events with travel time fitting root mean square errors greater than the error threshold are removed to obtain a preliminary located earthquake catalog; the preliminary located earthquake catalog includes several earthquake events and the corresponding earthquake occurrence time, earthquake occurrence location, P-wave phase arrival time at each seismic station and S-wave phase arrival time at each seismic station.

[0043] like Figure 5 As shown, in this embodiment, a seismic phase file is input in the format of P and S phase travel time, which is the relative time relative to the reference time (0 o'clock) with millisecond accuracy; the regional velocity structure model is determined, and noise and near-seismic tomography are performed to obtain the velocity structure. The regional velocity structure model is determined mainly based on the P-wave velocity of the exploration well data, supplemented by noise and near-seismic tomography; Taup is used to calculate the travel time table based on the regional structure velocity model, and the grid is set to 1 / 10-1 / 5 of the station spacing; earthquake correlation is performed using the REAL method to obtain earthquake correlation information; Hypoinverse is used to perform preliminary earthquake positioning, obtain an earthquake catalog, and remove earthquake events with large errors.

[0044] The precisely located earthquake catalog is specifically: The earthquake catalog of the initial location is used as the input for earthquake precise location. Two methods are used to precisely locate each earthquake event: The first method is to mesh the regional crustal velocity structure model to obtain several model units. The theoretical travel time of P waves and S waves from each seismic station to each model unit node is calculated using the ray tracing method. Using the ray tracing method, we construct several earthquake source trajectories with arrival time constraints and arrival time difference constraints based on the P-wave phase arrival time, S-wave phase arrival time, P-wave theoretical travel time, and S-wave theoretical travel time of each seismic station in the initially located earthquake catalog. Calculate the intersection point set of the earthquake source trajectories according to the earthquake source trajectories; Determine a distribution area V of the set of intersection points of the earthquake source trajectories used for positioning according to the set of intersection points of the earthquake source trajectories; Based on the geometric centroid of the set of earthquake source trajectory intersection points in the distribution area V, the observed arrival time residual of each earthquake source trajectory intersection point in the set of earthquake source trajectory intersection points is detected, and the observed arrival time residual is subtracted from the average value of the observed arrival time residual to obtain the observed arrival time difference value of each earthquake source trajectory intersection point. The earthquake source trajectory intersection point with an observed arrival time difference value greater than the difference threshold is marked as an abnormal earthquake source trajectory intersection point; Abnormal source trajectory intersection points are removed from the source trajectory intersection point set to obtain a remaining source trajectory intersection point set, and the geometric center of gravity of the remaining source trajectory intersection point set is calculated. The geometric center of gravity of the remaining source trajectory intersection point set is used as the source position, and the source position positioning error is obtained based on the earthquake occurrence position in the initially located earthquake catalog; Second method: The initially located earthquake catalog is input into the regional crust velocity structure model, and the hypoDD earthquake event precise location method is used to accurately locate the earthquake. Based on the earthquake occurrence location in the initially located earthquake catalog, the source location error is obtained. The earthquake events whose source location errors of the two methods differ by no more than the location error threshold are retained to obtain a set of precisely located earthquake events. The set of precisely located earthquake events is divided into several groups based on the shale gas depth, and then precisely located using the hypoDD earthquake event precise positioning method to determine the earthquake occurrence location of each precisely located earthquake event. The initially located earthquake catalog is updated based on the earthquake occurrence location of each precisely located earthquake event to obtain a precisely located earthquake catalog.

[0045] like Figure 6-Figure 8As shown, in this embodiment, the preliminary positioning catalog and seismic phase data are used as inputs for precise positioning; a seismic positioning graphic method based on ray tracing technology is used for precise earthquake positioning, the velocity model of the study area is gridded to obtain multiple model units, and the ray tracing method is used to calculate the theoretical travel time of the P wave or / and S wave of the seismic wave from each seismic station to each model unit node in the study area; the ray tracing method is used to construct multiple fine source trajectories with arrival time constraints or / and arrival time difference constraints based on all the observed arrival times and the theoretical travel times; the source trajectory intersection point set is calculated based on the multiple source trajectories; the distribution area V of the source trajectory intersection point set used for positioning is determined in the study area; based on the geometric center of gravity of the source trajectory intersection point set in the distribution area V, the source trajectory intersection point set is detected from the source trajectory intersection point set. Abnormal source trajectory intersection points are obtained; abnormal source trajectory intersection points are removed from the source trajectory intersection point set, and the geometric center of gravity of the remaining source trajectory intersection points in the source trajectory intersection point set within the distribution area V is calculated, and the calculated geometric center of gravity of the remaining source trajectory intersection points in the source trajectory intersection point set within the distribution area V is used as the source position, and the earthquake precise positioning error is given; the hypoDD method is used to precisely locate the earthquake, and the regional structure is constructed. The Hypoinverse is input to perform preliminary earthquake positioning results and phase data, and the earthquake is precisely positioned in different areas to give the positioning error; the positioning results and errors of the two methods are compared, and the data with higher consistency are separated, and the data with larger positioning errors and larger differences are classified; the region and depth are grouped according to the regional velocity structure and divided into three categories, the induced earthquake prone layer within 3 km thickness and the three layers above and below are divided into zones, and the hypoDD is used for secondary precise positioning to improve the positioning accuracy.

[0046] The spatial positioning and visualization of seismic activity includes: Data acquisition operation: Obtain structural map data, well location and well layer data, two-dimensional fault line data and topographic data; Data preprocessing operations: Convert the construction map data into XYZ format for spatial positioning and elevation mapping, and filter the elevation values; Perform depth and layer mapping on well locations and well layer data to ensure that the well information corresponds to the stratigraphic layers and structural sections in the structural map data; The two-dimensional fault line data is converted into three-dimensional coordinates according to the coordinates, depth and dip angle to form the fault XYZ control point data; Normalize the data of the precisely located earthquake catalog to form XYZ point data of the earthquake location according to latitude, longitude and depth; 3D reconstruction operation: Using fault XYZ control point data as input for the 3D modeling platform, the 3D reconstruction of underground structures is achieved by combining fault spatial expansion, structural surface regularization, and well data correction technology. Specifically, the fault XYZ control point data is interpolated using a neighboring grid triangulation algorithm to construct a fine fault surface; regular gridding technology is used to smooth the structural sections in the structural map data; and the stratigraphic depth and layer position are corrected using well location and well layer position data to improve the spatial accuracy of the structural section, thereby forming a complete regional 3D geological structure model. Seismic data integration operations: The earthquake location XYZ point data is integrated into the regional three-dimensional geological structure model in the form of point data to achieve spatial positioning and visualization of earthquake activities.

[0047] like Figure 9 As shown, in this embodiment, the fault and structure data preparation includes: Structural map data: geological reflection layer maps and contour maps reflecting structural features such as faults and folds in the region; Well location and well layer data: Provides information such as well coordinates, well depth, layer distribution, etc. for structural surface correction and stratigraphic comparison; Geological model data: covers information such as fault distribution and fault surface characteristics, providing a basis for fault feature description; Topographic data: Regional DEM data is used to reflect surface undulations and structural trends; Earthquake catalog data: covers precisely located earthquake information, providing data support for analyzing the relationship between faults and seismic activity.

[0048] The structural map data is converted into an XYZ format that can be used for spatial positioning and elevation mapping, and the elevation values ​​are reasonably filtered; the well data is mapped between depth and layer to ensure the accurate correspondence between well information and structural surfaces; the two-dimensional fault line data is numerically processed to form standard parameters that describe the spatial direction and extension length of the fault; and the earthquake catalog data is normalized to form a unified point data set.

[0049] The converted fault XYZ control point data serves as input to the 3D modeling platform. By combining fault spatial expansion, structural surface regularization, and well data correction techniques, 3D reconstruction of the subsurface structure is achieved. Using 2D fault line data, a fault plane is constructed using a near-grid triangulation interpolation algorithm. Structural surfaces are constructed based on geological data from each layer, and the model is smoothed using regular gridding techniques. Well data is used to correct stratum depth and horizon position, improving the spatial accuracy of the structural surface, ultimately forming a complete 3D regional geological structure model.

[0050] The XYZ point data of the earthquake precise positioning catalog are entered. The specific method includes: integrating the pre-processed earthquake catalog data into the three-dimensional model in the form of point data to realize the spatial positioning and visualization of earthquake activities.

[0051] The analysis of temporal trends in seismic activity, b-value changes, fault activation, and seismic risk is specifically as follows: Analyze the temporal trend of earthquake activity: Based on the spatial overlap of fault distribution characteristics in the regional 3D geological structure model and earthquake catalog data, seismological methods are used to analyze the temporal and spatial positional relationship of earthquakes through spatial distance calculation and projection. The distance between earthquake locations and faults in different time periods is calculated, and the number of earthquake events in the fault neighborhood is determined to obtain the temporal trend of earthquake activity. Analyze b-value changes: Based on the spatial overlap of fault distribution characteristics in the regional 3D geological structure model and earthquake catalog data, seismological methods are used to calculate the b-value relationship between earthquake magnitude and earthquake frequency in different time periods to determine the b-value changes; Analyze fault activation: Count the number of earthquakes detected in the year before shale gas development as the spatial background seismicity, and count the number of earthquakes in the fault neighborhood during the monitoring period. When the number of earthquakes in the fault neighborhood during the monitoring period is greater than twice the spatial background seismicity, it is determined that the fault has been activated. Analyze earthquake hazard: Obtain the fault scale and calculate the b-value of the relationship between the magnitude and earthquake frequency after industrial activities. When the b-value of the relationship between the magnitude and earthquake frequency after industrial activities is greater than the b-value threshold and the fault scale is greater than the scale threshold, earthquake hazard control is carried out.

[0052] The distance between the earthquake location and the fault is ;in, is the distance from the earthquake location to the point at the same depth on the fault plane; The fault direction.

[0053] The calculation of the b value includes: calculating the cumulative number of earthquake events of different magnitudes from the set minimum magnitude to the set maximum magnitude at intervals of 0.1 to obtain the total number of earthquakes, and then calculating the logarithm of the total number of earthquakes:

[0054] in, is the total number of earthquakes; is the minimum magnitude set; is the magnitude index; For the Magnitude range[ , ) number of earthquakes within; is the interval index within the magnitude range; is the maximum magnitude set; Constantly changing , and calculate the corresponding , get As the independent variable, is the fitted straight line of the dependent variable, and the slope of the fitted straight line is taken as the b value.

[0055] In this embodiment, based on the spatial overlap between the fault distribution characteristics of the model and the earthquake catalog data, seismological methods are used to analyze the temporal and spatial position relationship of earthquakes, calculate the distance between the earthquake location and the fault within a specified time period, and consider the fault scale and earthquake location accuracy. Earthquakes are limited to those within 200 meters of the fault to occur on that fault; The seismicity in the fault area in the two years before industrial activities were observed was analyzed to obtain the temporal and spatial background characteristics of seismic activity and to calculate the b-value of the relationship between earthquake magnitude and earthquake frequency. The b-value was calculated by calculating the cumulative number of earthquakes of different magnitudes from ML-3 to 5 earthquakes at intervals of 0.1, then calculating the logarithm of the number of earthquakes, and using a straight line with a slope to fit the change of the logarithm of the total number of earthquakes of different magnitudes to the maximum with the magnitude of the earthquake. The slope when the mean difference from the logarithmic value fitting was the smallest was the calculated b-value.

[0056] Compared with the background seismic activity, when the seismic activity near the fault increases exponentially relative to the background activity, it can be considered a clear sign of fault activation and can be judged as the fault activation.

[0057] Earthquake hazard analysis is conducted based on the b value after industrial activities and the fault scale. Generally, the b value of natural earthquakes is less than 1.0. According to the GR law and the statistical results of the b value in southern Sichuan and Zhaotong, if b>1.2, it is considered that the seismic activity is strong. If the b value fitting line intersects the horizontal axis of magnitude, the microseismic magnitude is the predicted maximum earthquake magnitude. If the b value is small, b<0.7, the intersection with the horizontal axis may have a larger magnitude, and there may be an increased risk of a larger earthquake. If the fault length is less than a few kilometers, it can be considered that the possibility of a larger earthquake is small, providing support for the judgment of earthquake hazard.

Claims

1. A water injection induced earthquake monitoring and fault activation detection system, characterized in that: It includes earthquake observation module, seismic phase data extraction module, earthquake preliminary positioning module, earthquake precise positioning module, fault and structure modeling module, and fault activation and earthquake hazard analysis module; The seismic observation module is used to obtain seismic monitoring data based on the seismic monitoring array deployed in the shale gas development area; The seismic phase data extraction module is used to identify microseisms based on seismic monitoring data to obtain final seismic phase data; The earthquake preliminary location module is used to perform earthquake event correlation and preliminary location based on the final seismic phase data to obtain a preliminary located earthquake catalog; The earthquake precise positioning module is used to perform secondary positioning based on the initially positioned earthquake catalog and the corresponding seismic phase data to obtain a precisely positioned earthquake catalog; The fault and structure modeling module is used to obtain and construct a three-dimensional geological structure model of the shale gas development area, and integrate the precisely located earthquake catalog into the three-dimensional geological structure model to perform spatial positioning and visualization of earthquake activities; The fault activation and earthquake hazard analysis module is used to analyze the temporal trend of earthquake activity, b-value changes, fault activation and earthquake hazard based on the three-dimensional geological structure model integrated from the earthquake catalog, combined with fault and earthquake data.

2. The water injection-induced earthquake monitoring and fault activation detection system according to claim 1 is characterized in that: The seismic monitoring array covers the well location area in the shale gas development area and the area within a kilometer near the well location area, and extends to c kilometers to the well location fault area. The distance between seismic stations is at least d kilometers.

3. The water injection-induced earthquake monitoring and fault activation detection system according to claim 1, characterized in that: The final seismic phase data is obtained by converting the seismic monitoring data into a day-based unit according to a standard seismic data format to obtain continuous earthquake identification input data; According to the continuous earthquake identification input data, two methods are used to identify the seismic phase: The first method: Based on the continuous earthquake recognition input data, the PhaseNet seismic phase selection algorithm is used to perform seismic phase recognition, and the first P-wave seismic phase recognition results and the first S-wave seismic phase recognition results are obtained; 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 second P-wave seismic phase identification results and the second S-wave seismic phase identification results: in, Classify the waveform to be measured into the seismic phase type The probability score of is the waveform to be measured; Phase type The cluster center of To calculate multiplication; Phase type The feature weights of Phase type ; Phase type The feature weights of is the total number of seismic phase types; is the recognition accuracy threshold corresponding to the target seismic phase; is the column normal form; 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 arrival time of the first-motion phase of the second S-wave phase identification result, the earthquake phase data in the area where the signal-to-noise ratio of the two methods is greater than the signal-to-noise ratio threshold and the travel time difference between the P-wave phase and the S-wave phase is not greater than the travel time difference threshold are retained as the final phase data; the final phase data includes the phase type and the arrival time of the first-motion phase.

4. The water injection-induced earthquake monitoring and fault activation detection system according to claim 1, characterized in that: The earthquake catalog obtained by preliminary positioning is specifically: Initialize regional crustal velocity structure models from well logging, seismic exploration, and near-seismic tomography; Based on the initialization results of the regional crustal velocity structure model, the Taup seismic wave travel time calculation software was used to calculate the P-wave initial motion phase travel time table and the S-wave initial motion phase travel time table; The grid spacing of the Taup seismic wave travel time calculation software is set to 1 / 10-1 / 5 of the station spacing; Based on the final seismic phase data, P-wave initial motion phase travel time table and S-wave initial motion phase travel time table, the REAL rapid earthquake correlation and positioning method is used to perform earthquake correlation and preliminary earthquake positioning, and obtain the P-wave phase information, S-wave phase information, earthquake occurrence events and earthquake occurrence locations belonging to the same earthquake event. Earthquake events with travel time fitting root mean square errors greater than the error threshold are removed to obtain a preliminary located earthquake catalog; the preliminary located earthquake catalog includes several earthquake events and the corresponding earthquake occurrence time, earthquake occurrence location, P-wave phase arrival time at each seismic station and S-wave phase arrival time at each seismic station.

5. The water injection-induced earthquake monitoring and fault activation detection system according to claim 1, characterized in that: The precisely located earthquake catalog is specifically: The earthquake catalog of the initial location is used as the input for earthquake precise location. Two methods are used to precisely locate each earthquake event: The first method is to mesh the regional crustal velocity structure model to obtain several model units. The theoretical travel time of P waves and S waves from each seismic station to each model unit node is calculated using the ray tracing method. Using the ray tracing method, we construct several earthquake source trajectories with arrival time constraints and arrival time difference constraints based on the P-wave phase arrival time, S-wave phase arrival time, P-wave theoretical travel time, and S-wave theoretical travel time of each seismic station in the initially located earthquake catalog. Calculate the intersection point set of the earthquake source trajectories according to the earthquake source trajectories; Determine a distribution area V of the set of intersection points of the earthquake source trajectories used for positioning according to the set of intersection points of the earthquake source trajectories; Based on the geometric centroid of the set of earthquake source trajectory intersection points in the distribution area V, the observed arrival time residual of each earthquake source trajectory intersection point in the set of earthquake source trajectory intersection points is detected, and the observed arrival time residual is subtracted from the average value of the observed arrival time residual to obtain the observed arrival time difference value of each earthquake source trajectory intersection point. The earthquake source trajectory intersection point with an observed arrival time difference value greater than the difference threshold is marked as an abnormal earthquake source trajectory intersection point; Abnormal source trajectory intersection points are removed from the source trajectory intersection point set to obtain a remaining source trajectory intersection point set, and the geometric center of gravity of the remaining source trajectory intersection point set is calculated. The geometric center of gravity of the remaining source trajectory intersection point set is used as the source position, and the source position positioning error is obtained based on the earthquake occurrence position in the initially located earthquake catalog; Second method: The initially located earthquake catalog is input into the regional crust velocity structure model, and the hypoDD earthquake event precise location method is used to accurately locate the earthquake. Based on the earthquake occurrence location in the initially located earthquake catalog, the source location error is obtained. The earthquake events whose source location errors of the two methods differ by no more than the location error threshold are retained to obtain a set of precisely located earthquake events. The set of precisely located earthquake events is divided into several groups based on the shale gas depth, and then precisely located using the hypoDD earthquake event precise positioning method to determine the earthquake occurrence location of each precisely located earthquake event. The initially located earthquake catalog is updated based on the earthquake occurrence location of each precisely located earthquake event to obtain a precisely located earthquake catalog.

6. The water injection-induced earthquake monitoring and fault activation detection system according to claim 1, characterized in that: The spatial positioning and visualization of seismic activity includes: Data acquisition operation: Obtain structural map data, well location and well layer data, two-dimensional fault line data and topographic data; Data preprocessing operations: Convert the construction map data into XYZ format for spatial positioning and elevation mapping, and filter the elevation values; Perform depth and layer mapping on well locations and well layer data to ensure that the well information corresponds to the stratigraphic layers and structural sections in the structural map data; The two-dimensional fault line data is converted into three-dimensional coordinates according to the coordinates, depth and dip angle to form the fault XYZ control point data; Normalize the data of the precisely located earthquake catalog to form XYZ point data of the earthquake location according to latitude, longitude and depth; 3D reconstruction operation: Using fault XYZ control point data as input for the 3D modeling platform, the 3D reconstruction of underground structures is achieved by combining fault spatial expansion, structural surface regularization, and well data correction technology. Specifically, the fault XYZ control point data is interpolated using a neighboring grid triangulation algorithm to construct a fine fault surface; regular gridding technology is used to smooth the structural sections in the structural map data; and the stratigraphic depth and layer position are corrected using well location and well layer position data to improve the spatial accuracy of the structural section, thereby forming a complete regional 3D geological structure model. Seismic data integration operations: The earthquake location XYZ point data is integrated into the regional three-dimensional geological structure model in the form of point data to achieve spatial positioning and visualization of earthquake activities.

7. The water injection-induced earthquake monitoring and fault activation detection system according to claim 1, characterized in that: The analysis of temporal trends in seismic activity, b-value changes, fault activation, and seismic risk is specifically as follows: Analyze the temporal trend of earthquake activity: Based on the spatial overlap of fault distribution characteristics in the regional 3D geological structure model and earthquake catalog data, seismological methods are used to analyze the temporal and spatial positional relationship of earthquakes through spatial distance calculation and projection. The distance between earthquake locations and faults in different time periods is calculated, and the number of earthquake events in the fault neighborhood is determined to obtain the temporal trend of earthquake activity. Analyze b-value changes: Based on the spatial overlap of fault distribution characteristics in the regional 3D geological structure model and earthquake catalog data, seismological methods are used to calculate the b-value relationship between earthquake magnitude and earthquake frequency in different time periods to determine the b-value changes; Analyze fault activation: Count the number of earthquakes detected in the year before shale gas development as the spatial background seismicity, and count the number of earthquakes in the fault neighborhood during the monitoring period. When the number of earthquakes in the fault neighborhood during the monitoring period is greater than twice the spatial background seismicity, it is determined that the fault has been activated. Analyze earthquake hazard: Obtain the fault scale and calculate the b-value of the relationship between the magnitude and earthquake frequency after industrial activities. When the b-value of the relationship between the magnitude and earthquake frequency after industrial activities is greater than the b-value threshold and the fault scale is greater than the scale threshold, earthquake hazard control is carried out.

8. The water injection-induced earthquake monitoring and fault activation detection system according to claim 7, characterized in that: The distance between the earthquake location and the fault is ;in, is the distance from the earthquake location to the point at the same depth on the fault plane; The fault direction.

9. The water injection-induced earthquake monitoring and fault activation detection system according to claim 7, characterized in that: The calculation of the b value includes: calculating the cumulative number of earthquake events of different magnitudes from the set minimum magnitude to the set maximum magnitude at intervals of 0.1 to obtain the total number of earthquakes, and then calculating the logarithm of the total number of earthquakes: in, is the total number of earthquakes; is the minimum magnitude set; is the magnitude index; For the Magnitude range[ , ) number of earthquakes within; is the interval index within the magnitude range; is the maximum magnitude set; Constantly changing , and calculate the corresponding , get As the independent variable, is the fitted straight line of the dependent variable, and the slope of the fitted straight line is taken as the b value.

Images