A three-dimensional reconstruction method of complex fault block geologic body based on joint constraint of covariance function and drift term

Abstract

A method for 3D reconstruction of complex fault-block geological bodies based on joint constraints of covariance function and drift term includes the following steps: acquiring basic stratigraphic data of the complex fault-block geological body and preprocessing the basic stratigraphic data; constructing an implicit potential field equation containing fault step terms and constructing a local signed distance field of the fault; based on the basic stratigraphic data and the implicit potential field equation, constructing a system of dual linear equations containing scalar constraints, gradient constraints, and drift term constraints; solving the system of dual linear equations to generate a full-space 3D scalar field, and then generating a 3D geological model of the complex fault-block geological body. The method provided by this invention can automate fault processing, has good dynamic update capability, improves interpolation accuracy, and significantly improves computational efficiency.

Description

Technical Field

[0001] This invention relates to the fields of geological engineering informatization and computer graphics, specifically to a method for three-dimensional reconstruction of complex fault-block geological bodies based on joint constraints of covariance function and drift term. Background Technology

[0002] In mining, tunnel engineering, and oil exploration, constructing accurate three-dimensional geological models is fundamental. Existing technologies for constructing three-dimensional geological models mainly fall into two categories: the first is explicit modeling, which relies on manually drawing stratigraphic boundaries section by section in CAD software and then connecting them; the second is implicit modeling, which mainly utilizes radial basis function (RBF) interpolation to generate isosurfaces.

[0003] Existing explicit modeling methods suffer from problems such as difficulty in handling complex topological structures and the need to rebuild the model when updating data. These methods are extremely inefficient and highly dependent on human experience.

[0004] Conventional implicit modeling methods require physically dividing the model into multiple independent blocks when dealing with faults, which can easily lead to the loss of mathematical representation of the continuity of strata on both sides of the fault. At the same time, non-physical oscillations (Runge) are prone to occur when borehole data is sparse. In addition, this method cannot automatically handle complex fault networks.

[0005] Existing methods for constructing 3D geological models are ineffective in scenarios such as complex fault block structures formed by multiple faults cutting each other, uneven distribution of exploration boreholes, sparse data, and the need for rapid model updates in response to new exploration data. Summary of the Invention

[0006] To address the problems of existing technologies in constructing 3D geological models that cannot automatically handle complex fault networks and unstable interpolation of sparse data, this invention provides a 3D reconstruction method for complex fault-block geological bodies based on joint constraints of covariance function and drift term. This method expresses faults mathematically, requires no mesh geometry cutting, is entirely data-driven, automatically processes sparse borehole data, and is dynamically updatable.

[0007] This invention provides a method for three-dimensional reconstruction of complex fault-block geological bodies based on joint constraints of covariance function and drift term, comprising the following steps:

[0008] The basic stratigraphic data of the complex fault-block geological body is obtained, and the basic stratigraphic data is preprocessed; the basic stratigraphic data includes borehole data, surface attitude data, and fault data; scalar constraint data is obtained based on the borehole data, and gradient constraint data is obtained based on the surface attitude data.

[0009] An implicit potential field equation including a fault step term is constructed, and the fault step term is constructed based on the acquired fault data; the fault step term is... ,in, This represents the total number of faults. For the first Vertical displacement of the fault plane; For the first The Heaviside step function of a fault is 1 when point x is located on the hanging wall and 0 when it is located on the footwall; specifically... , Let x be the sign distance from point x to the k-th fault plane;

[0010] Based on the aforementioned basic stratigraphic data and the aforementioned implicit potential field equations, a dual linear equation system is constructed.

[0011] Solving the dual linear equations generates a full-space three-dimensional scalar potential field, which in turn generates a three-dimensional geological model of the complex fault-block geological body.

[0012] Furthermore, the implicit potential field equation containing the fault step term is specifically as follows:

[0013] ;

[0014] in, The equations are implicit potential field equations;

[0015] Here, is the radial basis kernel function term; where, where is the weight coefficient of the i-th radial basis function; N is the total number of data constraint points; Here, r is the radial basis function kernel, and r is the Euclidean distance. , Let x be the three-dimensional spatial position vector of point x. Let i be the spatial location vector of the i-th data point;

[0016] This is a polynomial drift term.

[0017] Furthermore, the radial basis kernel function Use the Matérn kernel function, Gaussian kernel function, or cubic kernel function.

[0018] Furthermore, polynomial drift term It is represented by a first-order or second-order polynomial.

[0019] Furthermore, the dual linear equation system specifically comprises:

[0020] ;

[0021] in, Let covariance be the covariance matrix between scalar constraint points. , , The number of scalar constraint points;

[0022] The scalar-gradient cross covariance matrix is... , , The number of gradient constraint points. Let be the normal vector of the i-th attitude observation point;

[0023] Let be the second-order covariance matrix between gradient constraint points. , ;

[0024] U and V are drift term constraint matrices; for first-order polynomial drift... , , ,in Let i be the three-dimensional coordinates of the i-th scalar constraint point; , ,in Let J be the three-dimensional coordinates of the j-th scalar constraint point;

[0025] This is a vector of weight coefficients; For Lagrange multipliers; For the Lagrange multiplier of the drift term; It is a scalar constraint vector; This is the gradient constraint vector.

[0026] Furthermore, different solution methods are used depending on the size of the matrix when solving the dual linear equation system.

[0027] The method provided by this invention can achieve the following beneficial effects:

[0028] (1) Automated fault processing. It can automatically process dozens of intersecting fault networks without manual intervention in the topology, and the accuracy of fault intersection processing can reach the meter level;

[0029] (2) It has good dynamic update capability. When new borehole data is added, the model can be updated simply by resolving the matrix, reducing the time from several hours to several seconds to several minutes, which can support real-time exploration data integration;

[0030] (3) Improved interpolation accuracy. Compared with the traditional radial basis function (RBF) method, the error at the boundary is reduced by more than 60%, and the formation continuity can still be maintained in sparse data areas (borehole spacing > 500m);

[0031] (4) The computational efficiency is significantly improved. The method provided by this invention takes about 30 seconds to solve for tens of thousands of data points (using GPU acceleration) and about 5 minutes to solve for hundreds of thousands of data points (iterative solution). Attached Figure Description

[0032] To gain a more complete understanding of the invention, reference will now be made to the following description taken in conjunction with the accompanying drawings, wherein:

[0033] Figure 1 This is a schematic diagram of the overall process of the method described in this invention. Detailed Implementation

[0034] To clearly illustrate the purpose, technical details, and effective applications of this invention, and to facilitate understanding and implementation by those skilled in the art, a further detailed description will be provided below in conjunction with the embodiments and accompanying drawings. Obviously, the embodiments described herein are for illustrative and explanatory purposes only and are not intended to limit the scope of the invention.

[0035] This invention provides a method for three-dimensional reconstruction of complex fault-block geological bodies based on joint constraints of covariance function and drift term. See [link to relevant documentation]. Figure 1 It mainly includes the following steps:

[0036] Step S01: Obtain basic stratigraphic data of the geological body and preprocess the basic stratigraphic data; obtain stratigraphic interface constraint data, i.e., scalar constraint data, based on borehole data, and obtain gradient constraint data based on surface attitude data.

[0037] When performing 3D modeling of the complex fault-block geological body, the first step is to obtain the basic stratigraphic data of the geological body and extract the key data. For example, it is necessary to obtain borehole data, including the location, depth, and lithology of the boreholes; coordinates, dip, and dip angle data of surface attitude measurement points; and fault data, including surface traces, dip, and dip angle.

[0038] After cleaning and preprocessing the acquired basic stratigraphic data, the borehole data is read, and the stratigraphic boundary point coordinates (x, y, z) and lithological labels are extracted; the surface occurrence data is converted into normal vectors.

[0039] Assume that the geological strata are represented by isosurfaces of a three-dimensional continuous scalar field Z(x), and apply stratigraphic interface constraints and surface attitude constraints to them.

[0040] Constraint points on the formation interface That is, the scalar constraint is satisfied:

[0041]

[0042] in, The reference isosurface value for the target stratigraphic interface.

[0043] For constraint points with attitude observation Its scalar field gradient direction should be parallel to the formation interface normal vector n; since the formation interface tangent vector t is perpendicular to the normal vector, it satisfies:

[0044] , that is

[0045] in, t is the spatial gradient vector of the potential field; t is the tangent vector of the formation interface.

[0046] Regarding the aforementioned attitude data, the data obtained from actual exploration measurements is the dip direction. And dip angle It needs to be converted into a normal vector n:

[0047] .

[0048] Step S02: Construct an implicit potential field equation containing fault step terms, and construct a fault step function based on the obtained fault data.

[0049] Implicit potential field equations are the core mathematical equations used in implicit modeling to describe the spatial distribution of geometric shapes such as strata. Essentially, they are constructed by using continuous scalar potential field functions to divide spatial regions with different potential field values, thereby characterizing the boundaries and distribution of strata.

[0050] Unlike traditional methods that physically cut the mesh, this invention mathematically expresses the fault by superimposing a heaviside step function into the potential field equation. The implicit potential field equation containing the fault step term is specifically as follows:

[0051]

[0052] in,

[0053] It is an implicit scalar potential field function.

[0054] This is the radial basis kernel function term, used to control the spatial correlation of local interpolation. Wherein, is the weight coefficient of the i-th radial basis function; N is the total number of data constraint points. The radial basis function can be a Matérn kernel, a Gaussian kernel, or a cubic kernel, used to control spatial correlation. The specific form of expression can be:

[0055] Cubic kernel function:

[0056] Matérn kernel function:

[0057] Gaussian kernel function:

[0058] Where r is the Euclidean distance. Let x be the three-dimensional spatial position vector of point x. is the three-dimensional spatial position vector of the i-th data point; l is the correlation length scale of the kernel function (controlling the range of influence); v is the smoothness parameter; To correct the Bessel function.

[0059] This is a polynomial drift term used to capture large-scale trends in strata, such as monoclines and folds. It is typically represented by a low-order (first or second order) polynomial. Its specific expression can be:

[0060] First-order drift (linear trend): ;

[0061] Second-order drift (secondary trend): ;

[0062] Where a0, a1...a9 are polynomial drift coefficients.

[0063] For fault step terms, where M is the total number of faults; The vertical displacement of the k-th fault; Let x be the Heaviside step function of the k-th fault, which is 1 when point x is located on the hanging wall and 0 when it is located on the footwall. This causes a numerical jump in the potential field at the fault plane. It allows strata to naturally shift without the need for geometric cutting; furthermore, it can handle complex situations where multiple faults intersect, and the handling of fault intersections is achieved through function multiplication. Automatically implemented, among which, Let be the Heaviside step function of the i-th fault. Let be the Heaviside step function of the j-th fault.

[0064] Specifically, Calculate according to the following formula:

[0065]

[0066] in, Let x be the signed distance from point x to the k-th fault plane; let the equation of the k-th fault plane be... Then the signed distance of point x is: , where a, b, c, and d are the coefficients of the fault plane equation.

[0067] In the implicit potential field equations that include fault step terms, the relevant data preprocessed in step S01 is input to fit the relevant parameters.

[0068] Step S03: Construct a dual linear equation system that includes formation interface constraints, gradient constraints, and drift term constraints.

[0069] Construct a bi-dual linear system of equations containing scalar constraints, gradient constraints, and drift term constraints, specifically as follows:

[0070]

[0071] in,

[0072] Let covariance be the covariance matrix between scalar constraint points. , , The number of scalar constraint points.

[0073] The scalar-gradient cross covariance matrix is... , , The number of gradient constraint points. Let be the normal vector of the i-th attitude observation point.

[0074] Let be the second-order covariance matrix between gradient constraint points. , This is the quadratic form of the Hessian matrix of the kernel function.

[0075] U and V are drift term constraint matrices; for first-order polynomial drift... , , ,in Let i be the three-dimensional coordinates of the i-th scalar constraint point; , ,in Let be the three-dimensional coordinates of the i-th scalar constraint point.

[0076] This is a vector of weight coefficients; Lagrange multipliers with gradient constraints; For the Lagrange multiplier of the drift term; This is a scalar observation vector, i.e., a formation interface constraint data vector; This is the gradient observation vector, i.e., the gradient constraint data vector.

[0077] The present invention employs the above-described method, which enables high-precision control under sparse data by simultaneously solving for scalar values ​​and gradient vectors in a matrix system.

[0078] Step S04: Solve for the weighting coefficients This generates a three-dimensional scalar field across the entire space.

[0079] For the matrix dimension of the aforementioned dual linear equation system, assume that... A scalar constraint point Given P gradient constraint points and P polynomial basis functions, the system size is: Choose a solution strategy based on the system data scale.

[0080] Preferably, when the matrix size is less than 1000, a direct method can be used to solve the problem, employing Cholesky decomposition to obtain the weight coefficients. Other parameters; when the matrix size is between 1000 and 10000, a direct method combined with sparse optimization can be used to solve the problem. The sparse Cholesky decomposition method is employed to obtain the weight coefficients. Other parameters; when the matrix size > 10000, an iterative method is used to solve the problem. For the iterative method, the Preconditioned Conjugate Gradient Method (PCG) is preferred, which uses incomplete LU decomposition as a precondition and uses GMRES iterative solution.

[0081] Step S05: Generate a three-dimensional geological model.

[0082] Based on the above steps, the potential field value is calculated, the background mesh is sampled, and isosurfaces are extracted using the Moving Cubes algorithm to finally generate a three-dimensional geological model and output a triangular mesh.

[0083] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be considered illustrative and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the description of the embodiments above. Therefore, all variations falling within the meaning and scope of equivalents of the claims are intended to be embraced within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims. Furthermore, it is clear that the word "comprising" does not exclude other units or steps, and the singular does not exclude the plural. Multiple units or devices recited in the system claims may also be implemented by a single unit or device in software or hardware. The terms "first," "second," etc., are used to indicate names and do not indicate any particular order.

Claims

1. A method for three-dimensional reconstruction of complex fault-block geological bodies based on joint constraints of covariance function and drift term, comprising the following steps: The basic stratigraphic data of the complex fault-block geological body is obtained, and the basic stratigraphic data is preprocessed; the basic stratigraphic data includes borehole data, surface attitude data, and fault data; scalar constraint data is obtained based on the borehole data, and gradient constraint data is obtained based on the surface attitude data. An implicit potential field equation including a fault step term is constructed, and the fault step term is constructed based on the acquired fault data; the fault step term is... ,in, This represents the total number of faults. For the first Vertical displacement of the fault plane; For the first The Heaviside step function of a fault is 1 when point x is located on the hanging wall and 0 when it is located on the footwall; specifically... , Let x be the sign distance from point x to the k-th fault plane; Based on the aforementioned basic stratigraphic data and the aforementioned implicit potential field equations, a dual linear equation system is constructed. Solving the dual linear equations generates a full-space three-dimensional scalar potential field, which in turn generates a three-dimensional geological model of the complex fault-block geological body.

2. The method according to claim 1, characterized in that: The implicit potential field equations containing fault step terms are specifically as follows: ； in, The equations are implicit potential field equations; Here, is the radial basis kernel function term; where, where is the weight coefficient of the i-th radial basis function; N is the total number of data constraint points; Here, r is the radial basis function kernel, and r is the Euclidean distance. , Let x be the three-dimensional spatial position vector of point x. Let i be the spatial location vector of the i-th data point; This is a polynomial drift term.

3. The method according to claim 2, characterized in that: Radial basis kernel function Use the Matérn kernel function, Gaussian kernel function, or cubic kernel function.

4. The method according to claim 2, characterized in that: Polynomial drift term It is represented by a first-order or second-order polynomial.

5. The method according to claim 2, characterized in that: The specific dual linear equation system is as follows: ； in, Let covariance be the covariance matrix between scalar constraint points. , , The number of scalar constraint points; The scalar-gradient cross covariance matrix is... , , The number of gradient constraint points. Let be the normal vector of the i-th attitude observation point; Let be the second-order covariance matrix between gradient constraint points. , ; U and V are drift term constraint matrices; for first-order polynomial drift... , , ,in Let i be the three-dimensional coordinates of the i-th scalar constraint point; , ,in Let J be the three-dimensional coordinates of the j-th scalar constraint point; This is a vector of weight coefficients; For Lagrange multipliers; For the Lagrange multiplier of the drift term; It is a scalar constraint vector; This is the gradient constraint vector.

6. The method according to claim 1, characterized in that: When solving the dual linear equation system, different solution methods are used depending on the size of the matrix.

Images