Quantitative prediction and integrated modeling method for four-dimensional stress field of oil and gas reservoir

Abstract

The invention relates to the field of oil and gas field exploration and development, in particular to an oil and gas reservoir four-dimensional stress field quantitative prediction and integrated modeling method. Establishing finite element geological geometric models in different periods by recovering a paleostructure form in a key structure transformation period; predicting three-dimensional distribution of rock mechanics parameters by combining a rock mechanics experiment and adopting a well-seismic combination method; burial depth changes of rocks in different periods are determined through ancient and modern structure fusion, three-dimensional distribution of mechanical parameters of the ancient rocks is predicted, geomechanical heterogeneous models of oil and gas reservoirs in different periods are established, and numerical simulation of four-dimensional stress fields in different periods is achieved; and a multi-period stress unit cyclic coupling algorithm is adopted to realize coupling of multi-period stress field prediction data into one model. The invention provides a quantitative prediction and integrated modeling method for a four-dimensional stress field of an oil and gas reservoir. The method has practical significance in multiple aspects of ancient stress field numerical simulation, tectonic fracture distribution prediction, 'geology and engineering dessert 'evaluation and the like.

Description

Technical Field

[0001] The invention relates to the field of oil and gas field exploration and development, and in particular to a method for quantitatively predicting and integrally modeling a four-dimensional stress field of an oil and gas reservoir. Background Art

[0002] In the past decade, the exploration and development of traditional oil and gas resources has reached a bottleneck, and unconventional, deep-ultra-deep oil and gas resources have gradually attracted widespread attention from the industry and academia. Among them, deep sand / mudstone interbedded reservoirs and marine buried hill reservoirs have become the key areas and hotspots of oil and gas geological research. Unlike foreign oil and gas reservoirs, my country's oil and gas reservoirs are mainly continental deposits, usually experienced multiple tectonic movements, and generally have the characteristics of complex stress evolution process, strong reservoir heterogeneity, and multi-scale structural fracture development. At present, the geomechanical modeling method is an effective means to realize the numerical simulation of the stress field in the study area. Reservoir geomechanical modeling has gone through three stages: two-dimensional geomechanical modeling, simple three-dimensional geomechanical modeling, and three-dimensional geomechanical heterogeneous modeling. Stage I, the two-dimensional geomechanical modeling stage (1968-2008), the model is mainly based on the curved shell model and the thin plate model, the fault is mainly line or surface unit, contains little geological information, and cannot realize the simulation of three-dimensional stress. There are theoretical defects in simulating the mechanism of fault and fracture genesis and the current geostress field. Phase II, simple 3D geomechanical modeling phase (2008-2017), the model basically considers the 3D structural undulation and the change of fault geometry, and can realize the numerical simulation of 3D stress field, and the applicability of modeling is greatly improved. However, the 3D geomechanical model built in this phase is usually different from the 3D model geological type in the reservoir model. The former is often a stratum equal thickness model, which does not consider the information such as stratum pinch-out and thickness change, and cannot fully realize the gradual and sudden changes of rock mechanical parameters in the model; at the same time, the meshing accuracy cannot match the reservoir modeling. Phase III, 3D heterogeneous geomechanical modeling phase (2017-present), Liu et al. (2017a, b, 2018a) systematically proposed a reservoir geomechanical heterogeneous modeling method for the first time at home and abroad, and developed a corresponding software system, which achieved the matching of the accuracy of geomechanical modeling and reservoir geological modeling, and greatly improved the applicability of reservoir geomechanical modeling. The core elements of heterogeneous modeling in geomechanical modeling include heterogeneous geological models, heterogeneous mechanical protection models, heterogeneous boundary conditions, and heterogeneous fracture conditions. However, there are three problems in geomechanical modeling: (1) Paleomorphology and paleofault distribution in different periods cannot be considered in paleostress simulation; (2) Four-dimensional stress field simulation does not consider the assignment principles of rock mechanical parameters and mechanical parameters within the fault zone; (3) The geometric forms of models in different periods vary greatly. How to perform four-dimensional stress evolution in a model directly determines whether the four-dimensional stress field prediction results are practical, which is also crucial for fracture prediction in different regions. Therefore, systematic research on quantitative prediction and integrated modeling methods of four-dimensional stress fields in oil and gas reservoirs has important theoretical and practical significance for realizing fine reservoir geomechanical modeling, interpreting the tectonic evolution process of superimposed basins, improving the accuracy of reservoir fracture prediction, and optimizing reservoir transformation plans. Summary of the invention

[0003] The present invention aims to solve the above problems and provides a method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs. The method can accurately analyze the evolution process of stress field of complex oil and gas reservoirs and realize the display of multi-period stress field in one model.

[0004] The technical solution of the present invention is: a quantitative prediction and integrated modeling method of four-dimensional stress field of oil and gas reservoirs, and the specific steps are as follows: the first step is to restore the ancient structural morphology during the key structural change period and establish finite element geometric models of different periods;

[0005] Combined with the period of fault activity, crack formation time, tectonic deformation evolution process and regional stress background in the study area, the key tectonic change period of the study area was determined, and a comprehensive paleo-geomorphological restoration method combining forward modeling and inversion was used to restore the structural ups and downs and fault distribution of the target strata in different periods, clarify the activity intensity and activity nature of the faults, and establish geological geometric models for different periods. In the finite element software, the stratum ups and downs and fault geometry of different periods were established, and the intersection relationship between strata and faults, and between faults was determined. The three-dimensional distribution of the ancient tectonic morphology point cloud and the current tectonic morphology point cloud in different periods was obtained, and the finite element geometric models of different periods were established.

[0006] The comprehensive restoration method of paleo-geomorphology combining forward modeling and inversion is suitable for the restoration of paleo-geomorphology in complex tectonic areas, and mainly includes two types: forward modeling method and inversion method. The forward modeling method refers to the use of discrete element software, with known main fault displacement and residual sediment thickness as constraints, by determining the main stratum lithology combination and mechanical properties, according to the equilibrium profile restoration results and the regional tectonic stress evolution background, clarifying the discrete element numerical simulation boundary condition parameters, and through repeated iterative calculations, clarifying the displacement field, stress field, strain field, fault plane stress and fault displacement distribution map in three directions, determining the size and range of the uplift area, and then obtaining the quantitative forward prediction results of tectonic uplift and subsidence in the study area. The inversion method refers to the establishment of two-dimensional and three-dimensional structural geometric models, and sequentially carrying out decompaction, fault restoration and layer flattening operations to obtain paleo-geomorphology maps; combining the current stratum undulation, residual fault morphology and the size of the uplift obtained by forward simulation, the paleo-geomorphology of the previous period is obtained; finally, after comparing and correcting the landforms determined by the forward modeling method and the inversion method, the comprehensive and fine restoration of the paleo-geomorphology of the study area is completed.

[0007] The second step is to combine the paleo-tectonic morphology during the key tectonic change period, determine the size of the grid unit in the finite element geometric model according to the seismic accuracy, and determine the displacement direction, trajectory and size of different nodes in the finite element geometric model.

[0008] Combining the ancient structural morphology of key tectonic change periods and the finite element geometric models of different periods, as well as the ancient structural stress directions and deformation characteristics of different periods, the displacement directions, trajectories and sizes of different nodes in the finite element geometric model are determined, and the displacement vectors of units in different periods are stored, laying the foundation for the subsequent four-dimensional stress integrated modeling.

[0009] How to determine the size of the characterization unit of fractured rock mass is the primary problem faced by reservoir geomechanical modeling. In the traditional geomechanical modeling grid unit division, predecessors often determine the size of the grid unit based on oil and gas exploration and development requirements or computer hardware and software conditions. However, whether the size of the grid unit is reasonable often lacks systematic and scientific analysis. In geomechanical modeling, especially in fractured rock mass, if the grid unit size is too small, the grid unit will not be able to truly reflect the fracture combination style and development degree at that location, and thus cannot accurately reflect the size of the rock mass mechanical parameters at that location, that is, at the same location, the size of the rock mass mechanical parameters will be affected by the size of the simulation unit. On the other hand, if the grid unit is too large, it will affect the accuracy of the subsequent numerical simulation, resulting in a smaller difference between adjacent grid units, resulting in a decrease in the practicality of the simulation results. The method of determining the size of the grid unit in the finite element geometric model based on seismic accuracy described in the patent of the present invention refers to determining the size of the finite element grid division based on the size of the quarter wavelength length of the post-stack seismic data.

[0010] The third step is to determine the confining pressure range of rock mechanics experiments based on the ancient burial depth ranges of reservoirs in different periods, carry out variable confining pressure rock mechanics experiments, shear tests, and Brazilian splitting experiments, determine the change model between burial depth and static rock mechanics parameters, and establish a rock dynamic and static mechanics parameter conversion model.

[0011] Rock mechanics parameters (mainly including rock Poisson's ratio, rock strength parameters, various elastic moduli, internal friction angle and cohesion, etc.) are important basic data for ancient and modern stress field simulation, wellbore stability analysis and fracture parameter prediction. Representative cores were selected and processed into rock samples with flat end faces, 2.5 cm in diameter and 5.0 cm in length using drilling rigs, slicers and other equipment; the rock mechanics experiment was calibrated on the dynamic mechanical parameters interpreted by well logging, and a rock dynamic-static mechanical parameter conversion model was established to determine the three-dimensional distribution of static mechanical parameters of intact rocks (without cracks) in the study area.

[0012] Select rock samples with undeveloped cracks, conduct triaxial compression experiments and synchronous acoustic emission monitoring to reveal the in-situ fracture characteristics and mechanical properties of different types of rocks; combine the ancient structural morphology of different periods to determine the range of burial depth changes experienced by rocks in the geological mechanics period, conduct triaxial mechanical experiments and shear tests of rocks under different confining pressure conditions, and test the corresponding tensile strength and internal friction angle, establish a mathematical relationship model between rock Young's modulus, Poisson's ratio, shear strength, density and burial depth, and determine the type of change model between burial depth and static rock mechanical parameters. Combined with the dynamic mechanical parameters calculated from the logging data, find out the static mechanical experimental results of rocks determined by the variable confining pressure rock mechanical experiments, shear experiments, and Brazilian splitting experiments at the corresponding positions, and establish a rock dynamic and static mechanical parameter conversion model.

[0013] The types of variation models between burial depth and static rock mechanical parameters are divided into five categories: monotonic linear, two-stage, exponential, power law, and composite. By determining the variation model between burial depth and static rock mechanical parameters, the foundation is laid for the three-dimensional evolution of rock mechanical parameters in the later stage.

[0014] The fourth step is to directly invert the ratio of P-wave to S-wave velocity based on pre-stack seismic data and the one-dimensional distribution of rock mechanical parameters calculated from a single well, using well-seismic constraints to obtain the three-dimensional distribution data of current rock mechanical parameters.

[0015] Based on pre-stack seismic inversion, the inversion method of P-wave velocity ratio and Poisson's ratio is constructed to achieve the stability inversion of reservoir elastic parameters, which improves the identification effect of oil and gas layers with poor quality of pre-stack seismic data and small reservoir parameter differences. Based on the simplified Aki-Richard pre-stack rock mechanics parameter inversion method, the time domain data volume of P-wave velocity is directly inverted, and the depth domain data volume of P-wave velocity is obtained by combining the seismic interface of fine structural interpretation. The rock mechanics parameter calculation formula is used to calculate the rock mechanics parameter data volume, and the rock mechanics parameter calculation results of logging data are combined with the pre-stack inversion data volume to invert the current three-dimensional distribution data volume of rock mechanics parameters.

[0016] The rock mechanics parameter calculation formula is:

[0017]

[0018]

[0019] In formulas (1) to (2), E d is the dynamic Young's modulus of rock, GPa; μ d is the dynamic Poisson's ratio of rock, dimensionless; ρ b is the rock density interpreted by well logging, g / cm 3 ; Δt p is the longitudinal wave time difference of the rock, μs / m; Δts is the shear wave time difference of rock, μs / m.

[0020] The Aki-Richard simplified pre-stack rock mechanics parameter inversion method is:

[0021]

[0022] In formula (3), R(θ) is the approximate equation of the longitudinal wave reflection coefficient, v p is the average longitudinal wave velocity, m / s; v s is the average shear wave velocity, m / s; ρ is the average density, g / cm 3 ; Δv p is the longitudinal wave velocity difference, m / s; Δv s is the shear wave velocity difference, m / s; Δρ is the density difference, g / cm 3 ; θ is the average value of the reflection angle and the transmission angle.

[0023] The fifth step combines the current three-dimensional distribution of rock mechanical parameters and uses the mathematical model between burial depth and static rock mechanical parameters to determine the fitting coefficients of rock mechanical parameters in the vertical direction and burial depth at different locations in the study area.

[0024] Using the current three-dimensional distribution of rock mechanical parameters, the relationship model between the burial depth and rock mechanical parameters of units at different locations is established. According to the size of the fitting coefficient between the two, the type of variation model between the burial depth and the static rock mechanical parameters is selected to determine the fitting coefficient of the rock mechanical parameters at different locations in the study area in the vertical direction and the burial depth. The variation model types between the burial depth and the static rock mechanical parameters include monotonic linear, two-stage, exponential, power law, and composite.

[0025] The sixth step is to integrate ancient and modern structures, determine the variation range of burial depth of reservoirs at different locations at different times, predict the three-dimensional evolution of ancient rock mechanical parameters at different times, and obtain the three-dimensional distribution data body of ancient rock mechanical parameters of reservoirs at different times.

[0026] By using the three-dimensional distribution of current rock mechanics parameters, combined with the variation range of burial depth of reservoirs at different locations at different periods and the fitting coefficient of rock mechanics parameters in the vertical direction and burial depth at different locations in the study area, and using the variation model between burial depth and static rock mechanics parameters for correction, the three-dimensional evolution of ancient rock mechanics parameters in different periods is predicted, and the three-dimensional distribution data body of ancient rock mechanics parameters in different periods is obtained.

[0027] The fusion of ancient and modern structures refers to combining the ancient structural morphology point cloud data of different periods with the present structural morphology point cloud data, taking the present structural point data as the basis, selecting the ancient structural morphology point cloud data that is closest to it, determining the structural fluctuation changes in different periods corresponding to different present structural points, realizing the fusion of ancient and modern structures, and determining the variation range of the burial depth of reservoirs at different positions at different periods, laying the foundation for determining the mechanical parameters of ancient rocks in different periods.

[0028] The seventh step is to reveal the magnitude of rock mechanical parameters of rock masses containing fractures through rock mechanical experiments on fractured samples, find out the mechanical parameters of rock masses with filled fractures and unfilled fractures, and establish an assignment model for the mechanical parameters of underground fractures and fault rock masses.

[0029] In-situ mechanical experiments were carried out using rock samples with filled fractures, unfilled fractures, and rock samples without fractures to determine the mechanical parameters of the rock samples. Samples with different angles between the axis of the rock sample cylinder and the fracture were drilled to establish the relationship between the rock mechanical parameters of filled faults, unfilled faults, and samples without fractures. The influence of the angle between the stress loading direction and the fracture on the rock mechanical parameters and rock fracture was determined, and an assignment model for the mechanical parameters of rock masses with filled fractures and unfilled fractures was established.

[0030] The eighth step is to clarify the rock mechanical parameter assignment algorithm at different positions of the fault based on the three-dimensional distribution data of ancient rock mechanical parameters of reservoirs in different periods, combined with the assignment model of mechanical parameters of rock masses with filled fractures and unfilled fractures, and establish the three-dimensional distribution data of mechanical parameters of reservoirs containing faults in different periods.

[0031] Reasonable assignment of mechanical parameters of faults and their vicinity has always been a difficult problem. Combined with the paleo-tectonic morphology of the fault activity period and the key tectonic change period in the study area, the geometric morphology and burial depth of the faults in the study area at different periods were determined, and the unit volume data closest to the faults were selected from the three-dimensional distribution data body of paleo-rock mechanical parameters. Combined with the assignment model of rock mechanical parameters of filled and unfilled fractures, the filling properties of the faults were determined by the geological-geomechanical-geophysical integrated method, and finally the rock mechanical parameter assignment algorithm at different positions of the fault was determined, that is, the assignment algorithm of filled or unfilled faults was determined, and then the three-dimensional distribution data body of reservoir mechanical parameters containing faults at different periods was established.

[0032] The use of the comprehensive geological-geomechanical-geophysical method to determine the filling properties of faults includes three aspects: geological method, geomechanical method, and geophysical method. Among them, the geological method refers to the inference of the filling properties of adjacent faults through the filling properties of fractures in drilling and logging data. The geomechanical method refers to determining the filling properties of adjacent faults through the law of geostress steering. When the direction of the maximum horizontal principal stress turns perpendicular to the fault strike, it is inferred that the interior of the fault is mostly filled with quartz or calcite. When the direction of the maximum horizontal principal stress turns parallel to the fault strike, it is inferred that the interior of the fault is mostly unfilled or filled with mud. The geophysical method refers to inferring the filling properties of adjacent faults through the sound velocity and density characteristics near the fault determined in the geophysical data, or the filling characteristics of the fractures.

[0033] The ninth step is to determine the size of the finite element unit body according to the seismic accuracy, and use the geological stratification-coordinate adjacency method to assign the three-dimensional rock mechanics parameter distribution data body containing fault zones in different periods to the finite element geometric models of different periods to establish the geomechanical heterogeneous model of oil and gas reservoirs in different periods.

[0034] The geological stratification-coordinate adjacency method is used to couple the three-dimensional distribution data of reservoir mechanics parameters containing faults at different periods, the three-dimensional distribution data of ancient rock mechanics parameters of reservoirs at different periods, and the three-dimensional distribution data of current rock mechanics parameters to the corresponding finite element geometric models of different periods, and establish geomechanical heterogeneous models of oil and gas reservoirs at different periods.

[0035] The geological stratification-coordinate adjacency method, under the constraints of geological stratification, assigns the two types of unit bodies in the finite element geometric model and the geological model according to the principle of proximity, and assigns the geological model unit body to the finite element geometric model according to the principle of closest distance, thereby realizing geomechanical heterogeneous modeling in different periods.

[0036] The tenth step is to carry out rock acoustic emission testing, combine the structural deformation shortening and ancient and modern stress types in different periods, determine the boundary load conditions of the geomechanical heterogeneous model of the oil and gas reservoir in different periods, predict the three-dimensional distribution of the stress field in different periods, and realize the quantitative prediction of the four-dimensional stress field.

[0037] When rock is deformed by force, stress is concentrated around existing or newly generated cracks in the rock, and the strain energy is high; when the external force increases to a certain size, microscopic yielding or deformation will occur at the crack defect site, and the crack will expand, thereby relaxing the stress, and part of the stored energy will be released in the form of elastic waves, which is the acoustic emission phenomenon. It can be seen that acoustic emission is generated during the process of rock stress and deformation, and the stress evolution state of underground rock can be understood based on the number, size, frequency and other parameters of rock acoustic emission.

[0038] By using the shortening of structural deformation, the intensity of fault activity, and the ancient and modern stress types in different periods, combined with the test results of acoustic emission experiments, the direction and magnitude of ground stress in different periods are determined, and the corresponding stress magnitude is applied to different model boundaries. Displacement constraints are applied to the bottom surface of the model, and the Z-directional displacement constraint is set to zero to ensure that the model does not fall freely, and finally the mechanical boundary load conditions of the model in different periods are determined.

[0039] Using the geomechanical heterogeneous models of oil and gas reservoirs at different periods and the mechanical boundary load conditions of the models at different periods, we carried out finite element numerical simulations of ancient and modern stress fields, predicted the three-dimensional distribution of stress fields at different periods, and realized quantitative prediction of four-dimensional stress fields. The eleventh step is to establish a multi-period stress unit body cyclic coupling algorithm, couple the models of different periods to the current geostress model, and realize the integrated modeling of four-dimensional stress field data.

[0040] Combining the different node displacement directions, trajectories and sizes in the finite element geometric model, a multi-period stress unit body cyclic coupling algorithm is established, and the three-dimensional distribution data of the stress field in different periods are coupled to the current geostress geomechanical model, so that the multi-period stress field prediction data can be displayed in one model, achieving the goal of integrated modeling of four-dimensional stress field data, and laying the foundation for the later prediction of reservoir fracture parameters and analysis of wellpoint stress field evolution.

[0041] The multi-period stress unit body cyclic coupling algorithm refers to using the displacement direction and size of different periods to determine the unit body number along the displacement direction and the displacement amount of different periods, combining the center position coordinates of unit bodies of different periods and shapes, superimposing the displacement amounts of different periods along the displacement direction, and obtaining the unit body coordinates after tectonic movement, and then using the geological stratification-coordinate adjacency method to couple the size and direction of stress and strain in different periods to the current model, and coupling the models of different periods to the current ground stress model, so as to realize the integrated modeling of four-dimensional stress field data.

[0042] The beneficial effects of the present invention are as follows: by restoring the ancient structural morphology during the key structural change period, finite element geological geometry modeling of different periods is realized; through the ancient structural morphology of different periods, according to the variable confining pressure rock mechanics experiment, shear test, and Brazilian splitting experiment, the change model between the burial depth and the static rock mechanics parameters is modeled; through pre-stack data driving, the three-dimensional distribution data body of the current rock mechanics parameters is inverted; through the fusion of ancient and modern structures, the burial depth changes of rocks in different periods are determined, and according to the three-dimensional evolution law of the ancient rock mechanics parameters in different periods, the three-dimensional distribution data body of the ancient rock mechanics parameters in different periods is obtained; through the rock mechanics experiment of cracked samples, the mechanical parameter assignment principle of filled faults and unfilled faults is determined; through the ancient and modern rocks The three-dimensional evolution of mechanical parameters can realize the establishment of the fault zone mechanical parameter assignment model and the establishment of the three-dimensional rock mechanical parameter distribution data body of the fault zone; the grid unit size is determined by seismic accuracy to realize the geomechanical heterogeneous modeling of different periods; the boundary conditions of the geomechanical heterogeneous model of the oil and gas reservoir in different periods are determined by acoustic emission testing and the shortening of structural deformation in different periods, and the four-dimensional stress field numerical simulation of different periods is realized; the displacement field of the geomechanical model in different periods is established through the strain, displacement and displacement direction of the unit body in different periods; the multi-period stress field prediction data is displayed in one model based on the multi-period stress unit body cyclic coupling algorithm through the displacement direction and size of different periods, and the goal of integrated modeling of four-dimensional stress field data is achieved. The patent of this invention proposes a quantitative prediction and integrated modeling method of the four-dimensional stress field of oil and gas reservoirs, which has high practical value, low prediction cost and strong operability. This method can have practical significance in many aspects such as numerical simulation of ancient and modern stress fields, quantitative prediction of structural fracture distribution and evaluation of "geological and engineering sweet spots" of complex oil and gas reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 The figure is a flow chart of a quantitative prediction and integrated modeling method for four-dimensional stress field of oil and gas reservoirs.

[0044] Figure 2 It is the structural outline of the study area and the comprehensive columnar diagram of the buried hill strata (according to Yi Jian et al., 2022).

[0045] Figure 3 This is the structural map of the buried mountain top in the Bozhong 19-6 structural area (according to Wang Wei, 2023).

[0046] Figure 4 This is the structural evolution process of Bozhong 19-6 (according to Wang Wei, 2023).

[0047] Figure 5 The results of crack staging evaluation.

[0048] Figure 6 (A) is the elevation map of the buried mountain top during the Yanshanian period; (B) is the elevation map of the buried mountain top during the Himalayan period.

[0049] Figure 7 (A) is the finite element geometric model of the Yanshanian buried hill; (B) is the finite element geometric model of the Himalayan buried hill.

[0050] Figure 8 This is the displacement cloud map of the early Indochina study area.

[0051] Fig. 9 The quality of prestack seismic data in the study area.

[0052] Fig.10 This is the stress-strain test curve of rock sample 12-1 in the study area.

[0053] Fig.11 The compressive strength of the sample changes with the confining pressure.

[0054] Fig.12 Molar circle regression analysis for samples.

[0055] Fig.13 It is the dynamic and static conversion formula of the study area.

[0056] Fig.14 One-dimensional distribution of rock mechanical parameters calculated for a single well.

[0057] Fig.15 is the planar distribution of rock mechanics parameters.

[0058] Fig.16 This is the profile distribution of mechanical parameters of different types of rocks.

[0059] Fig.17 It is the relationship model between the burial depth of unit bodies at different locations and rock mechanical parameters.

[0060] Fig.18 It is a mathematical model between burial depth and static rock mechanical parameters.

[0061] Fig.19 Schematic diagram of the tectonic evolution of the ancient buried hill.

[0062] Fig. 20 It is a three-dimensional distribution data body of reservoir paleo-rock mechanical parameters in different periods.

[0063] Fig.21 These are the test results of rock mechanical parameters for unfilled fractures.

[0064] Fig. 22 Test results of rock mechanical parameters for filling cracks.

[0065] Fig.23 (A) is the geological geometric model of the late Indosinian fault in the study area; (B) is the finite element mesh division of the fault in the study area.

[0066] Fig.24It is a geomechanical heterogeneous model of oil and gas reservoirs at different periods.

[0067] Fig.25 This is a graph showing the relationship between the cumulative event counts of acoustic emission at 7-0° and time for the sample.

[0068] Fig.26 This is a graph showing the relationship between the cumulative event counts of acoustic emission from the sample at 7-45° and time.

[0069] Fig. 27 This is a graph showing the relationship between the cumulative event counts of acoustic emission from the sample at 7-90° and time.

[0070] Fig.28 This is the relationship between the cumulative event count and time of vertical acoustic emission of sample 7.

[0071] Fig.29 These are the stress field simulation results at different periods.

[0072] Fig.30 Integrated modeling of four-dimensional stress field data. DETAILED DESCRIPTION

[0073] The specific implementation of the present invention is described below in conjunction with the accompanying drawings:

[0074] The present invention patent takes the buried hill reservoir in the XY block of the BZ gas field in the Bozhong Sag as an example to illustrate the specific implementation process of the present invention. The Bozhong Sag is the largest hydrocarbon-rich sag in the oil field of the Bohai Bay Basin. The basement buried hill is composed of multiple strata such as the Archean, Lower Paleozoic, Upper Paleozoic and Mesozoic. The Archean buried hill has complex lithology, including mixed metamorphic rocks, mixed gneiss, gneiss and other lithologies ( Figure 2). Archean metamorphic buried hills are an important area for deep oil and gas exploration in the oil fields of the Bohai Bay Basin (Xue Yongan et al., 2018; Ye et al., 2022). In 2007 and 2012, large Archean buried hill oil fields such as Jinzhou 25-1 South and Penglai 9-1 were discovered respectively; the large Archean metamorphic buried hill condensate gas field of Bozhong 19-6 gas field discovered in 2019 marked the birth of another 100 billion cubic meters of large gas field in the Bohai Bay Basin, opening up a new situation for buried hill natural gas exploration in the oil-rich basin of the Bohai Bay; in 2021, the Bozhong 13-200 million ton Archean buried hill oil and gas field was discovered in the adjacent structure, further demonstrating the broad exploration prospects of deep Archean metamorphic buried hills in the oil fields of the Bohai Bay Basin (Zhou Xinhuai et al., 2017; Xue Yongan et al., 2018, 2021; Ye et al., 2022). Studies have shown that reservoir conditions are the key to determining whether Archean metamorphic rocks in the Bozhong Sag can form large buried-hill oil and gas reservoirs. Metamorphic rocks themselves are relatively dense and do not develop primary pores (Liu et al., 2021; Wang Qingbin et al., 2022). Studies on typical bedrock buried-hill oil and gas reservoirs at home and abroad have shown that metamorphic buried-hill reservoirs are formed by later transformations such as weathering and leaching, structural fractures, and hydrothermal dissolution, and structural fractures and dissolution pores are the main reservoir spaces (Li Xiongyan and Qin Ruibao, 2022; Wang et al., 2022). Fine characterization and prediction of the distribution of stress field evolution laws in buried hills are the basis for exploring the formation conditions and development laws of favorable metamorphic buried hill reservoirs in the oil fields of the Bohai Bay Basin (Fan Hongjun et al., 2021).

[0075] The first step is to restore the ancient structural morphology during the key structural change period and establish finite element geometric models for different periods;

[0076] Combined with the fault activity period, crack formation time, and structural deformation evolution process in the study area ( Figure 3-Figure 5 ), determine the key tectonic change period of the study area, and use a comprehensive paleo-geomorphological restoration method combining forward modeling and inversion to restore the structural fluctuations and fault distribution of the target strata in different periods ( Figure 6 ) and establish geological geometric models for different periods ( Figure 7 ); In the finite element software, the stratum fluctuations and fault geometry of different periods are established, and the intersection relationship between strata and faults is determined, the three-dimensional distribution of the point cloud of the ancient tectonic morphology and the point cloud of the current tectonic morphology of different periods is obtained, and the finite element geometric models of different periods are established.

[0077] The comprehensive restoration method of paleo-geomorphology combining forward modeling and inversion is suitable for the restoration of paleo-geomorphology in complex tectonic areas, and mainly includes two types: forward modeling method and inversion method. The forward modeling method refers to the use of discrete element software, with known main fault displacement and residual sediment thickness as constraints, by determining the main stratum lithology combination and mechanical properties, according to the equilibrium profile restoration results and the regional tectonic stress evolution background, clarifying the discrete element numerical simulation boundary condition parameters, and through repeated iterative calculations, clarifying the displacement field, stress field, strain field, fault plane stress and fault displacement distribution map in three directions, determining the size and range of the uplift area, and then obtaining the quantitative forward prediction results of tectonic uplift and subsidence in the study area. The inversion method refers to the establishment of two-dimensional and three-dimensional structural geometric models, and sequentially carrying out decompaction, fault restoration and layer flattening operations to obtain paleo-geomorphology maps; combining the current stratum undulation, residual fault morphology and the size of the uplift obtained by forward simulation, the paleo-geomorphology of the previous period is obtained; finally, after comparing and correcting the landforms determined by the forward modeling method and the inversion method, the comprehensive and fine restoration of the paleo-geomorphology of the study area is completed.

[0078] The second step is to combine the ancient structural morphology during the key structural change period, determine the size of the grid unit in the finite element geometric model according to the seismic accuracy, and determine the displacement direction, trajectory and size of different nodes in the finite element geometric model;

[0079] Combined with the ancient structural morphology of the key structural change period and the finite element geometric model of different periods, combined with the structural ancient stress direction and deformation characteristics of different periods, the displacement direction, trajectory and size of different nodes in the finite element geometric model are determined, and the displacement vectors of the units in different periods are stored ( Figure 8 ).

[0080] Determining the size of the grid unit in the finite element geometric model based on seismic accuracy means determining the size of the finite element grid division based on the size of the quarter wavelength length of the post-stack seismic data. Figure 8 As shown, the basic data conditions currently available for the target area gas field are: the seismic data are post-stack seismic data collected in 2019 and processed by PSDM in 2020, with a main frequency of 23 Hz and a quarter wavelength of approximately 30 m; therefore, the size of the grid unit in the finite element geometric model is determined to be 30 meters.

[0081] The third step is to determine the confining pressure range of rock mechanics experiments based on the ancient burial depth intervals of reservoirs in different periods, carry out variable confining pressure rock mechanics experiments, shear tests, and Brazilian splitting experiments, determine the variation model between burial depth and static rock mechanics parameters, and establish a rock dynamic and static mechanics parameter conversion model;

[0082] Combined with the ancient structural morphology of different periods, the range of burial depth changes experienced by rocks in the geomechanical period was determined, triaxial mechanical tests and shear tests of rocks under different confining pressure conditions were carried out, and the corresponding tensile strength and internal friction angle ( Figure 10-12 , Table 1), and establish a mathematical relationship model between rock Young's modulus, Poisson's ratio, shear strength, density, burial depth, and confining pressure ( Fig.11 ), determine the type of change model between burial depth and static rock mechanical parameters; combine the dynamic mechanical parameters calculated by logging data, find out the rock static mechanical experimental results determined by the variable confining pressure rock mechanics experiment, shear test, and Brazilian splitting experiment at the corresponding position, and establish the rock dynamic and static mechanical parameter conversion model ( Fig.13 ).

[0083] The types of variation models between the burial depth and the static rock mechanical parameters are divided into five categories: monotonic linear, two-stage, exponential, power law, and composite.

[0084] Table 1 Triaxial shear test results of rock samples in group 1

[0085]

[0086] Based on the multi-stage triaxial compression test, the strength parameters under confining pressure close to the original field stress are obtained, and the mathematical model of compressive strength and static Young's modulus is established. The relationship between tensile strength and static Young's modulus is determined through the Brazilian splitting test, and the rock strength calculation model of the study area is established ( Fig.11 , Fig.12 ).

[0087] The dynamic elastic parameters of rocks (Young's modulus, Poisson's ratio) are calculated using logging curves. The static elastic parameters obtained from rock mechanics experiments are selected to correspond to the dynamic elastic parameters at depth, and the two are intersected to obtain the dynamic-static conversion model of elastic parameters in the study area ( Fig.13 ).

[0088] The fourth step is to directly invert the ratio of P-wave and S-wave velocities based on pre-stack seismic data and the one-dimensional distribution of rock mechanical parameters calculated from a single well, using well-seismic constraints, and obtain the three-dimensional distribution data of current rock mechanical parameters.

[0089] The rock mechanics parameter calculation formula is:

[0090]

[0091]

[0092] In formulas (1) to (2), E d is the dynamic Young's modulus of rock, GPa; μ d is the dynamic Poisson's ratio of rock, dimensionless; ρ b is the rock density interpreted by well logging, g / cm 3 ; Δt p is the longitudinal wave time difference of the rock, μs / m; Δt sis the shear wave time difference of the rock, μs / m; Formulas (1)-(2) are used to calculate the one-dimensional distribution of rock mechanical parameters of a single well ( Fig.14 )

[0093] Based on the simplified Aki-Richard pre-stack rock mechanics parameter inversion method, the time domain data volume of P- and S-wave velocities is directly inverted, and the depth domain data volume of P- and S-wave velocities is obtained by combining the seismic interface of fine structural interpretation. The rock mechanics parameter data volume is calculated using the rock mechanics parameter calculation formula, and the rock mechanics parameter calculation results from logging data are combined with the pre-stack inversion data volume to invert the current three-dimensional distribution data volume of rock mechanics parameters.

[0094] The Aki-Richard simplified pre-stack rock mechanics parameter inversion method is:

[0095]

[0096] In formula (3), R(θ) is the approximate equation of the longitudinal wave reflection coefficient, v p is the average longitudinal wave velocity, m / s; v s is the average shear wave velocity, m / s; ρ is the average density, g / cm 3 ; Δv p is the longitudinal wave velocity difference, m / s; Δv s is the shear wave velocity difference, m / s; Δρ is the density difference, g / cm 3 ; θ is the average of the reflection angle and the transmission angle, °; Formula (3) is used to calculate the one-dimensional distribution of rock mechanical parameters of a single well ( Figure 15-16 ).

[0097] The fifth step is to combine the current three-dimensional distribution of rock mechanical parameters and use the mathematical model between burial depth and static rock mechanical parameters to determine the fitting coefficients of rock mechanical parameters in the vertical direction and burial depth at different locations in the study area;

[0098] Using the current three-dimensional distribution of rock mechanical parameters, the relationship model between the burial depth and rock mechanical parameters of different unit bodies is established. Fig.17 ), the optimal model type for the change between burial depth and static rock mechanical parameters ( Fig.18 ), determine the fitting coefficients of rock mechanical parameters at different locations in the study area in the vertical direction and the burial depth; the types of variation models between the burial depth and the static rock mechanical parameters include monotonic linear, two-stage, exponential, power law, and composite.

[0099] Based on the buried hill structural evolution process and the law of rock mechanical properties changing with depth, the vertical change rate of rock mechanical properties of buried hill reservoirs ( Fig.18 ).

[0100] The sixth step is to integrate ancient and modern structures, determine the depth variation of reservoirs at different locations at different times, predict the three-dimensional evolution of ancient rock mechanical parameters at different times, and obtain the three-dimensional distribution data of ancient rock mechanical parameters of reservoirs at different times;

[0101] Using the current three-dimensional distribution of rock mechanics parameters, combined with the depth variation of reservoirs at different locations at different times ( Fig.19 ) and the fitting coefficients of rock mechanical parameters in different locations in the study area in the vertical direction and burial depth, and the variation model between burial depth and static rock mechanical parameters is used for correction to predict the three-dimensional evolution of ancient rock mechanical parameters in different periods, and the three-dimensional distribution data of ancient rock mechanical parameters in different periods are obtained ( Fig. 20 ).

[0102] Based on the clarification of the fluctuation law of rock mechanical properties and the main controlling factors of rock mechanical properties in different structural parts, the distribution characteristics of three-dimensional rock mechanical parameters of buried hill reservoirs in different tectonic periods were clarified ( Fig. 20 ).

[0103] The fusion of ancient and modern structures refers to combining the ancient structural morphology point cloud data of different periods with the present structural morphology point cloud data, taking the present structural point data as the basis, selecting the ancient structural morphology point cloud data that is closest to it, determining the structural fluctuation changes in different periods corresponding to different present structural points, realizing the fusion of ancient and modern structures, and determining the variation range of the burial depth of reservoirs at different positions at different periods, laying the foundation for determining the mechanical parameters of ancient rocks in different periods.

[0104] The seventh step is to reveal the size of rock mechanical parameters of rock mass with fractures through rock mechanical experiments on fractured samples, find out the mechanical parameters of rock mass with filled fractures and unfilled fractures, and establish an assignment model for mechanical parameters of underground fractures and fault rock mass;

[0105] In-situ mechanical experiments were carried out using rock samples with filled fractures, unfilled fractures, and rock samples without fractures to determine the mechanical parameters of the rock samples ( Fig.21 ); drill samples with different angles between the axis of the rock sample cylinder and the crack, establish the relationship between the rock mechanical parameters of the filled fault, unfilled fault and non-cracked samples, and determine the influence of the angle between the stress loading direction and the crack on the rock mechanical parameters and rock fracture ( Fig. 22 ), establish the assignment model of mechanical parameters of rock mass with filled fractures and unfilled fractures;

[0106] Unfilled cracks affect the compressive strength and elastic parameters of rocks. When unfilled cracks develop, the rock strength decreases by about 50MPa. Filled cracks significantly affect the compressive strength of rocks. When the force crack angle (stress and crack angle) is 60°, the rock strength is the lowest and cracks are most likely to develop; when the force crack angle is 0°, that is, when the stress direction is parallel to the crack direction, it is most difficult for filled cracks to produce new fractures.

[0107] The eighth step is to clarify the rock mechanical parameter assignment algorithm at different positions of the fault based on the three-dimensional distribution data of the mechanical parameters of the ancient rocks in the reservoirs at different periods, combined with the mechanical parameter assignment model of the rock mass with filled fractures and unfilled fractures, and establish the three-dimensional distribution data of the mechanical parameters of the reservoirs with faults at different periods;

[0108] Combined with the paleo-tectonic morphology of the fault activity period and key tectonic change period in the study area, the geometric morphology and burial depth of the fault in different periods in the study area were determined ( Fig.23 ), select the unit volume data closest to the fault from the three-dimensional distribution data of ancient rock mechanical parameters, and combine the assignment model of rock mechanical parameters of filled fractures and unfilled fractures, and use the geological-geomechanical-geophysical integrated method to determine the filling properties of the fault, and finally determine the rock mechanical parameter assignment algorithm at different positions of the fault, that is, determine the assignment algorithm of filled or unfilled faults, and then establish the three-dimensional distribution data of reservoir mechanical parameters containing faults in different periods, and realize fault grid division ( Fig.23 ).

[0109] The above-mentioned comprehensive geological-geomechanical-geophysical method for determining the filling properties of faults includes three aspects: geological method, geomechanical method, and geophysical method. Among them, the geological method refers to inferring the filling properties of adjacent faults through the filling properties of cracks in drilling and logging data; the geomechanical method refers to determining the filling properties of adjacent faults through the law of ground stress turning. When the direction of the maximum horizontal principal stress turns perpendicular to the fault strike, it is inferred that the interior of the fault is mostly filled with quartz or calcite; when the direction of the maximum horizontal principal stress turns parallel to the fault strike, it is inferred that the interior of the fault is mostly unfilled or filled with mud; the geophysical method refers to inferring the filling properties of adjacent faults through the sound velocity and density characteristics near the fault determined in geophysical data, or the filling characteristics of the cracks.

[0110] The ninth step is to determine the size of the finite element unit according to the seismic accuracy, and use the geological stratification-coordinate adjacent method to assign the three-dimensional rock mechanics parameter distribution data of the fault zone in different periods to the finite element geometric model of different periods to establish the geomechanical heterogeneous model of the oil and gas reservoir in different periods;

[0111] The geological stratification-coordinate adjacency method is used to couple the three-dimensional distribution data of reservoir mechanics parameters containing faults at different periods, the three-dimensional distribution data of ancient rock mechanics parameters at different periods, and the three-dimensional distribution data of current rock mechanics parameters to the corresponding finite element geometric models of different periods, and establish the geomechanical heterogeneous model of oil and gas reservoirs at different periods ( Fig.24 ).

[0112] The geological stratification-coordinate adjacency method, under the constraints of geological stratification, assigns the two types of unit bodies in the finite element geometric model and the geological model according to the principle of proximity, and assigns the geological model unit body to the finite element geometric model according to the principle of closest distance, thereby realizing geomechanical heterogeneous modeling in different periods.

[0113] The tenth step is to carry out rock acoustic emission testing, combine the structural deformation shortening amount and ancient and modern stress types in different periods, determine the boundary load conditions of the geomechanical heterogeneous model of the oil and gas reservoir in different periods, predict the three-dimensional distribution of the stress field in different periods, and realize the quantitative prediction of the four-dimensional stress field;

[0114] Using the shortening of structural deformation in different periods, the intensity of fault activity and the ancient and modern stress types, combined with the test results of acoustic emission experiments ( Figure 25-28 ), determine the direction and magnitude of the ground stress at different periods (Table 2-Table 5), and apply the corresponding stress magnitude at different model boundaries; apply displacement constraints on the bottom surface of the model, set the Z-directional displacement constraint to zero to ensure that the model does not fall freely, and finally determine the mechanical boundary load conditions of the model at different periods.

[0115] Using geomechanical heterogeneous models of oil and gas reservoirs at different times and mechanical boundary load conditions of models at different times, finite element numerical simulations of ancient and modern stress fields are carried out to predict the three-dimensional distribution of stress fields at different times and realize quantitative prediction of four-dimensional stress fields ( Fig.29 ).

[0116] Table 2 Maximum and minimum Kaiser point stress values

[0117] Num Sample group number Confining pressure / MPa Maximum stress value (MPa) Minimum stress value (MPa) Intermediate stress value (MPa) 1 7 65 110.6 88.6 107.1 2 7-S1 70 110.2 72.9 90.4 3 14 56 105.0 69.7 98.6

[0118] Table 3 Vertical Kaiser point stress values

[0119] Num Sample group number Confining pressure MPa Kaiser point stress value of vertical sample (MPa) 1 7 65 96.1 2 7-S1 70 105.0 3 14 56 98.9

[0120] Table 4 Maximum and minimum horizontal ground stresses and ground stress gradients

[0121]

[0122] Table 5 Vertical geostress and geostress gradient

[0123] Num Sample group number Well depth(m) Vertical geostress (MPa) Vertical geostress gradient (MPa / 100) 1 7 4532.8 96.1 2.12 2 7-S1 4679.2 105.0 2.22 3 14 4373.0 98.9 2.26

[0124] The eleventh step is to establish a multi-period stress unit cycle coupling algorithm, couple the models of different periods to the current ground stress model, and realize the integrated modeling of four-dimensional stress field data;

[0125] Combined with the different node displacement directions, trajectories and sizes in the finite element geometric model, a multi-period stress unit body cyclic coupling algorithm is established to couple the three-dimensional distribution data of stress fields in different periods to the current geomechanical model of ground stress, so that the multi-period stress field prediction data can be displayed in one model, achieving the goal of integrated modeling of four-dimensional stress field data ( Fig.30 ), laying a foundation for the later prediction of reservoir fracture parameters and analysis of well point stress field evolution.

[0126] The multi-period stress unit body cyclic coupling algorithm refers to using the displacement direction and size of different periods to determine the unit body number along the displacement direction and the displacement amount of different periods, combining the center position coordinates of unit bodies of different periods and shapes, superimposing the displacement amounts of different periods along the displacement direction, and obtaining the unit body coordinates after tectonic movement, and then using the geological stratification-coordinate adjacency method to couple the size and direction of stress and strain in different periods to the current model, and coupling the models of different periods to the current ground stress model, so as to realize the integrated modeling of four-dimensional stress field data.

[0127] The present invention is described above by way of examples, but the present invention is not limited to the above specific embodiments, and any changes or modifications made based on the present invention belong to the scope of protection claimed by the present invention.

Claims

1. A quantitative prediction and integrated modeling method for four-dimensional stress field of oil and gas reservoirs, the implementation steps are as follows: The first step is to restore the ancient structural morphology during the key structural change period and establish finite element geometric models for different periods; Combined with the fault activity period, crack formation time, and structural deformation evolution process in the study area, the key structural change period of the study area was determined, and a comprehensive paleo-geomorphological restoration method combining forward modeling and inversion was used to restore the structural fluctuations and fault distribution of the target strata in different periods, and establish geological geometric models of different periods; in the finite element software, the stratum fluctuations and fault geometric forms of different periods were established, and the intersection relationship between strata and faults was determined, and the three-dimensional distribution of the ancient structural form point cloud and the current structural form point cloud of different periods was obtained, and the finite element geometric models of different periods were established; The second step is to combine the ancient structural morphology during the key structural change period, determine the size of the grid unit in the finite element geometric model according to the seismic accuracy, and determine the displacement direction, trajectory and size of different nodes in the finite element geometric model; Combined with the ancient structural morphology of the key structural change period and the finite element geometric model of different periods, combined with the structural ancient stress direction and deformation characteristics of different periods, the displacement direction, trajectory and size of different nodes in the finite element geometric model are determined, and the displacement vectors of the units in different periods are stored; The third step is to determine the confining pressure range of rock mechanics experiments based on the ancient burial depth intervals of reservoirs in different periods, carry out variable confining pressure rock mechanics experiments, shear tests, and Brazilian splitting experiments, determine the variation model between burial depth and static rock mechanics parameters, and establish a rock dynamic and static mechanics parameter conversion model; Combined with the ancient structural morphology of different periods, the range of burial depth changes experienced by rocks in the geological mechanical period is determined, triaxial mechanical experiments and shear tests of rocks under different confining pressure conditions are carried out, and the corresponding tensile strength and internal friction angle are tested. A mathematical relationship model between the rock Young's modulus, Poisson's ratio, shear strength, density and burial depth is established to determine the type of change model between burial depth and static rock mechanical parameters; combined with the dynamic mechanical parameters calculated from the logging data, the static mechanical test results of rocks determined by the variable confining pressure rock mechanical experiments, shear tests, and Brazilian splitting tests at the corresponding locations are found, and a rock dynamic and static mechanical parameter conversion model is established; The fourth step is to directly invert the ratio of P-wave and S-wave velocities based on pre-stack seismic data and the one-dimensional distribution of rock mechanical parameters calculated from a single well, using well-seismic constraints, and obtain the three-dimensional distribution data of current rock mechanical parameters. Based on the simplified pre-stack rock mechanics parameter inversion method of Aki-Richard, the time domain data volume of P- and S-wave velocities is directly inverted, and the depth domain data volume of P- and S-wave velocities is obtained by combining the seismic interface of fine structural interpretation. The rock mechanics parameter data volume is calculated using the rock mechanics parameter calculation formula, and the rock mechanics parameter calculation results of logging data are combined with the pre-stack inversion data volume to invert the current three-dimensional distribution data volume of rock mechanics parameters. The fifth step is to combine the current three-dimensional distribution of rock mechanical parameters and use the mathematical model between burial depth and static rock mechanical parameters to determine the fitting coefficients of rock mechanical parameters in the vertical direction and burial depth at different locations in the study area; Using the current three-dimensional distribution of rock mechanical parameters, the relationship model between the burial depth and rock mechanical parameters of the unit bodies at different locations is established respectively. According to the size of the fitting coefficient between the two, the type of variation model between the burial depth and the static rock mechanical parameters is selected, and the fitting coefficient of the rock mechanical parameters at different locations in the study area in the vertical direction and the burial depth is determined; the variation model types between the burial depth and the static rock mechanical parameters include five types: monotonic linear, two-stage, exponential, power law, and composite; The sixth step is to integrate ancient and modern structures, determine the depth variation of reservoirs at different locations at different times, predict the three-dimensional evolution of ancient rock mechanical parameters at different times, and obtain the three-dimensional distribution data of ancient rock mechanical parameters of reservoirs at different times; Using the three-dimensional distribution of current rock mechanical parameters, combined with the variation of burial depth of reservoirs at different locations at different times and the fitting coefficient of rock mechanical parameters in the vertical direction and burial depth at different locations in the study area, and using the variation model between burial depth and static rock mechanical parameters for correction, the three-dimensional evolution of ancient rock mechanical parameters at different times is predicted, and the three-dimensional distribution data body of ancient rock mechanical parameters at different times is obtained; The seventh step is to reveal the size of rock mechanical parameters of rock mass with fractures through rock mechanical experiments on fractured samples, find out the mechanical parameters of rock mass with filled fractures and unfilled fractures, and establish an assignment model for mechanical parameters of underground fractures and fault rock mass; In-situ mechanical experiments were carried out using rock samples with filled fractures, unfilled fractures, and rock samples without fractures to determine the size of the mechanical parameters of the rock samples; samples with different angles between the axis of the rock sample cylinder and the fracture were drilled respectively to establish the relationship between the rock mechanical parameters of the filled fault, unfilled fault, and samples without fractures, and to determine the influence of the angle between the stress loading direction and the fracture on the rock mechanical parameters and rock fracture, and to establish the assignment model of the mechanical parameters of the rock mass with filled fractures and unfilled fractures; The eighth step is to clarify the rock mechanical parameter assignment algorithm at different positions of the fault based on the three-dimensional distribution data of the mechanical parameters of the ancient rocks in the reservoirs at different periods, combined with the mechanical parameter assignment model of the rock mass with filled fractures and unfilled fractures, and establish the three-dimensional distribution data of the mechanical parameters of the reservoirs with faults at different periods; Combined with the paleo-tectonic morphology of the fault activity period and the key tectonic change period in the study area, the geometric morphology and burial depth of the faults in the study area at different periods are determined, and the unit volume data closest to the fault is selected from the three-dimensional distribution data body of paleo-rock mechanical parameters. Combined with the assignment model of rock mechanical parameters of filled fractures and unfilled fractures, the filling properties of the fault are determined by the geological-geomechanical-geophysical integrated method, and finally the rock mechanical parameter assignment algorithm at different positions of the fault is determined, that is, the assignment algorithm of filled or unfilled faults is determined, and then the three-dimensional distribution data body of reservoir mechanical parameters containing faults at different periods is established; The ninth step is to determine the size of the finite element unit according to the seismic accuracy, and use the geological stratification-coordinate adjacent method to assign the three-dimensional rock mechanics parameter distribution data of the fault zone in different periods to the finite element geometric model of different periods to establish the geomechanical heterogeneous model of the oil and gas reservoir in different periods; The geological stratification-coordinate adjacency method is used to couple the three-dimensional distribution data of reservoir mechanics parameters containing faults at different periods, the three-dimensional distribution data of ancient rock mechanics parameters at different periods, and the three-dimensional distribution data of current rock mechanics parameters to the corresponding finite element geometric models at different periods, and establish geomechanical heterogeneous models of oil and gas reservoirs at different periods; The tenth step is to carry out rock acoustic emission testing, combine the structural deformation shortening amount and ancient and modern stress types in different periods, determine the boundary load conditions of the geomechanical heterogeneous model of the oil and gas reservoir in different periods, predict the three-dimensional distribution of the stress field in different periods, and realize the quantitative prediction of the four-dimensional stress field; By using the shortening of structural deformation, the intensity of fault activity and the ancient and modern stress types in different periods, combined with the test results of acoustic emission experiments, the direction and magnitude of ground stress in different periods are determined, and the corresponding stress magnitude is applied to different model boundaries; displacement constraints are applied to the bottom surface of the model, and the Z-directional displacement constraint is set to zero to ensure that the model does not fall freely, and finally the mechanical boundary load conditions of the model in different periods are determined; Using geomechanical heterogeneous models of oil and gas reservoirs at different times and mechanical boundary load conditions of models at different times, finite element numerical simulations of ancient and modern stress fields are carried out to predict the three-dimensional distribution of stress fields at different times and realize quantitative prediction of four-dimensional stress fields. The eleventh step is to establish a multi-period stress unit cycle coupling algorithm, couple the models of different periods to the current ground stress model, and realize the integrated modeling of four-dimensional stress field data; Combining the different node displacement directions, trajectories and sizes in the finite element geometric model, a multi-period stress unit body cyclic coupling algorithm is established, and the three-dimensional distribution data of the stress field in different periods are coupled to the current geostress geomechanical model, so that the multi-period stress field prediction data can be displayed in one model, achieving the goal of integrated modeling of four-dimensional stress field data, and laying the foundation for the later prediction of reservoir fracture parameters and analysis of wellpoint stress field evolution.

2. A method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The comprehensive paleo-geomorphology restoration method combining forward modeling and inversion described in the first step is applicable to the restoration of paleo-geomorphology in complex tectonic areas, and mainly includes two types: forward modeling and inversion methods. The forward modeling method refers to the use of discrete element software, with known main fault displacement and residual sediment thickness as constraints, by determining the main stratum lithology combination and mechanical properties, according to the equilibrium profile restoration results and the regional tectonic stress evolution background, clarifying the discrete element numerical simulation boundary condition parameters, and through repeated iterative calculations, clarifying the displacement field, stress field, strain field, fault surface stress field in three directions. The distribution map of fault displacement is used to determine the size and range of the uplift area, and then the quantitative forward prediction results of the tectonic uplift and subsidence in the study area are obtained; the inversion method refers to the establishment of two-dimensional and three-dimensional structural geometric models, and the decompaction, fault recovery and layer flattening operations are carried out in sequence to obtain the paleo-geomorphological map; the paleo-geomorphology of the previous period is obtained by combining the current stratigraphic undulations, residual fault morphology and the uplift amount obtained by forward simulation; finally, after comparing and correcting the landforms determined by the forward method and the inversion method, the comprehensive and detailed restoration of the paleo-geomorphology of the study area is completed.

3. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: Determining the size of the grid unit in the finite element geometric model based on the seismic accuracy described in the second step means determining the size of the finite element grid division based on the size of the quarter wavelength length of the post-stack seismic data.

4. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The types of models for the change between burial depth and static rock mechanical parameters described in the third step are divided into five categories: monotonic linear, two-stage, exponential, power-law, and composite.

5. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The rock mechanics parameter calculation formula described in the fourth step is: In formulas (1) to (2), E d is the dynamic Young's modulus of rock, GPa; μ d is the dynamic Poisson's ratio of rock, dimensionless; ρ b is the rock density interpreted by well logging, g / cm 3 ; Δt p is the longitudinal wave time difference of the rock, μs / m; Δt s is the shear wave time difference of rock, μs / m.

6. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The Aki-Richard simplified pre-stack rock mechanics parameter inversion method described in the fourth step refers to: In formula (3), R(θ) is the approximate equation of the longitudinal wave reflection coefficient, v p is the average longitudinal wave velocity, m / s; v s is the average shear wave velocity, m / s; ρ is the average density, g / cm 3 ; Δv p is the longitudinal wave velocity difference, m / s; Δv s is the shear wave velocity difference, m / s; Δρ is the density difference, g / cm 3 ; θ is the average value of the reflection angle and the transmission angle.

7. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The fusion of ancient and modern structures described in the sixth step refers to combining the ancient structural morphology point cloud data of different periods with the present structural morphology point cloud data, taking the present structural point data as the basis, selecting the ancient structural morphology point cloud data that is closest to it, determining the structural fluctuation changes in different periods corresponding to different present structural points, realizing the fusion of ancient and modern structures, and determining the variation range of the burial depth of reservoirs at different locations at different periods, laying the foundation for determining the mechanical parameters of ancient rocks in different periods.

8. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The eighth step describes the use of a comprehensive geological-geomechanical-geophysical method to determine the filling properties of faults, which includes three aspects: geological method, geomechanical method, and geophysical method. The geological method refers to inferring the filling properties of adjacent faults through the filling properties of fractures in drilling and logging data; the geomechanical method refers to determining the filling properties of adjacent faults through the law of in-situ stress steering. When the direction of the maximum horizontal principal stress turns perpendicular to the fault strike, it is inferred that the interior of the fault is mostly filled with quartz or calcite; when the direction of the maximum horizontal principal stress turns parallel to the fault strike, it is inferred that the interior of the fault is mostly unfilled or filled with mud; the geophysical method refers to inferring the filling properties of adjacent faults through the sound velocity and density characteristics near the fault determined in geophysical data, or the filling characteristics of fractures.

9. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The geological stratification-coordinate adjacency method described in the ninth step, under the constraints of geological stratification, assigns the two types of units in the finite element geometric model and the geological model according to the principle of proximity, and assigns the geological model unit body to the finite element geometric model according to the principle of closest distance, thereby realizing geomechanical heterogeneous modeling in different periods.

10. The method for quantitative prediction and integrated modeling of four-dimensional stress field of oil and gas reservoirs according to claim 1, characterized in that: The multi-period stress unit cyclic coupling algorithm described in the tenth step refers to using the displacement direction and size of different periods to determine the unit number along the displacement direction and the displacement amount of different periods, combining the center position coordinates of units of different periods and shapes, and superimposing the displacement amounts of different periods along the displacement direction to obtain the unit coordinates after tectonic movement. Then, the geological stratification-coordinate adjacency method is used to couple the size and direction of stress and strain in different periods to the current model, and the models of different periods are coupled to the current ground stress model to realize the integrated modeling of four-dimensional stress field data.

Images