Three-dimensional geologic model modeling method and device

Abstract

The invention discloses a three-dimensional geologic model modeling method and device. The three-dimensional geological model modeling method comprises the following steps: acquiring to-be-established geological region data and actual geological boundary data; generating an initial stratum interface grid according to the to-be-established geological region data; correcting the initial stratigraphic interface grid to obtain a corrected stratigraphic interface grid; carrying out geological rule constraint on the corrected stratigraphic interface grid so as to form simple three-dimensional stratigraphic interface model data; correcting the simple three-dimensional stratigraphic interface model data according to the seismic profile information so as to obtain a composite geologic structure model; generating a three-dimensional fault plane model according to the composite geologic structure model; and generating a composite geological structure model including a fault network according to the three-dimensional fault plane model and the composite geological structure model. Through the three-dimensional geological model modeling method, the deviation between the model and an actual geological boundary can be remarkably reduced.

Description

Technical Field

[0001] The present application relates to the field of geological modeling technology, and in particular to a three-dimensional geological modeling method and a three-dimensional geological modeling device. Background Art

[0002] In geological exploration and resource evaluation, three-dimensional geological models are core tools for characterizing underground structures, analyzing geological patterns, and guiding engineering practices. Traditional three-dimensional geological modeling methods rely primarily on single data sources (such as borehole data or seismic profiles) for trend surface fitting or implicit modeling, which suffers from the following drawbacks:

[0003] Distortion of geological details: Relying solely on discrete control points (such as drill hole data) for global fitting makes it difficult to accurately depict local complex structures such as fault traces and lithologic mutation zones, resulting in centimeter to meter deviations between the model and actual geological boundaries (such as fault edges and unconformities).

[0004] Insufficient data fusion capabilities: A collaborative constraint mechanism for multi-source data (such as seismic profiles, gravity anomalies, and tectonic stress fields) has not been established, resulting in limited model expression of hidden structures (such as fold axis traces and salt dome boundaries). Summary of the Invention

[0005] The object of the present invention is to provide a three-dimensional geological modeling method to solve at least one of the above technical problems.

[0006] One aspect of the present invention provides a three-dimensional geological modeling method, the three-dimensional geological modeling method comprising:

[0007] Obtain data on the geological area to be established and actual geological boundary data;

[0008] Generate an initial stratigraphic interface grid based on the geological area data to be established;

[0009] Modifying the initial stratum interface grid to obtain a modified stratum interface grid;

[0010] The modified stratigraphic interface grid is constrained by geological rules to form simple three-dimensional stratigraphic interface model data;

[0011] The simple three-dimensional stratum interface model data is modified according to the seismic profile information to obtain a composite geological structure model;

[0012] Generate a three-dimensional fault surface model based on the composite geological structure model;

[0013] A composite geological structure model including a fault network is generated according to the three-dimensional fault surface model and the composite geological structure model.

[0014] Optionally, the geological area data to be established includes discrete stratigraphic data, stratigraphic boundary control points, gravity anomaly data, magnetic anomaly data, geological outcrop profile data, stratigraphic age data and lithologic data;

[0015] Generating an initial stratigraphic interface grid based on the geological area data to be established includes:

[0016] Gridding the geological area data to be established to form multiple grid data points;

[0017] Perform first-layer downsampling on each grid data point to obtain the information of each first-layer data point;

[0018] Performing second-layer downsampling on each first-layer data point information to obtain each second-layer data point information;

[0019] Performing third-layer downsampling on each second-layer data point information to obtain each third-layer data point information;

[0020] Performing fourth-layer downsampling on each third-layer data point information to obtain each fourth-layer data point information;

[0021] The geological rule coding information is associated with the first layer data point information, the second layer data point information, the third layer data point information and the fourth layer data point information respectively;

[0022] Using a quintic polynomial model to roughly fit the fourth layer data point information and the geological rule coding information associated with the fourth layer data point information, thereby obtaining the fourth layer polynomial coefficients;

[0023] Using a quadratic polynomial model to fit each first layer data point information and the geological rule coding information associated with each layer to obtain the first layer polynomial coefficients;

[0024] Using a quadratic polynomial model to fit each second layer data point information and the geological rule coding information associated with each layer to obtain the second layer polynomial coefficients;

[0025] Using a quadratic polynomial model to fit each third layer data point information and the geological rule coding information associated with each layer to obtain the third layer polynomial coefficients;

[0026] Calculate the construction complexity of each first-layer data point information, each second-layer data point information, and each third-layer data point information respectively;

[0027] When the construction complexity of any one of the first-layer data point information exceeds a preset threshold, the first-layer data point information exceeding the preset threshold is locally refitted through a quintic polynomial model to obtain a corrected first-layer polynomial coefficient;

[0028] When the construction complexity of any one of the second-layer data point information exceeds a preset threshold, the second-layer data point information exceeding the preset threshold is locally refitted through a quintic polynomial model to obtain a corrected second-layer polynomial coefficient;

[0029] When the construction complexity of any information in the third-layer data points exceeds a preset threshold, the third-layer data points exceeding the preset threshold are locally refitted through a quintic polynomial model to obtain a corrected third-layer polynomial coefficient;

[0030] A joint probability model is constructed based on the obtained corrected first-layer polynomial coefficients, corrected second-layer polynomial coefficients, corrected third-layer polynomial coefficients, fourth-layer polynomial coefficients, and various geological rule coding information;

[0031] Obtain attribute probability distribution of stratigraphic age and lithologic attributes;

[0032] Regular three-dimensional grid data is generated according to the optimized trend surface parameters. The regular three-dimensional grid data and the attribute probability distribution of the stratigraphic age and lithologic attributes constitute an initial stratigraphic interface grid.

[0033] Optionally, the actual geological boundary data include fault traces, unconformity surface contours, and measured point data of lithologic contact zones interpreted from seismic profiles; the geological area data to be established further include tectonic stress field data, Bouguer gravity anomaly data, and magnetic susceptibility anomaly data;

[0034] The step of correcting the initial stratum interface grid to obtain a corrected stratum interface grid includes:

[0035] Based on the initial stratigraphic interface grid, the actual geological boundary data is mapped to the initial stratigraphic interface grid surface through cubic spline interpolation, thereby obtaining the initial fused actual geological boundary grid model;

[0036] A parameterized constraint model is established based on the actual geological boundary data, thereby correcting the grid node coordinates in the initial grid model integrated with the actual geological boundary, thereby obtaining a final grid model integrated with the actual geological boundary;

[0037] The key geological boundary constraints are interpolated on the final grid model integrated with the actual geological boundary through the geological rule library to obtain the grid model after interpolation constraints;

[0038] The conflicting area geometry is corrected on the interpolation constrained grid model to obtain the corrected stratum interface grid.

[0039] Optionally, performing key geological boundary constraint interpolation on the final fused actual geological boundary grid model through the geological rule library to obtain the grid model after interpolation constraint includes:

[0040] Control points are arranged along the fault traces of the final fused actual geological boundary grid model, and a continuous fault boundary is generated using a cubic B-spline curve, thereby forming an initial fault geometry framework in the final fused actual geological boundary grid model;

[0041] Based on the initial fault geometric framework and the final fused actual geological boundary grid model, outside the fault influence domain, the fault boundary generated by the B-spline is used as a constraint, and a regional trend surface is generated through RBF smooth transition, thereby forming the fault boundary and the regional trend surface on the final fused actual geological boundary grid model that forms the initial fault geometric framework;

[0042] Establishing a lithologic probability field based on the final fusion actual geological boundary grid model and the measured point data of the lithologic contact zone in the actual geological boundary data, defining single point potential energy and adjacent node potential energy, thereby obtaining a lithologic probability field model;

[0043] Based on the lithologic probability field model and the measured point data of the lithologic contact zone, gradient constrained interpolation is performed along the normal direction of the lithologic mutation zone to obtain the morphology of the lithologic mutation zone;

[0044] Generate stress field constraints based on tectonic stress field data on the final fused actual geological boundary grid model that forms fault boundaries and regional trend surfaces;

[0045] On the grid model of actual geological boundaries that is finally integrated with the generated stress field constraints, density field constraints are generated according to the Bouguer gravity anomaly data, and density anomalies are inverted through the potential field forward modeling formula to generate density anomaly fields.

[0046] Perform 3D voxelized conflict detection on the fault boundaries and regional trend surfaces, lithologic mutation zone morphology, stress field constraints, and density field constraints on the final fusion grid model of the actual geological boundary to obtain the conflict area;

[0047] Obtain the First Mediation Rules;

[0048] Each conflicting area is processed according to the first mediation rule, thereby obtaining a grid model after interpolation constraints.

[0049] Optionally, performing geometric correction on the conflicting area of ​​the interpolation-constrained grid model to obtain a corrected stratum interface grid includes:

[0050] Calculate the curvature fractal dimension and maximum principal stress of each grid cell in the grid model after interpolation constraint respectively;

[0051] Determine whether the curvature fractal dimension meets the curvature subdivision requirements or whether the maximum principal stress meets the principal stress subdivision requirements. If either one meets the requirements, subdivide the grid cells that meet the requirements to obtain a conformally corrected subdivided grid.

[0052] The morphology of the subdivided grid after iterative conformal correction is minimized by the conjugate gradient method, thereby obtaining the grid morphology after potential field optimization.

[0053] The stress field-driven mesh deformation of the mesh shape after potential field optimization is calculated by the finite element method, thereby obtaining the mesh shape after stress field response;

[0054] Perform 3D voxelized conflict detection on the grid shape after stress field response, identify voxels occupied by multiple geological bodies at the same time, and obtain conflict areas for geometric correction;

[0055] Obtain the second mediation rules;

[0056] The conflicting regions for geometric correction are processed by a second mediation rule to obtain a corrected bed interface grid.

[0057] Optionally, constraining the modified stratum interface grid according to geological rules to form simple three-dimensional stratum interface model data includes:

[0058] Obtaining a preset geological rule library, the geological rule library including stratum contact relationship, fault cutting sequence information, lithologic mutation zone constraint information, tectonic stress field constraint information and density field constraint information, drilling data constraint information, and trend surface analysis information;

[0059] Generate a 3D grid model integrating geological rules based on the preset geological rule library and the revised stratigraphic interface grid;

[0060] Generate a 3D grid model with key geological boundary constraints based on the 3D grid model integrating geological rules and the preset geological rule library;

[0061] Obtain a three-dimensional grid model of coupled stress and density fields based on a three-dimensional grid model with key geological boundary constraints and a preset geological rule library;

[0062] Conflict detection and mediation are performed on the three-dimensional grid model of the coupled stress field and density field to obtain simple three-dimensional stratum interface model data.

[0063] Optionally, the seismic profile information includes information such as fault traces, fold axis traces, and intrusion body boundaries;

[0064] The method of modifying the simple three-dimensional stratum interface model data according to the seismic profile information to obtain the composite geological structure model includes:

[0065] Mapping seismic profile information onto a simple three-dimensional stratigraphic interface model;

[0066] According to the geometric characteristics of seismic profile information, the simple three-dimensional stratum interface model is geometrically transformed to make it better match the seismic information;

[0067] The seismic attributes from the seismic profile information are integrated with the attributes of a simple 3D stratigraphic interface model to optimize the morphology and attributes of the concealed structure.

[0068] The non-layered geological body modeling is carried out by integrating the simple three-dimensional stratum interface model with the seismic profile information to obtain the composite geological structure model.

[0069] Optionally, the seismic profile information further includes a structured fault database, wherein the structured fault database includes geometric parameters and attributes of each fault;

[0070] Generating a three-dimensional fault surface model based on the composite geological structure model data includes:

[0071] On the composite geological structure model, using the fault trace as a control line, evenly distributed control points are generated along the fault strike;

[0072] The parabola interpolation algorithm is used to generate a smooth and continuous three-dimensional fault plane with the control line as the skeleton.

[0073] For complex faults, segmented fitting is performed in combination with the sudden change points of formation attitude to ensure that the fault plane and the formation interface intersect reasonably;

[0074] Map fault attributes to fault plane grid nodes to generate a three-dimensional fault plane model with each attribute attached;

[0075] Interactively verifying each three-dimensional fault surface model with attached attributes with the composite geological structure model, thereby obtaining each interactively verified three-dimensional fault surface model;

[0076] Define fault priorities for each interactively validated 3D fault surface model based on fault activity phase and scale;

[0077] According to the priority order, the Boolean difference operation is performed on each interactively verified 3D fault surface model to achieve spatial cutting between faults and obtain a 3D geological structure model containing the fault network;

[0078] Laplace smoothing is performed on the fault surface edge of the three-dimensional geological structure model containing the fault network to eliminate the jagged shape, thereby obtaining the final three-dimensional fault surface model.

[0079] Optionally, generating a composite geological structure model including a fault network based on the three-dimensional fault plane model and the composite geological structure model includes:

[0080] Aligning the three-dimensional fault surface model and the composite geological structure model in the same coordinate system, thereby obtaining the aligned three-dimensional fault surface model and composite geological structure model;

[0081] Generate 3D fault plane models with priority labels and topology rules;

[0082] Fault cutting is achieved through Boolean difference operation, thus generating a composite geological structure model containing fault cutting traces;

[0083] The fault surface edge of the composite geological structure model containing fault cutting traces is smoothed to obtain a geometrically optimized composite geological structure model containing a fault network;

[0084] The present application also provides a three-dimensional geological modeling device, which includes:

[0085] A data acquisition module, wherein the data acquisition module is used to acquire the geological area data to be established and the actual geological boundary data;

[0086] An initial stratum interface grid generation module, wherein the initial stratum interface grid generation module is used to generate an initial stratum interface grid according to the geological area data to be established;

[0087] a correction module, the correction module being used to correct the initial formation interface grid, thereby obtaining a corrected formation interface grid;

[0088] A geological rule constraint module, which is used to perform geological rule constraints on the modified stratum interface grid, thereby forming simple three-dimensional stratum interface model data;

[0089] A seismic profile correction module, which is used to correct the simple three-dimensional stratum interface model data according to the seismic profile information, thereby obtaining a composite geological structure model;

[0090] A three-dimensional fault generation module, wherein the three-dimensional fault generation module is used to generate a three-dimensional fault surface model according to the composite geological structure model;

[0091] A composite geological structure model generation module for a fault network, wherein the composite geological structure model generation module for a fault network is used to generate a composite geological structure model including a fault network according to the three-dimensional fault surface model and the composite geological structure model.

[0092] Beneficial effects

[0093] The 3D geological modeling method proposed in this application achieves dynamic adaptation from discrete data to a continuous model by coupling multiple parameter constraints, including seismic profile fault traces, tectonic stress fields, Bouguer gravity anomalies, and magnetic susceptibility anomalies. This 3D geological modeling method can significantly reduce the deviation between the model and actual geological boundaries. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 It is a flowchart of a three-dimensional geological modeling method according to an embodiment of the present application. DETAILED DESCRIPTION

[0095] In order to make the purpose, technical solutions and advantages of the implementation of this application clearer, the technical solutions in the embodiments of this application will be described in more detail below in conjunction with the drawings in the embodiments of this application. In the drawings, the same or similar reference numerals throughout represent the same or similar elements or elements with the same or similar functions. The described embodiments are part of the embodiments of this application, not all of the embodiments. The embodiments described below with reference to the drawings are exemplary and are intended to be used to explain this application, and should not be understood as limitations on this application. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application. The embodiments of this application are described in detail below in conjunction with the drawings.

[0096] like Figure 1 The three-dimensional geological model building method shown includes:

[0097] Obtain data on the geological area to be established and actual geological boundary data;

[0098] Generate an initial stratigraphic interface grid based on the geological area data to be established;

[0099] Modifying the initial stratum interface grid to obtain a modified stratum interface grid;

[0100] The modified stratigraphic interface grid is constrained by geological rules to form simple three-dimensional stratigraphic interface model data;

[0101] The simple three-dimensional stratum interface model data is modified according to the seismic profile information to obtain a composite geological structure model;

[0102] Generate a three-dimensional fault surface model based on the composite geological structure model;

[0103] A composite geological structure model including a fault network is generated according to the three-dimensional fault surface model and the composite geological structure model.

[0104] In this embodiment, the geological regional data to be established include discrete stratigraphic data (strike, dip, inclination), stratigraphic boundary control points (drilling depth), gravity anomaly data, magnetic anomaly data, geological outcrop profile data, stratigraphic age data, and lithologic data;

[0105] Generating an initial stratigraphic interface grid based on the geological area data to be established includes:

[0106] Gridding the geological area data to be established to form multiple grid data points;

[0107] Perform first-layer downsampling on each grid data point to obtain the information of each first-layer data point;

[0108] Performing second-layer downsampling on each first-layer data point information to obtain each second-layer data point information;

[0109] Performing third-layer downsampling on each second-layer data point information to obtain each third-layer data point information;

[0110] Performing fourth-layer downsampling on each third-layer data point information to obtain each fourth-layer data point information;

[0111] The geological rule coding information is associated with the first layer data point information, the second layer data point information, the third layer data point information and the fourth layer data point information respectively;

[0112] Using a quintic polynomial model to roughly fit the fourth layer data point information and the geological rule coding information associated with the fourth layer data point information, thereby obtaining the fourth layer polynomial coefficients;

[0113] Using a quadratic polynomial model to fit each first layer data point information and the geological rule coding information associated with each layer to obtain the first layer polynomial coefficients;

[0114] Using a quadratic polynomial model to fit each second layer data point information and the geological rule coding information associated with each layer to obtain the second layer polynomial coefficients;

[0115] Using a quadratic polynomial model to fit each third layer data point information and the geological rule coding information associated with each layer to obtain the third layer polynomial coefficients;

[0116] Calculate the construction complexity of each first-layer data point information, each second-layer data point information, and each third-layer data point information respectively;

[0117] When the construction complexity of any one of the first-layer data point information exceeds a preset threshold, the first-layer data point information exceeding the preset threshold is locally refitted through a quintic polynomial model to obtain a corrected first-layer polynomial coefficient;

[0118] When the construction complexity of any one of the second-layer data point information exceeds a preset threshold, the second-layer data point information exceeding the preset threshold is locally refitted through a quintic polynomial model to obtain a corrected second-layer polynomial coefficient;

[0119] When the construction complexity of any information in the third-layer data points exceeds a preset threshold, the third-layer data points exceeding the preset threshold are locally refitted through a quintic polynomial model to obtain a corrected third-layer polynomial coefficient;

[0120] A joint probability model is constructed based on the obtained corrected first-layer polynomial coefficients, corrected second-layer polynomial coefficients, corrected third-layer polynomial coefficients, fourth-layer polynomial coefficients, and various geological rule coding information;

[0121] Obtain attribute probability distribution of stratigraphic age and lithologic attributes;

[0122] Regular three-dimensional grid data is generated according to the optimized trend surface parameters. The regular three-dimensional grid data and the attribute probability distribution of the stratigraphic age and lithologic attributes constitute an initial stratigraphic interface grid.

[0123] In this embodiment, the geological area data to be established are first preprocessed. Specifically, the ordinary kriging method is used to perform density enhancement on discrete stratigraphic data (spacing > 5 km) to generate a spatially continuous occurrence field (dip ∠φ(x, y), dip angle ∠α(x, y)). Then, anomalies are eliminated. Specifically, occurrence mutation points are filtered based on the 3σ principle, and occurrence data that conform to the regional tectonic background are retained. Finally, coordinate unification is performed, and all data are converted to the CGCS2000 coordinate system, and the elevation is converted to the ellipsoid.

[0124] In this embodiment, a four-layer data pyramid is constructed to downsample layer by layer to generate data layers of different resolutions. This technology can capture geological features at different scales, from macroscopic structural trends to microscopic geological details, providing multi-scale data support for subsequent geological modeling.

[0125] In this embodiment, the four-layer data pyramid can be constructed in the following manner:

[0126] Each grid data point is downsampled to the first level to obtain the information of each first-level data point. Specifically, this level is the bottom layer of the data pyramid and contains the original geological data with the highest resolution. In this embodiment, the data resolution of this level may be set to Δx = 10m and Δy = 10m to capture the most subtle geological changes.

[0127] Starting from the original data layer, lower resolution data layers are generated layer by layer through downsampling techniques (such as mean downsampling, maximum downsampling, etc.).

[0128] For example, the remaining layers can take the following parameters:

[0129] The second layer is generated by downsampling, with the resolution reduced to Δx = 20m and Δy = 20m. This layer begins to capture slightly larger-scale geological features.

[0130] Layer 3: Downsampling continues, reducing the resolution to Δx = 40m and Δy = 40m. This layer is used to capture geological structural trends over a larger area.

[0131] The fourth layer is the final downsampling layer, with the resolution reduced to Δx = 80m and Δy = 80m. This layer is mainly used to capture regional macroscopic geological features.

[0132] In this embodiment, the data point information of each layer needs to be associated with geological rule encoding information (for example, the constraint that “the earlier formed layer is below”).

[0133] In this embodiment, a quintic polynomial model is used to roughly fit the fourth layer of data point information, while the second-order polynomial model can be used for fine adjustment in other layers.

[0134] In this embodiment, the construction complexity can be obtained by the following formula:

[0135] C represents structural complexity, which describes the complexity of the stratigraphic interface. When structural complexity is high, a higher-order polynomial model may be required for fitting. Z represents the trend surface prediction result (the fitted elevation value), which represents the spatial distribution of the stratigraphic interface. x and y represent spatial coordinates, used to locate points on the stratigraphic interface.

[0136] Denotes partial derivatives, which describe the rate of change of the bed interface in the x and y directions. By calculating these partial derivatives, the gradient information of the bed interface can be obtained, and then the structural complexity can be calculated. In this embodiment, when the structural complexity C obtained from any data point information is greater than τ (a threshold automatically determined by the Otsu algorithm), the quintic polynomial model is switched.

[0137] In this embodiment, regardless of whether it is a quadratic polynomial model or a quintic polynomial model, the Levenberg-Marquardt algorithm is used to solve the parameters.

[0138] In this embodiment, the joint probability model is as follows:

[0139]

[0140] Where p(T,L,Z|x) represents the joint probability distribution of stratigraphic age data (T), lithologic data (L), and trend surface prediction results (Z) at a given spatial location x.

[0141] T stands for stratigraphic age data, which describes the time information of stratum formation or deposition and is an important basis for dividing and comparing strata in geology.

[0142] L stands for lithologic data, which describes the type, composition, structure and other characteristics of rocks. It is the basic data used in geology to analyze and interpret geological phenomena.

[0143] Z represents the trend surface prediction result, that is, the elevation value after fitting, which represents the distribution pattern of the stratum interface in space.

[0144] x represents the spatial position or coordinate, which is used to locate points on the stratigraphic interface and is an important parameter in geology that describes the spatial distribution of geological phenomena.

[0145] It represents the predicted or expected value of Z, that is, the best estimate of Z given the spatial location x and other relevant variables (such as T and L).

[0146] σ Z Represents the predicted value The standard deviation of the predictions is a measure of the uncertainty of the forecast. A smaller standard deviation indicates a more reliable forecast; a larger standard deviation indicates greater uncertainty in the forecast.

[0147] p(T|x) represents the probability distribution of stratigraphic age data T given a given spatial location x. This conditional probability describes the distribution of stratigraphic age data at different spatial locations. It is typically solved using a spatially varying kriging model, which accounts for the correlation between data points in space and produces more accurate and reliable interpolation results.

[0148] p(L|T,x) represents the probability distribution of lithologic data L given a given spatial location x and stratigraphic age data T. This conditional probability describes the distribution of lithologic data under different spatial locations and stratigraphic ages.

[0149] In this embodiment, Gibbs sampling is used for joint inference to obtain the attribute probability distribution of stratigraphic age and lithologic attributes. Specifically, the following steps can be used:

[0150] Fixed T, L, updated Z (trend surface correction)

[0151] Fixed Z, L, updated T (stratigraphic age inversion)

[0152] Fix Z, T and update L (lithology probability propagation).

[0153] In this embodiment, the matrix operations in the trend surface fitting process are calculated in a parallel computing manner, including matrix-vector multiplication and Cholesky decomposition to solve the normal equations.

[0154] Specifically, by performing matrix operations in parallel through CUDA kernel functions, the computational domain is divided into tiles, and each tile can perform calculations independently, which helps to make full use of computing resources and improve computing efficiency.

[0155] While tiles are performing independent calculations, boundary conditions are synchronized through MPI to ensure data consistency between tiles and the accuracy of the overall calculation results.

[0156] Generate regular 3D grid data based on the results of parallel computing and save it in .grd format.

[0157] In this embodiment, the actual geological boundary data include fault traces, unconformity surface contours, and measured point data of lithologic contact zones in seismic profiles; the geological region data to be established further include tectonic stress field data, Bouguer gravity anomaly data, and magnetic susceptibility anomaly data;

[0158] In this embodiment, correcting the initial stratum interface grid to obtain a corrected stratum interface grid includes:

[0159] Correcting the initial stratum interface grid to obtain a corrected stratum interface grid includes:

[0160] Based on the initial stratigraphic interface grid, the actual geological boundary data is mapped to the initial stratigraphic interface grid surface through cubic spline interpolation, thereby obtaining the initial fused actual geological boundary grid model;

[0161] A parameterized constraint model is established based on the actual geological boundary data, thereby correcting the grid node coordinates in the initial grid model integrated with the actual geological boundary, thereby obtaining a final grid model integrated with the actual geological boundary;

[0162] The key geological boundary constraints are interpolated on the final grid model integrated with the actual geological boundary through the geological rule library to obtain the grid model after interpolation constraints;

[0163] The conflicting area geometry is corrected on the interpolation constrained grid model to obtain the corrected stratum interface grid.

[0164] In this embodiment, the initial stratum interface grid is used as a base, and the actual geological boundary data is mapped to the initial stratum interface grid surface through cubic spline interpolation, thereby obtaining an initial fused actual geological boundary grid model, which includes:

[0165] Through the data credibility assessment model, dynamic weights are assigned to borehole data, seismic profile data, etc. The weighted least squares method is used to fit the trend surface, integrating the tectonic stress field and potential field data to obtain the initial fusion of the actual geological boundary grid model;

[0166] In this embodiment, the data credibility evaluation model adopts the following formula:

[0167]

[0168] Where W represents the weighted average result of data credibility assessment. This value is used to quantitatively evaluate the credibility of the fused data and is an important indicator of the final output of the model; λ borehole Represents the weight parameter of the well data. Represents the wellbore data index. seismic Represents the weight parameter of earthquake data. Represents earthquake data indicators.

[0169] The method of establishing a parameterized constraint model based on actual geological boundary data, thereby correcting the grid node coordinates in the initial grid model that integrates the actual geological boundary, and thus obtaining the final integrated actual geological boundary grid model, includes: based on the initial grid model that integrates multi-physics field constraints, topologically simplifying the fault traces using actual geological boundary data (such as fault traces and measured point data of lithologic contact zones), and retaining key turning points; establishing a lithologic mutation zone buffer zone model, and generating the final integrated actual geological boundary grid model (based on the initial integrated actual geological boundary grid model, adding geometric constraints of the fault traces and lithologic mutation zones);

[0170] In this embodiment, the key geological boundary constraints are interpolated on the final fused actual geological boundary grid model through the geological rule library, thereby obtaining the grid model after interpolation constraints, including:

[0171] Arrange control points along the fault traces of the final fused actual geological boundary grid model, and use cubic B-spline curves to generate continuous fault boundaries, thereby forming an initial fault geometry framework in the final fused actual geological boundary grid model (the initial fault geometry framework includes B-spline curve control points and interpolation results);

[0172] Based on the initial fault geometric framework and the final fused actual geological boundary grid model, outside the fault influence domain, with the fault boundary generated by B-spline as a constraint, a regional trend surface is generated through RBF smooth transition, thereby forming the fault boundary and regional trend surface on the final fused actual geological boundary grid model that forms the initial fault geometric framework (the radial basis function (RBF) interpolation result includes the fault boundary coordinates and the regional trend surface equation);

[0173] Using a Markov random field (MRF) model, a lithologic probability field is established based on the final fusion actual geological boundary grid model and the measured point data of the lithologic contact zone in the actual geological boundary data, and single-point potential energy and adjacent node potential energy are defined, thereby obtaining a lithologic probability field model (including lithologic probability and potential energy parameters of each node);

[0174] Based on the lithologic probability field model and the measured point data of the lithologic contact zone, gradient constrained interpolation is performed along the normal direction of the lithologic mutation zone to obtain the morphology of the lithologic mutation zone (the coordinates and properties of the lithologic interface after interpolation);

[0175] On the final fusion actual geological boundary grid model that forms the fault boundary and regional trend surface, stress field constraints (a set of grid node displacement vectors) are generated based on the tectonic stress field data. Specifically, the maximum principal stress direction is mapped to the grid node displacement vector and used as the interpolation constraint.

[0176] On the grid model of the actual geological boundary that is finally integrated with the generated stress field constraints, density field constraints (density anomaly field data) are generated based on the Bouguer gravity anomaly data. Specifically, density anomalies are inverted using the potential field forward modeling formula to generate a density anomaly field.

[0177] Perform 3D voxelized conflict detection on the fault boundaries and regional trend surfaces, lithologic mutation zone morphology, stress field constraints, and density field constraints on the final grid model that integrates the actual geological boundaries. Specifically, the grid model is voxelized and voxels simultaneously occupied by multiple geological bodies are detected to obtain conflicting areas.

[0178] Obtain a first mediation rule. In this embodiment, the first mediation rule can be obtained using the following method:

[0179] Define the conflict membership μ(x), which ranges from [0,1]. For example:

[0180] Where d is the Euclidean distance from the conflicting voxel to the nearest legal position, and k is the attenuation coefficient.

[0181] The first mediation rule can be set as needed, for example, if μ(x)>0.8, the stress field constraint is enforced.

[0182] If 0.5<μ(x)≤0.8, start the weighted average correction (weight = 0.6 stress field, 0.4 trend surface).

[0183] If μ(x)≤0.5, the original trend surface shape is retained.

[0184] Each conflicting area is processed according to the first mediation rule, thereby obtaining a grid model after interpolation constraints.

[0185] In this embodiment, performing geometric correction on the conflicting area of ​​the interpolation constrained grid model to obtain a corrected stratum interface grid includes:

[0186] Calculate the curvature fractal dimension and maximum principal stress of each grid cell in the grid model after interpolation constraint respectively;

[0187] Determine whether the curvature fractal dimension meets the curvature subdivision requirements or whether the maximum principal stress meets the principal stress subdivision requirements. If either one meets the requirements, subdivide the grid cells that meet the requirements to obtain a conformally corrected subdivided grid.

[0188] The morphology of the subdivided grid after iterative conformal correction is minimized by the conjugate gradient method, thereby obtaining the grid morphology after potential field optimization.

[0189] The stress field-driven mesh deformation of the mesh shape after potential field optimization is calculated by the finite element method, thereby obtaining the mesh shape after stress field response;

[0190] Perform 3D voxelized conflict detection on the grid shape after stress field response, identify voxels occupied by multiple geological bodies at the same time, and obtain conflict areas for geometric correction;

[0191] Obtain the second mediation rules;

[0192] The conflicting regions for geometric correction are processed by a second mediation rule to obtain a corrected bed interface grid.

[0193] In this embodiment, the modified stratum interface grid is constrained by geological rules to form a simple three-dimensional stratum interface model data including:

[0194] Obtaining a preset geological rule library, the geological rule library including stratum contact relationship, fault cutting sequence information, lithologic mutation zone constraint information, tectonic stress field constraint information and density field constraint information, drilling data constraint information, and trend surface analysis information;

[0195] A 3D grid model incorporating geological rules is generated based on a pre-set geological rule library and a modified stratigraphic interface grid. Specifically, radial basis function interpolation is used with borehole data constraints as hard constraints to forcibly correct grid node elevations, ensuring that the top and bottom depths and lithology of the stratigraphic layers revealed by the boreholes are accurately reflected in the model.

[0196] The radial basis function interpolation method uses the trend surface analysis results as soft constraints and generates a continuous three-dimensional interface through the RBF function, so that the overall trend of the model is consistent with the regional stratigraphic trend, while smoothly transitioning the regional geological characteristics and reducing the impact of data noise.

[0197] According to the contact relationship of the strata, the grid node elevation is forced to be adjusted to ensure that the old strata are located below the new strata; for the unconformity surface, the "truncation-overlap" algorithm is used to correct the stratum interface morphology and eliminate the stratum time-diachrony phenomenon;

[0198] According to the fault cutting sequence in the geological rule library, the elevation of the stratum interfaces on both sides of the fault is adjusted to ensure that the main fault preferentially cuts the secondary fault and that the fault activity period matches the stratum deposition time sequence. For the conjugate fault system, the voxel carving algorithm is used to ensure the topological relationship between the faults is reasonable and avoid unreasonable fault intersections, thereby obtaining a 3D grid model that integrates geological rules.

[0199] A 3D grid model with key geological boundary constraints is generated based on the 3D grid model that integrates geological rules and a preset geological rule library. Specifically, control points are placed along the fault traces in the actual geological boundary data in the 3D grid model that integrates geological rules, and a continuous fault boundary is generated using a cubic B-spline curve to ensure that the fault morphology is consistent with the actual interpretation.

[0200] Within the fault influence domain, using B-spline curves as constraints, RBF interpolation is used to generate the stratum interfaces on both sides of the fault, maintaining the spatial consistency between the fault plane and the stratum interface, and ensuring that the fault cutting relationship conforms to the fault cutting sequence in the geological rule library;

[0201] Based on the measured point data of the lithologic contact zone in the actual geological boundary data and the lithologic mutation zone constraint rules in the geological rule library, a lithologic probability field is established, the single point potential energy and the adjacent node potential energy are defined, and the possibility of lithologic mutation is quantified.

[0202] Gradient constrained interpolation is performed along the normal direction of the lithologic abrupt zone to generate a smooth lithologic interface transition zone, avoid the random diffusion of the lithologic abrupt zone, and ensure that the lithologic contact relationship conforms to the stratigraphic contact relationship in the geological rule library, thereby obtaining a 3D grid model with key geological boundary constraints.

[0203] A three-dimensional grid model of the coupled stress field and density field is obtained based on the three-dimensional grid model with key geological boundary constraints and a preset geological rule library. Specifically, the three-dimensional grid model with key geological boundary constraints is processed as follows to obtain the three-dimensional grid model of the coupled stress field and density field:

[0204] The tectonic stress field data is mapped into grid node displacement vectors. Based on the tectonic stress field constraint rules in the geological rule base, the grid morphology is iteratively optimized using the conjugate gradient method to minimize the stress field-driven deformation residual and ensure that the grid deformation conforms to the tectonic stress field distribution and the fault dislocation relationship rules in the geological rule base.

[0205] The Bouguer gravity anomaly data is inverted through the potential field forward modeling formula to generate the density anomaly field and quantify the formation density distribution;

[0206] Taking density anomaly as an additional constraint, according to the density field constraint rules in the geological rule library, the stratigraphic interface elevation is corrected to match the gravity anomaly distribution, avoiding model distortion caused by density anomaly and ensuring that the overall density distribution of the model is consistent with geological cognition;

[0207] Detect and mediate conflicts between the three-dimensional grid models of the coupled stress field and density field to obtain simple three-dimensional stratum interface model data;

[0208] Detect and mediate conflicts in 3D mesh models of coupled stress and density fields to obtain simple 3D stratigraphic interface model data. For example, the mesh model is voxelized to detect voxels simultaneously occupied by multiple geological bodies (e.g., faults, lithologic bodies), identify conflicting areas, and quantify the degree of conflict (e.g., fault cutting priority, consistency of normal orientation of lithologic contact zones).

[0209] In this embodiment, multiple mediation rules can also be set. For example, one mediation rule may prioritize preserving the spatial footprint of high-priority geological bodies (e.g., primary faults) and adjusting the interface morphology of low-priority geological bodies (e.g., secondary faults) to eliminate unreasonable fault segmentation. Another mediation rule may also be used to calculate the stress field response of the geometrically corrected conflicting regions using the finite element method to ensure that the mesh deformation conforms to the tectonic stress field distribution and avoid introducing new model distortion during the mediation process.

[0210] In this embodiment, the seismic profile information includes information such as fault traces, fold axis traces, and intrusion boundaries;

[0211] The method of modifying the simple three-dimensional stratum interface model data according to the seismic profile information to obtain the composite geological structure model includes:

[0212] Map the seismic profile information onto a simple 3D stratigraphic model. This step ensures the accurate location of the seismic information in 3D space.

[0213] According to the geometric characteristics of seismic profile information, the simple three-dimensional stratum interface model is geometrically transformed (such as rotation, translation, scaling, etc.) to make it better match the seismic information;

[0214] The seismic attributes from the seismic profile information are integrated with the attributes of a simple 3D stratigraphic interface model to optimize the morphology and attributes of the concealed structure.

[0215] The non-layered geological body modeling is carried out by integrating the simple three-dimensional stratum interface model with the seismic profile information to obtain the composite geological structure model.

[0216] In this embodiment, the non-layered geological body modeling is performed by fusing a simple three-dimensional stratum interface model with seismic profile information to obtain a composite geological structure model, including:

[0217] Contour drawing: For non-layered geological bodies (such as folds, salt domes, intrusions, etc.) revealed in seismic profiles, use 3D mesh editing tools to manually draw their contours.

[0218] 3D entity generation: Through lofting and other technologies, 3D entity models of these non-layered geological bodies are generated according to the outlined contour lines.

[0219] Model integration: Integrate the generated 3D solid model with the simple 3D stratum interface model to ensure their correct spatial positions and relationships.

[0220] Morphological optimization: Combined with seismic attribute data, the morphology of hidden structures (such as unexposed stratigraphic interfaces and hidden faults) is optimized. The modeling accuracy of hidden structures is improved by adjusting model parameters and applying geological constraints.

[0221] Attribute refinement: Based on seismic attribute information, the attributes of the hidden structure are refined, such as lithology, density, velocity, etc., to make it more consistent with the actual geological conditions.

[0222] In this embodiment, the seismic profile information further includes a structural fault database, which includes geometric parameters (trace coordinates, dip, fault throw) and attributes (properties, activity period) of each fault;

[0223] Generating a three-dimensional fault surface model based on the composite geological structure model data includes:

[0224] On the composite geological structure model, the fault trace is used as the control line, and control points are evenly distributed along the fault strike (the spacing is dynamically adjusted according to the complexity of the fault);

[0225] The parabola interpolation algorithm is used to generate a smooth and continuous three-dimensional fault plane with the control line as the skeleton.

[0226] For complex faults (such as shovel-shaped faults), segmented fitting is performed in combination with the sudden change points of the formation attitude to ensure that the fault plane and the formation interface intersect reasonably;

[0227] Map fault attributes (such as fault throw and activity period) to fault plane grid nodes to generate a three-dimensional fault plane model with various attached attributes;

[0228] Interactively verifying each three-dimensional fault surface model with attached attributes with the composite geological structure model, thereby obtaining each interactively verified three-dimensional fault surface model;

[0229] In this embodiment, interactive verification is performed on each of the three-dimensional fault surface models with the attached attributes with the composite geological structure model, thereby obtaining each interactively verified three-dimensional fault surface model, including:

[0230] Spatial intersection analysis: Detect the intersection line between the fault plane and the stratum interface to verify whether the fault cutting relationship conforms to geological rules (such as the main fault preferentially cuts the secondary fault).

[0231] Morphological optimization: For unreasonable intersection areas (such as fault planes penetrating strata that should not be cut), local corrections are made by adjusting fault plane control points or stratum interface nodes.

[0232] Conflict mediation: Identify voxels occupied by both fault planes and stratigraphic interfaces, and resolve conflicts based on pre-set mediation rules (e.g., prioritizing fault activity phases).

[0233] After completing the above three steps, the interactively verified three-dimensional fault surface model can be obtained.

[0234] Define fault priorities for each interactively validated 3D fault plane model based on the fault activity phase and scale (e.g., primary faults prioritize secondary faults).

[0235] Boolean subtraction operations are performed on the three-dimensional fault surface models after interactive verification in order of priority to achieve spatial cutting between faults and obtain a three-dimensional geological structure model containing the fault network. For conjugate faults or intersecting faults, a voxel carving algorithm is used to accurately depict the interaction relationship between faults (such as fault displacement transmission);

[0236] Laplace smoothing is performed on the fault surface edge of the three-dimensional geological structure model containing the fault network to eliminate the jagged shape, thereby obtaining the final three-dimensional fault surface model.

[0237] In this embodiment, each data in the structured fault database can be obtained in the following manner:

[0238] Fault trace extraction:

[0239] Fault trace coordinates are extracted from seismic profile interpretation data, and the geometric characteristics of each fault (such as strike, dip, and inclination) are recorded.

[0240] Spatial cluster analysis of fault traces was performed to identify fault combination patterns (e.g., graben, horst, and conjugate faults).

[0241] Drilling data verification:

[0242] Combined with drilling fault point data, the consistency of the spatial position and attributes (such as fault distance and activity period) of the fault trace is verified.

[0243] Conflicting data (e.g., fault traces that do not match the missing stratigraphic layers revealed by drilling) are flagged for subsequent reconciliation rule processing.

[0244] In this embodiment, generating a composite geological structure model including a fault network based on the three-dimensional fault plane model and the composite geological structure model includes:

[0245] Aligning the three-dimensional fault surface model and the composite geological structure model in the same coordinate system, thereby obtaining the aligned three-dimensional fault surface model and composite geological structure model;

[0246] Generate a 3D fault plane model with priority labels and topological rules. In this example, priority is defined based on fault activity phase and scale (e.g., throw and extension), with primary faults taking precedence over secondary faults. For conjugate fault systems, the interaction relationships between faults (e.g., displacement transmission direction) are labeled.

[0247] Fault cutting is achieved through a Boolean difference operation, thereby generating a composite geological structure model containing fault cutting traces. The method of achieving fault cutting through a Boolean difference operation to generate a composite geological structure model containing fault cutting traces includes: traversing the faults from high to low priority: for each fault, using a Boolean difference operation to cut the strata and non-stratified geological bodies in the composite geological structure model. The fragmented geological bodies produced by the fault cutting and their attributes (such as the lithology of the cut strata) are recorded. For conjugate faults or intersecting faults, a voxel carving algorithm is used to detect voxel occupancy conflicts in the fault intersection area. Voxel attribution is adjusted according to the fault displacement vector to ensure consistent displacement transfer logic between faults.

[0248] Smoothing the fault surface edges of a composite geological structure model containing fault cuts to obtain a geometrically optimized composite geological structure model containing a fault network includes: performing Laplace smoothing on the fault surface edges: iteratively optimizing the coordinates of the fault surface grid vertices to eliminate jagged shapes; maintaining the geometric characteristics of the fault trace (e.g., fault dip and strike); and constraining the smoothing process with stratigraphic strike data to ensure that the smoothed fault surface is consistent with the trend of the adjacent stratigraphic interface.

[0249] This application has the following advantages:

[0250] By coupling multiple parameter constraints, such as seismic profile fault traces, tectonic stress fields, Bouguer gravity anomalies, and magnetic susceptibility anomalies, dynamic adaptation from discrete data to a continuous model is achieved. For example, during the grid generation phase, cubic spline interpolation is used to map actual geological boundaries (such as fault traces and lithologic abrupt zones) onto the initial grid. Combined with RBF smoothing, this improves local accuracy to sub-meter levels, significantly reducing the deviation between the model and actual geological boundaries.

[0251] By converting geological rules such as stratigraphic contact relationships and fault cutting sequences into grid deformation constraints (such as stress field displacement vectors and density anomaly fields), geological contradictions can be dynamically corrected during the modeling process. For example, the conjugate gradient method is used to iteratively optimize the grid morphology to minimize stress-driven deformation residuals, ensuring that fault slip relationships are consistent with the regional stress field and avoiding the time-diachronic phenomenon of "old strata covering new strata" in traditional methods.

[0252] Based on a structured fault database, fault priorities and topological rules are defined, and a voxel carving algorithm is used to accurately characterize the displacement transfer relationship of conjugate faults. For example, for complex faults (such as shovel-shaped faults), segmented fitting is performed based on the sudden change points of formation attitude to ensure that the fault plane and the formation interface intersect properly. For conjugate fault systems, a voxel carving algorithm is used to detect voxel occupancy conflicts in fault intersection areas and adjust voxel ownership based on the fault displacement vector to avoid unreasonable fault intersections.

[0253] A conflict detection and quantitative mediation process is established for mesh subdivision, stress field optimization, and density field constraints. For example, mesh subdivision requirements are determined using the curvature fractal dimension and principal stress subdivision index. Membership functions are defined for conflicting areas and mediation rules are implemented, reducing manual intervention and improving overall model quality.

[0254] The 3D mesh editing tool is used to manually outline the contours of non-stratified geological bodies (such as folds, salt domes, and intrusions), and then combined with lofting techniques to generate a 3D solid model. For example, the fold axis trace revealed by a seismic profile is outlined using the 3D mesh editing tool, and then a 3D solid fold is generated through lofting. This is then integrated with the stratigraphic interface model to ensure the correct spatial position.

[0255] Deeply integrate the geological rule library with the modeling process to achieve rule-driven model optimization. For example, during the grid generation phase, using drillhole data as hard constraints and trend surfaces as soft constraints, a continuous 3D interface is generated through RBF interpolation. During the fault modeling phase, priorities are defined based on the stage and scale of fault activity to ensure that primary faults prioritize cutting secondary faults.

[0256] Bidirectional data-model adaptation is achieved through cubic spline interpolation and RBF smoothing, improving local accuracy while maintaining global trend consistency. For example, key geological boundaries such as fault edges and lithologic abrupt zones are interpolated with high precision, while RBF interpolation ensures that the overall model trend is consistent with regional stratigraphic trends.

[0257] Algorithms automate processes such as data fusion, rule constraints, and conflict resolution, reducing manual intervention. For example, during the mesh subdivision stage, the curvature fractal dimension and principal stress subdivision indicators are used to automatically determine subdivision requirements, avoiding the limitations of traditional methods that rely on manual experience.

[0258] Although the present invention has been described in detail above using general descriptions and specific embodiments, it will be apparent to those skilled in the art that modifications and improvements may be made based on the present invention. Therefore, such modifications and improvements, which do not depart from the spirit of the present invention, are intended to be within the scope of protection claimed herein.

Claims

1. A three-dimensional geological modeling method, characterized in that: The three-dimensional geological modeling method comprises: Obtain data on the geological area to be established and actual geological boundary data; Generate an initial stratigraphic interface grid based on the geological area data to be established; Modifying the initial stratum interface grid to obtain a modified stratum interface grid; The modified stratigraphic interface grid is constrained by geological rules to form simple three-dimensional stratigraphic interface model data; The simple three-dimensional stratum interface model data is modified according to the seismic profile information to obtain a composite geological structure model; Generate a three-dimensional fault surface model based on the composite geological structure model; A composite geological structure model including a fault network is generated according to the three-dimensional fault surface model and the composite geological structure model.

2. The three-dimensional geological modeling method according to claim 1, wherein: The geological regional data to be established include discrete stratigraphic data, stratigraphic boundary control points, gravity anomaly data, magnetic anomaly data, geological outcrop profile data, stratigraphic age data and lithologic data; Generating an initial stratigraphic interface grid based on the geological area data to be established includes: Gridding the geological area data to be established to form multiple grid data points; Perform first-layer downsampling on each grid data point to obtain the information of each first-layer data point; Performing second-layer downsampling on each first-layer data point information to obtain each second-layer data point information; Performing third-layer downsampling on each second-layer data point information to obtain each third-layer data point information; Performing fourth-layer downsampling on each third-layer data point information to obtain each fourth-layer data point information; The geological rule coding information is associated with the first layer data point information, the second layer data point information, the third layer data point information and the fourth layer data point information respectively; Using a quintic polynomial model to roughly fit the fourth layer data point information and the geological rule coding information associated with the fourth layer data point information, thereby obtaining the fourth layer polynomial coefficients; Using a quadratic polynomial model to fit each first layer data point information and the geological rule coding information associated with each layer to obtain the first layer polynomial coefficients; Using a quadratic polynomial model to fit each second layer data point information and the geological rule coding information associated with each layer to obtain the second layer polynomial coefficients; Using a quadratic polynomial model to fit each third layer data point information and the geological rule coding information associated with each layer to obtain the third layer polynomial coefficients; Calculate the construction complexity of each first-layer data point information, each second-layer data point information, and each third-layer data point information respectively; When the construction complexity of any one of the first-layer data point information exceeds a preset threshold, the first-layer data point information exceeding the preset threshold is locally refitted through a quintic polynomial model to obtain a corrected first-layer polynomial coefficient; When the construction complexity of any one of the second-layer data point information exceeds a preset threshold, the second-layer data point information exceeding the preset threshold is locally refitted through a quintic polynomial model to obtain a corrected second-layer polynomial coefficient; When the construction complexity of any information in the third-layer data points exceeds a preset threshold, the third-layer data points exceeding the preset threshold are locally refitted through a quintic polynomial model to obtain a corrected third-layer polynomial coefficient; A joint probability model is constructed based on the obtained corrected first-layer polynomial coefficients, corrected second-layer polynomial coefficients, corrected third-layer polynomial coefficients, fourth-layer polynomial coefficients, and various geological rule coding information; Obtain attribute probability distribution of stratigraphic age and lithologic attributes; Regular three-dimensional grid data is generated according to the optimized trend surface parameters. The regular three-dimensional grid data and the attribute probability distribution of the stratigraphic age and lithologic attributes constitute an initial stratigraphic interface grid.

3. The three-dimensional geological modeling method according to claim 2, wherein: The actual geological boundary data include fault traces, unconformity surface contours, and measured point data of lithologic contact zones interpreted from seismic profiles; the geological area data to be established further include tectonic stress field data, Bouguer gravity anomaly data, and magnetic susceptibility anomaly data; The step of correcting the initial stratum interface grid to obtain a corrected stratum interface grid includes: Based on the initial stratigraphic interface grid, the actual geological boundary data is mapped to the initial stratigraphic interface grid surface through cubic spline interpolation, thereby obtaining the initial fused actual geological boundary grid model; A parameterized constraint model is established based on the actual geological boundary data, thereby correcting the grid node coordinates in the initial grid model integrated with the actual geological boundary, thereby obtaining a final grid model integrated with the actual geological boundary; The key geological boundary constraints are interpolated on the final grid model integrated with the actual geological boundary through the geological rule library to obtain the grid model after interpolation constraints; The conflicting area geometry is corrected on the interpolation constrained grid model to obtain the corrected stratum interface grid.

4. The three-dimensional geological modeling method according to claim 3, wherein: The method of performing key geological boundary constraint interpolation on the final fused actual geological boundary grid model through the geological rule library to obtain the grid model after interpolation constraint includes: Control points are arranged along the fault traces of the final fused actual geological boundary grid model, and a continuous fault boundary is generated using a cubic B-spline curve, thereby forming an initial fault geometry framework in the final fused actual geological boundary grid model; Based on the initial fault geometric framework and the final fused actual geological boundary grid model, outside the fault influence domain, the fault boundary generated by the B-spline is used as a constraint, and a regional trend surface is generated through RBF smooth transition, thereby forming the fault boundary and the regional trend surface on the final fused actual geological boundary grid model that forms the initial fault geometric framework; Establishing a lithologic probability field based on the final fusion actual geological boundary grid model and the measured point data of the lithologic contact zone in the actual geological boundary data, defining single point potential energy and adjacent node potential energy, thereby obtaining a lithologic probability field model; Based on the lithologic probability field model and the measured point data of the lithologic contact zone, gradient constrained interpolation is performed along the normal direction of the lithologic mutation zone to obtain the morphology of the lithologic mutation zone; Generate stress field constraints based on tectonic stress field data on the final fused actual geological boundary grid model that forms fault boundaries and regional trend surfaces; On the grid model of actual geological boundaries that is finally integrated with the generated stress field constraints, density field constraints are generated according to the Bouguer gravity anomaly data, and density anomalies are inverted through the potential field forward modeling formula to generate density anomaly fields. Perform 3D voxelized conflict detection on the fault boundaries and regional trend surfaces, lithologic mutation zone morphology, stress field constraints, and density field constraints on the final fusion grid model of the actual geological boundary to obtain the conflict area; Obtain the First Mediation Rules; Each conflicting area is processed according to the first mediation rule, thereby obtaining a grid model after interpolation constraints.

5. The three-dimensional geological modeling method according to claim 4, wherein: The step of geometrically correcting the conflicting area of ​​the interpolation-constrained grid model to obtain a corrected stratum interface grid includes: Calculate the curvature fractal dimension and maximum principal stress of each grid cell in the grid model after interpolation constraint respectively; Determine whether the curvature fractal dimension meets the curvature subdivision requirements or whether the maximum principal stress meets the principal stress subdivision requirements. If either one meets the requirements, subdivide the grid cells that meet the requirements to obtain a conformally corrected subdivided grid. The morphology of the subdivided grid after iterative conformal correction is minimized by the conjugate gradient method, thereby obtaining the grid morphology after potential field optimization. The stress field-driven mesh deformation of the mesh shape after potential field optimization is calculated by the finite element method, thereby obtaining the mesh shape after stress field response; Perform 3D voxelized conflict detection on the grid shape after stress field response, identify voxels occupied by multiple geological bodies at the same time, and obtain conflict areas for geometric correction; Obtain the second mediation rules; The conflicting regions for geometric correction are processed by a second mediation rule to obtain a corrected bed interface grid.

6. The three-dimensional geological modeling method according to claim 5, wherein: The step of constraining the modified stratum interface grid with geological rules to form simple three-dimensional stratum interface model data includes: Obtaining a preset geological rule library, the geological rule library including stratum contact relationship, fault cutting sequence information, lithologic mutation zone constraint information, tectonic stress field constraint information and density field constraint information, drilling data constraint information, and trend surface analysis information; Generate a 3D grid model integrating geological rules based on the preset geological rule library and the revised stratigraphic interface grid; Generate a 3D grid model with key geological boundary constraints based on the 3D grid model integrating geological rules and the preset geological rule library; Obtain a three-dimensional grid model of coupled stress and density fields based on a three-dimensional grid model with key geological boundary constraints and a preset geological rule library; Conflict detection and mediation are performed on the three-dimensional grid model of the coupled stress field and density field to obtain simple three-dimensional stratum interface model data.

7. The three-dimensional geological modeling method according to claim 6, wherein: The seismic profile information includes fault traces, fold axis traces, intrusion body boundaries and other information; The method of modifying the simple three-dimensional stratum interface model data according to the seismic profile information to obtain the composite geological structure model includes: Mapping seismic profile information onto a simple three-dimensional stratigraphic interface model; According to the geometric characteristics of seismic profile information, the simple three-dimensional stratum interface model is geometrically transformed to make it better match the seismic information; The seismic attributes from the seismic profile information are integrated with the attributes of a simple 3D stratigraphic interface model to optimize the morphology and attributes of the concealed structure. The non-layered geological body modeling is carried out by integrating the simple three-dimensional stratum interface model with the seismic profile information to obtain the composite geological structure model.

8. The three-dimensional geological modeling method according to claim 7, wherein: The seismic profile information further includes a structural fault database, wherein the structural fault database includes geometric parameters and attributes of each fault; Generating a three-dimensional fault surface model based on the composite geological structure model data includes: On the composite geological structure model, using the fault trace as a control line, evenly distributed control points are generated along the fault strike; The parabola interpolation algorithm is used to generate a smooth and continuous three-dimensional fault plane with the control line as the skeleton. For complex faults, segmented fitting is performed in combination with the sudden change points of formation attitude to ensure that the fault plane and the formation interface intersect reasonably; Map fault attributes to fault plane grid nodes to generate a three-dimensional fault plane model with each attribute attached; Interactively verifying each three-dimensional fault surface model with attached attributes with the composite geological structure model, thereby obtaining each interactively verified three-dimensional fault surface model; Define fault priorities for each interactively validated 3D fault surface model based on fault activity phase and scale; According to the priority order, the Boolean difference operation is performed on each interactively verified 3D fault surface model to achieve spatial cutting between faults and obtain a 3D geological structure model containing the fault network; Laplace smoothing is performed on the fault surface edge of the three-dimensional geological structure model containing the fault network to eliminate the jagged shape, thereby obtaining the final three-dimensional fault surface model.

9. The three-dimensional geological modeling method according to claim 8, wherein: Generating a composite geological structure model including a fault network based on the three-dimensional fault plane model and the composite geological structure model includes: Aligning the three-dimensional fault surface model and the composite geological structure model in the same coordinate system, thereby obtaining the aligned three-dimensional fault surface model and composite geological structure model; Generate 3D fault plane models with priority labels and topology rules; Fault cutting is achieved through Boolean difference operation, thus generating a composite geological structure model containing fault cutting traces; The fault surface edge of the composite geological structure model containing fault cutting traces is smoothed to obtain a geometrically optimized composite geological structure model containing a fault network; 10. A three-dimensional geological modeling device, characterized in that: The three-dimensional geological modeling device comprises: A data acquisition module, wherein the data acquisition module is used to acquire the geological area data to be established and the actual geological boundary data; An initial stratum interface grid generation module, wherein the initial stratum interface grid generation module is used to generate an initial stratum interface grid according to the geological area data to be established; a correction module, the correction module being used to correct the initial formation interface grid, thereby obtaining a corrected formation interface grid; A geological rule constraint module, which is used to perform geological rule constraints on the modified stratum interface grid, thereby forming simple three-dimensional stratum interface model data; A seismic profile correction module, which is used to correct the simple three-dimensional stratum interface model data according to the seismic profile information, thereby obtaining a composite geological structure model; A three-dimensional fault generation module, wherein the three-dimensional fault generation module is used to generate a three-dimensional fault surface model according to the composite geological structure model; A composite geological structure model generation module for a fault network is used to generate a composite geological structure model including a fault network based on the three-dimensional fault surface model and the composite geological structure model.

Images