A method and system for 3D modeling of ore bodies based on multi-source exploration geological data

Abstract

This invention relates to the field of geological exploration and 3D geological modeling technology, specifically a method and system for 3D modeling of ore bodies based on multi-source exploration geological data. The method includes: acquiring multi-source spatial discrete data and spatial continuous gradient data within a preset 3D spatial range, and constructing a spatial voxel mesh based on the multi-source spatial discrete data and spatial continuous gradient data; calculating the spatial metric tensor matrix of each voxel node in the spatial voxel mesh to generate a spatial topological barrier field; using the source point of the target object as the initial zero-order jump set to generate an initial symbolic distance field; extracting the zero isosurface of the initial symbolic distance field to obtain the Euclidean deviation residual; if the Euclidean deviation residual exceeds a preset residual threshold, continuing until the Euclidean deviation residual does not exceed the residual threshold; when the Euclidean deviation residual does not exceed the residual threshold, generating a 3D triangular mesh of the target object. This invention achieves simultaneous continuous generation of the ore body surface and fault blocking.

Description

Technical Field

[0001] This invention relates to the field of geological exploration and three-dimensional geological modeling technology, specifically to a method and system for three-dimensional modeling of ore bodies based on multi-source exploration geological data. Background Technology

[0002] In deep ore body exploration and 3D modeling scenarios, the geological environment often exhibits characteristics of being significantly controlled by faults and sparsely revealed by boreholes, while also generating multi-source spatial discrete data, continuous gradient data, and topological fault surface data.

[0003] To construct 3D models of ore bodies, existing methods often employ direct distance field propagation or explicit geometric cutting techniques based on discrete source points. While these methods can generate preliminary ore body surfaces, they lack a unified fusion computing framework when faced with heterogeneous geological data with inconsistent semantics. On the one hand, relying solely on discrete source point evolution can lead to incorrect connectivity of the model across faults along the shortest Euclidean path. On the other hand, traditional explicit cutting can cause topological destruction and 3D mesh collapse, and it is difficult to perform closed-loop correction based on validation data, resulting in the final extracted manifold mesh failing to accurately reflect the tectonically constrained mineralization extension trend.

[0004] Therefore, how to unify multi-source exploration geological data and fault topological constraints into a computable implicit field framework, realize the continuous generation of grid surfaces, automatic fault blocking and error iterative correction, and thus improve the topological integrity and morphological accuracy of ore body 3D modeling, has become an urgent technical problem to be solved. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for three-dimensional modeling of ore bodies based on multi-source exploration geological data, which solves the following technical problems: avoiding the topological damage and mesh collapse caused by traditional explicit cutting, and effectively preventing cross-fault erroneous connectivity, thereby achieving the simultaneous completion of continuous generation of ore body surface, fault blocking and iterative correction.

[0006] The objective of this invention can be achieved through the following technical solutions: A method for 3D modeling of ore bodies based on multi-source exploration geological data includes: Acquire multi-source spatial discrete data and spatial continuous gradient data within a preset three-dimensional spatial range, wherein the multi-source spatial discrete data includes source points of the target object and discrete point data of the verification group, and construct a spatial voxel mesh based on the multi-source spatial discrete data and the spatial continuous gradient data; Calculate the spatial metric tensor matrix of each voxel node in the spatial voxel grid, and modify the eigenvalues ​​of the spatial metric tensor matrix of the voxel node at the position corresponding to the topological fracture surface data in the spatial voxel grid according to the preset topological fracture surface data, so as to generate a spatial topological obstacle field. The target object source point is used as the initial zero-step set to initialize the symbolic distance value for each voxel node. The symbolic distance value of each voxel node is updated by solving anisotropic partial differential equations in the spatial topological barrier field to generate the initial symbolic distance field. Extract the zero isosurface of the initial symbol distance field, calculate the nearest Euclidean distance between the points on the zero isosurface and the preset verification group discrete point data, and use the average or maximum value of the calculated nearest Euclidean distance as the Euclidean deviation residual. If the Euclidean deviation residual exceeds a preset residual threshold, the eigenvalues ​​of the spatial metric tensor matrix are corrected to adjust the anisotropy ratio, and the step of updating the symbolic distance values ​​of each voxel node by solving the anisotropic partial differential equation in the spatial topological obstacle field is re-executed to generate the current symbolic distance field. The step of extracting the zero isosurface of the current symbolic distance field is then executed until the Euclidean deviation residual does not exceed the residual threshold or the number of iterations reaches a preset upper limit. The initial symbolic distance field or the current symbolic distance field that satisfies the above conditions is then used as the target symbolic distance field. Manifold extraction is performed on the target symbol distance field to generate a three-dimensional triangular mesh of the target object.

[0007] As a further aspect of the present invention, the step of calculating the spatial metric tensor matrix of each voxel node in the spatial voxel mesh, and modifying the eigenvalues ​​of the corresponding spatial metric tensor matrix according to the preset topological fracture surface data includes: The main direction of the spatial metric tensor matrix of each voxel node is determined by using the gradient direction of the spatial continuous gradient data. Locate the corresponding position of the topological fracture surface data in the spatial voxel mesh; Calculate the normal direction of the topological fracture surface data at the corresponding position, and modify the eigenvalue of the spatial metric tensor matrix of the voxel node at the corresponding position in the direction parallel to the normal direction to a preset blocking eigenvalue, wherein the blocking eigenvalue is zero.

[0008] As a further aspect of the present invention, the steps of using the source point of the target object as the initial zero-order jump set, initializing the symbolic distance value for each voxel node, and updating the symbolic distance value of each voxel node by solving anisotropic partial differential equations in the spatial topological obstacle field to generate the initial symbolic distance field include: The source point of the target object is assigned the initial zero-order jump set, and the spatial voxel mesh is divided into a narrow band region and a background region according to the spatial distance from the initial zero-order jump set. The anisotropic partial differential equation is an anisotropic functional equation, and the parallel fast travel algorithm is used to solve the anisotropic functional equation in the narrow band region. The symbolic distance values ​​of each voxel node are updated along the fastest path defined by the anisotropic equation in the spatial topological barrier field until all voxel nodes have been updated, thereby generating the initial symbolic distance field.

[0009] As a further aspect of the present invention, the step of extracting the manifold from the target symbol distance field to generate a three-dimensional triangular mesh of the target object specifically includes: The moving cube algorithm is used to extract the isosurface of the target symbol range field; A watertight and non-self-intersecting 3D manifold structure is generated, and the 3D manifold structure is output as a 3D triangular mesh of the target object to a preset 3D rendering engine.

[0010] As a further aspect of the present invention, the multi-source spatial discrete data includes borehole point cloud data and geological boundary data, the spatial continuous gradient data includes geophysical inversion data, and the topological fault surface data includes geological fault plane data.

[0011] A 3D modeling system for ore bodies based on multi-source exploration geological data includes: The mesh construction module is used to acquire multi-source spatial discrete data and spatial continuous gradient data within a preset three-dimensional spatial range. The multi-source spatial discrete data includes source points of the target object and discrete point data of the verification group. A spatial voxel mesh is constructed based on the multi-source spatial discrete data and the spatial continuous gradient data. The tensor field generation module is used to calculate the spatial metric tensor matrix of each voxel node in the spatial voxel grid, and modify the eigenvalues ​​of the corresponding spatial metric tensor matrix according to the preset topological fracture surface data to generate a spatial topological barrier field. The distance field evolution module is used to take the source point of the target object as the initial zero-step set, initialize the symbolic distance value for each voxel node, and update the symbolic distance value of each voxel node by solving anisotropic partial differential equations in the spatial topological obstacle field to generate the initial symbolic distance field. The residual calculation module is used to extract the zero isosurface of the initial symbol distance field, calculate the nearest Euclidean distance between the points on the zero isosurface and the preset verification group discrete point data, and use the average or maximum value of the calculated nearest Euclidean distance as the Euclidean deviation residual. The manifold extraction module is used to, if the Euclidean deviation residual exceeds a preset residual threshold, correct the eigenvalues ​​of the spatial metric tensor matrix to adjust the anisotropy ratio, and re-execute the step of updating the symbolic distance values ​​of each voxel node in the spatial topological obstacle field by solving the anisotropic partial differential equation to generate the current symbolic distance field, and return to execute the step of extracting the zero isosurface of the current symbolic distance field until the Euclidean deviation residual does not exceed the residual threshold, and take the initial symbolic distance field or the current symbolic distance field when it does not exceed the residual threshold as the target symbolic distance field; When the Euclidean deviation residual does not exceed the residual threshold, manifold extraction is performed on the target symbolic distance field to generate a three-dimensional triangular mesh of the target object.

[0012] As a further aspect of the present invention, the tensor field generation module includes: The main direction determination unit is used to determine the main direction of the spatial metric tensor matrix of each voxel node using the gradient direction of the spatial continuous gradient data. A location unit is used to locate the corresponding position of the topological fracture surface data in the spatial voxel grid; The eigenvalue modification unit is used to calculate the normal direction of the topological fracture surface data at the corresponding position, and modify the eigenvalue of the spatial metric tensor matrix of the voxel node at the corresponding position in the direction parallel to the normal direction to a preset blocking eigenvalue, wherein the blocking eigenvalue is zero.

[0013] As a further aspect of the present invention, the distance field evolution module includes: A meshing unit is used to assign the source point of the target object to the initial zero-order set, and divide the spatial voxel mesh into a narrow band region and a background region according to the spatial distance from the initial zero-order set. The equation solving unit is used to solve the anisotropic equation within the narrow band region using a parallel fast travel algorithm. The anisotropic partial differential equation is the anisotropic equation. The distance update unit is used to update the symbolic distance value of each voxel node along the fastest path defined by the anisotropic equation in the spatial topological barrier field until all voxel nodes have been updated, so as to generate the initial symbolic distance field.

[0014] As a further aspect of the present invention, the manifold extraction module includes: The isosurface extraction unit is used to extract isosurfaces from the target symbol range field using the moving cube algorithm. The mesh output unit is used to generate a watertight and non-self-intersecting three-dimensional manifold structure, and output the three-dimensional manifold structure as a three-dimensional triangular mesh of the target object to a preset three-dimensional rendering engine.

[0015] As a further aspect of the present invention, the multi-source spatial discrete data includes borehole point cloud data and geological boundary data, the spatial continuous gradient data includes geophysical inversion data, and the topological fault surface data includes geological fault plane data.

[0016] The beneficial effects of this invention are: 1. This invention constructs a spatial voxel grid and calculates a spatial metric tensor matrix, using topological fault surface data as constraints, and modifies the tensor eigenvalues ​​in specific directions to zero to form a spatial topological barrier field. This mechanism unifies discrete and continuous data into an implicit field framework, effectively preventing cross-fault erroneous connectivity during distance field evolution, and realizing the simultaneous completion of continuous generation of the ore body surface and fault blocking.

[0017] 2. This invention introduces preset verification group discrete point data to calculate Euclidean deviation residuals. When the residuals exceed the threshold, the tensor matrix eigenvalues ​​are iteratively corrected to dynamically adjust the anisotropy ratio. This closed-loop correction mechanism can adaptively approximate the real mineralization extension trend constrained by structure, greatly improving the accuracy and objectivity of the model boundary in deep sparse borehole scenarios.

[0018] 3. In solving the anisotropic partial differential equation, this invention divides the voxel mesh into a narrow band region and a background region based on spatial distance, and uses a parallel fast traversal algorithm to solve the problem within the narrow band region. This strategy efficiently obtains the distance field that conforms to geological orientation in a large-scale space with reduced computational resource consumption, effectively avoiding the unnecessary computation cycle and memory consumption caused by global computation.

[0019] 4. After the residual meets the standard, the present invention uses the moving cube algorithm to extract the manifold from the symbolic distance field, directly generating a watertight and non-self-intersecting three-dimensional manifold structure. This conversion from implicit field to explicit mesh avoids the topological destruction and mesh collapse problems that are easily caused by traditional explicit geometric cutting. The output mesh can be seamlessly connected to the subsequent three-dimensional rendering and storage calculation process. Attached Figure Description

[0020] The invention will now be further described with reference to the accompanying drawings.

[0021] Figure 1 A flowchart illustrating a method for three-dimensional modeling of ore bodies based on multi-source exploration geological data, provided in this application embodiment; Figure 2 This is a schematic diagram of a module of a three-dimensional ore body modeling system based on multi-source exploration geological data provided in an embodiment of this application. Detailed Implementation

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

[0023] Please see Figure 1 A method for 3D modeling of ore bodies based on multi-source exploration geological data includes: acquiring multi-source spatial discrete data and spatial continuous gradient data within a preset 3D spatial range, wherein the multi-source spatial discrete data includes source points of the target object and discrete point data of the verification group, and constructing a spatial voxel grid based on the multi-source spatial discrete data and spatial continuous gradient data; calculating the spatial metric tensor matrix of each voxel node in the spatial voxel grid, and modifying the eigenvalues ​​of the spatial metric tensor matrix of the voxel node at the position corresponding to the topological fault surface data in the spatial voxel grid according to the preset topological fault surface data, so as to generate a spatial topological barrier field;

[0024] The target object source point is used as the initial zero-step set to initialize the symbolic distance value of each voxel node. The symbolic distance value of each voxel node is updated by solving the anisotropic partial differential equation in the spatial topological barrier field to generate the initial symbolic distance field. Extract the zero isosurface of the initial symbol distance field, calculate the nearest Euclidean distance between the points on the zero isosurface and the preset verification group discrete point data, and use the average or maximum value of the calculated nearest Euclidean distance as the Euclidean deviation residual. If the Euclidean deviation residual exceeds the preset residual threshold, the eigenvalues ​​of the spatial metric tensor matrix are corrected to adjust the anisotropy ratio, and the step of updating the symbolic distance values ​​of each voxel node by solving the anisotropic partial differential equation in the spatial topological obstacle field is re-executed to generate the current symbolic distance field. The step of extracting the zero isosurface of the current symbolic distance field is then executed until the Euclidean deviation residual does not exceed the residual threshold or the number of iterations reaches the preset upper limit. The initial symbolic distance field or the current symbolic distance field that meets the above conditions is taken as the target symbolic distance field. Manifold extraction is performed on the target symbolic distance field to generate a 3D triangular mesh of the target object.

[0025] This embodiment provides a 3D modeling mechanism for ore bodies based on multi-source exploration geological data. Specifically, this mechanism is designed for deep ore body scenarios that are clearly controlled by faults, have sparse borehole exposure but continuous geophysical coverage. It first constructs a spatial voxel grid, then forms a tectonically constrained spatial metric tensor field in the grid, evolves the symbolic distance field starting from the mineralization exposure point, and uses verification data to form a closed-loop correction. Finally, it outputs a 3D triangular mesh that can be directly entered into the 3D rendering or reserve calculation process. Collect multi-source data within a preset three-dimensional spatial range; this spatial range can be a cuboid bounding box of a deep prospecting target area in a mining area, such as a modeling area of ​​1200 meters east-west, 800 meters north-south, and 600 meters vertical; the received data includes at least two types: one is multi-source spatial discrete data, such as borehole mineralization points, geological boundaries formed by logging, trench control points, etc. Another type is spatially continuous gradient data, such as geophysical inversion bodies formed after interpolation; after the system performs unified coordinate transformation and scale normalization on various types of data, it establishes a spatial voxel grid according to a preset resolution; if the grid resolution is set to 20 meters, the aforementioned region can be discretized into 60×40×30 voxel units, and each voxel node is used to carry subsequent distance values ​​and tensor information. After the voxel grid is established, a spatial metric tensor matrix is ​​calculated for each voxel node. The tensor matrix is ​​used to express the spatial propagation weights in different directions. Specifically, along the stratigraphic distribution direction, the distance field has a propagation speed higher than the first preset tensor eigenvalue; perpendicular to the mineralization extension trend direction, the propagation speed is lower than the second preset tensor eigenvalue; if a fault plane is encountered, the propagation eigenvalue across the fault direction is modified to the preset blocking eigenvalue. Taking a simplified 3×3 matrix as an example: at node V1, if the main extension direction is east-west, then the matrix of this node can be represented as having a larger weight in the east-west direction and smaller weights in the other two directions. At node V2 near the fault, if the fault normal direction is close to the north-south direction, the eigenvalue corresponding to the north-south direction will be compressed to near zero or zero, thus forming a topological barrier at that location in that direction; after such processing, the entire space forms a spatial topological barrier field. The source point of the target object in the discrete data is used as the initial zero-step set to initialize the symbolic distance value for all voxel nodes. Here, the source point of the target object can be the midpoint of the mineralized section in a known mineralized borehole, or it can be a set of multiple mineralized contact points. During initialization, the distance value at the source point is set to 0, and the neighboring nodes are assigned positive, negative, or non-negative distance values ​​according to the initial geometric distance. Subsequently, anisotropic partial differential equations are solved in the aforementioned spatial topological barrier field, so that the sign distance value diffuses along the direction that is more geologically expansive, rather than penetrating along the fault blocking direction; The evolutionary calculation process is illustrated below: Assuming there are adjacent nodes A, B, and C, with A as the source point, there is no fault between A and B, and the distance between B and C is cut off by a fault. After the update, the distance between A and B increases smoothly, while the distance between C will not be rapidly shortened due to the presence of A and B. Instead, the zero isosurface stops or turns near B, rather than directly connecting to C through the fault. After obtaining the initial symbolic distance field, its zero isosurface is extracted to obtain the current ore body candidate surface; then the nearest Euclidean distance is calculated between this surface and the preset verification group discrete point data; the verification group discrete point data preferably uses another part of the borehole mineralization points, logging interface control points or subsequent supplementary exploration points that did not participate in the source point initialization, in order to ensure the objectivity of the verification; The verification group consists of three points: P1, P2, and P3. The closest distance from P1 to the zero isosurface is 8 meters, P2 is 11 meters, and P3 is 6 meters. Therefore, the average deviation residual can be calculated to be 8.33 meters, or the maximum value of 11 meters can be used as the evaluation index. If the residual exceeds the preset threshold, such as 10 meters, it indicates that the current anisotropy ratio is insufficient to approximate the actual ore body boundary, and the tensor matrix eigenvalues ​​of the corresponding region need to be corrected. Specifically, the system sets a basic adjustment step size such as 0.1, and calculates a dynamic correction coefficient based on the difference between the Euclidean deviation residual and the residual threshold at the verification point; Among them, the basic adjustment step size is a dimensionless parameter that is preset to control the convergence speed of the tensor matrix eigenvalue iteration; the residual threshold is an upper limit of the allowable error that is preset based on the mining area exploration accuracy requirements or the spatial voxel grid resolution. If a verification point is determined to be where the range field does not extend, then within its mapped local voxel neighborhood, the eigenvalue of the main extension direction is multiplied by (1 + base adjustment step size × the proportion exceeding the threshold), and the eigenvalue of the vertical direction is reduced accordingly to enhance the propagation accessibility in that direction; conversely, if overextension occurs, the reverse reduction is performed. For example, in regions where the Euclidean deviation residual exceeds the residual threshold, the eigenvalues ​​corresponding to the main extension direction are increased by a preset step size, while the eigenvalues ​​in the vertical diffusion direction are compressed to make the propagation more closely match the actual mineralization extension trend. After correction, the distance field evolution and zero isosurface extraction are re-executed until the residual does not exceed the threshold. The meaning of the anisotropic partial differential equation differs significantly in different steps: in the step of generating the initial symbolic distance field, the significance of solving the equation lies in establishing a global initial morphological benchmark with geological directionality based on the initially constructed tensor field; while in the iterative step of re-executing the solution after the residual exceeds the threshold, its meaning changes to providing an error-driven dynamic deformation mechanism, that is, using the spatial metric tensor matrix after correcting the eigenvalues ​​to guide the distance field diffusion in the local area, forcing the zero isosurface to undergo local deformation to fit the true verification group boundary. When the residual meets the requirements, manifold extraction is performed on the current symbolic distance field to form a three-dimensional triangular mesh of the target object; this triangular mesh can be further used for subsequent steps such as three-dimensional display, ore body volume calculation, and overlay analysis with tunnel engineering models. The mechanism for handling abnormal situations is as follows: during the data preparation stage, if some discrete data has missing coordinates or the hole depth is not closed, the data in question will first enter the repair queue and will not directly participate in the source point initialization. If the resolution of continuous gradient data is lower than the set grid resolution threshold, resampling and smoothing are performed first to avoid amplification of directional noise at the grid level. If certain regions lack both discrete control and continuous gradient support, the tensor matrix of these regions degenerates into an isotropic default matrix to avoid unfounded directional diffusion. If the residual still exceeds the threshold after the maximum number of iterations, the current model is output and high-uncertainty sub-regions are marked for priority deployment in subsequent resurveys. During the modeling process of a deep exploration area of ​​a skarn-type copper deposit in Southwest China, 12 boreholes have revealed copper mineralization. The midpoint of the mineralized section of 9 boreholes was used as the initial source point, and the other 3 boreholes were used as the verification group. Two near-NE-trending fault planes have been interpreted in the mining area, and electrical resistivity tomography (EDT) inversion bodies have been obtained. After unifying coordinates, the system constructs a three-dimensional voxel mesh, compresses the eigenvalues ​​of the normal direction of nodes near the fault plane to zero, and causes the sign distance field to expand along the interlayer and contact zone directions, while automatically stopping or bypassing at the fault; after the first round of extraction, the maximum deviation of the verification group was 18 meters; After the second round of increasing the anisotropy ratio in the interlayer direction in local areas, the maximum deviation was reduced to 9 meters, meeting the 10-meter threshold, and a closed orebody grid was output. The purpose of this step is to unify the originally semantically inconsistent discrete high-precision data, continuous low-frequency data and fault topological constraints into a computable implicit field framework, so as to achieve the simultaneous completion of continuous generation of ore body surface, fault blocking and iterative correction, and avoid the topological destruction and mesh collapse problems caused by traditional explicit cutting. To ensure the uniformity of technical solutions and the consistency of terminology, the eigenvalues ​​of the modified spatial metric tensor matrix are specifically governed by the following structured rules: the residual ratio is obtained based on the ratio of the portion of the Euclidean deviation residual that exceeds the residual threshold to the residual threshold. Multiply the residual ratio by the base adjustment step size to obtain the dynamic correction step size; if the range field does not extend to the verification point, increase the eigenvalue of the main extension direction by the proportion of the dynamic correction step size in the local voxel neighborhood, and correspondingly reduce the eigenvalues ​​of the other two vertical directions to increase the anisotropy ratio. If overextension occurs, decrease the eigenvalues ​​in the main extension direction and increase the eigenvalues ​​in the other two directions to reduce the anisotropy ratio. In this process, the mineralization exposure points, mineralization points, and midpoints of mineralization sections mentioned in the text are all strictly and uniformly mapped to the target object source points in the embodiments throughout the entire implementation process, ensuring the consistency of data reference throughout the text; through this structured rule, the calculation logic is clarified and the opaque propagation process is avoided.

[0026] In a preferred embodiment of the present invention, the steps of calculating the spatial metric tensor matrix of each voxel node in the spatial voxel grid and modifying the eigenvalues ​​of the corresponding spatial metric tensor matrix according to the preset topological fracture surface data include: determining the principal direction of the spatial metric tensor matrix of each voxel node using the gradient direction of the spatial continuous gradient data; locating the corresponding position of the topological fracture surface data in the spatial voxel grid; calculating the normal direction of the topological fracture surface data at the corresponding position; and modifying the eigenvalues ​​of the spatial metric tensor matrix of the voxel node located at the corresponding position in the direction parallel to the normal direction to preset blocking eigenvalues, wherein the blocking eigenvalues ​​are zero.

[0027] The embodiment provides a mechanism for generating a spatial metric tensor matrix and writing fracture constraints; specifically, in the aforementioned overall process, if the distance field is propagated only based on discrete source points, although a continuous surface can be generated, in the region where the spatial Euclidean distance on both sides of the fault is less than a preset distance threshold, the model may still incorrectly connect along the shortest Euclidean path, resulting in the problem of incorrect cross-fault connectivity. Therefore, this embodiment further introduces a processing method that determines the main direction based on continuous gradient data and performs eigenvalue blocking based on the normal direction of the fracture surface, so that spatial propagation conforms to the real geological topology; First, the gradient direction of the continuous gradient data in space is used to determine the main direction of each node. The continuous gradient data can be understood as a numerical change trend given at each point in three-dimensional space. The system reads the neighborhood window of each voxel node, such as the 3×3×3 neighborhood centered on the node, and obtains the change in the three directions by difference, and then combines them to form the local gradient direction. For regions with relatively stable gradients, directions perpendicular or parallel to the gradient can be mapped to the tensor principal axis directions to reflect the characteristic that mineralized bodies are more likely to extend along certain directions. If the changes in the three directions at node N1 are approximately 2 in the east-west direction, 1 in the north-south direction, and 0.2 in the vertical direction, the system can determine that there is a clear directionality near the node and construct a tensor ellipsoid with the principal axis biased towards the east-west direction. The main direction alone is not enough to solve the fault cutting problem; therefore, the system locates the corresponding position of the topological fault surface data in the voxel grid; the location can be achieved by projecting the fault surface triangle, wireframe or interpolated surface onto the voxel grid and marking the voxel nodes or voxel elements that are crossed by the fault surface. To reduce omissions, it is preferable to expand the fault plane along the normal direction by a preset width tolerance zone, such as half a voxel side length, so that nodes near the fault are also included in the constraint range. Next, the normal direction of the fracture surface at the corresponding position is calculated, and the eigenvalues ​​parallel to this direction are modified to the blocking eigenvalue of zero; the example calculation is as follows: the original tensor eigenvalues ​​of a certain node are set to 4, 2, and 1, which correspond to the propagation capabilities along the layer, along the strike, and vertically, respectively; If the fault normal direction is close to the second principal axis, the modified eigenvalues ​​of the node can become 4, 0, and 1, which means that the node no longer allows the range field to continue propagating in the direction of the second principal axis; similar modifications to adjacent nodes on the fault zone together form a continuous barrier, thereby naturally breaking the zero isosurface. The mechanism for handling abnormal situations is as follows: if the continuous gradient data in certain regions has large noise, causing the principal direction to frequently jump between adjacent nodes, then the gradient direction field is first locally smoothed, and then the tensor principal axis is calculated; if there is still no stable direction after smoothing, it degenerates into an isotropic tensor. If the fault surface data intersect or have local self-intersection, the system can prioritize sorting the fault surfaces according to the reliability of geological interpretation, and adopt the maximum blocking priority strategy for overlapping areas. That is, as long as any fault surface determines that the direction should be blocked, the feature value of that direction is compressed to zero. If the fracture surface passes through a voxel node but the normal direction is difficult to estimate stably, the weighted average of the normals of adjacent triangular pieces is used as an alternative. In the aforementioned skarn-type copper deposit in Southwest China, electrical resistivity inversion revealed a stable low-resistivity to high-polarizability gradient transition near the contact zone. After calculating the local gradient direction at each voxel node, the system found that the mineralization mainly extends along the near-NE-trending contact zone. Meanwhile, the geological interpretation results yielded two fault planes, F1 and F2. The system projected F1 and F2 onto the voxel grid and set the eigenvalues ​​to zero in the direction corresponding to the fault normal. In this way, even if the mineralization source point located on the eastern block of F1 is very close to the Euclidean distance of the control point on the western block, it will not be directly connected in the symbolic distance field. The purpose of this step is to simultaneously incorporate directional extension and rigid fault shearing into a unified tensor framework, thereby achieving the transformation from geometric proximity to geological topological proximity. In addition, it should be clearly pointed out that the topological fracture surface data in this embodiment and the faults and fault planes mentioned in the text strictly refer to the same geological entity that provides the blocking effect, so as to follow the consistency of terminology; in terms of the calculation rules for realizing the modification of eigenvalues, the principal direction of the spatial metric tensor matrix is ​​represented by a set of orthogonal eigenvectors, and each eigenvector corresponds to the propagation principal axis of each dimension of space; Modifying the eigenvalues ​​parallel to the normal direction to zero blocks the eigenvalues. The internal mechanism is to force the diagonal elements corresponding to the propagation weights of the eigenvector to be cleared to zero. Through this extreme structural decomposition method, the internal calculation logic of the tensor barrier field generation is clarified, ensuring the logical integrity and execution consistency of the algorithm scheme.

[0028] In a preferred embodiment of the present invention, the steps of taking the source point of the target object as the initial zero-order jump set, initializing the symbolic distance value of each voxel node, and updating the symbolic distance value of each voxel node by solving anisotropic partial differential equations in the spatial topological obstacle field to generate the initial symbolic distance field include: assigning the source point of the target object as the initial zero-order jump set, and dividing the spatial voxel mesh into a narrow band region and a background region according to the spatial distance from the initial zero-order jump set; The anisotropic partial differential equation is an anisotropic equation. The parallel fast traversal algorithm is used to solve the anisotropic equation in a narrow region. The symbolic distance value of each voxel node is updated along the fastest path defined by the anisotropic equation in the spatial topological barrier field until all voxel nodes are updated to generate the initial symbolic distance field.

[0029] This embodiment provides a mechanism for initialization and parallel evolution of the symbolic distance field. Specifically, after the tensor barrier field has been constructed, if a global uniform solution is directly performed on all voxel nodes, two problems will arise: first, the computational load is large and the contribution is limited in regions far from the mineralization control point; Secondly, in the large-volume mining area model, global iteration will cause the computation time and memory usage to exceed the preset resource limit; therefore, this embodiment uses narrow-band region partitioning, parallel fast travel solution and fastest path update to make the generation of the distance field both accurate and efficient. First, assign the source points of the target object to the initial zero-order set; if the mineralization source points are multiple discrete points, these points can be regarded as a set, and the distance value of each point in the set is set to 0; partition the voxel mesh according to the spatial distance from the set. The region whose spatial distance from the source point of the target object is less than or equal to the first preset distance threshold and is most likely to affect the zero contour surface morphology is defined as the narrow band region, and the remaining region whose spatial distance from the source point of the target object is greater than the first preset distance threshold is defined as the background region. If we set the voxel side length to 20 meters and the narrowband radius to 120 meters, then nodes that are no more than 6 voxels away from any source point are classified into the narrowband region. Through this classification, the system prioritizes fine-grained solution within the narrowband, while the background region can be updated coarsely or delayed. Within the narrow band region, the system employs a parallel fast traversal algorithm to solve the anisotropic equations; the system classifies node states into three types: determined, pending update, and unvisited; initially, all source node neighborhoods are added to the pending update set. The system selects nodes to be updated in ascending order of their current distance values, and calculates the new distance value of the node by combining the distance values ​​of the nodes around it that have been determined and the local tensor matrix. In this process, the local tensor matrix is ​​introduced as a weighting coefficient of the anisotropic equation. The specific derivation is as follows: Using the distance values ​​of the determined neighboring nodes, the distance gradient of the current node is estimated through the finite difference method, requiring that the inner product of this gradient vector and the local tensor matrix satisfy the evolution constant constraint; its specific anisotropic partial differential equation is expressed as follows: ,in Let be the spatial gradient vector of the distance value to be solved. The tensor matrix is ​​a local spatial metric. In the specific solution, the gradient components of the partial differential equation are discretized by combining the parallel fast traversal algorithm and the upwind finite difference scheme, so as to obtain the new distance value of the node to be updated. Since the tensor eigenvalues ​​of the main extension direction are large, the corresponding gradient components are penalized the least in the calculation, so the required distance increment in this direction is smaller. If the new value is smaller, it is written back and its adjacent nodes are added to the set to be updated; in parallel implementation, the nodes to be updated in different sub-blocks can be processed synchronously, and then a distance consistency merge is performed at the block boundary. Example calculation: There are nodes A, B, and C around the source point S. A is along the main extension direction, B is perpendicular to the main extension direction, and C is blocked by the fault normal. After the update, the distance value of A may change from infinity to 1.0, B to 2.5, while C remains at its maximum value or does not participate in the propagation. In this way, the same geometric distance will result in different propagation results under different tensor directions. After completing the iteration of the narrowband region, the system continues to update the symbolic distance values ​​of all nodes outward along the fastest path defined by the spatial topological barrier field; the fastest path here is not a Euclidean straight line, but the minimum cost path after combining the tensor matrix and the fracture blockage. The system classifies node states into three types: determined, pending update, and unvisited; once all nodes are updated, the complete initial symbolic distance field is obtained. The handling mechanism for abnormal situations is as follows: if the number of source points is less than the preset number threshold, resulting in the narrowband area coverage being less than the preset range threshold and unable to form a stable zero isosurface, the system will automatically expand the narrowband radius or supplement reliable points near the geological boundary as auxiliary source points. If a region is completely closed by a fault and cannot be reached from the existing source point, the distance value of that region will remain at a maximum value and will not be misidentified as a continuous ore body during subsequent zero isosurface extraction. If there is an update conflict between parallel blocks, the results are merged according to the principle of prioritizing the smaller distance value; if two candidate values ​​are the same, the path from the main extension direction is retained first to reduce jagged propagation. In the aforementioned mining area, the midpoints of the nine ore-bearing boreholes were set as zero-order jump sets; the voxel side length was 20 meters, and the narrow band radius was 100 meters. Since mineralization is controlled by the interlayer contact zone, the node update speed in the contact zone direction is significantly faster than that in the vertical direction. Near the F2 fault, the nodes marked as blocking directions maintain high distance values. After parallel rapid travel, the initial symbol distance field formed a continuous mineralization channel in the contact zone, but split into two unconnected sub-regions on both sides of the fault. The purpose of this step is to efficiently obtain the distance field that conforms to geological orientation in a large-scale voxel space with low computational resource consumption, thereby laying a stable foundation for subsequent zero isosurface extraction. Meanwhile, in order to clarify the specific calculation rules of the evolution process of the above anisotropic partial differential equations and to ensure the consistency of terminology definitions, the evolution constant constraint mentioned here is specifically manifested as follows: at any voxel node, the norm value of the spatial gradient vector of the distance value to be solved under the weighting effect of the spatial metric tensor matrix is ​​always equal to a preset constant such as 1. In this constraint, the preset constant represents the normalized value of the ideal expansion speed of the evolution interface, which is usually assigned a value of 1 to ensure the convergence, numerical stability and consistency of physical dimensions in the solution of partial differential equations. This deterministic calculation rule ensures that the propagation of distance values ​​is limited by local preset eigenvalues: the larger the eigenvalue, the smaller the distance increment accumulated in the calculation, and thus it can occupy a higher priority access order in the min-heap structure of the parallel fast traversal algorithm. Furthermore, in this embodiment, the narrowband region and the background region constitute a set of mutually exclusive and completely enclosing the entire spatial voxel mesh, eliminating undefined states in symbolic computation and ensuring strict consistency between global data flow and the physical meaning of the algorithm.

[0030] In a preferred embodiment of the present invention, the step of extracting the manifold of the target symbolic distance field to generate a three-dimensional triangular mesh of the target object specifically includes: The moving cube algorithm is used to extract the isosurface of the target symbol range field; Generate a watertight and non-self-intersecting 3D manifold structure, and output the 3D manifold structure as a 3D triangular mesh of the target object to the preset 3D rendering engine.

[0031] This embodiment provides a manifold extraction mechanism from the symbolic distance field to the explicit triangular mesh. Specifically, after the aforementioned iterative correction is completed, if only the voxel-level distance field is retained, although it is sufficient for internal judgment, it is not convenient for geological interpreters to conduct intuitive browsing, and it is also difficult to uniformly overlay with external models such as mine shafts, mining areas, and reserve boundaries. Therefore, this embodiment extracts the zero isosurface using the moving cube algorithm and further generates a watertight and non-self-intersecting three-dimensional manifold structure; The system iterates through each cube cell in the current symbol distance field and reads the distance values ​​of its 8 vertices. If there are both greater than zero and less than zero values ​​among the 8 vertices, it means that the zero isosurface passes through the cube cell, and the corresponding triangle needs to be generated by looking up the table based on the vertex symbol combination. An exemplary deduction can be made: the distance values ​​of the 8 vertices of a cube are respectively If a zero isosurface exists inside, the system interpolates the edges according to the standard topological template to find the intersection point. For example, if the distances between the two ends of an edge are -0.5 and 0.5, the intersection point is located at the midpoint of the edge; if they are -0.2 and 0.8, the intersection point is closer to the -0.2 end. After all cubes are processed, a complete set of triangular pieces is obtained. To ensure that the final mesh can be used for engineering analysis, a manifold check is also required. The system performs edge sharing consistency check, normal consistency check, and local self-intersection check on the generated triangular patch set. If an edge is used by only one triangular patch, it indicates that there is a crack, and the system can use the isosurface topology patch of adjacent voxels. If an edge is shared by more than two triangles, it indicates the existence of a non-manifold structure. The system splits it into multiple legal regions through local reconstruction. For adjacent triangles with inconsistent normals, the normal orientation is uniformly adjusted according to the gradient direction of the distance field. After processing, a watertight and self-intersecting 3D manifold structure is obtained and output to the preset 3D rendering engine. In handling abnormal situations, if the gradient of the distance field in a local area changes between adjacent nodes, it may cause the interior angle of the extracted triangular piece to be less than the preset angle threshold. In this case, smoothing with a preset intensity can be performed without changing the topology of the zero isosurface. However, the smoothing radius should be limited to one or two voxels to avoid re-smoothing the fault fault features. If the thickness of some thin-layered ore bodies is close to the side length of a voxel, the system can temporarily use sub-voxel interpolation or local subdivision of the grid for extraction to avoid missing thin-layer features. If multiple disconnected grid bodies are found before output, they can be grouped and identified by volume, source point affiliation, or spatial location to distinguish the main ore body from isolated small bodies. In the aforementioned mining area model, when the maximum residual drops to 9 meters, the system extracts isosurfaces from the current symbolic distance field. Since the distance fields on both sides of the F1 fault have been blocked, the grid generated by the moving cube appears as two closed ore bodies that are staggered from each other, rather than a continuous large block that has been cut later. Geologists can directly view fault fault relationships in the rendering engine and further overlay mining engineering wireframes and reserve zoning boundaries; The purpose of this step is to stably transform the implicit field results into an engineering-usable explicit triangular mesh, thereby enabling visualization analysis, subsequent calculations, and multi-model integration.

[0032] In a preferred embodiment of the present invention, the multi-source spatial discrete data includes borehole point cloud data and geological boundary data, the spatial continuous gradient data includes geophysical inversion data, and the topological fault surface data includes geological fault plane data.

[0033] This embodiment provides a multi-source data type mapping mechanism. Specifically, in the aforementioned modeling process, although it is helpful to abstractly describe discrete data, continuous data, and fractured data, it is still necessary to clarify the typical data sources and their roles in modeling when deploying them in practice. This embodiment further provides a set of specific data combinations suitable for ore body modeling in order to realize data access and engineering implementation. Multi-source spatial discrete data preferably includes borehole point cloud data and geological boundary data; borehole point cloud data can be obtained by converting borehole trajectory, sampling location in borehole, and mineralization start and end depth, where each discrete point has three-dimensional coordinates and can also carry grade, lithology or mineralization identifiers; Geological boundary data can come from tunnel logging lines, surface geological boundaries, contact zone interpretation lines, etc. These linear data can supplement geometric control in voxel grids; spatial continuous gradient data preferably include geophysical inversion data, such as resistivity, polarizability, magnetic susceptibility or density inversion data. This type of data has a large coverage and strong continuity, making it suitable as a background field for determining the main direction and the strength of anisotropy; the topological fault surface data preferably includes geological fault plane data, which can be derived from geological interpretation profiles, three-dimensional seismic interpretation surfaces or integrated geological modeling results, and is used to write rigid cut constraints. To illustrate how various types of data work together, an exemplary simulation can be performed: Suppose that there are mineralized points D1, D2, and D3 in the borehole point cloud, there is a contact line L1 in the geological boundary, the geophysical inversion body shows a continuous high gradient zone near L1, and the fault plane F1 cuts through this high gradient zone. The system can first set D1, D2, and D3 as source points, then determine the main direction around the contact zone by the inversion body, and set the local blocking direction to zero by F1; the distance field formed in this way is constrained by discrete mineralized points and guided by the continuous background field, and will not cross faults to connect incorrectly. In handling anomalies, if the spatial distribution of borehole point clouds is too unilateral, geological boundaries can be used as supplementary constraints to improve model stability; if there are local null values ​​in the geophysical inversion body, adjacent voxel interpolation or low confidence markers can be used instead of being directly used as the basis for strong principal direction. If the fault plane data only has two-dimensional profile lines and no complete three-dimensional surface, an initial fault plane can be formed by connecting the profiles first, and then blocking constraints can be written; if a certain type of data is temporarily missing, the system can still run on the existing data, but the model confidence level of that area should be reduced accordingly. In the aforementioned mining area project, the borehole point cloud came from 12 deep boreholes, the geological boundary came from the surface contact zone interpretation line and two tunnel logging lines, the continuous gradient data was obtained from the three-dimensional resistivity volume obtained by electrical inversion, and the fault plane data came from the F1 and F2 fault surfaces formed by comprehensive logging and profile interpretation. The system uniformly projects the above data to the same mining area coordinate system before entering the modeling process. The purpose of this step is to clarify the division of responsibilities among multi-source geological data within a unified modeling framework, thereby achieving standardized data access and consistency in subsequent calculations.

[0034] Please see Figure 2 A three-dimensional modeling system for ore bodies based on multi-source exploration geological data includes: a mesh construction module, used to acquire multi-source spatial discrete data and spatial continuous gradient data within a preset three-dimensional spatial range, wherein the multi-source spatial discrete data includes source points of the target object and discrete point data of the verification group, and constructs a spatial voxel mesh based on the multi-source spatial discrete data and spatial continuous gradient data; The tensor field generation module is used to calculate the spatial metric tensor matrix of each voxel node in the spatial voxel mesh, and modify the eigenvalues ​​of the corresponding spatial metric tensor matrix according to the preset topological fracture surface data to generate a spatial topological barrier field. The distance field evolution module is used to initialize the symbolic distance value of each voxel node by taking the source point of the target object as the initial zero-order jump set, and update the symbolic distance value of each voxel node by solving anisotropic partial differential equations in the spatial topological barrier field to generate the initial symbolic distance field. The residual calculation module is used to extract the zero isosurface of the initial symbol distance field, calculate the nearest Euclidean distance between the points on the zero isosurface and the preset verification group discrete point data, and use the average or maximum value of the calculated nearest Euclidean distance as the Euclidean deviation residual. The manifold extraction module is used to correct the eigenvalues ​​of the spatial metric tensor matrix to adjust the anisotropy ratio if the Euclidean deviation residual exceeds the preset residual threshold. It then re-executes the step of updating the symbolic distance values ​​of each voxel node by solving the anisotropic partial differential equation in the spatial topological obstacle field to generate the current symbolic distance field. Finally, it returns to the step of extracting the zero isosurface of the current symbolic distance field until the Euclidean deviation residual does not exceed the residual threshold. The initial symbolic distance field or the current symbolic distance field when the residual threshold is not exceeded is then used as the target symbolic distance field. When the Euclidean deviation residual does not exceed the residual threshold, manifold extraction is performed on the target symbolic distance field to generate a three-dimensional triangular mesh of the target object; This embodiment provides a software system or a hardware-software hybrid system for implementing the aforementioned method; specifically, the system can be deployed on a graphics workstation, server cluster, or dedicated equipment with parallel computing capabilities in a mining technology center, and completes a complete closed loop from data access to 3D mesh output in a modular manner; The grid construction module is responsible for data reading, format unification, coordinate transformation, and spatial voxel grid creation; its input can be connected to borehole databases, geological interpretation result libraries, and geophysical inversion result files. Its output generates a voxel node table with uniform numbering; the tensor field generation module reads the voxel node table and topological fault surface data, writes the metric tensor matrix at each node position, and performs blocking feature value modification on nodes near the fault, outputting a spatial topological barrier field. The distance field evolution module receives the target object source point set and the spatial topological obstacle field, initializes the symbolic distance value, and generates the initial or current symbolic distance field through a parallel solver; the residual calculation module extracts the zero isosurface from the distance field, calculates the nearest Euclidean distance by comparing it with the discrete point data of the verification group, and outputs the residual index. The manifold extraction module has a dual role: when the residual exceeds the threshold, it feeds back the correction parameters to the tensor field generation module and triggers re-evolution; when the residual meets the threshold, it directly performs manifold extraction on the current distance field and outputs a 3D triangular mesh. The following is an example of the data flow and coordination relationship between the various modules of the system: The mesh construction module outputs the node set G1; the tensor field generation module writes the tensor field T1 onto G1; the distance field evolution module generates the distance field D1 based on the source point set S1 and T1; the residual calculation module calculates the residual R1 from D1. If R1 is greater than the threshold, for example, R1 = 14 meters and the threshold is 10 meters, the manifold extraction module will not output the mesh for the time being, but will instead send back a set of local correction parameters C1, for example, increasing the principal direction eigenvalue by 20% in region Z1; The tensor field generation module generates T2 based on this, the distance field evolution module generates D2, and the residual calculation module obtains R2=8 meters; at this time, the manifold extraction module performs the final mesh output. In handling abnormal situations, if any module has missing input, the system can return a list of missing items. For example, if the residual calculation module does not receive the validation group data, the system can output only the initial model, but will mark the status as uncorrected. If the mesh building module detects that the input data coordinate system is inconsistent, it will pause the execution of subsequent modules to prevent the generation of pseudo-models on incorrect coordinate references. If the matrix output by the tensor field generation module contains singular or illegal values, it will automatically revert to the previous round of valid parameters and mark the abnormal region; if the manifold extraction module detects that the extraction result is empty, it means that the current zero isosurface has not formed a closed body. At this time, the system can prompt to expand the narrowband range, increase the source point, or reduce the local blocking strength. In the aforementioned mining project, the system was deployed on a workstation with graphics computing capabilities. Geologists first imported the borehole database, electrical resistivity inversion volume, and fault interpretation surface. The mesh construction module completed the three-dimensional discretization. The tensor field generation module wrote the direction and blocking information in the contact zone and fault neighborhood. The distance field evolution module generates the initial ore body; the residual calculation module uses subsequent exploration boreholes to check for deviations; after two rounds of correction, the manifold extraction module outputs the final ore body triangular mesh, which can be directly called by the mining 3D platform; The purpose of this mechanism is to break down the complex modeling process into collaborative, traceable, and iterative functional modules, thereby achieving stable deployment and automated operation in the engineering environment.

[0035] In a preferred embodiment of the present invention, the tensor field generation module includes: a main direction determination unit, used to determine the main direction of the spatial metric tensor matrix of each voxel node using the gradient direction of the spatial continuous gradient data; and a position positioning unit, used to locate the corresponding position of the topological fracture surface data in the spatial voxel grid. The eigenvalue modification unit is used to calculate the normal direction of the topological fracture surface data at the corresponding position, and modify the eigenvalue of the spatial metric tensor matrix of the voxel node at the corresponding position in the direction parallel to the normal to a preset blocking eigenvalue, wherein the blocking eigenvalue is zero.

[0036] This embodiment provides a refined structure for a tensor field generation module. Specifically, in a systematic implementation, if the generation of tensor fields is entirely handled by a single processing unit, two problems may arise in engineering applications: First, the estimation of continuous gradient directions differs from the accuracy requirements for fault location, and mixed processing is not conducive to debugging. Secondly, if the fault normal blocking is calculated incorrectly, it will directly affect the propagation of all subsequent distance fields; therefore, in this embodiment, the tensor field generation module is further divided into a main direction determination unit, a position positioning unit, and an eigenvalue modification unit. The principal direction determination unit extracts the changing trend around each node from continuous gradient data and outputs the principal direction result of the node; its input is the inverted volume or other continuum data, and the output can be a combination of the direction vector and the initial eigenvalue of each node. The location unit is responsible for mapping the fracture surface data to the voxel grid coordinates and outputting the set of nodes in the fault-affected area; the eigenvalue modification unit reads the fault normal direction from these node sets and sets the eigenvalues ​​in the tensor matrix that are parallel to the normal to zero, thereby forming a blockage. An exemplary deduction can be performed: the main direction determination unit outputs the main axis direction as vector U1 at node Q1; the position positioning unit determines that Q1 belongs to the F1 fault influence area; the feature value modification unit then reads the normal direction N1 of F1 at Q1. If N1 and U1 do not coincide, then blockage is performed on the axis corresponding to N1, while retaining high accessibility in the direction of U1. Thus, the final tensor at Q1 no longer simply propagates down the gradient, but propagates along the main geological direction and is prohibited from propagating across faults. In handling abnormal situations, if the output direction of the main direction determination unit is unstable, the position positioning unit and the feature value modification unit can still work independently to ensure that the fault blocking can at least take effect. If the location unit finds that the fault plane is tangent to the mesh boundary, causing some nodes to be marked repeatedly, then the nodes closer to the fault plane will be written first; if the eigenvalue modification unit detects that the angle between the normal direction and the existing principal axis direction is less than the preset angle threshold, which may cause tensor degradation to exceed the allowable range, then only one or two layers of neighborhood can be set to zero, and a preset non-zero small value can be used in the outer layers as a numerical stability transition. In the aforementioned mining area modeling system, the main direction determination unit identifies the contact zone extension direction based on the resistivity inversion volume, the location positioning unit maps the F1 fault surface to the voxel space, and the eigenvalue modification unit writes zero eigenvalues ​​in the fault normal direction. After the three are executed in sequence, the subsequent distance field evolution responds to the fault significantly better than the case that only depends on the gradient direction. The purpose of this mechanism is to improve the controllability and maintainability of tensor field generation by splitting the unit, thereby achieving precise superposition of directional guidance and fault blocking.

[0037] In a preferred embodiment of the present invention, the distance field evolution module includes: a mesh division unit, used to assign the source point of the target object to an initial zero-order jump set, and divide the spatial voxel mesh into a narrow band region and a background region according to the spatial distance from the initial zero-order jump set; The equation solving unit is used to solve anisotropic equations in a narrow region using a parallel fast-margin algorithm. The anisotropic partial differential equations are anisotropic equations. The distance update unit is used to update the symbolic distance value of each voxel node along the fastest path defined by the anisotropic equation in the spatial topological barrier field until all voxel nodes have been updated, so as to generate the initial symbolic distance field.

[0038] This embodiment provides a refined structure for the distance field evolution module; specifically, in system operation, if the distance field evolution exists only as a coarse-grained function as a whole, it is difficult to control the narrow band range, solution efficiency, and propagation path anomalies separately; Especially when the mining area expands or the number of source points increases, the computational resources required for solving the problem will increase non-linearly; therefore, in this embodiment, the distance field evolution module is further divided into a grid generation unit, an equation solving unit, and a distance update unit. The meshing unit is responsible for dividing the voxel space into narrowband regions and background regions based on the source point location; its output is not a simple binary label, but can also include distance hierarchy, such as the first ring, second ring, and third ring neighborhood, for priority sorting of solutions; The equation solving unit performs parallel fast movement within a narrow region to generate locally reliable distance values; the distance update unit then uses the solution results to propagate to a larger area and completes the global update according to the fastest path principle. An exemplary deduction is as follows: Assume that the source point set contains only two points, S1 and S2; the mesh generation unit marks the nodes that are no more than 5 voxels away from S1 and S2 as narrow bands, and the rest as background; The equation solving unit first calculates the distance values ​​of nodes A, B, and C in the narrow band as 0.8, 1.3, and unreachable, respectively. The unreachable state is caused by a fault blocking the path from A to C. The distance update unit continues to propagate outward based on this. If the background node D can be reached by detouring through A, it is updated to 2.6. If node E can only be reached through the blocked path, it is kept at its maximum value. In this way, a global distance field matching the geological obstacles can be formed. In handling abnormal situations, if the narrow band defined by the mesh generation unit is too narrow, resulting in insufficient accuracy at the edge of the zero isosurface, the system can automatically expand the narrow band according to the residual distribution of the previous round. If the narrowband is too wide and exceeds the memory budget, it is split into multiple subbands for solving. If the equation solving unit converges slowly in a local region, the distance update unit can temporarily freeze that region, process other stable regions first, and then backfill. If multiple equivalent candidates appear for the fastest path, the path that crosses fewer blocking areas and has a smoother path curvature is selected first. In the aforementioned mining area, as the number of subsequent exploration boreholes increased, the number of source points increased from 9 to 15. The system automatically expanded the narrow band range by dividing the grid into units, but did not allow all voxels to enter the high-precision solution. Instead, it only performed a fine solution near the contact zone where the source points were dense and the residuals were large. This maintained the modeling accuracy and avoided the extra time consumption caused by global solution. The purpose of this mechanism is to improve the solution efficiency and numerical stability of large-scale mining area models by focusing computational resources on key areas through the division of labor within the distance field evolution internal units.

[0039] In a preferred embodiment of the present invention, the manifold extraction module includes: The isosurface extraction unit is used to extract isosurfaces from the target symbol range field using the moving cube algorithm. The mesh output unit is used to generate a watertight and non-self-intersecting 3D manifold structure, and outputs the 3D manifold structure as a 3D triangular mesh of the target object to the preset 3D rendering engine.

[0040] This embodiment provides a refined structure for a manifold extraction module; specifically, in the aforementioned system closed loop, if isosurface extraction and mesh output are placed in the same processing unit, the source of the problem may be easily obscured between the two stages. For example, mesh cracks may originate from isosurface interpolation or from topology finishing before output; therefore, in this embodiment, the manifold extraction module is further divided into isosurface extraction units and mesh output units. The isosurface extraction unit is only responsible for extracting the zero isosurface from the current symbolic distance field and outputting the original triangular patch set; the mesh output unit is responsible for manifolding, validity checking and rendering interface adaptation of the original triangular patch set; in this way, geometry generation and engineering output can be separated. An exemplary deduction is as follows: The isosurface extraction unit generates triangular pieces T1, T2, and T3 in the local voxel block B1; after checking, the mesh output unit finds that there is a boundary gap between T2 and T4 in the adjacent block, so it merges T2 and T4 by reconstructing through shared vertices, and at the same time unifies the normals, and finally outputs a closed mesh block M1; multiple mesh blocks are then spliced ​​together to form a complete ore body mesh; In handling abnormal situations, if the number of local patches generated by the isosurface extraction unit is too small, it means that the distance field change in the region is insufficient to form a clear surface. In this case, the block can be temporarily not output and wait for the next round of correction. If the mesh output unit finds that the watertightness cannot be repaired by local patching, it can backtrack and indicate that there are abnormal fluctuations in the preceding distance field, requiring re-extraction or local recalculation; if the rendering engine requires a specific format, such as vertex normals, material groups or sub-mesh numbers, the mesh output unit will generate them during output. In the aforementioned mining project, the isosurface extraction unit first generates original triangular patch sets for the hanging wall and footwall of the fault, respectively. The mesh output unit detects a thin-film crack in a certain part of the hanging wall mesh because the interpolation points at the local voxel boundaries are repeated. After merging the shared vertices, two non-self-intersecting and closed ore body meshes are obtained and successfully loaded into the mining 3D platform for display. The purpose of this mechanism is to decouple zero isosurface extraction from engineering-grade mesh output, thereby improving problem localization capabilities and the usability of the final model.

[0041] In a preferred embodiment of the present invention, the multi-source spatial discrete data received by the system includes borehole point cloud data and geological boundary data, the spatial continuous gradient data includes geophysical inversion body data, and the topological fault surface data includes geological fault plane data.

[0042] This embodiment provides a system input data configuration scheme. Specifically, during the system deployment phase, different mining areas may have data with different levels of completeness. If there is a lack of input type constraints, it is easy to cause confusion of data roles when modules are called, such as mistaking fault interpretation lines as continuous fields or mistaking inversion volumes as discrete source points. Therefore, this embodiment defines the typical input components of the system as borehole point cloud data, geological boundary data, geophysical inversion body data, and geological fault plane data; After the borehole point cloud data is connected to the grid construction module, it mainly serves two purposes: providing source points and discrete verification. Geological boundary data, like borehole point cloud data, is a discrete constraint, but it is more inclined to supplement the shape and contour. Geophysical inversion volume data is connected to the tensor field generation module, which mainly serves to provide spatial continuous gradients. Geological fault plane data is integrated into the positioning and eigenvalue modification stage of the tensor field generation module, mainly undertaking the function of topological cutting; in this way, all types of data are given a clear responsibility from the moment they enter the system, avoiding the repeated misuse of the same data in different modules; An exemplary simulation can be performed: During system initialization, four input channels are established, denoted as C1, C2, C3, and C4 respectively; C1 receives borehole point clouds, C2 receives geological boundaries, C3 receives geophysical inversion bodies, and C4 receives fault planes; if a project only provides C1 and C3, the system can still model, but the output results will not have the rigid fault cutting effect. If C4 is present, a topological barrier field can be formed; if C2 is missing, the edge of the ore body may rely more on continuous gradient guidance in areas of discrete control sparseness. In handling abnormal situations, if the input data type does not match the expectation, such as mistakenly sending a two-dimensional profile image directly into the cross-sectional channel, the system should perform format verification during the import stage and prompt that it needs to be converted into a three-dimensional surface or point set first. If the data time versions of multiple channels are inconsistent, for example, the borehole database has been updated but the fault plane is still the old version, the system can prompt a version conflict to prevent the mixing of old and new interpretations; if a certain type of data is completely missing, the system will automatically reduce the function mode based on the remaining data and write the corresponding confidence label in the output model. When the aforementioned mining area modeling system went online, technicians imported mineralized point clouds from 12 boreholes from the database, contact zone boundaries from the geological interpretation platform, three-dimensional resistivity inversion bodies from geophysical software, and F1 and F2 fault planes from the structural interpretation results library. After the system completes registration through the established channels, the aforementioned modules can automatically identify the purpose of this data and successfully complete the 3D modeling of the ore body. The purpose of this mechanism is to achieve standardized access and stable operation of the modeling process by clearly defining the types and processing roles of the system's input data.

[0043] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for three-dimensional modeling of ore bodies based on multi-source exploration geological data, characterized in that, include: Acquire borehole point cloud data, geological boundary data, geophysical inversion body data and geological fault plane data within a preset three-dimensional spatial range. The borehole point cloud data includes the source point of the target object and the discrete point data of the verification group that did not participate in the source point initialization. Construct a spatial voxel grid based on the borehole point cloud data, the geological boundary data and the geophysical inversion body data. For each voxel node, the geophysical inversion volume data is differentially processed within a 3×3×3 neighborhood centered on that voxel node to obtain the changes in three spatial directions. Local gradient directions are then formed based on these changes, and the principal direction of the spatial metric tensor matrix of that voxel node is determined using these local gradient directions. When the gradient directions of adjacent voxel nodes frequently change, local smoothing is performed on the gradient direction field. If a stable direction cannot be obtained after smoothing, the spatial metric tensor matrix of the corresponding voxel node is set to an isotropic tensor. The geological fault plane data is projected onto the spatial voxel grid, and the geological fault plane data is expanded along the normal direction by a tolerance zone of half the voxel side length to mark the voxel nodes that the geological fault plane data passes through and is close to. When the normal direction of the geological fault plane is difficult to estimate stably, the weighted average of the normals of the triangular patches of adjacent fault planes is used as the normal direction at that location. Within one or two neighborhoods of the geological fault plane, the eigenvalues ​​parallel to the normal direction in the corresponding spatial metric tensor matrix are modified to zero, and the eigenvalues ​​parallel to the normal direction in the outermost layer are modified to preset non-zero small values ​​as a numerically stable transition, thereby generating a spatial topological barrier field. The target object source point is used as the initial zero-step set. The symbolic distance value at the target object source point is set to zero, and the neighboring voxel nodes are assigned positive or negative or non-negative distance values ​​according to the initial geometric distance. Based on the spatial distance to the source point of the target object, the spatial voxel mesh is divided into a narrow band region and a background region. In the narrow band region, the anisotropic equation is solved using the upwind finite difference scheme and the parallel fast travel algorithm. The symbol distance value of each voxel node is updated along the fastest path defined by the anisotropic equation in the spatial topological obstacle field to generate an initial symbol distance field. The expression for the anisotropic equation is: ; in Let be the spatial gradient vector of the distance value to be solved. Let be the local spatial metric tensor matrix; the norm of the spatial gradient vector of the distance value to be solved under the weighting effect of the spatial metric tensor matrix is ​​equal to 1; Voxel units are identified in the initial symbolic distance field where 8 vertices simultaneously have distance values ​​greater than zero and less than zero. Based on the vertex symbolic combination of the voxel units, zero-value intersection points are obtained by interpolation on the edges between vertices with different symbolic distance values, and zero isosurfaces are extracted based on the zero-value intersection points. Calculate the nearest Euclidean distance between the discrete points of each verification group and the zero isosurface, and use the average or maximum value of the nearest Euclidean distance as the Euclidean deviation residual; When the Euclidean deviation residual exceeds a preset residual threshold, the residual ratio is obtained based on the ratio of the portion of the Euclidean deviation residual that exceeds the residual threshold to the residual threshold, and the residual ratio is multiplied by a preset basic adjustment step size to obtain the dynamic correction step size. For cases where the symbolic distance field is determined not to extend to the corresponding validation set discrete point, the eigenvalues ​​in the main extension direction are increased according to the dynamic correction step size within the local voxel neighborhood mapped to the validation set discrete point, and the eigenvalues ​​in the other two perpendicular directions are correspondingly decreased. For cases where the symbolic distance field is determined to be overextended, the eigenvalues ​​in the main extension direction are decreased according to the dynamic correction step size within the corresponding local voxel neighborhood, and the eigenvalues ​​in the other two perpendicular directions are increased. The anisotropic equation is re-solved based on the corrected spatial metric tensor matrix to generate the current symbolic distance field. The extraction of the zero isosurface and the calculation of the Euclidean deviation residual are re-executed until the Euclidean deviation residual does not exceed the residual threshold or the number of iterations reaches the preset upper limit. The initial symbolic distance field or the current symbolic distance field that meets the above conditions is taken as the target symbolic distance field. Manifold extraction is performed on the target symbol distance field to generate a three-dimensional triangular mesh of the target object.

2. The method for three-dimensional modeling of ore bodies based on multi-source exploration geological data according to claim 1, characterized in that, The steps of extracting the manifold from the target symbol range field to generate a 3D triangular mesh of the target object specifically include: The moving cube algorithm is used to extract the isosurface of the target symbol range field; A watertight and non-self-intersecting 3D manifold structure is generated, and the 3D manifold structure is output as a 3D triangular mesh of the target object to a preset 3D rendering engine.

Images