Method for depicting breaking joint body in breaking joint body development area

Abstract

The invention relates to the technical field of shale gas geology, in particular to a method for depicting a breaking joint body in a breaking joint body development area, and the method comprises the steps: firstly, generating an initial breaking joint body probability body based on interpretable machine learning fusion seismic attributes, and building a three-dimensional space distribution model containing large and medium faults after manual interactive interpretation and well seismic calibration; on the basis, large and medium-sized faults serve as a deterministic skeleton, a fracture development rule is predicted in combination with paleotectonic stress field simulation and a rock fracture criterion, a three-dimensional discrete fracture network model constrained by geomechanics is constructed and then coupled with a matrix model, and a model set is generated by defining an uncertainty parameter space. And screening out an equivalent model group through historical matching, performing prediction verification and iterative optimization by using recent production data, and finally obtaining a three-dimensional space distribution model of the high-precision breaking joint body. According to the method, closed loop of static identification and dynamic verification is realized, and the depicting precision and the geological reliability are remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of shale gas geology, and more particularly to a method for characterizing fractures in fracture development zones. Background Technology

[0002] Fractured zones, as important locations for hydrocarbon accumulation, possess a complex network system of faults and associated fractures, which are key geological elements controlling hydrocarbon seepage and accumulation. As the focus of oil and gas exploration and development gradually shifts to deep, unconventional, and other complex reservoirs, a precise characterization of the spatial distribution and internal connectivity of fractured zones has become a core prerequisite for achieving efficient and stable oil and gas reservoir development. Existing technologies typically rely on comprehensive identification and modeling methods such as seismic attribute analysis, imaging logging interpretation, and geostatistical simulation.

[0003] However, existing technologies have significant limitations in addressing the isochronous comparison and dynamic verification of the spatial distribution of fracture bodies. Traditional methods largely rely on matching seismic data with static well data. While the models established can maintain accuracy near the well points, they exhibit extreme ambiguity in inter-well regions, especially when dealing with highly heterogeneous fracture body systems. Furthermore, the lack of effective constraints and verification from dynamic production data results in insufficient predictive ability for fluid flow behavior, making it difficult to accurately guide key measures such as drainage and water shut-off in the later stages of development. Summary of the Invention

[0004] To overcome the above deficiencies, this invention provides a method for characterizing fracture bodies in the fracture body development zone, aiming to improve the significant limitations of existing technologies in solving the problems of isochronous comparison and dynamic verification of fracture body spatial distribution.

[0005] In a first aspect, the present invention provides the following technical solution: a method for characterizing a fracture body in a fracture development zone, comprising the following steps: Acquire 3D seismic data and wellbore imaging logging data of the target work area, process the 3D seismic data, extract multiple seismic attribute volumes, and use interpretable machine learning algorithms to fuse the multiple seismic attribute volumes to generate an initial fracture probability volume that can be interactively interpreted by human experts. Human expert interaction is introduced to interpret and correct the initial fracture probability volume, and the wellbore imaging logging data is combined for calibration to output a calibrated three-dimensional spatial distribution model of the fracture volume, in which large and medium-sized faults have been identified. Based on the calibrated three-dimensional spatial distribution model of the fracture body, the large and medium-sized faults identified therein are used as a deterministic framework. At the same time, the paleotectonic stress field is restored based on the regional tectonic evolution history, and the rock fracture criterion is used to perform geomechanical simulation to predict the development zone, orientation and density of fractures. Using the results of the geomechanical simulation as constraints, non-deterministic fractures are generated outside the deterministic framework to jointly construct a three-dimensional discrete fracture network model. Obtain the reservoir matrix attribute model pre-established for the target work area, couple the constructed three-dimensional discrete fracture network model with the reservoir matrix attribute model to construct an integrated three-dimensional geological model of fracture body-matrix, and define the uncertainty parameters and their value ranges in the integrated three-dimensional geological model. An experimental design method is used to generate an integrated model set consisting of multiple parameters. Fluid numerical simulation is performed on each model in the integrated model set, and the simulated production dynamic data is automatically compared with the actual historical production data of the oilfield. A set of equivalent models whose historical matching error meets the preset conditions is selected. Production forecasting is performed using the selected set of equivalent models, and the forecast results are compared and verified with recent production data that has not been matched with historical data, until a final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements is obtained.

[0006] Preferably, the extraction process for the multiple seismic attribute volumes includes: The acquired three-dimensional seismic data is subjected to amplitude compensation and frequency broadening processing to enhance the effective signal and form a preprocessed three-dimensional seismic data volume. The preprocessed 3D seismic data volume is subjected to construction-guided filtering to suppress noise and preserve the boundary features of the fracture body, thereby forming a 3D seismic data volume with a high signal-to-noise ratio. Based on the high signal-to-noise ratio three-dimensional seismic data volume, multiple seismic attribute volumes for identifying fracture bodies are extracted. These multiple seismic attribute volumes include at least: a coherent attribute volume for characterizing stratigraphic discontinuities, a curvature attribute volume for characterizing structural curvature changes, and an ant-tracking attribute volume for enhancing fracture profiles.

[0007] Preferably, the process for generating the initial fracture probability body includes: Based on the aforementioned multiple seismic attribute bodies, at the drilled well points, according to the developed and undeveloped sections of the fracture body identified by the wellbore imaging logging data, corresponding sample data are extracted from the multiple seismic attribute bodies to form a labeled training sample set. The training sample set is input into the interpretable machine learning algorithm for training, thereby generating a fracture body recognition model. Subsequently, the fracture body recognition model is used to perform global scanning and fusion calculation of the various seismic attribute bodies of the target work area, and output an initial fracture body probability body. While generating the initial fault probability volume, the interpretable machine learning algorithm simultaneously calculates and outputs the feature importance ranking of the various seismic attribute volumes for the prediction results. The feature importance ranking is used to assist human experts in understanding the geological origin of the initial fault probability volume and to interactively interpret it.

[0008] Preferably, the output process of the calibrated three-dimensional spatial distribution model of the fracture body includes: Based on the initial fracture probability volume and its corresponding feature importance ranking, artificial experts identify and track fracture structures in three-dimensional space, and eliminate false anomalies in the initial fracture probability volume according to geological knowledge, enhance and supplement ambiguous or missing fracture structures, and generate a preliminary interpretation result consisting of a series of three-dimensional spatial point, line and surface data. The preliminary interpretation results are precisely compared with the wellbore imaging logging data at the well trajectory position. When the preliminary interpretation results match the wellbore data, they are confirmed. When there is a contradiction, the preliminary interpretation results are corrected based on the wellbore imaging logging data to generate a fracture body interpretation result after well control correction. After completing the well-seismic calibration of the entire area, the interpretation results of the fracture body after well control correction are integrated in three-dimensional space and converted into a complete three-dimensional data volume model. In this model, large and medium-sized faults with a scale reaching the preset threshold value are clearly identified, and finally the calibrated three-dimensional spatial distribution model of the fracture body is output.

[0009] Preferably, the prediction process for the development zone, orientation, and density of the cracks includes: Based on the calibrated three-dimensional spatial distribution model of the fracture body, the spatial geometry and attitude of the identified large and medium-sized faults are extracted and used as an immovable and unmodifiable deterministic skeleton to be directly constructed in the initial three-dimensional discrete fracture network model. Based on the regional tectonic evolution history, key tectonic movement periods are determined, and paleotectonic restoration technology is used to reconstruct the paleotectonic framework and paleoburial depth of the target strata during the key hydrocarbon accumulation period. Then, the paleotectonic stress field of that period is inverted through rock mechanics parameters to obtain the magnitude and direction of paleostress. The deterministic framework, as a discontinuous interface in the paleotectonic framework, participates in the stress field simulation. Using the restored paleotectonic stress field as input, and combined with the rock mechanics parameters of the target strata, the rock fracture potential value of each grid point is calculated in a three-dimensional spatial grid based on the Coulomb fracture criterion or the Griffith criterion. Based on the calculated rock fracture potential value, the development zone, orientation, and density of fractures are quantitatively predicted. Specifically, the high value area of ​​the rock fracture potential value is predicted as the dominant fracture development zone, the direction of the maximum principal stress of paleostress at each grid point is predicted as the dominant orientation of the fracture, and the spatial distribution of the rock fracture potential value is mapped to a normalization function to quantitatively calculate the relative density distribution of fractures.

[0010] Preferably, the construction process of the three-dimensional discrete crack network model includes: The dominant crack development zone, the dominant crack orientation, and the relative crack density distribution are respectively converted into a crack density field and a crack strike-dip field in three-dimensional space as constraint parameter fields. Based on the crack density field and the crack strike-dip field, a large number of discrete nondeterministic cracks are generated in three-dimensional space using a stochastic simulation method. The spatial location, development scale, orientation and density of the nondeterministic cracks are strictly controlled by the constraint parameter field. The deterministic skeleton and the non-deterministic cracks are fused and their geometric topological relationships are checked in three-dimensional space to ensure that the spatial contact relationships between cracks of different scales are reasonable, and finally a complete three-dimensional discrete crack network model is constructed together.

[0011] Preferably, the construction process of the integrated three-dimensional geological model of the fracture body and matrix includes: Obtain a pre-established reservoir matrix attribute model for the target work area, wherein the reservoir matrix attribute model includes a three-dimensional attribute field of matrix porosity and permeability; The three-dimensional discrete fracture network model is spatially superimposed and coupled with the reservoir matrix property model; Based on the geometric parameters and conductivity of each crack in the three-dimensional discrete crack network model, the permeability properties of the matrix mesh through which the crack passes are calculated, thereby transforming the background permeability field of the matrix into an equivalent permeability field that can characterize the interaction between the fracture body and the matrix. The transformed equivalent permeability field is integrated with the matrix porosity field to form a dual-media model that simultaneously characterizes the high-speed seepage channels of the fracture body and the matrix storage space, namely the fracture-matrix integrated three-dimensional geological model.

[0012] Preferably, the definition process for the uncertainty parameters and their value ranges in the integrated three-dimensional geological model includes: Based on the construction process of the three-dimensional discrete crack network model and the calculation principle of the equivalent permeability field, key uncertainty parameters that significantly affect the fluid simulation results of the model are identified. Based on the core test data of the target work area, the fracture aperture statistics interpreted by imaging logging, and the dynamic permeability retrieved from the well test or production history data, a predetermined value range is determined for each of the identified key uncertainty parameters. Each of the identified key uncertainty parameters is combined with its corresponding value range to form a multi-dimensional parameter space, which is used to characterize all geological uncertainties of the integrated three-dimensional geological model before historical matching.

[0013] Preferably, the screening process for the set of equivalent models includes: Using Latin hypercube sampling or factorial analysis experimental design methods, a large number of samples are taken in the multi-dimensional parameter space. Each set of specific parameters obtained from each sampling constitutes a specific implementation of the integrated three-dimensional geological model. Through multiple sampling, an integrated model set containing hundreds to thousands of different parameter implementations is generated. Each model implementation in the integrated model set is sequentially loaded into a pre-made fluid numerical simulator. Under the same well location, production system and simulation duration settings, batch numerical simulation calculations are performed to obtain the simulated production dynamic data corresponding to each model implementation. The simulated production dynamic data of each model is automatically compared with the actual historical production data of the oilfield, and the historical matching error is calculated. A historical matching error threshold is set as the preset condition. All models with historical matching errors lower than the threshold are implemented and selected as the group of equivalent models.

[0014] Preferably, the process for obtaining the final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements includes: Using the selected set of equivalent models, production forecasts for future periods are performed in the fluid numerical simulator. For each equivalent model, predicted production dynamics data are obtained. The forecast results from all the equivalent models are then combined to form a comprehensive forecast curve covering a range of uncertainties. Obtain the actual production data newly generated after the end of the historical matching period, and use it as the recent production data that did not participate in the historical matching. Compare the comprehensive prediction curve band with the recent production data. If the recent production data falls within the range of the comprehensive prediction curve band, the model is deemed to have passed the verification. When the verification is passed, the model with the smallest historical matching error is selected from the set of equivalent models. The integrated three-dimensional geological model corresponding to it and the three-dimensional discrete fracture network model that forms the basis of its construction are back-determined as the most reliable model in terms of geological understanding. The model is then output as the final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements. If the verification fails, it indicates that there is a systematic bias in the geological understanding of the set of equivalent models. At this time, based on the difference between the prediction results and the recent production data, the cause of the geological understanding bias is analyzed, and the process returns to the step of introducing human expert interaction to interpret and correct the initial fracture probability volume, or to the step of generating non-deterministic fractures outside the deterministic framework using the results of the geomechanical simulation as constraints. After correcting the geological parameters and modeling scheme, the subsequent steps are re-executed until the verification passes.

[0015] The present invention has the following beneficial effects: 1. In this invention, the progressive process of interpretable machine learning for initial identification, human expert interaction for correction, and well seismic data for precise calibration effectively overcomes the black box limitations and multiple solutions of a single intelligent algorithm. It deeply integrates the efficiency of artificial intelligence with the cognitive experience of geological experts, ensuring that the final output three-dimensional spatial distribution model of the fracture body conforms to both data-driven laws and geological laws.

[0016] 2. In this invention, by introducing paleotectonic stress field simulation and rock fracture criteria, the traditional discrete fracture network modeling method, which is mainly based on statistics, is transformed into a deterministic-nondeterministic hybrid modeling driven by geomechanical mechanisms. The fracture system generated by this method not only has a reasonable geometric shape, but also has a clear geological origin, thus enabling more accurate prediction of fracture development in unexplored areas and reducing exploration risks.

[0017] 3. In this invention, by selecting multiple equivalent models that can match historical production data and using them for production prediction and independent verification, the model group is eliminated. This process not only quantifies the uncertainty of geological understanding, but also forms an iterative optimization closed loop of dynamic data-driven correction of geological understanding, ensuring that the final model not only has a good historical fit, but also has a good future prediction capability, providing an extremely reliable geological basis for development adjustment decisions. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the method for depicting the fracture development zone of the fracture body proposed in this invention. Detailed Implementation

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

[0020] In a first embodiment of the present invention, the present invention provides a method for characterizing the fracture body in the fracture development region, such as... Figure 1 As shown, it includes the following steps: Acquire 3D seismic data and wellbore imaging logging data of the target work area, process the 3D seismic data, extract multiple seismic attribute volumes, and use interpretable machine learning algorithms to fuse multiple seismic attribute volumes to generate an initial fracture probability volume that can be interactively interpreted by human experts.

[0021] Furthermore, the extraction process for multiple seismic attribute volumes includes: The acquired 3D seismic data is subjected to amplitude compensation and frequency broadening processing to enhance the effective signal and form a preprocessed 3D seismic data volume. Construction-guided filtering is applied to the preprocessed 3D seismic data volume to suppress noise and preserve the boundary features of the fracture volume, forming a 3D seismic data volume with a high signal-to-noise ratio. Based on high signal-to-noise ratio 3D seismic data volumes, multiple seismic attribute volumes are extracted for identifying fracture bodies. These multiple seismic attribute volumes include at least: coherent attribute volumes for characterizing stratigraphic discontinuities, curvature attribute volumes for characterizing structural curvature changes, and ant-tracking attribute volumes for enhancing fracture profiles.

[0022] Furthermore, the process for generating the initial fracture probability volume includes: Based on multiple seismic attribute bodies, at the drilled well points, according to the developed and undeveloped sections of the fracture body identified by the wellbore imaging logging data, corresponding sample data are extracted from multiple seismic attribute bodies to form a labeled training sample set. The training sample set is input into an interpretable machine learning algorithm for training, thereby generating a fracture body recognition model. Subsequently, the fracture body recognition model is used to perform a global scan and fusion calculation of multiple seismic attribute volumes in the target work area, and output an initial fracture body probability volume. While generating the initial fault probability volume, the interpretable machine learning algorithm simultaneously calculates and outputs the feature importance ranking of various seismic attribute volumes for the prediction results. The feature importance ranking is used to assist human experts in understanding the geological origin of the initial fault probability volume and to interactively interpret it.

[0023] Specifically, firstly, the original 3D seismic data volume collected from the target work area is subjected to amplitude preservation processing, mainly including amplitude compensation and frequency broadening, in order to restore the deep signal intensity and improve the vertical resolution. Subsequently, the preprocessed data volume is processed using a structurally guided filtering algorithm. This algorithm can effectively suppress random noise while protecting the boundary features of geological bodies such as faults and fractures, ultimately obtaining a 3D seismic data volume with a high signal-to-noise ratio. Based on this high signal-to-noise ratio data volume, various seismic attributes sensitive to fracture bodies are extracted. In this embodiment, the core extracted attributes include coherence attribute volume, curvature attribute volume, and ant-tracking attribute volume. Among them, the coherence attribute volume is used to characterize the discontinuity of the strata and is calculated using the third-generation C3 coherence algorithm. Its formula can be expressed as: ; in, To analyze the amplitude of the i-th seismic track within the time window, The coherence value represents the average amplitude of all participating traces; the lower the coherence value, the stronger the formation discontinuity. Curvature attribute volumes are used to characterize the curvature of formations; fractures often develop in areas of high curvature. Maximum curvature is calculated based on seismic horizons or data volumes. Ant tracking attribute body is an edge detection technology based on ant colony optimization algorithm, which can automatically enhance and connect the contour of the fracture system and generate fracture skeleton data body with better continuity; At drilled well points, using wellbore imaging logging (such as FMI / UBI) interpretation results as the gold standard, clearly identified fractured sections are labeled as positive samples (label: 1), while dense, fracture-free sections are labeled as negative samples (label: 0). Subsequently, along the well trajectory, attribute values ​​for corresponding depth points are extracted from various seismic attribute volumes such as coherence, curvature, and ant-like structures, forming a feature-label paired training sample set. ,in It is a feature vector containing multiple attribute values. These are the corresponding sample labels, and n represents the total number of samples; In this implementation, random forest is preferred as an interpretable machine learning algorithm. The sample set D is input into the random forest algorithm for training, generating a fracture recognition model. The random forest constructs multiple decision trees and integrates them, and its prediction probability for sample X is... The formula for calculating the probability (i.e., the probability that the location is a fracture) can be simplified to the average of the votes from all decision tree predictions: ; in, The total number of decision trees, The prediction result for the t-th tree is 1 or 0. As an indicator function, after training, the model is used to perform a global scan of the 3D seismic attribute volume of the entire work area, and the probability value of each grid point belonging to the fault zone is calculated. This generates a spatially continuous initial fracture probability volume; While generating the probability volume, the Random Forest algorithm simultaneously outputs a ranking of feature importance. This ranking is typically calculated based on metrics such as the reduction in Gini impurity or the decrease in average precision, and its calculation formula can be expressed as: ; in, This represents the j-th feature (such as a coherent attribute). It is a feature The sum of Gini impurity reductions at all nodes of the t-th decision tree, ranked by the importance of this feature, directly informs geological experts of the respective contributions of attributes such as coherence, curvature, and ant-body in the model's decision-making process. For example, if the ant-body has the highest importance, it means that this attribute contributes the most to identifying fracture bodies in this work area. This provides key evidence and entry points for experts to conduct subsequent manual interpretation, enabling them to understand the model's decision logic and make targeted corrections to the probability volume.

[0024] Human expert interaction is introduced to interpret and correct the initial fracture probability volume, and calibration is performed in combination with wellbore imaging logging data to output a calibrated three-dimensional spatial distribution model of the fracture volume. Large and medium-sized faults have been identified in the calibrated three-dimensional spatial distribution model of the fracture volume.

[0025] Furthermore, the output process of the calibrated three-dimensional spatial distribution model of the fracture body includes: Based on the initial probability volume of the fracture body and its corresponding feature importance ranking, human experts identify and track the fracture body structure in three-dimensional space. They also remove false anomalies in the initial fracture body probability volume based on geological knowledge, and enhance and supplement the blurred or missing fracture body structure to generate a preliminary interpretation result consisting of a series of three-dimensional spatial point, line and surface data. The preliminary interpretation results are precisely compared with the wellbore imaging logging data at the well trajectory location. When the preliminary interpretation results match the wellbore data, they are confirmed. When there is a contradiction, the preliminary interpretation results are corrected based on the wellbore imaging logging data to generate a fracture body interpretation result that has been corrected by well control. After completing the well-seismic calibration of the entire area, the interpretation results of the fracture bodies after well control correction are integrated in three-dimensional space and converted into a complete three-dimensional data volume model. In this model, large and medium-sized faults that reach the preset threshold value are clearly identified, and finally, the calibrated three-dimensional spatial distribution model of the fracture bodies is output.

[0026] Specifically, geological experts load an initial fault probability volume and its corresponding feature importance ranking into professional 3D geological interpretation software. First, they analyze the feature importance ranking to determine the key seismic attributes (such as coherence or ant-like structures) that dominate fault identification, thereby understanding the decision-making logic of machine learning and using it as a core reference. Subsequently, experts preliminarily extract candidate fault regions by setting a dynamic probability threshold. This threshold can be adjusted according to the actual situation of the work area. Its function can be described as follows: for any grid point voxel(x,y,z) in the probability volume, if its probability value P(x,y,z) satisfies... If so, it is initially selected as a candidate point for the fracture body. It can be set to 0.5 to 0.7; Based on this, experts implement core manual interventions, including the removal of false anomalies and enhancement and completion. The removal of false anomalies involves identifying and manually deleting high-probability anomalies that are inconsistent with geological laws, caused by channel sand body boundaries, special lithological bodies, or noise, by combining geological knowledge such as sedimentary facies diagrams and tectonic trends. The enhancement and completion involves manually tracing, connecting, and completing the drawing for intermittent or ambiguous fault responses of probability bodies, based on regional tectonic styles (such as imbricate and feather structures) and stress field directions, to ensure the rationality of the fault system and the integrity of its spatial distribution. The final output of this step is a preliminary interpretation result in vector format, consisting of a series of three-dimensional spatial point (such as fault intersection points), line (such as fault lines and fracture centerlines), and surface (such as fault plane triangular mesh) data. At each well trajectory location with imaging logging, the interpreted fracture features (such as fault planes) are checked one by one to see if they match the densely fractured zones, fault gouge, or structural stress concentration sections identified on the imaging logging. When the vertical depth of the interpreted fracture features at the well trajectory is within the allowable error range of the imaging logging response depth... When the distance is consistent within ±10 meters, that is... If the following contradictions occur, corrections should be made using well data as the absolute benchmark: False positive (multiple interpretations): The interpreted fracture body crosses the well, but the corresponding depth segment in the imaging log shows no response. In this case, the interpretation element at that location should be deleted. False negative (missed interpretations): The imaging log shows a significant fracture zone or breakpoint. But There is no corresponding explanation within the scope. In this case, A new breakpoint is added at the well point, and the interpretation section is corrected or added accordingly. This process generates a well-controlled corrected fracture body interpretation result with hard data guaranteeing accuracy at the well point. The interpretation results of fracture bodies after well control correction and calibration of all well points in the entire area are integrated in three-dimensional space. Using the gridding algorithm of geological modeling software (such as convergent interpolation or indicator kriging), the discrete point, line, and surface vector data are converted into a continuous three-dimensional gridded data volume model M(x,y,z), where each grid is assigned a value of 1 (representing the fracture body) or 0 (representing the matrix). In this model M, large and medium-sized faults need to be clearly identified. This invention calculates the extension length (L) of a single fault and the average fault displacement along the fault strike. To quantitatively identify faults, a preset threshold value is set, and faults that simultaneously meet the following conditions are identified as large or medium-sized faults: ; in, and For example, a threshold set according to the scale of the work area structure. =1000 meters, =20 meters; The final output, a calibrated three-dimensional spatial distribution model of the fracture body, is a high-precision model that combines the objective laws of machine learning, the cognitive experience of geological experts, and the hard data constraints of well points, providing a reliable structural framework for subsequent geomechanical analysis and fluid simulation.

[0027] Based on the calibrated three-dimensional spatial distribution model of the fracture body, the large and medium-sized faults identified in it are used as the deterministic framework. At the same time, the paleotectonic stress field is restored based on the regional tectonic evolution history, and the rock fracture criterion is used to conduct geomechanical simulation to predict the development zone, orientation and density of fractures.

[0028] Furthermore, the prediction process for the development zone, orientation, and density of cracks includes: Based on the calibrated three-dimensional spatial distribution model of fracture bodies, the spatial geometry and attitude of the identified large and medium-sized faults are extracted and used as an immovable and unmodifiable deterministic skeleton to be directly constructed in the initial three-dimensional discrete fracture network model. Based on the regional tectonic evolution history, key tectonic movement periods are identified, and paleotectonic restoration technology is used to reconstruct the paleotectonic framework and paleoburial depth of the target strata during the key hydrocarbon accumulation period. Then, the paleotectonic stress field of this period is inverted through rock mechanics parameters to obtain the magnitude and direction of paleostress. Among them, the deterministic framework, as the discontinuous interface in the paleotectonic framework, participates in the stress field simulation. Using the restored paleotectonic stress field as input, and combined with the rock mechanical parameters of the target strata, the rock fracture potential value of each grid point is calculated in a three-dimensional spatial grid based on the Coulomb fracture criterion or the Griffith criterion. Based on the calculated rock fracture potential value, the development zone, orientation, and density of fractures are quantitatively predicted. Specifically, the high value area of ​​the rock fracture potential value is predicted as the dominant fracture development zone, the direction of the maximum principal stress of paleostress at each grid point is predicted as the dominant orientation of the fracture, and the spatial distribution of the rock fracture potential value is mapped with the normalization function to quantitatively calculate the relative density distribution of fractures.

[0029] Specifically, based on a calibrated three-dimensional spatial distribution model of fracture bodies, spatial geometric elements of large and medium-sized faults are extracted, including their fault location, strike, dip angle, and plunging angle. These faults are used as a deterministic framework and directly imported into three-dimensional discrete fracture network modeling software. In subsequent modeling processes, the spatial location and geometry of this deterministic framework will be fixed as unmodifiable constraints. Based on the research results of regional tectonic evolution history, the key tectonic movement periods affecting the work area are determined. Using paleotectonic restoration techniques (such as software based on the principle of equilibrium profiles), the overlying strata are stripped back to reconstruct the paleotectonic framework and paleoburial depth of the target strata during the key hydrocarbon accumulation period. Based on this, numerical simulation methods (such as the finite element method) are used to invert the paleotectonic stress field of this period. In the simulation, the established deterministic skeleton is used as the internal friction discontinuity interface in the model for calculation. Its mechanical properties are described by the Coulomb friction criterion and can be expressed as: ; in, The shear stress on the cross-section, For fault cohesion, The normal stress (usually compressive stress) acting on the cross-section is used to obtain the three principal stresses of the target stratum during the paleotectonic period through inversion. ,and The magnitude and direction of the paleotectonic stress field; Using the restored paleotectonic stress field as input, combined with the target formation rock mechanical parameters (such as Young's modulus E, Poisson's ratio ν, cohesion) obtained through core experiments and well logging data, internal friction angle In a three-dimensional spatial grid, the rock fracture potential value at each grid point is calculated based on the Coulomb fracture criterion. The fracture value CF of the Coulomb fracture criterion can be calculated by the following formula: ; When CF≥1, it indicates that the rock has undergone shear fracture. The larger the CF value, the easier it is for the rock to fracture under the current stress state, that is, the higher the fracture potential. In this embodiment, the calculated CF value is used as the rock fracture potential value to form a three-dimensional CF(x,y,z) data volume. Based on the calculated 3D crack flow (CF) data, three key parameters of the crack are quantitatively predicted, including the crack development zone, crack orientation, and crack density. The crack development zone is defined by predicting areas with high CF values ​​as dominant crack development zones, achieved by setting a CF threshold (e.g., ...). ), which will satisfy CF(x,y,z)≥ The grid area is delineated as a favorable target area for crack development, with the crack orientation as shown. The direction of the maximum principal stress at each grid point in the paleostress field is determined. The direction of the crack is predicted to be the dominant orientation for cracks (especially shear cracks) to originate at that location, because the crack orientation is usually perpendicular to the minimum principal stress. And the tilt angle and intermediate principal stress Correspondingly, the crack density is such that the spatial distribution of CF values ​​is mapped to the relative density distribution of cracks through a normalization function. An exemplary linear mapping formula is as follows: ; in, and These are the maximum and minimum values ​​in the entire 3D CF data volume, respectively. It is a dimensionless value between 0 and 1, whose magnitude represents the relative density of crack development at that point compared to other areas within the work area. The closer the value is to 1, the denser the cracks.

[0030] Using the results of geomechanical simulation as constraints, nondeterministic fractures are generated outside the deterministic framework to jointly construct a three-dimensional discrete fracture network model.

[0031] Furthermore, the construction process of the three-dimensional discrete crack network model includes: The dominant crack development zone, the dominant crack orientation, and the relative crack density distribution are respectively converted into a crack density field and a crack strike-dip field in three-dimensional space as constraint parameter fields. Based on the crack density field and crack orientation-dip field, a large number of discrete nondeterministic cracks are generated in three-dimensional space using a stochastic simulation method. The spatial location, development scale, orientation and density of the nondeterministic cracks are strictly controlled by the constraint parameter field. The deterministic skeleton and the non-deterministic cracks are fused and their geometric topological relationships are checked in three-dimensional space to ensure that the spatial contact relationships between cracks of different scales are reasonable, and finally a complete three-dimensional discrete crack network model is constructed together.

[0032] Specifically, the fracture parameters predicted by geomechanical simulation are converted into a three-dimensional constrained parameter field that can be used to control stochastic simulations, wherein the predicted relative density distribution of fractures ( Converted to the surface density field commonly used in engineering (P32, unit: m) 2 / m 3 The conversion formula can be expressed as: P32 ; in, and The upper and lower limits of fracture surface density (e.g., based on core observations and imaging logging statistics in the work area) are determined. =1.0m 2 / m 3 , =0.1m 2 / m 3 The resulting P32(x,y,z) field directly controls the density of crack development at various points in space during subsequent simulations, and predicts the dominant crack orientation (i.e., the maximum principal stress). The direction is converted into the crack orientation and dip angle. For a grid point, the crack orientation α and dip angle β can be calculated from the stress direction vector of that point, ultimately forming a three-dimensional orientation field α(x,y,z) and dip angle field β(x,y,z), which serve as the dominant directional constraint controlling the crack orientation. Based on the aforementioned constraint parameter field, a large number of nondeterministic cracks are generated in three-dimensional space using stochastic simulation methods (such as stochastic modeling algorithms based on Poisson point processes). The simulation process is strictly controlled. Specifically, the spatial location is controlled such that the probability distribution of crack center points in space is proportional to the P32 field; that is, in regions with high P32 values, crack center points are more densely packed. The orientation (azimuth) is controlled such that the strike and dip angle of the generated cracks are no longer completely random, but rather follow a Gaussian distribution (or von Mises distribution) with means α(x,y,z) and β(x,y,z). Its probability density function can be expressed as: ; in, This represents the direction or dip angle of the crack to be generated. The dominant strike α or dominant dip angle β represents the location, and σ is the standard deviation that controls the degree of dispersion of the attitude (which can be set according to geological uncertainty, for example, σ=15°). For size control, the radius (r) of the crack is usually assumed to follow a power-law distribution, with the probability density function being: ; in, The power-law exponent is C, which is the normalization constant. The size distribution of cracks is controlled by this function. The imported deterministic framework (large and medium-sized faults) and randomly generated non-deterministic fractures (small and micro fractures) are fused in three-dimensional space. Specifically, a geometric algorithm is used to automatically handle the termination of fractures based on their intersection relationships. For example, a rule is set that low-level fractures terminate when they encounter high-level fractures. Typically, non-deterministic fractures terminate when they encounter the deterministic fault framework, thus forming a dense fracture zone near the fault. Then, the spatial contact relationship between all fracture elements is checked and ensured to avoid non-geological suspended fractures or abnormal intersection relationships. Unreasonable connections are manually or automatically corrected. Finally, a complete three-dimensional discrete fracture network model with reasonable hierarchical and spatial configuration relationships, including large faults to micro fractures, is constructed. This model not only has realistic geometric morphology but also has a solid geomechanical genetic basis, providing a key carrier for subsequent reservoir property calculations and fluid flow simulations.

[0033] Obtain the pre-established reservoir matrix property model of the target work area, couple the constructed three-dimensional discrete fracture network model with the reservoir matrix property model to construct an integrated three-dimensional geological model of fracture body-matrix, and define the uncertainty parameters and their value range in the integrated three-dimensional geological model.

[0034] Furthermore, the construction process of the integrated three-dimensional geological model of the fracture body and matrix includes: Obtain the pre-established reservoir matrix property model of the target work area. The reservoir matrix property model includes a three-dimensional property field of matrix porosity and permeability. Spatial superposition and attribute coupling are performed between the three-dimensional discrete fracture network model and the reservoir matrix attribute model; Based on the geometric parameters and conductivity of each fracture in the three-dimensional discrete fracture network model, the permeability properties of the matrix mesh through which the fracture passes are calculated, thereby transforming the background permeability field of the matrix into an equivalent permeability field that can characterize the interaction between the fracture body and the matrix. The transformed equivalent permeability field is integrated with the matrix porosity field to form a dual-media model that simultaneously characterizes the high-speed seepage channels of the fracture body and the matrix storage space in terms of properties, namely, the fracture-matrix integrated three-dimensional geological model.

[0035] Furthermore, the definition process for uncertainty parameters and their value ranges in the integrated three-dimensional geological model includes: Based on the construction process of the three-dimensional discrete crack network model and the calculation principle of the equivalent permeability field, key uncertainty parameters that significantly affect the fluid simulation results of the model are identified. Based on the core test data of the target work area, the fracture aperture statistics interpreted by imaging logging, and the dynamic permeability inverted from the well test or production history data, a predetermined value range is determined for each of the identified key uncertainty parameters. Each key uncertainty parameter is combined with its corresponding value range to form a multi-dimensional parameter space, which is used to characterize all geological uncertainties of the integrated three-dimensional geological model before historical matching.

[0036] Specifically, a reservoir matrix property model pre-established using geostatistical methods is obtained for the target work area. This model consists of two three-dimensional data volumes, specifically the matrix porosity field. and matrix permeability field Simultaneously, the three-dimensional discrete crack network model constructed in the aforementioned steps is loaded. Before coupling, it must be ensured that the two models are in the same three-dimensional spatial coordinate system and mesh system. For each matrix mesh traversed by the crack, its equivalent permeability tensor The calculation is crucial. This embodiment uses the Oda method or its improved version for approximate calculation. The principle of the Oda method is to accumulate the contribution of fractures to the mesh permeability tensor. For a given matrix mesh, the increment of its equivalent permeability tensor is... The contribution from all crack segments passing through the mesh is superimposed and can be expressed as: ; in, It is the volume of the matrix mesh. This represents summing over the i-th crack segment that passes through the grid. It is a dimensionless correction factor for the conductivity of the crack, and is an empirical parameter related to the degree of crack filling. It is the hydraulic aperture of the i-th fracture segment within the grid. It is a unit tensor. It is the unit normal vector of the i-th crack segment. This formula represents the tensor product (exterior product). A tensor was constructed with its principal direction perpendicular to the crack surface to characterize the strong anisotropic permeability caused by the crack. Finally, the equivalent permeability tensor of this mesh was determined. The sum of matrix permeability and the contribution of all cracks, i.e.: ; By traversing all the grids, each grid's (x,y,z) are replaced with the calculated values. (x,y,z), thus generating a new equivalent permeability field; The newly generated equivalent permeability field (x,y,z) Original matrix porosity field For integration, regarding porosity, it is generally considered that the storage capacity of cracks is much smaller than that of the matrix. Therefore, in this embodiment, the porosity field of the integrated model still mainly adopts... Ultimately, this results in a structure that simultaneously contains attributes from... and The matrix pore storage and seepage system and the fracture high-velocity conduction network explicitly characterized by the DFN model, this integrated model with dual media properties, is the fracture-matrix integrated three-dimensional geological model. Based on the principles of DFN modeling and equivalent calculation, the key uncertainty parameters that have the most significant impact on fluid simulation results were identified. These parameters include at least the fracture equivalent permeability multiplier (FPS). ), fracture hydraulic aperture (a) and matrix absolute permeability ( ), where the fracture equivalent permeability multiplier ( ) represents a global multiplier used to correct for potential systematic biases in the Oda method. The final equivalent permeability field of the fracture used for simulation is: The hydraulic aperture of the fracture (a) is a key parameter in the DFN model, directly affecting... The calculation of the matrix absolute permeability has the greatest uncertainty. This can be explained as follows: the matrix background permeability itself also has interpretable and modeling uncertainties. Based on hard data from multiple sources, a reasonable range of values ​​for the above parameters is determined. The fracture hydraulic aperture (a) is implemented as follows: based on the statistical results of fracture static aperture interpreted by imaging logging, its distribution range is determined. For example, the P10 to P90 quantiles are taken as the lower limit. and upper limit ; Crack equivalent permeability multiplier ( The implementation involves calibrating the ratio of the system permeability obtained through well test interpretation to the initial equivalent permeability calculated by the model. The range of this ratio can be set as follows: ; Absolute permeability ( The implementation is as follows: its uncertainty range is directly derived from the standard deviation given in the geostatistical modeling. Its value range can be set to ; Combining the above parameters with their range of values ​​forms a multidimensional parameter space Ω, where any point p in this space represents a specific combination of parameters, which can be expressed as: ]; This parameter space Ω quantitatively characterizes all the key geological uncertainties contained in the integrated three-dimensional geological model prior to historical matching.

[0037] An experimental design method is used to generate an integrated model set consisting of multiple parameters. Fluid numerical simulation is performed on each model in the integrated model set, and the simulated production dynamic data is automatically compared with the actual historical production data of the oilfield. A set of equivalent models whose historical matching error meets the preset conditions is selected.

[0038] Furthermore, the screening process for a set of equivalent models includes: Using Latin hypercube sampling or factorial analysis experimental design methods, a large number of samples are taken in a multi-dimensional parameter space. Each set of specific parameters obtained from each sampling constitutes a specific realization of the integrated three-dimensional geological model. Through multiple sampling, an integrated model set containing hundreds to thousands of different parameter realizations is generated. Each model implementation in the integrated model set is loaded sequentially into a pre-made fluid numerical simulator. Under the same well location, production regime, and simulation duration settings, batch numerical simulation calculations are performed to obtain the simulated production dynamic data corresponding to each model implementation. The simulated production dynamic data of each model is automatically compared with the actual historical production data of the oilfield, and the historical matching error is calculated. Set a historical matching error threshold as a preset condition, implement all models whose historical matching errors are lower than the threshold, and select them as a group of equivalent models.

[0039] Specifically, this embodiment employs the Latin hypercube sampling (LHS) experimental design method to efficiently sample within the defined multidimensional parameter space Ω. The advantage of LHS is that it can comprehensively and unbiasedly cover the entire parameter space with fewer sampling attempts. Assuming the parameter space Ω contains M uncertain parameters and N model implementations are planned, the LHS sampling steps are as follows: [The text then abruptly shifts to a different topic:] ...the domain of each parameter... The system divides the parameter into N equally spaced intervals. For each parameter, a value is randomly selected from each of the N intervals. These M randomly selected values ​​are then combined to ensure that each parameter is selected only once in each interval, resulting in N M-dimensional sampling points. Each sampling yields a specific set of parameter combinations. That is, a specific implementation of an integrated three-dimensional geological model. A unified model set is generated through N samplings (N is typically hundreds to thousands). ; The generated model set Each model implementation is sequentially loaded into a pre-built fluid numerical simulator (such as Eclipse, CMG, or other commercial software). The simulator needs to be pre-configured with simulation conditions that are completely consistent with the actual production history. Specifically, the well location and trajectory are consistent with the real oilfield, the production regime includes constant oil production and constant bottom hole flowing pressure, and the simulation duration covers the entire historical production period that needs to be matched. Subsequently, batch numerical simulation calculations are performed to obtain the results for each model implementation. The corresponding simulated production dynamic data typically includes, but is not limited to, bottom hole flowing pressure. Daily oil production Daily gas production and daily water production ; The simulated production dynamics data from each model are automatically compared with the actual historical production data of the oilfield. To do this, a quantified historical matching error function needs to be defined. This is used to measure how well the i-th model fits the historical data. A comprehensive error function is usually the weighted sum of squared errors of each dynamic parameter. For the i-th model, its total error is... The calculation is as follows: ; in, The pressure matching error of the i-th model is calculated as follows: ; Let be the production matching error of the i-th model, calculated as follows: ; and These are weighting coefficients used to balance the contributions of pressure and yield data to the total error, and are typically set based on data quality and importance (e.g., =0.5, =0.5); Set a historical matching error threshold As a prerequisite, the historical matching error of all N models is... Models below this threshold are filtered out, i.e., satisfying the following: ; These selected models constitute a set of equivalent models, whose geological parameters ( Although they are different, their dynamic responses can be matched well with historical production data. Therefore, from the perspective of fluid flow, they are all equivalent interpretations of the actual underground conditions. This set of equivalent models represents all possible geological conditions under the existing data, providing a probabilistic basis for the next step of prediction and decision-making.

[0040] Production forecasting is performed using a set of equivalent models selected from the pool, and the forecast results are compared and verified with recent production data that has not been matched with historical data, until a final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements is obtained.

[0041] Furthermore, the process for obtaining the final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements includes: Using a selected set of equivalent models, production forecasts for future periods are performed in a fluid numerical simulator. The predicted production dynamics data for each equivalent model are obtained, and the forecast results from the set of equivalent models are combined to form a comprehensive forecast curve covering a range of uncertainties. Obtain the actual production data newly generated after the end of the historical matching period, as the recent production data that did not participate in the historical matching, and compare the comprehensive prediction curve band with the recent production data. If the recent production data falls within the range of the comprehensive prediction curve band, the model is deemed to have passed the validation. When the verification is successful, the model with the smallest historical matching error is selected from a set of equivalent models. The integrated three-dimensional geological model corresponding to it and the three-dimensional discrete fracture network model that forms the basis of its construction are back-determined as the most reliable model in terms of geological understanding. The model is then output as the final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements. If the verification fails, it indicates that there is a systematic bias in the geological understanding of a set of equivalent models. At this time, based on the difference between the prediction results and the latest production data, the cause of the geological understanding bias is analyzed, and the process returns to the step of introducing human expert interaction to interpret and correct the initial fracture probability volume, or to the step of generating non-deterministic fractures outside the deterministic framework using the results of geomechanical simulation as constraints. After correcting the geological parameters and modeling scheme, the subsequent steps are re-executed until the verification passes.

[0042] Specifically, a set of equivalent models selected using historical matching: (Where M is the number of equivalent models), the simulation for the prediction period is performed in a pre-built fluid numerical simulator. The prediction period is set after the end of the historical production period. The simulation adopts a production strategy consistent with the actual oilfield's future plan (such as constant pressure production). For each equivalent model... This allows us to obtain a series of predictive production dynamics data during the forecast period, such as the cumulative gas production forecast curve. ; By synthesizing the prediction results of all equivalent models, a comprehensive prediction curve is generated to quantify the prediction range caused by geological uncertainties, for any prediction time point. Calculate the statistical distribution of all equivalent model predictions at that moment. The prediction curve band at that moment is then defined by an upper and lower limit. A commonly used definition is to take the P10 and P90 quantiles, where the upper limit of the prediction curve is expressed as: ; The lower limit of the prediction curve is expressed as: ; thus, This forms a comprehensive forecast curve band from the present to a certain point in the future, acquiring newly generated, recent actual production data that has not participated in any model calibration after the historical matching period ends. The recent data is compared and verified with the predicted curve band. The verification condition is that, within a certain verification period, the vast majority of recent data points should fall within the predicted curve band, which can be quantified as follows: ; For the vast majority If this condition is met, the model is considered to have passed validation. If the validation passes, select historical matching errors from a set of equivalent models. The smallest model ,Right now: ; The best model The corresponding integrated three-dimensional geological model, the three-dimensional discrete fracture network model that constructs the integrated model, and the three-dimensional spatial distribution model of the fracture body that is artificially calibrated and serves as the basis of the DFN model are back-determined to be the most reliable model in terms of geological understanding under the current data and cognitive level. Finally, this three-dimensional spatial distribution model of the fracture body is output as the final model that meets the accuracy requirements and is used to guide oilfield development decisions. When the verification fails, it indicates a systematic bias in the current geological understanding represented by the equivalent model set. In this case, it is necessary to analyze the difference pattern between the prediction results and recent data, such as whether it is a systematic overestimation or underestimation. Based on the difference analysis results, return to the key steps in the aforementioned process for correction. Specifically, if the bias stems from inaccurate understanding of the spatial distribution of the fracture body, return to the step of introducing human expert interaction to interpret and correct the initial fracture body probability volume, and adjust the interpretation scheme of faults and fractures. If the bias stems from unreasonable parameters such as fracture density or aperture, return to the step of generating non-deterministic fractures outside the deterministic framework using the results of geomechanical simulation as constraints, and adjust the constraints of geomechanical parameters or DFN modeling. After correcting the geological understanding and modeling scheme, re-execute all modeling, coupling, history matching, and prediction verification processes from this step onwards, forming an iterative optimization closed loop until the model's prediction results are independently verified by recent data, thereby finally obtaining a realistic and reliable fracture body model.

[0043] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for depicting fracture bodies in the fracture development zone, characterized in that, Includes the following steps: Acquire 3D seismic data and wellbore imaging logging data of the target work area, process the 3D seismic data, extract multiple seismic attribute volumes, and use interpretable machine learning algorithms to fuse the multiple seismic attribute volumes to generate an initial fracture probability volume that can be interactively interpreted by human experts. Human expert interaction is introduced to interpret and correct the initial fracture probability volume, and the wellbore imaging logging data is combined for calibration to output a calibrated three-dimensional spatial distribution model of the fracture volume, in which large and medium-sized faults have been identified. Based on the calibrated three-dimensional spatial distribution model of the fracture body, the large and medium-sized faults identified therein are used as a deterministic framework. At the same time, the paleotectonic stress field is restored based on the regional tectonic evolution history, and the rock fracture criterion is used to perform geomechanical simulation to predict the development zone, orientation and density of fractures. Using the results of the geomechanical simulation as constraints, non-deterministic fractures are generated outside the deterministic framework to jointly construct a three-dimensional discrete fracture network model. Obtain the reservoir matrix attribute model pre-established for the target work area, couple the constructed three-dimensional discrete fracture network model with the reservoir matrix attribute model to construct an integrated three-dimensional geological model of fracture body-matrix, and define the uncertainty parameters and their value ranges in the integrated three-dimensional geological model. An experimental design method is used to generate an integrated model set consisting of multiple parameters. Fluid numerical simulation is performed on each model in the integrated model set, and the simulated production dynamic data is automatically compared with the actual historical production data of the oilfield. A set of equivalent models whose historical matching error meets the preset conditions is selected. Production forecasting is performed using the selected set of equivalent models, and the forecast results are compared and verified with recent production data that has not been matched with historical data, until a final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements is obtained.

2. The method for depicting the fracture development zone of a fracture body according to claim 1, characterized in that, The extraction process for the various seismic attribute volumes includes: The acquired three-dimensional seismic data is subjected to amplitude compensation and frequency broadening processing to enhance the effective signal and form a preprocessed three-dimensional seismic data volume. The preprocessed 3D seismic data volume is subjected to construction-guided filtering to suppress noise and preserve the boundary features of the fracture body, thereby forming a 3D seismic data volume with a high signal-to-noise ratio. Based on the high signal-to-noise ratio three-dimensional seismic data volume, multiple seismic attribute volumes for identifying fracture bodies are extracted. These multiple seismic attribute volumes include at least: a coherent attribute volume for characterizing stratigraphic discontinuities, a curvature attribute volume for characterizing structural curvature changes, and an ant-tracking attribute volume for enhancing fracture profiles.

3. The method for depicting the fracture development zone of a fracture body according to claim 1, characterized in that, The process for generating the initial fracture probability body includes: Based on the aforementioned multiple seismic attribute bodies, at the drilled well points, according to the developed and undeveloped sections of the fracture body identified by the wellbore imaging logging data, corresponding sample data are extracted from the multiple seismic attribute bodies to form a labeled training sample set. The training sample set is input into the interpretable machine learning algorithm for training, thereby generating a fracture body recognition model. Subsequently, the fracture body recognition model is used to perform global scanning and fusion calculation of the various seismic attribute bodies of the target work area, and output an initial fracture body probability body. While generating the initial fault probability volume, the interpretable machine learning algorithm simultaneously calculates and outputs the feature importance ranking of the various seismic attribute volumes for the prediction results. The feature importance ranking is used to assist human experts in understanding the geological origin of the initial fault probability volume and to interactively interpret it.

4. The method for depicting the fracture development zone of a fracture body according to claim 3, characterized in that, The output process of the calibrated three-dimensional spatial distribution model of the fracture body includes: Based on the initial fracture probability volume and its corresponding feature importance ranking, artificial experts identify and track fracture structures in three-dimensional space, and eliminate false anomalies in the initial fracture probability volume according to geological knowledge, enhance and supplement ambiguous or missing fracture structures, and generate a preliminary interpretation result consisting of a series of three-dimensional spatial point, line and surface data. The preliminary interpretation results are precisely compared with the wellbore imaging logging data at the well trajectory position. When the preliminary interpretation results match the wellbore data, they are confirmed. When there is a contradiction, the preliminary interpretation results are corrected based on the wellbore imaging logging data to generate a fracture body interpretation result after well control correction. After completing the well-seismic calibration of the entire area, the interpretation results of the fracture body after well control correction are integrated in three-dimensional space and converted into a complete three-dimensional data volume model. In this model, large and medium-sized faults with a scale reaching the preset threshold value are clearly identified, and finally the calibrated three-dimensional spatial distribution model of the fracture body is output.

5. The method for depicting the fracture development zone of a fracture body according to claim 1, characterized in that, The prediction process for the development zone, orientation, and density of the cracks includes: Based on the calibrated three-dimensional spatial distribution model of the fracture body, the spatial geometry and attitude of the identified large and medium-sized faults are extracted and used as an immovable and unmodifiable deterministic skeleton to be directly constructed in the initial three-dimensional discrete fracture network model. Based on the regional tectonic evolution history, key tectonic movement periods are determined, and paleotectonic restoration technology is used to reconstruct the paleotectonic framework and paleoburial depth of the target strata during the key hydrocarbon accumulation period. Then, the paleotectonic stress field of that period is inverted through rock mechanics parameters to obtain the magnitude and direction of paleostress. The deterministic framework, as a discontinuous interface in the paleotectonic framework, participates in the stress field simulation. Using the restored paleotectonic stress field as input, and combined with the rock mechanics parameters of the target strata, the rock fracture potential value of each grid point is calculated in a three-dimensional spatial grid based on the Coulomb fracture criterion or the Griffith criterion. Based on the calculated rock fracture potential value, the development zone, orientation, and density of fractures are quantitatively predicted. Specifically, the high value area of ​​the rock fracture potential value is predicted as the dominant fracture development zone, the direction of the maximum principal stress of paleostress at each grid point is predicted as the dominant orientation of the fracture, and the spatial distribution of the rock fracture potential value is mapped to a normalization function to quantitatively calculate the relative density distribution of fractures.

6. The method for depicting the fracture development zone of a fracture body according to claim 5, characterized in that, The construction process of the three-dimensional discrete crack network model includes: The dominant crack development zone, the dominant crack orientation, and the relative crack density distribution are respectively converted into a crack density field and a crack strike-dip field in three-dimensional space as constraint parameter fields. Based on the crack density field and the crack strike-dip field, a large number of discrete nondeterministic cracks are generated in three-dimensional space using a stochastic simulation method. The spatial location, development scale, orientation and density of the nondeterministic cracks are strictly controlled by the constraint parameter field. The deterministic skeleton and the non-deterministic cracks are fused and their geometric topological relationships are checked in three-dimensional space to ensure that the spatial contact relationships between cracks of different scales are reasonable, and finally a complete three-dimensional discrete crack network model is constructed together.

7. The method for depicting the fracture development zone of a fracture body according to claim 6, characterized in that, The construction process of the integrated three-dimensional geological model of the fracture body and matrix includes: Obtain a pre-established reservoir matrix attribute model for the target work area, wherein the reservoir matrix attribute model includes a three-dimensional attribute field of matrix porosity and permeability; The three-dimensional discrete fracture network model is spatially superimposed and coupled with the reservoir matrix property model; Based on the geometric parameters and conductivity of each crack in the three-dimensional discrete crack network model, the permeability properties of the matrix mesh through which the crack passes are calculated, thereby transforming the background permeability field of the matrix into an equivalent permeability field that can characterize the interaction between the fracture body and the matrix. The transformed equivalent permeability field is integrated with the matrix porosity field to form a dual-media model that simultaneously characterizes the high-speed seepage channels of the fracture body and the matrix storage space, namely the fracture-matrix integrated three-dimensional geological model.

8. The method for depicting the fracture development zone of a fracture body according to claim 1, characterized in that, The process for defining the uncertainty parameters and their value ranges in the integrated three-dimensional geological model includes: Based on the construction process of the three-dimensional discrete crack network model and the calculation principle of the equivalent permeability field, key uncertainty parameters that significantly affect the fluid simulation results of the model are identified. Based on the core test data of the target work area, the fracture aperture statistics interpreted by imaging logging, and the dynamic permeability retrieved from the well test or production history data, a predetermined value range is determined for each of the identified key uncertainty parameters. Each of the identified key uncertainty parameters is combined with its corresponding value range to form a multi-dimensional parameter space, which is used to characterize all geological uncertainties of the integrated three-dimensional geological model before historical matching.

9. The method for depicting the fracture development zone of a fracture body according to claim 8, characterized in that, The screening process for the set of equivalent models includes: Using Latin hypercube sampling or factorial analysis experimental design methods, a large number of samples are taken in the multi-dimensional parameter space. Each set of specific parameters obtained from each sampling constitutes a specific implementation of the integrated three-dimensional geological model. Through multiple sampling, an integrated model set containing hundreds to thousands of different parameter implementations is generated. Each model implementation in the integrated model set is sequentially loaded into a pre-made fluid numerical simulator. Under the same well location, production system and simulation duration settings, batch numerical simulation calculations are performed to obtain the simulated production dynamic data corresponding to each model implementation. The simulated production dynamic data of each model is automatically compared with the actual historical production data of the oilfield, and the historical matching error is calculated. A historical matching error threshold is set as the preset condition. All models with historical matching errors lower than the threshold are implemented and selected as the group of equivalent models.

10. The method for depicting the fracture development zone of a fracture body according to claim 1, characterized in that, The process for obtaining the final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements includes: Using the selected set of equivalent models, production forecasts for future periods are performed in the fluid numerical simulator. For each equivalent model, predicted production dynamics data are obtained. The forecast results from all the equivalent models are then combined to form a comprehensive forecast curve covering a range of uncertainties. Obtain the actual production data newly generated after the end of the historical matching period, and use it as the recent production data that did not participate in the historical matching. Compare the comprehensive prediction curve band with the recent production data. If the recent production data falls within the range of the comprehensive prediction curve band, the model is deemed to have passed the verification. When the verification is passed, the model with the smallest historical matching error is selected from the set of equivalent models. The integrated three-dimensional geological model corresponding to it and the three-dimensional discrete fracture network model that forms the basis of its construction are back-determined as the most reliable model in terms of geological understanding. The model is then output as the final three-dimensional spatial distribution model of the fracture body that meets the accuracy requirements. If the verification fails, it indicates that there is a systematic bias in the geological understanding of the set of equivalent models. At this time, based on the difference between the prediction results and the recent production data, the cause of the geological understanding bias is analyzed, and the process returns to the step of introducing human expert interaction to interpret and correct the initial fracture probability volume, or to the step of generating non-deterministic fractures outside the deterministic framework using the results of the geomechanical simulation as constraints. After correcting the geological parameters and modeling scheme, the subsequent steps are re-executed until the verification passes.

Images