A fine three-dimensional geological modeling method based on geophysical three-dimensional forward and inversion algorithm

Abstract

This invention relates to a refined 3D geological modeling method based on geophysical 3D forward and inverse algorithms, belonging to the field of geological exploration. The method includes: acquiring multi-source data and constructing an initial model; jointly constraining a 3D inversion algorithm to obtain a 3D physical property model; creating a lithology classifier and generating a 3D lithology probability volume based on the 3D physical property model, followed by geological structure boundary enhancement and geological classification based on the probability volume; extracting isosurfaces from the 3D lithology probability volume, generating preliminary geological interfaces, performing structural processing, and then coupling attribute models to output the model. This invention directly constructs a spatially continuous model through 3D joint inversion, avoiding interpolation errors between survey lines; forced matching of borehole data + trend control of geophysical data balances local accuracy and global rationality; the entire process from inversion to lithology classification is algorithm-driven, reducing manual intervention; and it is suitable for high-precision scenarios such as urban underground space, mineral exploration, and geological hazard assessment.

Description

Technical Field

[0001] This invention relates to the field of geological exploration, and in particular to a refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms. Background Technology

[0002] With the continuous development of the construction and application of 3D geophysical engineering in China, refined geological models, as the foundation of Digital China, are becoming increasingly important. Traditional geological modeling generally relies on drilling data, but drilling data is point-based, and modeling between boreholes uses inter-hole difference methods, resulting in poor accuracy, especially when boreholes are sparse or the distance between boreholes is large, leading to severe model distortion. Existing technologies can use lithological stratification data and physical property data obtained from drilling as constraints, repeatedly performing two-dimensional forward and inverse iterative calculations to continuously reduce the difference between the surface physical response and geophysical measured data, thereby obtaining an optimal subsurface model that satisfies linear unbiased estimation and improving the local accuracy of geological modeling. However, to ensure exploration resolution, geophysical exploration generally involves linear measurements, not isometric measurements. Therefore, it can only obtain a locally optimal subsurface model along the direction of the geophysical survey lines. If inter-line modeling is used in the areas between survey lines, the model will also face the problem of distortion. Summary of the Invention

[0003] The purpose of this invention is to overcome the shortcomings of the prior art and provide a refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms, thus solving the deficiencies of the prior art.

[0004] The objective of this invention is achieved through the following technical solution: a refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms, the method comprising:

[0005] Inversion calculation steps: acquire multi-source data and construct an initial model, and use a joint constrained 3D inversion algorithm to obtain a 3D physical property model;

[0006] Geological classification steps: Create a lithology classifier, generate a 3D lithology probabilistic volume based on the 3D physical property model, then enhance the geological structure boundary, and perform geological classification based on the probabilistic volume;

[0007] Model construction steps: Extract isosurfaces from the 3D lithology probability volume, generate preliminary geological interfaces, perform structural processing, and then output the model after coupling with the attribute model.

[0008] The geological classification based on probabilistic volume includes the following steps: stratigraphic boundary delineation, trap boundary definition, boundary trimming, discontinuous boundary connection, cavity detection and repair, and non-boundary point geological classification.

[0009] The steps for defining the closed interface include: tracing closed interfaces that do not penetrate the borehole; when a closed interface intersects with or is contained by an existing interface marked as a boundary point, the closed interface is removed; when the number of grid points contained in a closed interface is less than a preset threshold, the closed interface is removed; the closed interface has polarity, which is determined by the average value of the physical property parameters inside and outside the closed interface; the closed interface with the average value of the parameters inside the closed interface is positive, and vice versa; only discontinuous interfaces with the same polarity can participate in the construction of continuous interfaces.

[0010] The interface trimming step includes: for interfaces that intersect starting from different borehole boundary points, if they are the same interface, trim the intersecting interface and retain the interface before the intersection; if they are different interfaces, compare the average physical property parameters of the area to be trimmed after the intersection with the area retained before the intersection, trim the interface with small differences, and retain the interface with large differences.

[0011] The discontinuous interface connection step includes: collecting discontinuous interfaces with the same boundary number, finding the shortest distance point pairs between each pair of interfaces, marking the boundary point, boundary number, and search status of the grid points through which the short distance point pair connection passes, and constructing a continuous interface network.

[0012] The process of acquiring multi-source data and constructing an initial model includes:

[0013] Acquire geophysical data, including linearly distributed electromagnetic, gravity, and seismic profile data;

[0014] Acquire drilling / trenching / logging data, including borehole / trench lithological stratification data, logging physical property data, and borehole / trench spatial coordinate data. Among these, the logging physical property data includes resistivity, density, and wave velocity.

[0015] The initial model is constructed by interpolating along the direction of the geophysical exploration line, using borehole lithological stratification as a hard constraint to generate a three-dimensional initial physical property field model.

[0016] The joint constraint 3D inversion algorithm yields a 3D physical property model including:

[0017] Set the objective function as ,in, Represents a vector of geophysical measured data. This represents the three-dimensional forward response, where m represents the three-dimensional matrix composed of physical property parameters. A priori model representing borehole constraints; Represents a weighted diagonal matrix of data. Indicates the model balance constraint weights. W represents the gradient constraint weights. m Let ∠m be the model weighting matrix, and ∠m be the gradient of parameter m.

[0018] Based on the current model parameters of the a-th iteration Calculate the forward response And calculate the residuals Solve for the Jacobian or Hessian matrix of the first or second order Taylor expansion of the objective function, as well as the gradient and step size, and update the model parameters for the (a+1)th iteration. ;

[0019] If the residual decrease rate is less than the threshold or the maximum number of iterations is reached, the final three-dimensional physical property model is output.

[0020] The creation of the lithology classifier includes:

[0021] Lithological points, located across multiple dimensions including physical properties, depth, and lithology, are used as vertices in an undirected graph. The edges between vertices are defined by the Gaussian distance between them. Points that are closer together have a larger Gaussian distance and a higher weight. Let x be the Gaussian distance between vertices p and q. p and x q Let be the eigenvectors of vertices p and q, respectively, and σ be the bandwidth parameter of the Gaussian kernel function.

[0022] The graph partitioning problem is transformed into the problem of dividing the graph into k non-interacting subsets and minimizing the sum of edge weights between each subset and all other subsets. The objective function is: A1, A2, ..., A k Let there be k mutually exclusive subsets. For subset With the remaining subsets The sum of edge weights between them;

[0023] Adjust the objective function to To prevent single-sample segmentation and to focus more on the weights of vertices within a subset, The sum of the weights of all vertices in the i-th subset is vol(A). The closer the vertex is to the other vertex, the greater its weight, and the smaller the objective function. i ) is a subset A i The sum of the weights of all vertices in the equation;

[0024] Introducing A j Indicator variables This indicates the indication of sample data i to subset j, where This represents the i-th node in the graph;

[0025] get H is the indicator variable h of all samples. i The matrix formed, T r For trace operation, T represents transpose;

[0026] At this point, the objective function is transformed into finding... The H matrix in the model is used to obtain the eigenvectors corresponding to the k smallest eigenvalues ​​of the Laplacian matrix L, thereby obtaining the lithology classifier.

[0027] The process of generating a three-dimensional lithology probability volume based on a three-dimensional physical property model includes:

[0028] The inverted three-dimensional physical property model is input into the lithology classifier, and the lithology probability distribution is calculated grid by grid.

[0029] Using borehole lithology as the calibration point, a three-dimensional lithology probability body with 100% matching of borehole lithology classification is obtained by correcting the probability field through Kriging co-simulation.

[0030] By adopting lithological probabilistic volume fusion modeling, and through normalization processing, physical property parameters with different dimensions and ranges are uniformly converted into standard probabilities, thereby avoiding the implicit weight bias caused by the difference in parameter values ​​and ensuring the objectivity of multi-source information fusion.

[0031] The step of inputting the inverted three-dimensional physical property model into the lithology classifier and calculating the lithology probability distribution grid by grid includes:

[0032] Set the t-th lithology L t The material property vector x follows a multivariate Gaussian distribution, i.e. μ t Indicates lithology L t The mean vector of physical properties, Σ t Represent the covariance matrix;

[0033] The inverted 3D physical property model mesh node values ​​m l Input to In the calculation, the lithology L of the l-th grid is obtained. t The probability is , P(L t Let be the prior probability of the t-th lithology. Let L be the posterior probability, i.e., the t-th lithology is L. t At that time, the l-th mesh of the three-dimensional physical property model is m l The probability estimate, P(L) s ) represents the prior probability of the s-th lithology;

[0034] Finally, N types of lithological probability vectors are generated for each grid [P(L1),P(L2),...,P(L...]. N ]], forming a three-dimensional lithology probability volume, where P(L1) is the first lithology prior probability, P(L N ) represents the prior probability of the Nth lithology.

[0035] The model construction steps specifically include the following:

[0036] Isosurfaces are extracted from lithological probability volumes. Stratigraphic topology networks are constructed based on stratigraphic interface division results in geological classification and borehole stratification data, and contact relationships are closed.

[0037] In identifying fault trajectories in physical discontinuities, implicit function interpolation is used to generate fault planes, spatial morphology is delineated based on high-resolution geophysical anomalies, and surfaces are generated using boreholes as control points.

[0038] The inverted physical property parameters are loaded as a continuous attribute field onto a 3D mesh to form an attribute model. The geological model generated by structural modeling and the attribute model are dynamically linked through a shared mesh topology to support parameter zoning statistics by lithology. Finally, the model is output.

[0039] The stratigraphic interface delineation steps include: starting from the borehole boundary point, searching in different directions using the weighted average probability of the lithology of the upper and lower layers, taking the grid point closest to the weighted average probability of the two layers as the boundary point, and obtaining the interface between the two lithologies. This process is repeated for different boundary points of different boreholes to obtain all interfaces within the survey area. For closed interfaces that do not penetrate the borehole, the corresponding interfaces are searched and marked using the same method based on the weighted average probability of the lithology of each interface in each borehole, ultimately obtaining all interfaces within the survey area. Simultaneously, the interface grid points are marked with boundary points, boundary numbers, and search status.

[0040] The hole detection and repair steps include: collecting all grid points on the same interface, calculating the Euclidean distance between any two grid points, and if the Euclidean distance is greater than the model space resolution but less than the maximum connection distance, it is judged as a hole. The boundary point, boundary number and search status of the grid points through which the connection line of this group of grid points passes are marked to complete the hole repair.

[0041] The non-boundary point geological classification step includes: assigning geological classification values ​​to other grid points of non-boundary points between adjacent boundaries; randomly selecting a non-boundary point; querying the geological classification of the six nearest directional boundary points; if the boundary numbers of at least two directional boundary points are different, then determining the geological classification of the non-boundary point based on the common geological classification of the boundary points with different boundary numbers.

[0042] This invention has the following advantages: a refined 3D geological modeling method based on geophysical 3D forward and inverse algorithms, which directly constructs a spatially continuous model through 3D joint inversion, avoiding interpolation errors between survey lines; forced matching of borehole data + trend control of geophysical data, taking into account both local accuracy and global rationality; algorithm-driven process from inversion to lithology classification, reducing manual intervention; applicable to high-precision scenarios such as urban underground space, mineral exploration, and geological hazard assessment. Attached Figure Description

[0043] Figure 1 This is a schematic diagram of the process of the present invention;

[0044] Figure 2 This is a flowchart illustrating the three-dimensional inversion calculation process.

[0045] Figure 3 A schematic diagram for modeling physical property-lithology mapping;

[0046] Figure 4 This is a schematic diagram of the polarity of the closed interface;

[0047] Figure 5 The diagrams illustrate the trimming of the interface, where (a) shows the interface after trimming the intersection, and (b) shows the interface between the two parts after trimming. Figure 1 (c) is a schematic diagram of the interface between the two trimming processes. Figure 2 . Detailed Implementation

[0048] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. The components of the embodiments of this application described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the detailed description of the embodiments of this application provided below with reference to the accompanying drawings is not intended to limit the scope of protection of the claimed application, but merely represents selected embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without inventive effort are within the scope of protection of this application. The present invention will be further described below with reference to the accompanying drawings.

[0049] This invention specifically relates to a refined 3D modeling method for geological models based on the fusion of geophysical and drilling data. Addressing the low accuracy of traditional borehole interpolation modeling and model distortion between geophysical survey lines, it proposes a point-line fusion 3D inversion technology framework. Using borehole lithological stratification and physical properties as hard constraints and geophysical profile data as trend control, a globally optimal subsurface physical property model is constructed through 3D regularized inversion iterative calculation. Based on physical property-lithology probability mapping and geological boundary reinforcement, a high-precision 3D lithological structure is generated. Combining topological network construction and attribute coupling, automated modeling of geological bodies with centimeter- to meter-level accuracy is achieved. This method overcomes the limitations of point-line data, significantly improves the accuracy of cross-survey line regional modeling, and achieves deep synergy between the local accuracy of borehole data and the regional coverage of geophysical data, providing a high-fidelity geological foundation for the construction of a realistic 3D China.

[0050] like Figure 1 and Figure 2As shown, the process includes three steps: inversion calculation, geological classification, and model building. The inversion calculation specifically includes the following:

[0051] 1. Acquire multi-source data and build an initial model;

[0052] Geophysical data: Linearly distributed electromagnetic / gravity / seismic profile data (such as resistivity, wave velocity, density field).

[0053] Drilling / Trenching / Logging Data: Borehole / trench lithological stratification, logging physical properties (density, wave velocity, resistivity), borehole / trench spatial coordinates; among which, lithological stratification data refers to borehole and trench geological logging data, mainly including lithological description, stratum depth, and other data.

[0054] Initial model construction: Using borehole lithological stratification as a hard constraint, a three-dimensional initial physical property field model (resistivity / velocity / density volume) is generated by interpolation along the direction of the geophysical probe.

[0055] 2. Jointly constrained 3D inversion algorithm;

[0056] The objective function is designed as follows:

[0057] ,

[0058] in, Represents a vector of geophysical measured data. This represents the three-dimensional forward response, where m represents the three-dimensional matrix composed of physical property parameters, and represents the model parameters to be solved. A priori model representing borehole constraints; This represents a weighted diagonal matrix of data. High-reliability data points (such as low-noise or near-drilling points) are assigned high weights, forcing priority fitting of these data points during inversion. Low-reliability data points (such as high-noise or far from the control area) are assigned low weights, allowing for larger fitting errors to avoid excessive model oscillations. It also serves as a data normalization mechanism, eliminating the influence of different physical dimensions (such as resistivity vs. seismic wave velocity), ensuring stable convergence of the inversion. m The model weighting matrix assigns high weights to model parameters at borehole locations, forcing a fit to those model parameters during inversion. The weights represent the model's balance constraints (which control inversion stability). ∠m represents the gradient constraint weight. By adjusting the model gradient regularization term, the model space is constrained, the multivariable nature of the model is reduced, and the ill-posedness of the inversion is resolved (highlighting the geological boundary). ∠m is the gradient of parameter m.

[0059] Based on the current model parameters of the a-th iteration Calculate the forward response and calculating residuals Solve for the Jacobian or Hessian matrix of the first or second order Taylor expansion of the objective function, as well as the gradient and step size, and update the model parameters for the (a+1)th iteration. If the residual decrease rate is less than the threshold or the maximum number of iterations is reached, the final three-dimensional physical property model is output.

[0060] Furthermore, the geological classification specifically includes the following:

[0061] 1. Physical property-lithology mapping modeling;

[0062] (1) Establish the statistical mapping relationship between the physical property parameters (such as resistivity / density / wave velocity, etc.), depth location and lithology of the borehole location;

[0063] (2) Use clustering methods such as Gaussian mixture model (GMM), random forest, K-means or spectral clustering to classify the lithological categories corresponding to the physical property range;

[0064] Furthermore, such as Figure 3 As shown, this is achieved through spectral clustering, which is a problem of segmenting undirected graphs. Lithological points, located across multiple dimensions including physical properties (resistivity / density / wave velocity), depth, and lithology, are considered as vertices of the undirected graph. The edges between vertices are represented by the Gaussian distance between them. S i,j Let x be the Gaussian distance between vertices p and q. p and x q These are the feature vectors of vertices p and q (such as resistivity, density, wave velocity, depth location, lithology, etc.). Points that are closer to each other have a larger Gaussian distance and a higher weight.

[0065] Therefore, the graph partitioning problem is transformed into the problem of dividing the graph into k disjoint subsets, and minimizing the sum of edge weights between each subset and all other subsets. The objective function is... as follows:

[0066] ,

[0067] A1, A2, …, A k Let there be k mutually exclusive subsets. For subset With the remaining subsets The sum of edge weights between them.

[0068] To prevent single-sample segmentation and to focus more on the weights of each vertex within a subset, the new objective function... The adjustments are as follows, where the denominator... The sum of the weights of all vertices in the i-th subset is given. The closer the vertex is to the other vertex, the greater its weight and the smaller the objective function.

[0069] ,

[0070] vol(A i ) is a subset A i The sum of the weights of all vertices in the equation, i.e., vol(A) i )= , where dv is the degree of vertex v (the sum of the weights of the connecting edges).

[0071] Introducing A j Indicator variable h j,i This indicates the indication of sample data i to subset j, specifically:

[0072] ,

[0073] The Laplacian matrix L = DS, where D is the degree matrix, with each element representing the number of lines connecting each vertex; it is a diagonal matrix. S is the adjacency matrix, with each element representing the weight between any two points, i.e., the Gaussian distance S between them. i,j It is also a diagonal matrix, which can be proven by the property of Laplace:

[0074] ,

[0075] Among them, h i It is an indicator variable h j,i The vector representation of h, where each element is h j,i This indicates that sample i corresponds to the indicator value of a different subset, and T represents the transpose.

[0076] Further results were obtained:

[0077] ,

[0078] At this point, the objective function is transformed into calculating the H matrix in the above equation, where H is the indicator variable h of all samples. i The matrix formed, T r The trace operation is used to transform the expression of Ncut into matrix form. According to the Rayleigh-Ritz method, it is only necessary to obtain the eigenvectors corresponding to the k smallest eigenvalues ​​of the Laplace matrix L, thus obtaining the lithology classifier. The value of k is a hyperparameter that needs to be tuned. The number of all borehole lithology types in the modeling area can be selected as the value of k. T represents transpose.

[0079] 2. Generation of three-dimensional lithological probabilistic volumes;

[0080] (1) Input the three-dimensional physical property model obtained by inversion into the lithology classifier, and calculate the lithology probability distribution P(L) grid by grid. t |m) ;

[0081] Taking Gaussian Mixture Model (GMM) as an example (Gaussian Mixture Model (GMM) and Random Forest natively support probabilistic outputs, while clustering methods such as K-means or spectral clustering are essentially distance partitions rather than probabilistic models, directly outputting discrete labels, requiring probability transformation through techniques such as distance transformation and resampling):

[0082] Assume the t-th lithology L t The material property vector x = [ρ, R, V] (ρ is density, R is resistivity, V is wave velocity) follows a multivariate Gaussian distribution:

[0083] ,

[0084] Where, μ t Indicates lithology L t The mean vector of physical properties (statistics from borehole data), Σ t This represents the covariance matrix (characterizing the correlation between physical property parameters).

[0085] The inverted 3D physical property model mesh node values ​​m l =[ρ l ,R l, V l Enter to In the calculation, the lithology L of the l-th grid is obtained. t The probability is , P(L t Let P(L) be the prior probability of the t-th lithology, and T denote the transpose. s ) represents the prior probability of the s-th lithology;

[0086] Finally, N types of lithological probability vectors are generated for each grid [P(L1),P(L2),...,P(L...]. N )], forming a three-dimensional lithological probability body.

[0087] (2) Using the borehole lithology as the calibration point, the probability volume obtained in the previous step is corrected by Kriging co-simulation to obtain a three-dimensional lithology probability volume with 100% matching of lithology classification at the borehole.

[0088] Specifically, it includes the following:

[0089] (a) Forced probability setting at the borehole location: If the borehole location is sandstone, then P(L1)=1.0, P(L2)=0,...P(L... N )=0, P(L1) is the first lithological prior probability, P(L N ) represents the prior probability of the Nth lithology.

[0090] (b) Constructing the lithology probability covariance matrix: Obtain borehole physical property data for different lithologies, including resistivity ρ and wave velocity V. Statistical analysis is performed on the physical property data for each lithology to obtain its mean vector μ. t Using borehole data, the covariance between each type of lithological physical property parameter is calculated, resulting in the covariance matrix Σ. t The elements Σ of the covariance matrix ij The formula Σij=Cov(ρ) is used. i V j The calculation reflects the linear correlation between physical property parameters. The physical property covariance matrices of each lithology are integrated to form a lithology probability covariance matrix. This matrix describes the correlation between physical property parameters of different lithologies, providing a basis for subsequent probability calculations and kriging interpolation.

[0091] (c) Solving for the Kriging weights:

[0092] ,

[0093] in, and The weighting coefficients fitted by the variogram function can be obtained by solving the co-kriging equations. To assist the physical property field at point x r The value of x0 is the spatial coordinate vector (X0, Y0, Z0) of the current grid point (unknown point) to be predicted. The position coordinates of the u-th borehole data point (X) u ,Y u Z u ), where u = 1, 2, ..., n (a total of n holes). Let X be the position coordinates (X) of the r-th geophysical auxiliary data point. r ,Y r Z r ), where r=1,2,...,m (a total of m physical property measurement points, usually from the mesh of the three-dimensional physical property model). For the lithology L at location x0 t The co-kriging prediction probability, such as the sandstone probability P=0.85 for the prediction point (100,200,50).

[0094] 3. Enhanced geological structural boundaries;

[0095] (1) At the point where the physical property gradient changes abruptly ( Apply boundary reinforcement and calculate the property gradient magnitude for each mesh:

[0096] .

[0097] in, The property gradient magnitude represents the rate of change of a property parameter (such as density, resistivity, wave velocity) at the l-th grid point, where ρ is density, R is resistivity, V is wave velocity, and ε is a threshold. If the value exceeds the threshold, boundary reinforcement is performed.

[0098] (2) Use anisotropic diffusion filtering to smooth data within homogeneous regions and sharpen lithological interfaces:

[0099] ,

[0100] If the diffusion coefficient c approaches 0 in regions with large gradients, then sharp boundaries are preserved; if the diffusion coefficient c approaches 1 in regions with small gradients, then noise in homogeneous regions is smoothed. P represents the lithology L at spatial location x=(X,Y,Z). i The attribution probability, where t represents pseudo-time (here, the number of iterations), is used to simulate the gradual evolution of the probability field P during the diffusion process. After this step, the boundary of the three-dimensional lithological probability volume will be enhanced, and the noise in the homogeneous region within the boundary will be smoothed. Let P be the gradient of parameter P.

[0101] (3) Integrate regional geological patterns (such as fault dip and stratigraphic attitude) to correct boundary geometry.

[0102] 4. Geological classification based on probabilistic volumes;

[0103] Even after the above processing, the three-dimensional lithological probability body still contains discontinuous strata such as bubble-like and clumpy structures in some areas, which does not match the actual geological conditions. Therefore, further processing is required. Homogeneous rock strata need to cross the discontinuous probability bodies such as bubble-like and clumpy structures to meet the actual geological conditions of continuous strata.

[0104] (1) Delineation of stratigraphic boundaries: Starting from the borehole boundary point, the probability weighted average of the upper and lower lithologies is used to search in different directions. The grid point closest to the probability weighted average of the two lithologies is taken as the boundary point to obtain the boundary between the two lithologies. The different boundary points of different boreholes are searched in turn to obtain all the boundaries in the survey area. For the closed boundary that does not pass through the borehole, the corresponding boundary is searched and marked in turn according to the lithology probability weighted average of each boundary of each borehole, in the same way as above, to obtain all the boundaries in the survey area. At the same time, the boundary point, boundary number, search status and other attribute information are marked on the boundary grid points.

[0105] (2) Delineation of Trap Interfaces: Due to factors such as measurement errors, method limitations, and environmental differences, the physical properties of continuous homogeneous rock strata differ, leading to the generation of discontinuous trap interfaces. Furthermore, these factors often result in blurred boundaries between the physical properties of different rock strata, although relative differences still exist. Therefore, determining the trap interface delineation rules and identifying the physical property differences inside and outside the trap interface is crucial for subsequent connection of discontinuous interfaces. For tracing trap interfaces that do not penetrate the borehole, the trap interface is removed when it intersects with or is contained within an existing interface marked as a boundary point. The trap interface is also removed when the number of grid points contained within it is less than a preset threshold. Additionally, trap interfaces have polarity, such as... Figure 4 As shown, the polarity is determined by the average value of the physical property parameters inside and outside the closed interface. The positive polarity is the one where the average value of the parameters inside the closed interface is greater than that outside the closed interface, and the negative polarity is the one where the average value of the parameters inside the closed interface is greater than that outside the closed interface. This provides a basis for whether or not the discontinuous interfaces are connected in step (4). Only discontinuous interfaces with the same polarity can participate in the construction of continuous interfaces.

[0106] (3) Trim the interface: such as Figure 5 As shown in the figure, (a) is a schematic diagram of the interface after trimming the intersection, specifically reflecting that 1 and 2 are the same interface from which two boreholes originate, and the interface after trimming the intersection; (b) is a schematic diagram of the interface between the two after trimming. Figure 1 Specifically, it reflects the different interfaces starting from boreholes 1 and 2. After the intersection, the average physical property parameters of region 3 are closer to those of region 4 before the intersection. The interface between the two is trimmed. (c) is a schematic diagram of the trimmed interface. Figure 2 Specifically, this reflects that regions 1 and 2 are different interfaces originating from two boreholes. After the intersection, the average physical property parameters of region 3 are closer to those of region 5 before the intersection, so the interface between them is trimmed. The traditional interface tracking algorithm, which stops when it encounters a search point, is affected by the borehole search order and has unreasonable results. The improved algorithm, which stops only when it encounters a boundary point, ensures the closure of the starting point of the same interface among multiple boreholes and the sufficiency of searching different interface termination states, but it introduces new problems such as abnormal interface intersections and redundant interface generation. Therefore, interface trimming is necessary. For interfaces intersecting from different borehole boundary points, if they are the same interface, the interface after the intersection is trimmed, and the interface before the intersection is retained. If they are different interfaces, the average physical property parameters of the area to be trimmed after the intersection are compared with those of the area retained before the intersection. The interface with smaller differences is trimmed, and the interface with larger differences is retained.

[0107] (4) Connecting discontinuous interfaces: Due to the existence of closed interfaces, discontinuous strata such as bubble-like or clumpy formations may appear locally, which does not conform to the actual geological conditions. Therefore, further processing is required. Homogeneous rock strata need to cross discontinuous probability bodies such as bubble-like or clumpy formations to meet the actual geological conditions. Collect discontinuous interfaces with the same boundary number, find the shortest distance point pairs between each pair of interfaces, and mark the boundary points, boundary numbers, search status, and other attribute information of the grid points through which the short distance point pairs are connected, and construct a continuous interface network.

[0108] (5) Hole detection and repair: Collect all grid points on the same interface, calculate the Euclidean distance between each pair of grid points. If the Euclidean distance is greater than the model space resolution and less than the maximum connection distance, it is judged as a hole. Mark the boundary point, boundary number, search status and other attribute information of the grid points through which the connection line of this group of grid points passes, and complete the hole repair.

[0109] (6) Geological Classification of Non-Boundary Points: Geological classification values ​​are assigned to other grid points at non-boundary points between adjacent interfaces. A non-boundary point is randomly selected, and the geological classifications of its six nearest directional boundary points are queried. If at least two directional boundary points have different boundary numbers, the geological classification of the non-boundary point can be determined based on the common geological classification of these boundary points. For the remaining closed or semi-closed non-boundary points that cannot be selected, their geological classification values ​​can be determined based on another existing geological classification method for boundary points on the interface.

[0110] Furthermore, model construction specifically includes the following:

[0111] 1. Structural modeling;

[0112] (1) Generation of topological relationships: Based on the stratigraphic interface division results in geological classification and borehole stratification data, a stratigraphic topological network is constructed, and contact relationships (conformity, unconformity, fault) are forcibly closed.

[0113] (2) Complex structural processing: Fault modeling: Identify fault trajectories in physical discontinuities and generate fault planes using implicit function (RBF) interpolation; Lens / thin layer processing: Delineate spatial morphology based on high-resolution geophysical anomalies and generate surfaces using boreholes as control points. After the above two steps, a geological model is finally formed.

[0114] 2. Attribute model coupling;

[0115] The inverted physical property parameters (density, resistivity) are loaded as continuous attribute fields onto the three-dimensional mesh;

[0116] The geological model and attribute model generated by structural modeling are dynamically linked through a shared grid topology, and parameters are statistically analyzed by lithology partition.

[0117] 3. Model output;

[0118] Multi-scale output: Supports the generation of LOD (Level of Detail) models, with exploration area scale (hundred-meter grid): regional tectonic framework; engineering scale (meter grid): lithological interfaces and attribute parameters.

[0119] The above description is merely a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the forms disclosed herein and should not be construed as excluding other embodiments. It can be used in various other combinations, modifications, and improvements, and can be altered within the scope of the concept described herein through the above teachings or related technologies or knowledge. Modifications and variations made by those skilled in the art that do not depart from the spirit and scope of the present invention should be within the protection scope of the appended claims.

Claims

1. A refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms, characterized in that: The method includes: Inversion calculation steps: acquire multi-source data and construct an initial model, and use a joint constrained 3D inversion algorithm to obtain a 3D physical property model; Geological classification steps: Create a lithology classifier, generate a 3D lithology probabilistic volume based on the 3D physical property model, then enhance the geological structure boundary, and perform geological classification based on the probabilistic volume; Model construction steps: Extract isosurfaces from the 3D lithology probability volume, generate preliminary geological interfaces, perform structural processing, and then output the model after coupling with the attribute model. The geological classification based on probabilistic volume includes the following steps: stratigraphic boundary delineation, trap boundary definition, boundary trimming, discontinuous boundary connection, cavity detection and repair, and non-boundary point geological classification. The steps for defining the closed interface include: tracing closed interfaces that do not penetrate the borehole; when a closed interface intersects with or is contained by an existing interface marked as a boundary point, the closed interface is removed; when the number of grid points contained in a closed interface is less than a preset threshold, the closed interface is removed; the closed interface has polarity, which is determined by the average value of the physical property parameters inside and outside the closed interface; the closed interface with the average value of the parameters inside the closed interface is positive, and vice versa; only discontinuous interfaces with the same polarity can participate in the construction of continuous interfaces. The interface trimming step includes: for interfaces that intersect starting from different borehole boundary points, if they are the same interface, trim the intersecting interface and retain the interface before the intersection; if they are different interfaces, compare the average physical property parameters of the area to be trimmed after the intersection with the area retained before the intersection, trim the interface with small differences, and retain the interface with large differences. The discontinuous interface connection step includes: collecting discontinuous interfaces with the same boundary number, finding the shortest distance point pairs between each pair of interfaces, marking the boundary point, boundary number, and search status of the grid points through which the short distance point pair connection passes, and constructing a continuous interface network.

2. The refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithm according to claim 1, characterized in that: The process of acquiring multi-source data and constructing an initial model includes: Acquire geophysical data, including linearly distributed electromagnetic, gravity, and seismic profile data; Acquire drilling / trenching / logging data, including borehole / trench lithological stratification data, logging physical property data, and borehole / trench spatial coordinate data. Among these, the logging physical property data includes resistivity, density, and wave velocity. The initial model is constructed by interpolating along the direction of the geophysical exploration line, using borehole lithological stratification as a hard constraint to generate a three-dimensional initial physical property field model.

3. The refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms according to claim 1, characterized in that: The joint constraint 3D inversion algorithm yields a 3D physical property model including: Set the objective function as ,in, Represents a vector of geophysical measured data. This represents the three-dimensional forward response, where m represents the three-dimensional matrix composed of physical property parameters. A priori model representing borehole constraints; Represents a weighted diagonal matrix of data. Indicates the model balance constraint weights. W represents the gradient constraint weights. m Let ∠m be the model weighting matrix, and ∠m be the gradient of parameter m. Based on the current model parameters of the a-th iteration Calculate the forward response And calculate the residuals Solve for the Jacobian or Hessian matrix of the first or second order Taylor expansion of the objective function, as well as the gradient and step size, and update the model parameters for the (a+1)th iteration. ; If the residual decrease rate is less than the threshold or the maximum number of iterations is reached, the final three-dimensional physical property model is output.

4. The refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithm according to claim 3, characterized in that: The creation of the lithology classifier includes: Lithological points, located across multiple dimensions including physical properties, depth, and lithology, are used as vertices in an undirected graph. The edges between vertices are defined by the Gaussian distance between them. Points that are closer together have a larger Gaussian distance and a higher weight. Let x be the Gaussian distance between vertices p and q. p and x q Let be the eigenvectors of vertices p and q, respectively, and σ be the bandwidth parameter of the Gaussian kernel function; The graph partitioning problem is transformed into the problem of dividing the graph into k non-interacting subsets and minimizing the sum of edge weights between each subset and all other subsets. The objective function is: A1, A2, ..., A k Let there be k mutually exclusive subsets. For subset With the remaining subsets The sum of edge weights between them; Adjust the objective function to ,in, The sum of the weights of all vertices in the i-th subset is vol(A). The closer the vertex is to the other vertex, the greater its weight, and the smaller the objective function. i ) is a subset A i The sum of the weights of all vertices in the equation; Introducing A j Indicator variables This indicates the indication of sample data i to subset j, where This represents the i-th node in the graph; get H is the indicator variable h of all samples. i The matrix formed, T r For trace operation, T represents transpose; At this point, the objective function is transformed into finding... The H matrix in the model is used to obtain the eigenvectors corresponding to the k smallest eigenvalues ​​of the Laplacian matrix L, thereby obtaining the lithology classifier.

5. A refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms according to claim 1, characterized in that: The process of generating a three-dimensional lithology probability volume based on a three-dimensional physical property model includes: The inverted three-dimensional physical property model is input into the lithology classifier, and the lithology probability distribution is calculated grid by grid. Using borehole lithology as the calibration point, a three-dimensional lithology probability body with 100% matching of borehole lithology classification is obtained by correcting the probability field through Kriging co-simulation. A lithological probabilistic volume fusion model was adopted, and through normalization, physical property parameters with different dimensions and ranges were uniformly converted into standard probabilities.

6. A refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms according to claim 5, characterized in that: The step of inputting the inverted three-dimensional physical property model into the lithology classifier and calculating the lithology probability distribution grid by grid includes: Set the t-th lithology L t The material property vector x follows a multivariate Gaussian distribution, i.e. μ t Indicates lithology L t The mean vector of physical properties, Σ t Let T denote the covariance matrix; The inverted 3D physical property model mesh node values ​​m l Input to In the calculation, the lithology L of the l-th grid is obtained. t The probability is , P(L t Let be the prior probability of the t-th lithology. Let L be the posterior probability, i.e., the t-th lithology is L. t At that time, the l-th mesh of the three-dimensional physical property model is m l The probability estimate, P(L) s ) represents the prior probability of the s-th lithology; Finally, N types of lithological probability vectors are generated for each grid [P(L1),P(L2),...,P(L...]. N ]], forming a three-dimensional lithology probability volume, where P(L1) is the first lithology prior probability, P(L N ) represents the prior probability of the Nth lithology.

7. A refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms according to claim 1, characterized in that: The model construction steps specifically include the following: Isosurfaces are extracted from lithological probability volumes. Stratigraphic topology networks are constructed based on stratigraphic interface division results in geological classification and borehole stratification data, and contact relationships are closed. In identifying fault trajectories in physical discontinuities, implicit function interpolation is used to generate fault planes, spatial morphology is delineated based on high-resolution geophysical anomalies, and surfaces are generated using boreholes as control points. The inverted physical property parameters are loaded as a continuous attribute field onto a 3D mesh to form an attribute model. The geological model generated by structural modeling and the attribute model are dynamically linked through a shared mesh topology, and finally the model is output.

8. A refined three-dimensional geological modeling method based on geophysical three-dimensional forward and inverse modeling algorithms according to claim 1, characterized in that: The stratigraphic interface delineation steps include: starting from the borehole boundary point, searching in different directions using the weighted average probability of the upper and lower lithologies, taking the grid point closest to the weighted average probability of the two lithologies as the boundary point, and obtaining the interface between the two lithologies. This process is repeated for different boundary points of different boreholes to obtain all interfaces within the survey area. For closed interfaces that do not pass through boreholes, the corresponding interfaces are searched and marked in the same way based on the weighted average probability of the lithologies of each interface in each borehole, thus obtaining all interfaces within the survey area. At the same time, the boundary points, boundary numbers, and search status are marked on the interface grid points. The hole detection and repair steps include: collecting all grid points on the same interface, calculating the Euclidean distance between any two grid points, and if the Euclidean distance is greater than the model space resolution but less than the maximum connection distance, it is judged as a hole. The boundary point, boundary number and search status of the grid points through which the connection line of this group of grid points passes are marked to complete the hole repair. The non-boundary point geological classification step includes: assigning geological classification values ​​to other grid points of non-boundary points between adjacent boundaries; randomly selecting a non-boundary point; querying the geological classification of the six nearest directional boundary points; if the boundary numbers of at least two directional boundary points are different, then determining the geological classification of the non-boundary point based on the common geological classification of the boundary points with different boundary numbers.

Images