Shale gas platform site selection method and device based on seismic risk assessment, equipment, medium and product

Abstract

This application discloses a method, apparatus, equipment, medium, and product for shale gas platform site selection based on seismic risk assessment, relating to the field of platform site selection. The method includes: obtaining the fault hazard of a grid based on the shear stress and normal stress of triangular facet elements; inverting the stress field based on focal mechanism data to obtain the in-situ stress field of the triangular facet elements; clustering the maximum principal stress directions of the triangular facet elements and updating the maximum principal stress directions based on the cluster centers to obtain the maximum principal stress directions of the grid; interpolating the in-situ stress fields of the triangular facet elements to obtain the in-situ stress field of the grid; obtaining the seismic susceptibility of the grid based on factors such as the grid's historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and maximum principal stress direction; and determining the optimal well placement location based on the grid's seismic susceptibility, fault hazard, and maximum principal stress direction. This application improves the rationality of shale gas platform site selection.

Description

Technical Field

[0001] This application relates to the field of platform site selection, and in particular to a method, apparatus, equipment, medium and product for shale gas platform site selection based on seismic risk assessment. Background Technology

[0002] With the continuous expansion of shale gas exploration and development, especially in seismically active areas, the risk of induced earthquakes has gradually become a key factor restricting the green and safe development of the shale gas industry. During intensive fracturing and long-term development, factors such as underground fault activation and stress disturbance may trigger medium-intensity earthquake events, threatening operational safety, infrastructure, and the stability of residential areas. Therefore, developing earthquake susceptibility assessment methods to guide platform site selection and well layout, and improving the rationality of shale gas platform site selection, is a pressing technical bottleneck that needs to be addressed in current shale gas development. Although some results of earthquake susceptibility assessment technology have been applied to the analysis of shale gas-induced earthquakes, there is no evidence that the assessment results have been used for platform site optimization, resulting in a lack of effective connection between earthquake assessment and engineering decision-making. Summary of the Invention

[0003] The purpose of this application is to provide a method, apparatus, equipment, medium, and product for shale gas platform site selection based on seismic risk assessment, which can guide the site selection of shale gas platforms based on seismic susceptibility assessment results and improve the rationality of shale gas platform site selection.

[0004] To achieve the above objectives, this application provides the following solution: In the first aspect, this application provides a shale gas platform site selection method based on seismic risk assessment, including: gridding the study area and determining the historical seismic frequency, historical seismic energy and historical seismic activity of each grid.

[0005] The slip tendency of each triangular element is calculated based on the shear stress and normal stress of each triangular element, and the fault hazard of each grid in the study area is obtained based on the slip tendency of each triangular element.

[0006] The stress field is inverted based on the focal mechanism data corresponding to each standard earthquake event to obtain the in-situ stress field of each triangular element in the study area; the standard earthquake events are historical earthquake events in the study area that simultaneously meet the study time range, study magnitude range, and study earthquake depth range.

[0007] Clustering is performed on the maximum principal stress directions in the in-situ stress field of all triangular facets to obtain clustering results. Then, the maximum principal stress directions in the in-situ stress field of each triangular facet are updated based on the cluster centers of the clustering results to obtain the updated maximum principal stress directions of each triangular facet.

[0008] The maximum principal stress direction of each mesh is obtained by updating the maximum principal stress direction of each triangular element.

[0009] Interpolation analysis is performed on the in-situ stress field of each triangular element to obtain the interpolated in-situ stress field of each triangular element. Then, the in-situ stress field of each grid is obtained based on the interpolated in-situ stress field of each triangular element.

[0010] For any given grid, the seismic susceptibility of the grid is determined by using the grid's historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and direction of maximum principal stress as factors.

[0011] The optimal well placement location is determined based on the seismic susceptibility, fault hazard, and direction of maximum principal stress for each grid.

[0012] Secondly, this application provides a shale gas platform site selection device based on seismic risk assessment, including: a historical data determination module, used to grid the study area and determine the historical seismic frequency, historical seismic energy and historical seismic activity of each grid.

[0013] The fault hazard determination module is used to calculate the slip trend of each triangular element based on the shear stress and normal stress of each triangular element, and to obtain the fault hazard of each grid in the study area based on the slip trend of each triangular element.

[0014] The in-situ stress field inversion module is used to invert the stress field based on the focal mechanism data corresponding to each standard earthquake event, and obtain the in-situ stress field of each triangular element in the study area; the standard earthquake events are historical earthquake events in the study area that simultaneously meet the study time range, study magnitude range, and study earthquake depth range.

[0015] The maximum principal stress direction update module performs clustering processing on the maximum principal stress directions in the in-situ stress field of all triangular facets to obtain clustering results, and updates the maximum principal stress directions in the in-situ stress field of each triangular facet based on the cluster centers of the clustering results, thus obtaining the updated maximum principal stress directions of each triangular facet.

[0016] The stress field direction determination module is used to obtain the maximum principal stress direction of each mesh based on the updated maximum principal stress direction of each triangular element.

[0017] The in-situ stress field interpolation module performs interpolation analysis on the in-situ stress field of each triangular element to obtain the interpolated in-situ stress field of each triangular element. Then, based on the interpolated in-situ stress field of each triangular element, the in-situ stress field of each grid is obtained.

[0018] The earthquake susceptibility determination module is used to determine the comprehensive evaluation vector for any given grid by using the grid's historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and direction of maximum principal stress as factors, thereby obtaining the seismic susceptibility of the grid.

[0019] The optimal well placement module is used to determine the optimal well placement location based on the seismic susceptibility, fault hazard, and direction of maximum principal stress of each grid.

[0020] Thirdly, this application provides a computer device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the above-described shale gas platform site selection method based on seismic risk assessment.

[0021] Fourthly, this application provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the aforementioned shale gas platform site selection method based on seismic risk assessment.

[0022] Fifthly, this application provides a computer program product, including a computer program that, when executed by a processor, implements the aforementioned shale gas platform site selection method based on seismic risk assessment.

[0023] According to the specific embodiments provided in this application, this application has the following technical effects: This application provides a shale gas platform site selection method, apparatus, equipment, medium, and product based on seismic risk assessment. Using the historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and maximum principal stress direction of the grid as factors, a fuzzy hierarchical analysis method is used to determine the comprehensive evaluation vector, thereby obtaining the seismic susceptibility of the grid. Based on the seismic susceptibility, fault hazard, and maximum principal stress direction of each grid, the optimal well location is determined. The seismic susceptibility assessment results can guide the shale gas platform site selection, improving the rationality of shale gas platform site selection. This addresses the lack of effective connection between seismic assessment and engineering decision-making, thus providing a scientific and quantitative basis for site selection decisions and reducing the risk of induced earthquakes. Attached Figure Description

[0024] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0025] Figure 1This is a flowchart illustrating a shale gas platform site selection method based on seismic risk assessment, provided as an embodiment of this application.

[0026] Figure 2 This is a schematic diagram illustrating a shale gas platform site selection method based on seismic risk assessment, provided as an embodiment of this application.

[0027] Figure 3 A flowchart of the fuzzy hierarchical analysis method provided in an embodiment of this application.

[0028] Figure 4 This is a schematic diagram of the functional modules of a shale gas platform site selection device based on seismic risk assessment, provided in an embodiment of this application.

[0029] Figure 5 This is a schematic diagram of the structure of a computer device provided in an embodiment of this application. Detailed Implementation

[0030] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0031] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0032] In one exemplary embodiment, such as Figure 1 and Figure 2 As shown, a shale gas platform site selection method based on seismic risk assessment is provided, including the following steps 101 to 108.

[0033] Step 101: Grid the study area and determine the historical earthquake frequency, historical earthquake energy, and historical earthquake activity for each grid.

[0034] Step 102: Calculate the slip trend of each triangular element based on the shear stress and normal stress of each triangular element, and obtain the fault hazard of each grid in the study area based on the slip trend of each triangular element.

[0035] Step 103: Invert the stress field based on the focal mechanism data corresponding to each standard earthquake event to obtain the in-situ stress field of each triangular element within the study area; the standard earthquake events are historical earthquake events within the study area that simultaneously meet the study time range, study magnitude range, and study earthquake depth range. Existing earthquake susceptibility assessment methods mostly use simplified parameters or regional averages, lacking high-precision stress field modeling methods, and are difficult to reflect stress differences under local plateaus. This application solves the problem of existing methods lacking high-precision stress field modeling methods and being unable to reflect stress differences under local plateaus through stress field inversion.

[0036] Step 104: Cluster the maximum principal stress directions in the in-situ stress field of all triangular facets to obtain clustering results, and update the maximum principal stress directions in the in-situ stress field of each triangular facet based on the cluster centers of the clustering results to obtain the updated maximum principal stress directions of each triangular facet.

[0037] Step 105: Obtain the maximum principal stress direction of each mesh based on the updated maximum principal stress direction of each triangular element.

[0038] Step 106: Perform interpolation analysis on the in-situ stress field of each triangular element to obtain the interpolated in-situ stress field of each triangular element, and then obtain the in-situ stress field of each grid based on the interpolated in-situ stress field of each triangular element.

[0039] Step 107: For any grid, using the grid's historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and direction of maximum principal stress as factors, the fuzzy hierarchical analysis method is used to determine the comprehensive evaluation vector, thereby obtaining the seismic susceptibility of the grid.

[0040] Step 108: Determine the optimal well locations based on the seismic susceptibility, fault hazard, and direction of maximum principal stress for each grid.

[0041] In practical applications, before gridding the study area and determining the historical earthquake frequency, historical earthquake energy, and historical seismic activity of each grid, the following steps are taken: importing an earthquake catalog table, whose header format is: longitude, latitude, year, month, day, hour, minute, second, magnitude, depth, which is then sorted into: longitude, latitude, depth, magnitude, date.

[0042] The earthquake data corresponding to historical earthquake events are entered into the earthquake catalog table. The entered earthquake catalog table is then deduplicated and missing items are removed, leaving only the normal earthquake data.

[0043] Based on the latitude and longitude range of the study area, the study time range, the study magnitude range, and the study earthquake depth range, earthquake events that meet the criteria in the table are selected as standard earthquake events. Finally, a table consisting of earthquake data corresponding to the standard earthquake events is output, which is called the standard earthquake catalog table.

[0044] Import the focal mechanism table. The focal mechanism table must contain the fields: strike, dip, and slip angle. Optional fields are: longitude, latitude, depth, magnitude, and date.

[0045] Fill the focal mechanism data corresponding to historical earthquake events into the focal mechanism table, and then perform a deduplication operation on the focal mechanism table of the filled data.

[0046] Using the previously set latitude and longitude range, study time range, study magnitude range, and study earthquake depth range as the study area filtering conditions, and combining optional fields, the focal mechanism data that meets the conditions are selected as the focal mechanism data corresponding to the standard earthquake events; finally, a table consisting of the focal mechanism data corresponding to each standard earthquake event is output, which is called the standard focal mechanism table.

[0047] In another exemplary embodiment of this application, determining the historical earthquake frequency, historical earthquake energy, and historical seismic activity of each grid specifically includes: determining the number of standard earthquake events corresponding to each grid to obtain the historical earthquake frequency of each grid. Specifically, earthquake frequency, also known as seismic intensity, refers to the number of earthquakes that occur in a certain region and within a certain time period. To describe the local earthquake frequency in a region, the study area is gridded, a surface vector mask is created for each grid, and the number of standard earthquake events occurring in each grid is counted to obtain the historical earthquake frequency of each grid. This method can calculate the earthquake frequency distribution within three planar distributions: longitude-latitude, longitude-depth, and latitude-depth, i.e., the earthquake frequency distribution of three planar distributions in three-dimensional space.

[0048] Based on the magnitude of each standard earthquake event, the seismic energy corresponding to each standard earthquake event is obtained, and the seismic energy corresponding to each standard earthquake event is interpolated to each grid to obtain the historical seismic energy of each grid.

[0049] The historical seismic activity of each grid is calculated based on the magnitude of each standard earthquake event within each grid.

[0050] In practical applications, the seismic energy corresponding to each standard earthquake event is obtained based on its magnitude. This seismic energy is then interpolated to each grid cell to obtain the historical seismic energy for each grid cell. Specifically, seismic energy refers to the energy released during rapid fault slippage in the event of an earthquake, generally including mechanical and thermal energy. Its magnitude is directly related to the earthquake's intensity and affected area. The magnitude of a standard earthquake event is converted into energy using the Gutenberg-Richter relation, as shown in the following formula: (1) in, and They are the first The seismic energy and magnitude of a standard earthquake event.

[0051] Seismic energy propagates outward through the underground medium. Referring to the seismic motion attenuation characteristics, seismic energy decreases exponentially with increasing distance (e). Therefore, the energy of the seismic event is interpolated into a grid across the entire study area, as shown in the following formula: (2) in, Indicates the first j Historical earthquake energy of each grid It is the rate of energy attenuation during an earthquake. It is the first The standard earthquake event and the first The distance of the first grid cell is specifically based on the distance of the second grid cell. The latitude and longitude of the first standard earthquake event, compared with the first... The latitude and longitude of each grid are calculated, N represents the total number of standard earthquake events in the standard earthquake catalog table, and exp represents an exponential function with the natural constant e as the base.

[0052] This method considers projecting seismic events within a certain spatial range onto a plane, filtering seismic events near the grid using a surface vector mask, and roughly considering the geodetic arc problem when calculating distances.

[0053] In practical applications, the historical seismic activity of each grid is calculated based on the magnitudes corresponding to standard earthquake events within each grid. Specifically, this includes: seismic activity primarily targeting earthquakes... a Value and b These are two key parameters used in the analysis, describing the seismic activity of the study area. a The value represents the overall level of earthquake frequency within a given magnitude range and is positively correlated with seismic activity; bThe value represents the slope between earthquake magnitude and frequency, describing the relationship between the frequency of larger magnitude earthquakes and smaller earthquakes. For any grid within the study area, the relationship between the a-value, b-value, and earthquake frequency is estimated according to the Gutenberg-Richard relation as follows: (3) in, The magnitude in the current grid is greater than The number of standard earthquake events, Represents the minimum integrity magnitude, calculated The algorithms include: Maximum Curvature (MAXC), Goodness-of-Fit (GFT), and B-value Stabilization (MBS).

[0054] First, calculate using the maximum curvature method, goodness-of-fit test method, or b-value stability method. Value. Among them, the maximum curvature method is: calculate the maximum cumulative number of earthquake events in the magnitude bin (the magnitude sequence is divided into several bins, and the bins contain a different number of earthquakes. The cumulative number refers to the number of earthquakes in the stacked bins), and take the corresponding bin magnitude as the minimum integrity magnitude.

[0055] The goodness-of-fit test method is as follows: the minimum integrity magnitude is estimated by comparing the observed magnitude-frequency distribution with the synthesized Gutenberg-Richard distribution. The calculation formula is as follows: (4) in, The cumulative residual between observed and estimated values ​​is defined as the sum of the absolute values ​​of the differences between the observed cumulative number of events and the theoretical cumulative number of events within each magnitude interval, normalized to the sum of the observed cumulative number of events. Based on this residual index, the goodness-of-fit level is defined as 1- R Its value ranges from 0 to 1, and the larger the value, the better the fitting effect of the magnitude-frequency relationship. and They represent the first and second elements within the research grid, respectively. The cumulative distribution of theoretical magnitude-frequency and corresponding observed seismic events for each magnitude interval. and The minimum and maximum values ​​of the search interval are defined respectively, and the magnitude range is discretized by using fixed magnitude intervals (forming magnitude intervals). The symbol indicates that the absolute value is taken.

[0056] Finally, a goodness-of-fit level range is set, and the magnitude that meets the confidence level condition is selected as the minimum integrity magnitude. That is, if the goodness-of-fit level is 1- R If the set horizontal range is met, the corresponding assumed magnitude is the minimum integrity magnitude. The specific operation involves finding the assumed magnitude and iteratively calculating the goodness-of-fit level until the set horizontal range is met.

[0057] The b-value stabilization method is as follows: assuming the actual minimum integrity magnitude is below the actual minimum integrity magnitude, as the cutoff magnitude approaches the minimum integrity magnitude, the low-magnitude deviation in the cumulative frequency-magnitude distribution will be gradually eliminated, and eventually the b-value will stabilize. The cutoff magnitude at this point is the minimum integrity magnitude, and its calculation formula is as follows: (5) in, It is the mean b-value calculated up to the magnitude sequence. This represents a series of b values ​​calculated up to the magnitude sequence. The value of b is uncertain, which is calculated by the maximum likelihood method (see formula (6)).

[0058] The above three methods calculate a more accurate minimum integrity magnitude, and then the maximum likelihood method or least squares method is used to calculate a more accurate historical seismic activity: b Value and a value.

[0059] 1) Maximum likelihood estimation method.

[0060] This method can calculate the value of b using an empirical formula, as well as the uncertainty of the value of b. The specific calculation formula is as follows: (6) in, It is the average magnitude. It is the earthquake magnitude packing interval. It is the minimum integrity magnitude. It is the number of standard earthquake events in the grid with a magnitude greater than the minimum integrity magnitude.

[0061] Substitute the calculated value of b into formula (3) to obtain the value of a.

[0062] 2) Least squares method.

[0063] Calculate using the least squares method b Value and a The value is obtained by fitting the magnitude and frequency distribution with a straight line (see formula (3)) and expressed as the slope. b Value, intercept representation a Value, its b The method for estimating the uncertainty of the value is consistent with the maximum likelihood estimation method.

[0064] In another exemplary embodiment of this application, the slip trend of each triangular element is calculated based on the shear stress and normal stress of each triangular element, and the fault hazard of each grid in the study area is obtained according to the slip trend of each triangular element. Specifically, this includes: for any triangular element in the study area, calculating the ratio of the shear stress to the normal stress of the triangular element to obtain the slip trend of the triangular element.

[0065] Using the distance between each grid and each triangular element as weights, the inverse distance weighted interpolation method is used to process the slip trend of each triangular element, thereby obtaining the fault hazard of each grid.

[0066] In practical applications, for any triangular element within the study area, the ratio of shear stress to normal stress of the triangular element is calculated to obtain the slip trend of the triangular element. Specifically, this includes: (1) fault data collection, processing and standardization.

[0067] Import the 3D fault sampling data corresponding to each fault plane. The 3D fault sampling data includes point or line coordinate information describing the fault morphology, and must include a sequence number, longitude, latitude, depth, and group number. This step aims to standardize the fault data format. For one or more fault planes, the entire plane or surface is usually not stored directly; instead, the sampling point data is stored in point cloud form. Initial point data only contains longitude, latitude, and depth, without distinguishing the fault plane to which the sampling point belongs. Therefore, it is necessary to group and sort the sampling points according to the research object. Line data has a clear structure, directly containing a sequence number, longitude, latitude, depth, and group number, which can be used to identify different fault planes. For example, using the group number information (0, 1, 2, ... 9; 0, 1, 2, ..., 6), fault segments can be distinguished as fault plane 1 and fault plane 2. Subsequently, based on the longitude, latitude, and depth range of the research area, the fault planes that meet the criteria (the fault planes corresponding to the research area) are selected. Each fault plane (represented by a group number) is divided into multiple fault segments, and each segment is assigned a serial number starting from 0 and incrementing. These sampling serial numbers allow for the identification and grouping of each fault segment. Subsequently, fault planes meeting the criteria are selected based on the longitude, latitude, and depth range of the study area. Therefore, the fields for 3D fault sampling data must include: serial number, longitude, latitude, depth, and group number. Each group of fault sampling points corresponds to one or more fault segments. As the research object, a corresponding triangular discrete unit is established for each fault sampling point to describe the 3D spatial distribution of the fault plane.

[0068] A fault element model is constructed using 3D sampled point clouds within a qualified fault plane. Each element consists of three triangular elements (fault elements, patches, or subfaults) composed of three sampled points, used to describe the 3D spatial distribution of the fault plane. Through this point cloud discretization method, the fault plane can be directly generated from the sampled points without first constructing an overall plane, resulting in an element model that reflects both the overall fault morphology and finely depicts local fault features.

[0069] (2) Calculation of fault slip trend.

[0070] After dividing the fault plane into triangular facets, the following describes how to calculate the slip tendency of the fault on a single triangular facet. The probability of fault slip depends on whether the balance between the driving force and friction on the triangular facet is broken. On a triangular facet in a fractured rock stratum, there exists a normal stress in the normal direction, a shear stress perpendicular to the normal stress, and a defined coefficient of friction. According to the Mohr-Coulomb law and Byerlee's law, the stability of the triangular facet reaches the failure criterion when the ratio of shear stress to normal stress on the triangular facet is greater than the coefficient of friction. However, the presence of fluid in the rock strata creates pore pressure that can counteract the normal stress on the triangular facet, as shown in the following formula: (7) in, and These are the effective normal stress and shear stress on the triangular facet element, respectively. It is pore pressure. It is the coefficient of friction of the triangular facet. This represents the normal stress on a triangular element.

[0071] Given the vertical stress, maximum horizontal stress, and minimum horizontal stress of the triangular element, the maximum principal stress in the Mohr circle can be calculated. Intermediate principal stress and minimum principal stress stress tensor If the vertical stress > the maximum horizontal stress > the minimum horizontal stress, then the vertical stress = Maximum horizontal stress = Minimum horizontal stress = This simply means sorting by size. Using the geometric parameters of the triangular elements, the normal direction of each triangular element can be calculated, and then determined based on the rotation matrix. Transform the stress tensor from a Cartesian reference frame (X, Y, Z) into a triangular facet reference frame (normal vector, orientation, tilt angle). The formula is as follows: (8) in, yes The axis of rotation between the unit vector of the direction and the plane normal vector. yes The angle between the unit vector of the direction and the plane normal vector. It is the identity matrix. yes The cross product matrix, yes The outer product, express The transpose of .

[0072] Finally, the calculated It is a 3x3 matrix, and the first row of the first column represents the normal stress in the normal direction of the triangular facets, which corresponds to the magnitude of the normal stress in formula (7). The second and third rows of the first column correspond to the two stress components in the tangential direction of the triangular element. The square root of their sum gives the magnitude of the shear stress in formula (7). Subsequently, by combining formulas (7) and (8), the current fluid pore pressure can be determined. and coefficient of friction Under certain conditions, the ratio of shear stress to effective normal stress in the stress tensor of the triangular facet element is calculated to determine the slip tendency of each triangular facet element.

[0073] In practical applications, after the study area is meshed, the attenuation value can be calculated based on the spatial distance between the mesh and the triangular elements, typically using the inverse ratio of distance as the attenuation coefficient. Simultaneously, the slip trend of the triangular elements is used as a weight for weighting. Specifically, an inverse distance-weighted interpolation method is used to interpolate the slip trend of each triangular element onto the mesh nodes, ultimately obtaining the fault hazard distribution for each mesh. Taking one mesh as an example: (9) in, Indicates the first Fault hazard of each grid, It is the first fault plane in the study area The slippage trend of each triangular element It is the first The triangular facet element and the first j The distance between grid cells is specifically calculated based on the latitude and longitude of the cell center and the latitude and longitude of the grid cells. N This indicates the number of triangular facets present in this fault.

[0074] In practical applications, the stress field is inverted based on the source mechanism data corresponding to each standard earthquake event to obtain the in-situ stress field of each triangular surface element in the study area, specifically including: (3) stress field inversion based on source mechanism.

[0075] The purpose of focal mechanism inversion of the stress field is to deduce the stress state within a region from seismic data. Assuming uniform structural stress within the region, the earthquake occurring on an existing fault, with the slip vector pointing towards the shear stress direction on the fault, the four parameters of the stress tensor can be determined using existing stress inversion algorithms: three defined principal stress directions (…). , and (azimuth and dip angles) and shape ratio The formula for shape ratio is as follows: (10) In the stress field inversion, expressions for the normal stress and shear stress on the triangular element were established, and based on the Wallace-Bott assumption, the direction of the shear stress and its relationship with the fault slip direction were determined. Consistent, and further assuming shear stress on the same activated triangular element. The seismic events in all studies remain consistent. Therefore, the triangular element normal... Sliding direction and stress tensor of quantity The relationship is as follows: (11) in, These are all components of the stress tensor, with the superscript T indicating transpose. It originates from the normal of the triangular facet. A 3×5 matrix: (12) in, n 1 , n 2 and n 3 These represent the three components of the triangular facet normal: horizontally eastward, horizontally northward, and vertically upward from the Earth's surface.

[0076] Next, adopt Generalized linear inversion under norms solves for : (13) in, express The generalized inverse matrix, because It is not a square matrix, so it needs to be solved using least squares. .

[0077] Ultimately, through Construct a 3x3 deviatoric stress tensor matrix, calculate the eigenvalues ​​of this matrix (a 1x3 vector), and then... Determine the direction of the stress field Substituting into formula (10) will allow you to calculate the shape ratio.

[0078] However, this inversion method, when the orientation of the triangular facet is unknown, has a 50% probability of selecting either the triangular facet or an auxiliary facet, and the triangular facet and auxiliary facet are orthogonal. Therefore, fault instability is introduced to improve the accuracy of selecting triangular facets during the inversion. The formula is as follows: (14) in, and These represent the normal stress and shear stress calculated on the current analysis plane, respectively. and The distribution represents the shear stress and effective normal stress of the optimal triangular element. It is the friction coefficient of the triangular facet element. Equation (12) shows that it is independent of the absolute stress value, and the fault instability can be determined from... Shape ratio The direction cosine of the angle between the defined principal stress axis and the triangular element (i.e. the component of the fault plane normal in the principal stress coordinate system, consistent with formula (12)). To evaluate this, we scale the stress tensor: (15) At this point, the initial value is calculated using formulas (11), (12), and (13). and Then, by scaling the initial value to the interval [-1, 1], we can obtain... And thus calculate the shape ratio This practice is consistent with... Substitute into formula (10) to calculate the shape ratio The principle is the same.

[0079] Therefore, for and : (16) Substituting formula (16) into formula (14), the expression for fault instability is rewritten as: (17) At this point, formula (17) becomes the final formula for calculating fault instability. Furthermore, based on the strike, dip, and slip angle of the imported focal mechanism, the normal to the triangular element is calculated and decomposed into three directions: horizontal east, horizontal north, and vertical upward from the surface. This decomposes the normal vector of the triangular element into its corresponding vector. coefficient of friction and shape ratio You can calculate first. and Thus, the fault instability can be determined. However, in actual calculations of fault instability, a friction coefficient is required, typically ranging from 0.2 to 0.8. Therefore, a mesh search method is used to iterate through multiple friction coefficients, calculate the fault instability for each coefficient, and select the friction coefficient that maximizes fault instability as the optimal value. Subsequently, the principal stress direction corresponding to this optimal friction coefficient is determined. Substitute into formula (10) to calculate the shape ratio This determines that the plane being analyzed is a triangular facet rather than an auxiliary facet, thus determining the orientation, tilt angle, and slip angle of the triangular facet. Simultaneously, it determines the directions of the three principal stress axes. Each is represented by its azimuth and inclination angle.

[0080] In practical applications, the clustering results are obtained by clustering the directions of the maximum principal stress in the in-situ stress field of all triangular facets. Specifically, the K-means clustering algorithm is used to cluster the directions of the maximum principal stress in the in-situ stress field of all triangular facets.

[0081] In another exemplary embodiment of this application, the maximum principal stress direction in the in-situ stress field of each triangular facet element is updated according to the cluster center of the clustering results to obtain the updated maximum principal stress direction of each triangular facet element. Specifically, this includes: for any cluster center, updating the maximum principal stress direction in the in-situ stress field of the target triangular facet element according to the cluster center to obtain the updated maximum principal stress direction in the in-situ stress field of each target triangular facet element; the target triangular facet element is the triangular facet element corresponding to the cluster center. Specifically, the maximum principal stress direction in the in-situ stress field of the target triangular facet element is changed to the cluster center.

[0082] In practical applications, the maximum principal stress direction of each mesh is obtained based on the updated maximum principal stress direction of each triangular element. Specifically, this involves dividing the entire study area into several meshes, calculating the distance between each mesh and each triangular element, selecting triangular elements within a certain range of their vicinity, and calculating the average of the maximum principal stress directions of these triangular elements, assigning this average to the mesh. If no triangular elements are found near the mesh, a null value is assigned, indicating that it does not possess a maximum principal stress direction.

[0083] In another exemplary embodiment of this application, the in-situ stress field of each triangular facet element is interpolated to obtain the in-situ stress field of each triangular facet element after interpolation. Specifically, this includes: sequentially using the Kriging interpolation method and the spatial variability interpolation method to interpolate the in-situ stress field of each triangular facet element to obtain the in-situ stress field of each triangular facet element after interpolation.

[0084] In practical applications, Kriging interpolation and spatial variability interpolation are used sequentially to analyze the in-situ stress field of each triangular facet element, yielding the interpolated in-situ stress field. Specifically, the in-situ stress field is an important physical quantity describing the stress state of underground rock masses, resulting from the combined effects of natural stress and external forces. Studying the in-situ stress field allows for a better understanding of crustal mechanical behavior and the prediction of geological hazards. Using Kriging interpolation to analyze the in-situ stress field is a commonly used and scientific method, considering spatial correlation and providing an estimate of the uncertainty of the predicted values.

[0085] In-situ stress field data typically include: longitude, latitude, depth, and maximum principal stress. Intermediate principal stress Minimum principal stress , Azimuth, Azimuth, Azimuth, tendency, tendency, The study considers factors such as tendency, measurement time, and measurement methods. First, it ensures that the in-situ stress field data covers the target interpolation region; otherwise, the interpolation region is reduced to eliminate uncertainties in uncovered areas. Then, it observes whether the stress values ​​in the interpolation region are depth-dependent. If a relationship exists, generalized kriging is used to determine the maximum principal stress. Intermediate principal stress Minimum principal stress , Azimuth, Azimuth, Azimuth, tendency, Tendency and Prefers interpolation; conversely, uses ordinary kriging for the maximum principal stress. Intermediate principal stress Minimum principal stress , Azimuth, Azimuth, Azimuth, tendency, Tendency and It tends to perform interpolation.

[0086] Based on the spatial distance between observation points and their stress values, where observation points refer to the sample locations with stress values, the spatial variability of the in-situ stress field at each observation point is calculated. The variability function includes spherical models, exponential models, and Gaussian models. The formula for the variability function is as follows: (18) The formulas above represent the variograms of the spherical model, exponential model, and Gaussian model, respectively. This is the theoretical variogram value, representing the distance. The variability between two points. and These are the baseline variation and the degree of variation, respectively. The range is used to describe the rate of decay of variability, and exp represents an exponential function with the natural constant e as the base.

[0087] The relationship between the above variogram and the actual observed values ​​is as follows: (19) in, and They are in the position and The observed values, The localized variable sample point focal distance is The number of point pairs.

[0088] The specific calculation process is as follows: after importing the stress data, first select formula (19) to calculate its experimental variation value. Note that formula (19) has the following characteristics: It is the experimental variation value, which quantifies the change in real data with distance. The spatial differences vary, therefore, a series of spatial difference points will be output. Next, substitute the spatial difference points into formula (18) to fit the theoretical variogram, noting the following here: The variance is the value of the fitted theoretical model. The parameters in formula (18) can be solved by least squares or nonlinear fitting to obtain the fitted theoretical variance function. .

[0089] Then, the unknown points (Interpolation point) Import theoretical variability function Thus, the Kriging equation is constructed: (20) in, It is the Kriging weight. It is a Lagrange multiplier. This indicates the number of known sample points involved in the interpolation. and These are known samples and subscript, This represents the unknown point to be interpolated. A total of [number] points were generated. Equations There are several unknowns. By constructing the covariance matrix between known sample points and calculating the covariance between known sample points and unknown points, a Kriging weighted matrix equation is established. By solving this system of linear equations, the contribution weight of each known data point to the target point can be obtained. The weights minimize the variance of the estimation error between the predicted and observed values ​​under unbiased constraints (see formula (20)). Finally, a weighted sum is applied to all known points to obtain the predicted value for that point, as shown in the following formula: (twenty one) in, It is the predicted value of regional change at a certain point. It is the first The weight coefficients of the known points It is the first Observations at known points The number of samples is known. These are the predicted values ​​for the points to be interpolated.

[0090] The formula for the error between the predicted value and the actual observed value at this point is as follows: (twenty two) Based on the error value, it is determined whether Kriging interpolation and spatial variability difference need to be performed again. Finally, after Kriging interpolation of the stress value or principal stress direction of the in-situ stress data, the interpolation results are compared with the stress values ​​and other data of the observation points that did not participate in the interpolation. This allows for the analysis and evaluation of the rock mass stability, induced earthquake risk, and impact of underground engineering activities in the study area.

[0091] In practical applications, the historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and direction of maximum principal stress of the grid are used as factors. A fuzzy hierarchical analysis method is employed to determine the comprehensive evaluation vector, thereby obtaining the seismic susceptibility of the grid, specifically including: This step will construct a fuzzy hierarchical model to quantitatively analyze the seismic susceptibility of shale gas extraction areas and classify their hazard levels. For example... Figure 3 As shown, it consists of the following 5 steps: (1) Divide the levels based on the main factors.

[0092] The index layers are defined as the primary and secondary components of stress accumulation, as well as the influence of historical earthquakes, denoted by U1 to U3 respectively. U1 includes the in-situ stress field (U... 11 ) and the direction of maximum principal stress (U 12 U2 includes historical earthquake frequencies (U...). 21 ) and historical earthquake energy (U 22 U3 includes historical seismic activity (U3). 31 ) and fault hazard (U32 The final evaluation conclusion uses three indicators: frequent small and medium earthquakes correspond to high-risk rating V1; frequent small earthquakes and few medium earthquakes correspond to medium-risk rating V2; and both small and medium earthquakes are rare, corresponding to low-risk rating V3.

[0093] (2) Establish the fuzzy membership matrix.

[0094] The core task of the fuzzy comprehensive evaluation method is to construct a membership matrix by establishing fuzzy relationships between each factor and the evaluation layer based on the characteristics of each factor. The evaluation layer is divided into V1-V3 levels, and the factor layer is represented by U... ij In other words, using R i Indicator layer U i The membership matrix has the following form: (twenty three) r 12 U i The degree to which the first factor belongs to the second comment has a range of (0,1).

[0095] A weighted average fuzzy evaluation model (·,⊕) is adopted to balance all factors according to their weights.

[0096] (twenty four) in, This represents the weight of the j-th comment. Indicates the first k The weights of each factor, where m represents the total number of factors and n represents the total number of comments.

[0097] The frequency normalization method was used to calculate membership degrees. This method first selects appropriate parameters based on the characteristics of each factor to determine its different values ​​for each comment. Then, the specific algorithm for the final score of the data in each grid cell is to use the frequency of each type of comment received by the surrounding nine squares as the three membership degrees of the central square. Ultimately, a membership degree matrix for all factors is generated in each small square. For example, if 1, 2, and 6 squares in the surrounding nine squares receive ratings V1, V2, and V3 respectively, then the score of the central square is 1 / 9, 2 / 9, and 6 / 9. The advantage of this calculation is that it weakens the subjective definition of factor boundaries while simultaneously normalizing the data. When determining the initial score for a factor in each grid cell, it is essential to thoroughly investigate the calculation method and classification criteria for that factor.

[0098] (3) Weight calculation based on AHP.

[0099] An expert scoring method was used to conduct three rounds of importance rating: once at the indicator level and once each at the two factor levels. Since it was difficult for the scoring experts to accurately express relative importance using numbers between 1 and 9, a three-scale method was employed to determine the judgment matrix. C kl This makes it easy for experts to provide an accurate and simple comparison of importance.

[0100] (25) Then calculate the importance ranking index of each factor. r k .

[0101] (26) Next, we calculate the pairwise comparison matrix under the new scale, where any comparison coefficient... Calculated this way: (27) in b kl =r max -r min .

[0102] use b kl Construct an n×n pairwise comparison matrix B, and then use the eigenvector corresponding to its largest eigenvalue to represent the corresponding weight vector: (28) Among them, here λ max for n Comparison matrix B The largest eigenvalue, ω for λ max The corresponding eigenvector.

[0103] After obtaining the eigenvectors, a consistency check needs to be performed, and the consistency ratio coefficient needs to be calculated. CRC : (29) RI The value can be obtained by querying the consistency check table. When the final calculated CRC is less than 0.1, the feature vector is considered to have passed the consistency check.

[0104] (4) Comprehensive evaluation vector.

[0105] use W 1 = ( ω11 ,ω 12 ,ω 13 ) and W 2 = ( ω 21 ,ω 22 ,ω 23 The following represents the weight vectors between factors in indicator layers U1 and U2, respectively. Next, the intermediate weight vectors between the two indicator layers are calculated: (30) Then, using n The intermediate evaluation matrix is ​​constructed from the weight vectors: (31) Finally, according to the multi-level fuzzy comprehensive evaluation model, using W= ( ω 1 , ω 2， ) represents the indicator layer U i The weight vectors between the indicators are used to perform a composite operation on the indicator layer to obtain the comprehensive evaluation vector, which is then: (32) Here, A is the final evaluation vector. ω 12 This indicates the weight of the second factor in indicator layer U1.

[0106] In another exemplary embodiment of this application, the optimal well placement location is determined based on the seismic susceptibility, fault hazard, and maximum principal stress direction of each grid. Specifically, this includes obtaining a seismic susceptibility score for each grid based on its seismic susceptibility.

[0107] The fault hazard distribution score for each grid is obtained based on the fault hazard of each grid within the study area.

[0108] The stress field orientation suitability score for each grid is obtained based on the direction of the maximum principal stress of each grid.

[0109] The optimal well locations are determined based on the seismic susceptibility score, fault hazard distribution score, and stress field orientation suitability score for each grid.

[0110] Fuzzy hierarchical analysis was used to quantitatively analyze multiple factors, including historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and the direction of maximum principal stress. Based on expert scoring, earthquake susceptibility and hazard assessments were performed. Finally, the spatial distribution characteristics of each factor were constructed into a standardized evaluation vector, which was then linearly weighted and summed with a normalized fuzzy weight vector to output a comprehensive earthquake hazard index. This index comprehensively reflects the potential seismic induced risk level of the study area under the background of seismic activity and geological structure. It can be used to identify high-risk grid areas and to provide avoidance or set engineering control boundaries during well site selection. Therefore, in practical applications, the seismic susceptibility of the grid is considered... The seismic susceptibility score of the grid is obtained, specifically by: according to the formula (33) Calculate the seismic susceptibility score of the grid. . The lower the better. This represents the normalization function.

[0111] Before assessing fault hazard, it is necessary to estimate the pore pressure and horizontal and vertical in-situ stress near the fault in the study area, and calculate the probability of fault instability based on fault slip theory. To more effectively reflect the spatial attenuation characteristics of the influence of triangular elements, an inverse distance weighting method is introduced to construct a fault hazard distribution map, which can capture the asymmetric control of triangular elements on surrounding grid cells—that is, the closer to the triangular element, the higher the slip risk. This method integrates the geometric shape and stress field characteristics of triangular elements while maintaining high spatial resolution, highlighting the criticality of triangular elements under existing stress states and reflecting their potential slip-induced seismic sensitivity under external disturbances (such as fracturing and water injection). This index is used in the spatial avoidance and engineering control strategy design during the platform well site selection stage, prioritizing the avoidance of areas near high-risk triangular elements, or developing supporting monitoring and mitigation measures to reduce the probability of seismic induction. Therefore, in practical applications, the fault hazard of each grid in the study area is considered. The fault hazard distribution score for each grid is obtained as follows: based on the formula... (34) Calculate the fault hazard distribution score of the grid. , The lower the better.

[0112] The direction of the maximum principal stress in the inverted region is determined by the focal mechanism, and K-means clustering is used to divide the spatial stress field distribution into main controlling sub-regions, extracting more locally representative principal stress direction features. Simultaneously, combined with the spatial interpolation results of the in-situ stress field, it is recommended that horizontal wells be designed with a layout perpendicular to the direction of the maximum principal stress to facilitate the propagation of fracturing fractures along the principal stress control direction, thereby improving fracturing efficiency. Therefore, in practical applications, the direction of the maximum principal stress in each grid should be considered. The stress field orientation suitability score for each grid is obtained by calculating the angle between the maximum principal stress direction and the candidate well section direction. Specifically, the stress field orientation suitability score is calculated by subtracting the angle between the maximum principal stress direction and the candidate well section direction.

[0113] Based on the suitability of the stress field direction of each grid The stress field orientation suitability score of each grid is obtained. The specific formula is as follows: (35), The closer to vertical, the better.

[0114] Both earthquake hazard and fault risk are treated linearly, while the stress field orientation suitability is treated with a quadratic term to reduce the influence of the direction of the maximum principal stress.

[0115] In another exemplary embodiment of this application, the optimal well placement location is determined based on the seismic susceptibility score, fault hazard distribution score, and stress field direction suitability score of each grid. Specifically, the seismic susceptibility score, fault hazard distribution score, and stress field direction suitability score corresponding to the grid are weighted and summed to obtain the suitability score for well placement in each grid. The specific formula is as follows: (36) in, , and yes , and The weights need to satisfy The recommended initial value is: 5, , .

[0116] The optimal well placement location is determined based on the suitability score of each grid's well locations. Specifically, if the final suitability score is not less than 0.7, the grid is considered a high-quality candidate site; if the suitability score is below 0.3, site selection is not recommended. The specific weights and suitability score thresholds are adjusted according to the actual situation.

[0117] This application constructs a regional gridded well placement suitability evaluation system based on multi-source geological and seismic information. Combining seismic susceptibility, fault hazard, and the direction of maximum principal stress, it enables the optimal placement of shale gas platforms. It can propose scientific planning for optimizing shale gas platform site selection from the perspectives of seismic susceptibility and seismic hazard.

[0118] Compared with existing technologies, the shale gas platform site selection method proposed in this application has the following main advantages: 1. High assessment accuracy.

[0119] Current technologies often rely solely on single factors such as historical earthquake frequency or fault distribution for earthquake risk assessment, neglecting tectonic stress environment and the interaction of multiple factors, resulting in low accuracy of assessment results. This application, however, has been approved. (1) By combining fault geometry information with regional stress state, the fault slip trend can be calculated, which improves the ability to quantify tectonic activity and induced risks.

[0120] (2) Based on the source mechanism, the real underground stress field is inverted and then interpolated to make the stress information more consistent with the actual underground structural conditions.

[0121] (3) The fuzzy hierarchical analysis method is adopted to integrate multiple factors and introduce expert judgment and membership function, which effectively overcomes the uncertainty of geological information and significantly improves the scientificity and accuracy of earthquake susceptibility assessment.

[0122] 2. It has strong applicability and allows for more reasonable platform location selection.

[0123] This application not only provides a complete static assessment system for seismic susceptibility, but also further applies the assessment results to platform site selection decisions, constructing a comprehensive well placement evaluation index system based on seismic susceptibility, fault hazard, and the direction of maximum principal stress. Through a gridded regional assessment method, more detailed platform site selection and well placement optimization can be achieved, overcoming the problem of "disconnect between assessment and engineering decision-making" in existing technologies, and improving the engineering usability and regional adaptability of the scheme.

[0124] 3. High processing efficiency and strong process integration.

[0125] Current methods suffer from low data standardization and fragmented processing workflows, hindering large-scale deployment and practical engineering applications. This application introduces standardized and normalized processing methods for seismic and fault data, facilitating the automatic reading, interpolation, and analysis of large-scale geological and seismic data. This addresses the integration and quantification challenges of multi-source uncertain geological information fusion processing, and promotes large-scale deployment and practical engineering applications. Furthermore, the slip trend calculation and stress field inversion employ modular mathematical modeling methods, facilitating integration into GIS platforms or other engineering support systems. This enables efficient calculation and automated analysis from data collection to platform site selection, significantly improving workflow efficiency.

[0126] In summary, through technological integration and methodological innovation, this application significantly outperforms related technologies in terms of assessment accuracy, platform site selection rationality, and processing efficiency, and is particularly suitable for seismic risk prediction and site selection planning in shale gas extraction areas.

[0127] Based on the same inventive concept, this application also provides a shale gas platform site selection device based on seismic risk assessment for implementing the above-mentioned shale gas platform site selection method based on seismic risk assessment. The solution provided by this device is similar to the solution described in the above method. Therefore, the specific limitations of one or more embodiments of the shale gas platform site selection device based on seismic risk assessment provided below can be found in the limitations of the shale gas platform site selection method based on seismic risk assessment described above, and will not be repeated here.

[0128] In one exemplary embodiment, such as Figure 4 As shown, a shale gas platform site selection device based on seismic risk assessment is provided, including: a historical data determination module, used to grid the study area and determine the historical seismic frequency, historical seismic energy and historical seismic activity of each grid.

[0129] The fault hazard determination module is used to calculate the slip trend of each triangular element based on the shear stress and normal stress of each triangular element, and to obtain the fault hazard of each grid in the study area based on the slip trend of each triangular element.

[0130] The in-situ stress field inversion module is used to invert the stress field based on the focal mechanism data corresponding to each standard earthquake event, and obtain the in-situ stress field of each triangular element in the study area; the standard earthquake events are historical earthquake events in the study area that simultaneously meet the study time range, study magnitude range, and study earthquake depth range.

[0131] The maximum principal stress direction update module performs clustering processing on the maximum principal stress directions in the in-situ stress field of all triangular facets to obtain clustering results, and updates the maximum principal stress directions in the in-situ stress field of each triangular facet based on the cluster centers of the clustering results, thus obtaining the updated maximum principal stress directions of each triangular facet.

[0132] The stress field direction determination module is used to obtain the maximum principal stress direction of each mesh based on the updated maximum principal stress direction of each triangular element.

[0133] The in-situ stress field interpolation module performs interpolation analysis on the in-situ stress field of each triangular element to obtain the interpolated in-situ stress field of each triangular element. Then, based on the interpolated in-situ stress field of each triangular element, the in-situ stress field of each grid is obtained.

[0134] The earthquake susceptibility determination module is used to determine the comprehensive evaluation vector for any given grid by using the grid's historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and direction of maximum principal stress as factors, thereby obtaining the seismic susceptibility of the grid.

[0135] The optimal well placement module is used to determine the optimal well placement location based on the seismic susceptibility, fault hazard, and direction of maximum principal stress of each grid.

[0136] In one exemplary embodiment, a computer device is provided, which may be a server or a terminal, and its internal structure diagram may be as follows. Figure 5 As shown, the computer device includes a processor, memory, input / output (I / O) interfaces, and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is also connected to the system bus via the I / O interfaces. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and a database. The internal memory provides the environment for the operating system and computer programs in the non-volatile storage media to run. The database stores shale gas platform site selection data. The I / O interfaces are used for information exchange between the processor and external devices. The communication interface is used for communication with external terminals via a network connection. When the computer program is executed by the processor, it implements a shale gas platform site selection method based on seismic risk assessment.

[0137] Those skilled in the art will understand that Figure 5 The structure shown is merely a block diagram of a portion of the structure related to the present application and does not constitute a limitation on the computer device to which the present application is applied. Specific computer devices may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0138] In one exemplary embodiment, a computer device is provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the above-described method embodiments.

[0139] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the above-described method embodiments.

[0140] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the above-described method embodiments.

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

[0142] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM).

[0143] The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc., and are not limited to these.

[0144] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0145] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A shale gas platform site selection method based on seismic risk assessment, characterized in that, The shale gas platform site selection method based on seismic risk assessment includes: The study area was gridded, and the historical earthquake frequency, historical earthquake energy, and historical earthquake activity of each grid were determined. The fault plane corresponding to the study area is divided into multiple triangular facets; The slip tendency of each triangular element is calculated based on the shear stress and normal stress of each triangular element, and the fault hazard of each grid in the study area is obtained based on the slip tendency of each triangular element. The stress field is inverted based on the focal mechanism data corresponding to each standard earthquake event to obtain the in-situ stress field of each triangular element in the study area; the standard earthquake events are historical earthquake events in the study area that simultaneously meet the study time range, study magnitude range, and study earthquake depth range. Clustering is performed on the maximum principal stress directions in the in-situ stress field of all triangular facets to obtain clustering results. The maximum principal stress directions in the in-situ stress field of each triangular facet are updated based on the cluster centers of the clustering results to obtain the updated maximum principal stress directions of each triangular facet. The maximum principal stress direction of each mesh is obtained based on the updated maximum principal stress direction of each triangular element; Interpolation analysis is performed on the in-situ stress field of each triangular facet to obtain the interpolated in-situ stress field of each triangular facet. Then, the in-situ stress field of each grid is obtained based on the interpolated in-situ stress field of each triangular facet. For any given grid, the comprehensive evaluation vector is determined using the grid's historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and direction of maximum principal stress, thus obtaining the grid's seismic susceptibility. The optimal well placement location is determined based on the seismic susceptibility, fault hazard, and direction of maximum principal stress for each grid.

2. The shale gas platform site selection method based on seismic risk assessment according to claim 1, characterized in that, Determine the historical earthquake frequency, historical earthquake energy, and historical earthquake activity for each grid, specifically including: The historical earthquake frequency of each grid is obtained by determining the number of standard earthquake events corresponding to each grid. Based on the magnitude of each standard earthquake event, the seismic energy corresponding to each standard earthquake event is obtained, and the seismic energy corresponding to each standard earthquake event is interpolated to each grid to obtain the historical seismic energy of each grid. The historical seismic activity of each grid is calculated based on the magnitude of each standard earthquake event within each grid.

3. The shale gas platform site selection method based on seismic risk assessment according to claim 1, characterized in that, The slip tendency of each triangular element is calculated based on its shear stress and normal stress. Based on this slip tendency, the fault hazard of each grid within the study area is obtained, specifically including: For any triangular element within the study area, calculate the ratio of shear stress to normal stress of the triangular element to obtain the slip tendency of the triangular element. Using the distance between each grid and each triangular element as weights, the inverse distance weighted interpolation method is used to process the slip trend of each triangular element, thereby obtaining the fault hazard of each grid.

4. The shale gas platform site selection method based on seismic risk assessment according to claim 1, characterized in that, Based on the cluster centers of the clustering results, the direction of the maximum principal stress in the in-situ stress field of each triangular element is updated to obtain the updated direction of the maximum principal stress for each triangular element, specifically including: For any cluster center, the updated maximum principal stress direction in the in-situ stress field of each target triangular element is obtained based on the maximum principal stress direction in the in-situ stress field of the updated target triangular element; the target triangular element is the triangular element corresponding to the cluster center.

5. The shale gas platform site selection method based on seismic risk assessment according to claim 1, characterized in that, Interpolation analysis is performed on the in-situ stress field of each triangular facet element to obtain the interpolated in-situ stress field of each triangular facet element, specifically including: The in-situ stress field of each triangular element was interpolated using Kriging interpolation and spatial variability interpolation in sequence to obtain the in-situ stress field of each triangular element after interpolation.

6. The shale gas platform site selection method based on seismic risk assessment according to claim 1, characterized in that, The optimal well placement is determined based on the seismic susceptibility, fault hazard, and direction of maximum principal stress for each grid, specifically including: Based on the seismic susceptibility of the grid, a seismic susceptibility score for the grid is obtained; The fault hazard distribution score for each grid is obtained based on the fault hazard of each grid within the study area; The stress field orientation suitability score for each grid is obtained based on the direction of the maximum principal stress for each grid. The optimal well locations are determined based on the seismic susceptibility score, fault hazard distribution score, and stress field orientation suitability score for each grid.

7. A shale gas platform site selection device based on seismic risk assessment, characterized in that, The shale gas platform site selection device based on seismic risk assessment includes: The historical data determination module is used to grid the study area and determine the historical earthquake frequency, historical earthquake energy, and historical earthquake activity of each grid. The fault hazard determination module is used to calculate the slip trend of each triangular element based on the shear stress and normal stress of each triangular element, and to obtain the fault hazard of each grid in the study area based on the slip trend of each triangular element. The in-situ stress field inversion module is used to invert the stress field based on the focal mechanism data corresponding to each standard earthquake event, and obtain the in-situ stress field of each triangular element in the study area; the standard earthquake events are historical earthquake events in the study area that simultaneously meet the study time range, study magnitude range, and study earthquake depth range. The maximum principal stress direction update module performs clustering processing on the maximum principal stress directions in the in-situ stress field of all triangular facets to obtain clustering results, and updates the maximum principal stress directions in the in-situ stress field of each triangular facet based on the cluster center of the clustering results to obtain the updated maximum principal stress directions of each triangular facet. The stress field direction determination module is used to obtain the maximum principal stress direction of each mesh based on the updated maximum principal stress direction of each triangular element. The in-situ stress field interpolation module performs interpolation analysis on the in-situ stress field of each triangular element to obtain the in-situ stress field of each triangular element after interpolation. Then, the in-situ stress field of each grid is obtained based on the in-situ stress field of each triangular element after interpolation. The earthquake susceptibility determination module is used to determine the comprehensive evaluation vector for any given grid by taking historical seismic activity, fault hazard, in-situ stress field, historical earthquake frequency, historical earthquake energy, and direction of maximum principal stress as factors, and thus obtain the earthquake susceptibility of the grid. The optimal well placement module is used to determine the optimal well placement location based on the seismic susceptibility, fault hazard, and direction of maximum principal stress of each grid.

8. A computer device, comprising: A memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor executes the computer program to implement the shale gas platform site selection method based on seismic risk assessment as described in any one of claims 1-6.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the shale gas platform site selection method based on seismic risk assessment as described in any one of claims 1-6.

10. A computer program product, comprising a computer program, characterized in that, When executed by a processor, the computer program implements the shale gas platform site selection method based on seismic risk assessment as described in any one of claims 1-6.

Images