A method, system, device and medium for automatically modeling ore body morphology based on three-dimensional topological reconstruction

Abstract

The present application relates to a kind of ore body form automatic modeling method, system, equipment and medium based on three-dimensional topological reconstruction, method includes: to discrete drilling data is spatialized, according to boundary grade threshold value generates internal point set, boundary point set and wall rock point set;With three kinds of point set respectively exert positive, zero, negative constraint, construct radial basis function implicit potential field and solve coefficient to generate potential function;Potential function is carried out isovolume voxel level tracking, generates initial triangular mesh;Initial triangular mesh is successively filtered to connected component noise, local encryption based on skeleton extraction is extracted in branch and pinch-out area, and potential field is directed hole filling, generates manifold triangular mesh;With manifold triangular mesh as spatial constraint domain, establish block model, based on grade data carries out spatial variation analysis and Kriging interpolation to generate block grade model, and weighted accumulation obtains ore body total volume, ore quantity and metal quantity, output three-dimensional entity model and resource reserve evaluation report.

Description

Technical Field

[0001] This invention relates to the field of electronic digital data processing technology, specifically to an automated modeling method, system, equipment, and medium for ore body morphology based on three-dimensional topological reconstruction. Background Technology

[0002] In solid mineral exploration and mining design, constructing a three-dimensional digital model of the ore body using discrete borehole data is fundamental for resource reserve estimation and mining plan development. Traditionally, geological engineers employ explicit modeling methods that combine manual profile drawing with contour line stitching. This involves projecting borehole sampling points onto multiple profiles, manually drawing the ore body boundary lines, and then fitting these two-dimensional contour lines to a surface in three-dimensional space to form a solid model of the ore body. In recent years, implicit surface modeling technology has been gradually introduced into this field. Its core idea is to treat the ore body boundary as an isosurface of a certain three-dimensional potential function. By using constraint values ​​at discrete sampling points (such as assigning positive values ​​to internal points and negative values ​​to external points), the potential function is solved inversely, thereby automatically extracting the ore body surface. This type of method can handle large-scale data and reduce manual interaction.

[0003] However, existing orebody modeling methods based on implicit surfaces still have significant drawbacks. On the one hand, the potential field constructed by traditional radial basis functions or kriging methods imposes equal-weight constraints on all training points, which can easily lead to oversmoothing or oscillations in sparse borehole regions, resulting in orebody boundary distortion. On the other hand, the extracted initial isosurface mesh often has topological defects such as non-manifold edges, fine noise debris, branching fractures, and voids. Especially for orebodies with complex morphologies such as branching and pinch-out recurrence, existing methods struggle to automatically maintain correct topological connectivity, requiring substantial manual correction. Furthermore, geometric morphology modeling and internal grade space interpolation are often separated into two independent modules. The solid model output by the former lacks strict spatial consistency constraints between the block model established by the latter, which can easily lead to errors in resource quantity calculation. Summary of the Invention

[0004] Based on this, the purpose of this invention is to provide an automated modeling method, system, equipment, and medium for ore body morphology based on three-dimensional topology reconstruction, which can automatically perform three-dimensional topology reconstruction and achieve integrated modeling of geometric shape and grade distribution from discrete borehole data.

[0005] The objective of this invention is achieved through the following solution:

[0006] In a first aspect, the present invention provides an automated modeling method for ore body morphology based on three-dimensional topological reconstruction, comprising the following steps:

[0007] S1: Spatialize the acquired discrete borehole data and classify the sampling points into internal points of the ore body, boundary crossing points, and surrounding rock points according to the boundary grade threshold, generating internal point sets, boundary point sets, and surrounding rock point sets;

[0008] S2: Apply positive constraints to points in the internal point set, zero constraints to points in the boundary point set, and negative constraints to points in the surrounding rock point set to construct a radial basis function implicit potential field. Solve for the radial basis function coefficients and drift term coefficients of the radial basis function implicit potential field to generate the potential function.

[0009] S3: Call the moving cube algorithm to perform voxel-level tracking of the potential function at the isosurface level. For each cube cell that intersects with the isosurface in the modeling space, determine the connection method of the triangular facets inside the cube cell by looking up the positive and negative combinations of the potential function values ​​at the eight vertices of the cube cell, and generate an initial triangular mesh that reflects the boundary of the ore body.

[0010] S4: Perform 3D topological reconstruction on the initial triangular mesh. This involves removing noise fragments through connected component analysis, refining and reconstructing branches and pinch-out regions locally through skeleton extraction, filling hole boundaries through gradient-guided iterative projection of potential functions, and reconstructing local potential fields and mesh patches multiple times to transform the non-manifold mesh into a manifold triangular mesh in which each edge is shared by at most two triangular patches, thus generating a manifold triangular mesh.

[0011] S5: Using a manifold triangular mesh as the spatial constraint domain, a block model is established within the spatial constraint domain. Based on the grade data within the internal point set and boundary point set, spatial variation analysis is performed to obtain the variation function parameters. The spatial interpolation method is called to estimate the grade of the center point of each block, generating a block grade model. Based on the block grade model, all valid blocks are traversed to perform a weighted summation of volume and grade, generating a three-dimensional solid model of the ore body and a resource reserve assessment report.

[0012] Secondly, the present invention provides an automated ore body morphology modeling system based on three-dimensional topological reconstruction, which is configured with the following modules:

[0013] The borehole data classification module is used to spatialize the acquired discrete borehole data. Based on the boundary grade threshold, the sampling points are classified into points inside the ore body, boundary crossing points, and surrounding rock points, generating internal point sets, boundary point sets, and surrounding rock point sets.

[0014] The implicit potential field construction module is used to construct a radial basis function implicit potential field by applying positive constraints to points in the internal point set, zero constraints to points in the boundary point set, and negative constraints to points in the surrounding rock point set, and to solve for the radial basis function coefficients and drift term coefficients of the radial basis function implicit potential field to generate the potential function.

[0015] The orebody mesh initialization module is used to call the moving cube algorithm to perform voxel-level tracking of the potential function at the isosurface level. For each cube cell that intersects with the isosurface in the modeling space, the connection method of the triangular facets inside the cube cell is determined by a lookup table based on the positive and negative combinations of the potential function values ​​at the eight vertices of the cube cell, thereby generating an initial triangular mesh that reflects the orebody boundary.

[0016] The mesh topology optimization module is used to perform three-dimensional topology reconstruction on the initial triangular mesh. It removes noise fragments through connected component analysis, performs local densification and reconstruction of branches and pinch-out regions through skeleton extraction, fills the hole boundaries through iterative projection guided by potential function gradient, and transforms the non-manifold mesh into a manifold triangular mesh in which each edge is shared by at most two triangular patches through multiple local potential field reconstructions and iterative reconnection of mesh patches, thus generating a manifold triangular mesh.

[0017] The orebody reserve modeling module is used to build block models within a spatially constrained domain using a manifold triangular mesh. Based on the grade data within the internal and boundary point sets, spatial variation analysis is performed to obtain variation function parameters. Spatial interpolation methods are then used to estimate the grade at the center point of each block, generating a block grade model. Finally, based on the block grade model, all valid blocks are traversed to perform a weighted summation of volume and grade, generating a three-dimensional solid model of the orebody and a resource reserve assessment report.

[0018] Thirdly, this application provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement any of the above-mentioned automated modeling methods for ore body morphology based on three-dimensional topological reconstruction.

[0019] Fourthly, this application provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements any of the above-mentioned automated modeling methods for ore body morphology based on three-dimensional topological reconstruction.

[0020] In summary, the automated orebody morphology modeling method based on 3D topology reconstruction provided in this application, based on the implicit potential field construction of radial basis functions and adaptive constraint configuration, can achieve stable fitting of the boundary of sparse borehole regions, avoiding boundary distortion caused by equal-weight interpolation in traditional methods, thereby improving the geometric fidelity of orebody morphology reconstruction. Secondly, through connected component noise filtering, skeleton extraction-guided branching and pinch-out region local refinement, and potential function gradient-driven hole filling, automated topology repair of non-manifold meshes can be achieved, enabling the final generated manifold triangular mesh to accurately maintain the topological continuity of complex orebodies such as branching and pinch-out reproduction, thereby eliminating… In addition to the tediousness and subjective errors of manual repair, this method uses manifold triangular meshes as spatial constraint domains to establish block models. By combining spatial variation analysis and kriging interpolation, it can achieve strict consistency between geometric morphology models and attribute grade models on a spatial basis, thus avoiding resource quantity calculation deviations caused by the separation of geometric entities and block grade models in traditional processes. Finally, through weighted accumulation of volume and grade, a three-dimensional entity model and resource reserve assessment report can be automatically output, achieving integrated processing from discrete borehole data to high-precision ore body digital models. This significantly improves modeling efficiency and assessment reliability in the exploration and mining design of complex ore deposits.

[0021] To better understand and implement this invention, the following detailed description is provided in conjunction with the accompanying drawings. Attached Figure Description

[0022] Figure 1 A flowchart illustrating an automated ore body morphology modeling method based on three-dimensional topological reconstruction, provided for an embodiment of this application;

[0023] Figure 2 This is a schematic diagram of the structure of an automated ore body morphology modeling system based on three-dimensional topology reconstruction, provided as another embodiment of this application. Detailed Implementation

[0024] To facilitate understanding of the present invention, a more complete description will be given below with reference to the accompanying drawings. Preferred embodiments of the invention are shown in the drawings. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a thorough and complete understanding of the disclosure of the invention.

[0025] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0026] In one embodiment, such as Figure 1 As shown, an automated ore body morphology modeling method based on 3D topological reconstruction is provided. This embodiment illustrates the application of this method to a terminal. It is understood that this method can also be applied to a server, and to a system including both a terminal and a server, and implemented through interaction between the terminal and the server. In this embodiment, the method includes the following steps:

[0027] S1: Spatialize the acquired discrete borehole data and classify the sampling points into internal points of the ore body, boundary crossing points, and surrounding rock points according to the boundary grade threshold, generating internal point sets, boundary point sets, and surrounding rock point sets.

[0028] Specifically, discrete borehole data collected during solid mineral exploration is acquired. This discrete borehole data includes at least the spatial coordinates of each borehole (X, Y, and Z three-dimensional coordinates, where the Z coordinate corresponds to the borehole depth, uniformly calibrated using a geodetic coordinate system), grade data of the borehole sampling points (such as the content of useful components), borehole number, and sampling depth. Since borehole data may contain coordinate deviations, missing values, and outliers during field exploration, spatial preprocessing is necessary. This preprocessing includes three steps: coordinate calibration, outlier removal, and data completion. Coordinate calibration uses a combination of GPS positioning calibration and comparison with known geological control points to uniformly convert all borehole data to the same spatial coordinate system (such as the National Geodetic Coordinate System 2000), eliminating coordinate deviations caused by different exploration stages and different acquisition equipment. Outlier removal uses a 3D model... The criteria involve statistical analysis of the grade data of the sampling points to remove outliers exceeding the mean ± 3 standard deviations, thus avoiding the impact of outliers on the accuracy of subsequent modeling. For missing sampling point data, linear interpolation of adjacent sampling points is used to supplement the missing data, ensuring the continuity and integrity of the borehole data.

[0029] The system determines boundary grade thresholds, which are based on the geological exploration report of the mining area, industry standards, and ore industrial utilization standards. The system classifies pre-processed sampling points according to these thresholds, dividing them into points within the ore body, boundary crossing points, and surrounding rock points. Sampling points whose grades meet the ore body determination criteria are grouped into an internal point set, which characterizes the spatial location and grade characteristics of the ore body's internal region. Adjacent sampling points within the same borehole whose grades meet both the ore body and surrounding rock determination criteria are grouped into a boundary point set, which constrains the spatial location of the ore body boundary. Sampling points whose grades meet the surrounding rock determination criteria are grouped into a surrounding rock point set, which characterizes the surrounding rock region outside the ore body, forming a reverse constraint with the internal points.

[0030] S2: Apply positive constraints to points in the internal point set, zero constraints to points in the boundary point set, and negative constraints to points in the surrounding rock point set to construct a radial basis function implicit potential field. Solve for the radial basis function coefficients and drift term coefficients of the radial basis function implicit potential field to generate the potential function.

[0031] Specifically, the system calls the generated internal point set, boundary point set, and surrounding rock point set, applying different potential constraints to the three types of point sets. The system applies positive constraints to all sampling points in the internal point set, setting the potential function value of these points to a positive number, indicating that the area is inside the ore body. The system applies zero constraints to all sampling points in the boundary point set, setting the potential function value of these points to zero, indicating that the area is on the boundary of the ore body. The system applies negative constraints to all sampling points in the surrounding rock point set, setting the potential function value of these points to a negative number, indicating that the area is outside the ore body. These three types of constraints form a spatial distribution gradient of the potential field, ensuring that the potential function can accurately distinguish between the ore body and the surrounding rock. The system selects radial basis functions to construct the implicit potential field; the type of radial basis function is determined based on the distribution characteristics of the borehole data in the mining area and the modeling accuracy requirements. The system constructs a system of linear equations for the potential function, with the potential function expression including radial basis function terms and drift terms, the drift terms being in polynomial form. The system solves the system of linear equations to obtain the coefficients of the radial basis functions and the coefficients of the drift terms. The system substitutes the obtained coefficients into the potential function expression to generate the potential function of the complete radial basis function implicit potential field. This potential function can reflect the relative positional relationship between any point in space and the ore body boundary. Regions with potential function values ​​greater than zero correspond to the interior of the ore body, regions with potential function values ​​equal to zero correspond to the ore body boundary, and regions with potential function values ​​less than zero correspond to the surrounding rock region.

[0032] S3: Call the moving cube algorithm to perform voxel-level tracing of the potential function at the isosurface level. For each cube cell that intersects with the isosurface in the modeling space, determine the connection method of the triangular facets inside the cube cell by looking up the positive and negative combinations of the potential function values ​​at the eight vertices of the cube cell, and generate an initial triangular mesh that reflects the boundary of the ore body.

[0033] Specifically, the system defines the modeling space based on the exploration area and borehole distribution of the mining area. The system discretizes the modeling space into cubic units, with the unit division criteria determined by balancing modeling accuracy and computational efficiency. The system calculates the potential function values ​​at the eight vertices of each cubic unit, using the potential functions generated by S2. The system determines whether each cubic unit intersects with a zero isosurface based on the positive and negative combinations of the potential function values ​​at its eight vertices. If a vertex has both positive and negative potential function values, then the cubic unit intersects with the zero isosurface. The system uses the moving cube algorithm to trace the isosurface for each cubic unit intersecting with the zero isosurface. The system pre-constructs a vertex positive / negative combination lookup table, storing the triangular facet connection methods corresponding to each vertex positive / negative combination.

[0034] The system matches the corresponding triangular facet connection method in a lookup table based on the positive and negative combinations of the eight vertices of the cube element. It then calculates the coordinates of the intersection points of the isosurface and the edges of the cube element using linear interpolation; these intersection coordinates are the vertex coordinates of the triangular facet. The system connects the intersection points according to the connection order specified in the lookup table to generate triangular facets. This process is repeated for all intersecting cube elements, and all generated triangular facets are aggregated to form an initial triangular mesh reflecting the boundary contour of the ore body. This initial triangular mesh can preliminarily characterize the three-dimensional morphology of the ore body; however, due to the discreteness of the borehole data and the small fluctuations in the potential field, the initial triangular mesh contains topological defects such as non-manifold edges, small noise fragments, branch fractures, and voids.

[0035] S4: Perform 3D topological reconstruction on the initial triangular mesh. This involves removing noise fragments through connected component analysis, refining and reconstructing branches and pinch-out regions locally through skeleton extraction, filling hole boundaries through gradient-guided iterative projection of potential functions, and reconstructing local potential fields and mesh patches multiple times. This transforms the non-manifold mesh into a manifold triangular mesh where each edge is shared by at most two triangular patches, thus generating a manifold triangular mesh.

[0036] Specifically, the system calls the generated initial triangular mesh and performs 3D topology reconstruction on it, gradually eliminating topological defects. The system uses a region growing method to perform connected component analysis on the initial triangular mesh, dividing it into multiple independent connected components based on the adjacency relationships of triangular faces. The system sets a criterion for determining connected components, removing noise fragments based on this criterion and retaining qualified connected components as the basis for subsequent topology reconstruction. The system uses a distance transform method to extract the central skeleton of the retained connected components. This central skeleton can characterize the branching trend, pinch-out location, and extension trend of the ore body. The system adds sampling points around key nodes of the skeleton, calculates the potential value of the new sampling points using the potential function in S2, and adds them to the boundary point set. Based on the supplemented boundary point set, the system reconstructs the local radial basis function potential field, performs local refinement and fitting on the triangular mesh of the branching and pinch-out regions, and repairs branch fracture defects.

[0037] Furthermore, the system identifies void regions in the initial triangular mesh, which are areas with unclosed boundaries. The system calculates the gradient direction of the potential function at each vertex on the void boundary, with the gradient direction pointing in the direction of increasing potential function value. The system iteratively projects the vertices on the void boundary along their gradient directions, calculating the potential function value of each vertex after projection until the vertex's potential function value reaches a set standard. This vertex is then used as a new vertex for void filling, and the above process is repeated until the void boundary is closed. The system identifies local regions containing non-manifold edges, extracts internal points, boundary points, and surrounding rock points within these regions, reconstructs the local radial basis function potential field, and adjusts the potential function distribution in these regions. Based on the adjusted local potential field, the system iteratively reconnects the triangular patches surrounding the non-manifold edges, adjusting the vertex connection relationships of the patches to ensure that each edge is shared by at most two triangular patches. The system repeats the above local reconstruction process until there are no topological defects in the entire mesh, generating a complete, continuous, and topologically correct manifold triangular mesh.

[0038] S5: Using a manifold triangular mesh as the spatial constraint domain, a block model is established within the spatial constraint domain. Based on the grade data within the internal point set and boundary point set, spatial variation analysis is performed to obtain the variation function parameters. The spatial interpolation method is called to estimate the grade of the center point of each block, generating a block grade model. Based on the block grade model, all valid blocks are traversed to perform a weighted summation of volume and grade, generating a three-dimensional solid model of the ore body and a resource reserve assessment report.

[0039] Specifically, the system invokes a manifold triangular mesh, using it as the spatial constraint domain. Within this domain, a block model is established, employing regular cubic blocks. The block division criteria are determined based on the required grade interpolation accuracy. The system uses a spatial pruning algorithm to remove blocks exceeding the manifold triangular mesh, retaining only valid blocks within the constraint domain to ensure spatial consistency between the block model and the ore body geometry. The system then calls grade data from internal and boundary point sets for spatial variation analysis. A spherical model is used as the variation function model. By calculating the experimental variation function of the grade data, variation function parameters are fitted. These parameters reflect the spatial distribution pattern of the grade data, providing a basis for subsequent spatial interpolation. Finally, the system uses a spatial interpolation method to estimate the grade at the center point of each block. This method is determined based on the distribution characteristics of the ore area's grade data.

[0040] The system uses grade data from internal and boundary point sets as sample points, and the potential function value generated by S2 as an auxiliary variable. Substituting this value into the spatial interpolation formula, it calculates the grade estimate for the center point of each valid block. This estimate is then used as the grade value for the corresponding block, generating a complete block grade model. The system iterates through all valid blocks, calculates the volume of each block, and weights and sums the volume of each block with its corresponding grade value to obtain the total resource volume of the ore body. Based on the block grade model, the system uses 3D rendering technology to generate a 3D solid model of the ore body. This 3D solid model is used to display the spatial morphology, grade distribution patterns, and boundary characteristics of the ore body. The system organizes various parameters from the modeling process, block model data, resource reserve calculation results, and error analysis data to generate a standardized resource reserve assessment report, providing data support for mining area resource reserve estimation and mining plan formulation.

[0041] In summary, the automated orebody morphology modeling method based on 3D topology reconstruction provided in this application, based on the implicit potential field construction of radial basis functions and adaptive constraint configuration, can achieve stable fitting of the boundary of sparse borehole regions, avoiding boundary distortion caused by equal-weight interpolation in traditional methods, thereby improving the geometric fidelity of orebody morphology reconstruction. Secondly, through connected component noise filtering, skeleton extraction-guided branching and pinch-out region local refinement, and potential function gradient-driven hole filling, automated topology repair of non-manifold meshes can be achieved, enabling the final generated manifold triangular mesh to accurately maintain the topological continuity of complex orebodies such as branching and pinch-out reproduction, thereby eliminating… In addition to the tediousness and subjective errors of manual repair, this method uses manifold triangular meshes as spatial constraint domains to establish block models. By combining spatial variation analysis and kriging interpolation, it can achieve strict consistency between geometric morphology models and attribute grade models on a spatial basis, thus avoiding resource quantity calculation deviations caused by the separation of geometric entities and block grade models in traditional processes. Finally, through weighted accumulation of volume and grade, a three-dimensional entity model and resource reserve assessment report can be automatically output, achieving integrated processing from discrete borehole data to high-precision ore body digital models. This significantly improves modeling efficiency and assessment reliability in the exploration and mining design of complex ore deposits.

[0042] In one embodiment, S1 of the automated ore body morphology modeling method based on three-dimensional topological reconstruction provided by the present invention specifically includes the following steps:

[0043] S11: Spatialize the acquired discrete borehole data, bind the spatial coordinates of each sampling point with the grade value to generate a three-dimensional spatial point set, and divide the sampling points into points inside the ore body and points in the surrounding rock based on the preset boundary grade threshold to generate preliminary classification results.

[0044] Specifically, the system acquires discrete borehole data, which includes the spatial coordinates of each borehole, the grade data of the borehole sampling points, the borehole number, and the sampling depth. The system performs spatialization preprocessing on the discrete borehole data, encompassing coordinate calibration, outlier removal, and data completion. Coordinate calibration transforms all borehole data to the same spatial coordinate system, eliminating coordinate deviations caused by different exploration stages and acquisition equipment. Outlier removal involves statistically analyzing the grade data of the sampling points to filter and remove outliers that do not conform to the overall distribution pattern. Data completion uses linear interpolation between adjacent sampling points to supplement missing sampling point data, ensuring the continuity and integrity of the borehole data. After spatialization preprocessing, the system binds the spatial coordinates of each sampling point to its corresponding grade value, generating a three-dimensional spatial point set containing the spatial location and grade information of all sampling points. The system calls a preset boundary grade threshold, which is determined by combining industry standards with the actual geological conditions of the mining area. Based on this boundary grade threshold, the system divides all sampling points in the three-dimensional spatial point set. Sampling points with grade values ​​higher than the boundary grade threshold are classified as points inside the ore body, and sampling points with grade values ​​lower than the boundary grade threshold are classified as surrounding rock points. This completes the preliminary classification of the sampling points and generates preliminary classification results, which are used for the accurate identification of boundary points.

[0045] S12: Perform boundary crossing detection processing on the preliminary classification results. Compare the relative relationship between the grade values ​​of adjacent sampling points and the boundary grade threshold segment by segment along the depth direction of each borehole. Identify the boundary points where the grade crosses the boundary grade threshold from high to low and from low to high, and generate a set of boundary positions to be interpolated.

[0046] Specifically, the system calls the generated preliminary classification results and performs boundary crossing detection processing on them. The core purpose is to identify the boundary positions that cross the boundary grade threshold in the borehole depth direction. The system processes all boreholes one by one according to their borehole numbers. For each borehole, the system extracts all sampling points corresponding to that borehole and sorts the sampling points according to the order of sampling depth. The system compares the sorted sampling points segment by segment along the depth direction of each borehole, comparing the relative relationship between the grade values ​​of two adjacent sampling points and the boundary grade threshold. The system identifies two types of boundary points through comparison: one type is where the grade value of the preceding sampling point is higher than the boundary grade threshold and the grade value of the following sampling point is lower than the boundary grade threshold, corresponding to the boundary point where the grade crosses the boundary grade threshold from high to low; the other type is where the grade value of the preceding sampling point is lower than the boundary grade threshold and the grade value of the following sampling point is higher than the boundary grade threshold, corresponding to the boundary point where the grade crosses the boundary grade threshold from low to high. The system records the location and grade relationship of adjacent sampling points corresponding to each type of boundary point, summarizes the spatial locations corresponding to all identified boundary points, and generates a set of boundary locations to be interpolated. This set contains all boundary location information that needs to be further calculated for precise coordinates.

[0047] S13: Perform linear interpolation on each boundary position in the set of boundary positions to be interpolated. Take the sampling points on both sides of the boundary position as high-grade points and low-grade points, respectively. Calculate the straight line equation between the high-grade points and low-grade points and substitute the boundary grade threshold to obtain the scaling factor. Use the scaling factor to perform a weighted average of the spatial coordinates of the high-grade points and the spatial coordinates of the low-grade points to obtain the accurate three-dimensional coordinates of the boundary crossing points, generate the boundary point set, and classify all internal points of the ore body into the internal point set and all surrounding rock points into the surrounding rock point set.

[0048] Specifically, the system calls the set of boundary locations to be interpolated and performs linear interpolation on each boundary location in the set. For each boundary location, the system extracts sampling points on both sides, designating sampling points whose grade values ​​meet the ore body identification criteria as high-grade points and those whose grade values ​​meet the surrounding rock identification criteria as low-grade points. The system obtains the spatial coordinates and corresponding grade values ​​of the high-grade and low-grade points, and establishes a linear equation based on their spatial coordinates. This linear equation characterizes the spatial relationship between the high-grade and low-grade points. The system substitutes the boundary grade threshold into the linear equation and calculates the scaling factor, which is used to determine the spatial position of the boundary crossing point between the high-grade and low-grade points. The system performs a weighted average calculation of the spatial coordinates of the high-grade and low-grade points using the scaling factor, strictly following the linear interpolation principle to ensure the accuracy of the calculation results. The system uses the coordinates obtained from the weighted average calculation as the precise three-dimensional coordinates of the boundary crossing point, and organizes all boundary crossing points to generate a boundary point set. The system also organizes the preliminary classification results generated by S11, assigning all internal points of the ore body to the internal point set and all surrounding rock points to the surrounding rock point set, ultimately forming the internal point set, boundary point set, and surrounding rock point set.

[0049] In one embodiment, S2 of the automated ore body morphology modeling method based on three-dimensional topological reconstruction provided by the present invention specifically includes the following steps:

[0050] S21: Configure constraint conditions for the internal point set, boundary point set, and surrounding rock point set. Based on the geometric constraint principle that the potential function value inside the ore body is positive, the potential function value at the boundary of the ore body is zero, and the potential function value outside the ore body is negative, configure a positive constraint value for each point in the internal point set, a zero constraint value for each point in the boundary point set, and a negative constraint value for each point in the surrounding rock point set, thereby generating a training data point set with constraint values.

[0051] Specifically, in this embodiment, based on the implicit expression principle of ore body geometry, a numerical distribution rule for the potential function in three-dimensional space is set. This rule requires the potential function in the internal region of the ore body to be positive, the potential function on the boundary surface of the ore body to be zero, and the potential function in the external region of the ore body to be negative. Based on this rule, the system assigns a scalar constraint value to each sampling point in the three point sets. For each point in the internal point set, the system configures the same positive constraint number. The absolute value of this positive constraint determines the amplitude level of the potential field inside the ore body, but does not affect the position of the zero isosurface. For each point in the boundary point set, the system configures a zero constraint. These points are forced to fall on the zero isosurface of the final potential function, thereby precisely controlling the boundary position of the ore body obtained through borehole interpolation. For each point in the surrounding rock point set, the system configures a negative constraint number. The absolute value of this negative constraint is equal to the absolute value of the positive constraint configured for the internal points, to ensure that the potential field is symmetrically distributed about the zero isosurface and to avoid overall potential field shift.

[0052] After configuring all constraint values, the system combines the 3D coordinates of each original sampling point with its corresponding constraint value into a training data unit. All training data units together constitute a training data point set with constraint information. This point set serves as input data for subsequent potential field construction steps. Each training point contains its spatial coordinates and the target value that the potential function at that location should satisfy. The system does not introduce any subjective adjustments during constraint value configuration; all internal points share the same positive constraint value, all surrounding rock points share the same negative constraint value, and the constraint values ​​at boundary points are uniformly set to zero, thus ensuring the objectivity and repeatability of the constraint configuration process.

[0053] S22: Perform implicit potential field construction processing on the training data point set with constraint values ​​using radial basis functions. Use compactly supported radial basis functions as interpolation kernel functions, assign a radial basis function to each training point, and weight and superimpose all radial basis functions to form the main part of the potential field. Combine it with the linear drift function to form a composite implicit potential field model, and generate a potential function expression containing unknown coefficients.

[0054] Specifically, the system selects the tightly supported Wendland kernel function as the basic form of the radial basis functions. This kernel function takes a non-zero value within the support radius and a zero value outside the support radius. The system assigns a radial basis function to each training point in the training data set, with the center of the radial basis function located at the spatial coordinates of the training point. The system introduces an adaptive weight coefficient multiplied by each radial basis function. This weight coefficient is calculated based on the spatial distance between the training point and its nearest neighbor borehole sampling point. Specifically, the system calculates the Euclidean distance between each training point and its nearest neighbor sampling point on its corresponding borehole trajectory. This distance reflects the training point's responsiveness to local sampling density. The system maps this distance to the adaptive weight coefficient using a sigmoid function, the form of which controls the sensitivity of the weight coefficient to distance variations. The system linearly superimposes all weighted radial basis functions to form the main body of the potential field.

[0055] Furthermore, the system adds a linear drift term to the linear superposition result. This drift term consists of a set of basis functions and coefficients to be solved. The basis functions include a constant term and three spatial coordinate component terms. The linear drift term is used to improve the behavior of the potential field in regions far from the training point, preventing abnormal oscillations outside the boundary. The system adds the sum of the weighted radial basis functions to the linear drift term to form the complete expression of the composite implicit potential field model. In this expression, the radial basis function coefficients, the shape and scale parameters in the adaptive weight coefficients, and the drift term coefficients are all unknowns to be solved. The system outputs a potential function expression containing unknown coefficients, which mathematically defines the mapping relationship between the potential function value and the three-dimensional spatial coordinates. Preferably, the expression of the composite implicit potential field model is:

[0056]

[0057]

[0058] in, This is the functional expression for the composite implicit potential model. Represents a three-dimensional spatial coordinate vector. To tightly support the Wendland kernel function, , Let J be the adaptive weight coefficient for the j-th training point. Let be the coefficients to be solved for the j-th radial basis function. For linear drift basis vectors, The vector of drift term coefficients. , Let be the Euclidean distance between the j-th training point and its nearest neighbor borehole sampling point. The scaling parameter is set based on the variance of the distribution of all training point locations. The shape parameter is used to adjust the tilt of the sigmoid curve.

[0059] S23: Solve the potential function expression using a system of linear equations. Establish a system of linear algebraic equations based on the condition that the potential function value at each training point is equal to the corresponding constraint value. Rearrange the coefficient matrix and right-hand vector of the system of equations into a symmetric saddle point system. Solve the saddle point system using the conjugate gradient method to obtain all radial basis function coefficients and drift term coefficients. Substitute the solved radial basis function coefficients and drift term coefficients back into the potential function expression to generate a potential function with a complete analytical expression.

[0060] Specifically, the system substitutes the potential function expression containing unknown coefficients into each training point in the training data point set generated in step S21. For the i-th training point, its spatial coordinates are... The system writes out the potential function expression. The expanded form is then compared with the constraint values ​​configured for the training points. The conditions are equal. This condition requires that the potential function strictly passes through the given constraint value at the training point location. The system applies this condition to all N training points one by one, obtaining... There are several linear equations. The unknowns in these equations are the coefficients of the radial basis functions. and drift term coefficient vector The system also introduces a set of additional constraints to ensure the stability of the solution. These additional constraints require that the weighted sum of all radial basis function coefficients be zero, and that the weighted sum of each radial basis function coefficient and its center point coordinates be zero.

[0061] The system will this The equations and additional constraints are combined to form a symmetric saddle point system. The coefficient matrix of the saddle point system is in block matrix form. The upper left block is the radial basis function kernel matrix, and each element of this matrix is ​​calculated from the adaptive weight coefficients and the kernel function. The upper right and lower left blocks are geometric moment matrices, composed of the coordinates of the training points. The lower right block is a zero matrix. The right-hand vector of the system is composed of the constraint values ​​of each training point, and zeros are filled in the positions corresponding to the additional constraints. The system calls a numerical linear algebra solver to solve the saddle point system. Due to the symmetry and sparsity of the coefficient matrix, the system uses the preconditional conjugate gradient method for iterative solution. This method does not require explicit storage of the complete matrix, but only matrix-vector multiplication operations. During the solution process, the system monitors the magnitude of the residual vector, and terminates the iteration when the residual magnitude drops below the preset relative tolerance. The system extracts the coefficients of each radial basis function from the solution results. and the drift term coefficient vector The system substitutes these solved coefficients into the potential function expression, replacing the original unknown symbols, to obtain a complete analytical expression. This expression can return a unique potential function value for any input three-dimensional spatial coordinates and strictly satisfies the constraints at all training points.

[0062] In one embodiment, S3 of the automated ore body morphology modeling method based on three-dimensional topological reconstruction provided by the present invention specifically includes the following steps:

[0063] S31: Perform spatial meshing processing on the potential function. Determine the bounding box of the modeling space based on the spatial distribution range of the borehole data. Divide the interior of the bounding box into uniform cubic cell grids according to the set grid resolution. Calculate the potential function values ​​at the eight vertices of each cubic cell to generate a voxel grid field with potential function values.

[0064] Specifically, the system obtains the complete analytical potential function. The system collects the three-dimensional spatial coordinates of all sampling points within the internal point set, boundary point set, and surrounding rock point set from step S1, and calculates the minimum and maximum values ​​of these coordinates in the X, Y, and Z axes. Using these minimum and maximum values ​​as boundaries, the system constructs an axis-aligned cuboid bounding box that completely covers the spatial distribution range of all sampling points and extends outwards in each direction by a distance correlated with the average borehole spacing to ensure that the ore body boundary does not exceed the bounding box's range.

[0065] The system sets a mesh resolution parameter, which determines the number of voxel mesh divisions in the three directions. Based on the bounding box's length along the X, Y, and Z axes and the set mesh resolution, the system calculates the voxel edge length in each direction. According to the calculated voxel edge lengths, the system divides the bounding box's interior into regular cubic cell meshes, each with the same dimensions. The system iterates through each cubic cell, extracting the 3D spatial coordinates of its eight vertices for the current cell. The system then calls the potential function. The coordinates of each vertex are calculated to obtain the potential function values ​​at the eight vertices. Preferably, the formula for calculating the potential function at the vertices is:

[0066]

[0067] in, Indicates that it is located in the grid index The coordinates of the vertex of the cube cell at that location. To obtain the coefficients of the m-th radial basis function, Let m be the adaptive weight coefficient for the m-th training point. To tightly support the Wendland kernel function, Let m be the spatial coordinates of the m-th training point. For linear drift basis vectors, This is the vector of drift term coefficients. The system associates and stores these potential function values ​​with the corresponding vertex coordinates and spatial position indices of the cube elements, forming a voxel mesh field data structure. This voxel mesh field records the potential function distribution at discrete points throughout the modeling space.

[0068] S32: Perform isosurface intersection detection on the voxel grid field. Traverse each cube cell and check the sign distribution of the potential function values ​​of the eight vertices of the cube cell. If all eight potential function values ​​have the same sign, the cube cell does not intersect with the ore body boundary and is skipped directly. If there are opposite signs, it is determined that the cube cell intersects with the isosurface and a set of cube cells to be generated for isosurfaces is generated.

[0069] Specifically, the system reads the voxel grid field, which contains the potential function values ​​at the eight vertices of each cubic unit. The system sets the target value of the isosurface to zero, and the ore body boundary corresponds to the surface with a zero potential function. The system traverses each cubic unit in spatial index order, extracting the pre-stored potential function values ​​at its eight vertices for the current cubic unit. The system checks the sign of the potential function value at each vertex one by one, using the following rule: if the potential function value is greater than or equal to zero, it is marked as positive; if the potential function value is less than zero, it is marked as negative. The system records the distribution of the eight signs. The system counts the number of positive and negative signs. If all eight signs are positive or all eight signs are negative, it indicates that the cubic unit is completely inside or completely outside the ore body, and the isosurface will not pass through the interior of the cubic unit. The system skips this cubic unit and does not include it in the subsequent processing flow. If both positive and negative signs exist among the eight signs, i.e., opposite signs occur, the system determines that the cubic unit intersects with the isosurface, and the isosurface must pass through the interior region of the cubic unit. The system records the spatial index of the cube element, along with the coordinates of its eight vertices and the potential function values, and adds it to a set of cube elements for which isosurfaces are to be generated. After the system completes traversal and detection of all cube elements, the set contains all cube elements that intersect with the zero isosurface.

[0070] Preferably, the determination condition for this step can be expressed as:

[0071]

[0072] in, and These are two distinct vertices within the same cubic unit. Let be a symbolic function, defined as: if If ≥0, the sign is positive; if A value less than 0 indicates a negative sign. A cube element is considered to intersect an isosurface when at least one pair of vertices have different signs. This condition is equivalent to the case where the product of the potential function values ​​of the eight vertices is less than or equal to zero and not all of them are zero. The system generates a set of cube elements for which isosurfaces are to be generated by performing the above judgment on all cube elements. ,in This represents the eight vertices of the cube element c.

[0073] S33: Perform isosurface tracing processing on each cube element in the cube element set. Based on the positive and negative distribution pattern of the potential function values ​​at the eight vertices of the cube element, obtain the corresponding triangular facet connection topology from the predefined connection method lookup table. Construct one or more triangular faces inside the cube element according to the found connection method. Determine the vertex positions of the triangular faces through linear interpolation of the potential function values ​​and output them as a continuous mesh structure to generate the initial triangular mesh.

[0074] Specifically, the system calls a pre-stored lookup table of connection methods. This table contains entries corresponding to the positive and negative distribution patterns of the potential function values ​​at the vertices of the cube element. Each entry defines the number of triangular faces to be generated under the corresponding pattern and the edge information of each triangular face's vertex. The system iterates through each cube element in the cube element set. For the current cube element, the system obtains the potential function values ​​of all its vertices. The system generates an eight-bit binary index value based on the positive and negative distribution of the vertex potential function values. Each bit of the index value corresponds to a vertex; a positive potential function value corresponds to a 1, and a negative or zero potential function value corresponds to a 0.

[0075] The system uses the index value of the current cube element to look up a join table and obtain a list of triangular faces that should be generated within the current cube element. This list contains the edge identifiers corresponding to each triangular face. For each triangular face in the list, the system sequentially determines the spatial coordinates of its three vertices. The system then locates the edge of the cube element to which the current vertex belongs, obtaining the coordinates and potential function values ​​of the two endpoints of the edge. The potential function values ​​of the two endpoints must satisfy the condition that they have opposite signs or one of them is zero. The system performs linear interpolation calculations targeting isosurfaces with zero potential function values. The linear interpolation formula is:

[0076]

[0077] in, and Let the coordinates be the coordinates of the two endpoints of the edge. and For the potential function values ​​at the corresponding endpoints, the interpolation parameters are... Depend on and The calculated value of t is within the interval. The system uses this formula to calculate the coordinates of the intersection points of the isosurface and the edges. The system repeats the above interpolation process for the three intersection points of each triangular facet to obtain the coordinates of the three intersection points to form a complete triangular facet.

[0078] The system adds all triangular faces generated within the current cube element to the global triangular face set. During this process, the system removes duplicate vertices, merging vertices with coordinate distances less than a set tolerance threshold into a single vertex, and updates the vertex indices of the triangular faces. After processing all cube elements, the system organizes the global triangular face set into a vertex list and a face index list, forming a continuous mesh data structure. The system outputs this mesh as the initial triangular mesh, which reflects the discrete triangular approximation of the zero isosurface of the potential function.

[0079] In one embodiment, S4 of the automated ore body morphology modeling method based on three-dimensional topological reconstruction provided by the present invention specifically includes the following steps:

[0080] S41: Perform connectivity component analysis and noise filtering on the initial triangular mesh. Construct a graph connectivity model based on the edge-sharing relationship of triangular facets. Divide the mesh into several independent connected subgraphs using breadth-first search. Use the signed tetrahedral volume method to estimate the three-dimensional volume value of each connected subgraph. Delete the connected subgraphs with volume values ​​lower than the preset noise fragment volume threshold from the mesh to generate a denoised triangular mesh.

[0081] Specifically, the system obtains an initial triangular mesh, which contains a set of vertices and a set of triangular faces. Using triangular faces as graph nodes and connecting triangles by sharing an edge, the system constructs an undirected connected graph model. The system performs a breadth-first search traversal of this graph, starting from any unvisited triangle and visiting all adjacent faces sharing the same edge, until no further expansion is possible. All traversed triangles form a connected subgraph. The system repeats this process until all triangles are assigned to a connected subgraph, thus dividing the initial triangular mesh into one or more independent connected subgraphs, each corresponding to a spatially separated mesh patch.

[0082] The system estimates the volume of each connected subgraph. For a given triangular mesh surface represented by a connected subgraph, the system first determines whether the surface is a closed manifold. If the surface is not closed, the system uses the signed tetrahedral volume method to calculate the approximate volume. It selects the 3D centroid of any triangular facet in the subgraph as a reference point, and constructs a tetrahedron with each triangular facet in the subgraph. The system calculates the signed volume of each tetrahedron and sums them. The absolute value of the sum is used as the estimated volume of the subgraph. The formula for calculating the signed tetrahedral volume is:

[0083]

[0084] Where m is the number of triangular faces contained in the connected subgraph. For the reference point coordinates, the centroid of the first triangular facet in the subgraph is usually taken. For the k-th triangular facet, the coordinates of its three vertices are as follows: , , , This represents the determinant of a matrix with three vectors as column vectors. The absolute value of the summation divided by six gives the approximate volume of the space enclosed by the subgraph. This formula assumes that the reference point is located inside the entity enclosed by the subgraph. If this condition is not met, the volume estimation will have a sign error, but the absolute value can still be used as a basis for comparing relative volumes.

[0085] The system presets a noise debris volume threshold, which is determined proportionally based on the overall scale of the modeling space and the average borehole spacing. This threshold is used to distinguish the main body of geologically significant ore bodies from small noise debris caused by data sparsity or isosurface extraction errors. The system compares the volume estimate of each connected subgraph with the noise debris volume threshold. If the volume estimate is lower than the threshold, the system determines the subgraph as a noise debris and removes all its triangular faces from the global mesh data structure. If the volume estimate is greater than or equal to the threshold, the system retains the subgraph. After filtering all subgraphs, the system outputs a denoised triangular mesh, which retains only the major connected components whose volume exceeds the threshold.

[0086] S42: Based on skeleton extraction technology, local densification of topologically complex regions of the denoised triangular mesh is performed. The topological skeleton of the mesh is represented by a Reeb graph using the gradient field of the potential function and the connectivity of the level set. All branch nodes and the topological structure of the branches in the topological skeleton are extracted. The branch compound and pinch-out reproduction regions in the ore body morphology are identified. The triangular patches in the branch compound and pinch-out reproduction regions are densified and reconstructed to generate a topologically enhanced triangular mesh.

[0087] Specifically, the system acquires the denoised triangular mesh and the complete potential function. The system uses Reeb graph technology to extract the topological skeleton of the mesh. The Reeb graph is defined based on the change in the connectivity of the level set of a monotonic function on the mesh. The system selects the potential function as the defining function of the Reeb graph. The potential function value changes monotonically within the ore body, and the change in the connectivity of its isosurfaces intersecting the mesh reflects topological events. The system projects each vertex of the denoised triangular mesh onto the potential function value axis and sorts the vertices according to the potential function value in ascending order. The system sets a series of isosurface levels within the potential function value range. For each level, the system calculates the potential function value equal to the intersection of the level line and the mesh at that level and analyzes the connectivity components of the intersection. The system calculates the gradient field of the potential function. The formula for calculating the gradient field is:

[0088]

[0089] in, Let be the gradient field of the potential function. In three-dimensional space coordinates, These are the partial derivatives of the potential function in the x, y, and z directions, respectively. When the potential function value changes continuously, the appearance, splitting, merging, or disappearance of connected components corresponds to nodes in the Reeb graph, and continuous intervals between adjacent nodes correspond to edges in the Reeb graph. The system constructs a complete Reeb graph representation based on these node-edge connections. The system analyzes the topological structure of the Reeb graph, identifying branch nodes (nodes with degree three or higher) and terminal nodes (nodes with degree one). Branch nodes correspond to branching and compounding positions in the ore body morphology, and terminal nodes correspond to pinch-out and reappearance positions.

[0090] The system extracts the coordinates of these nodes in 3D space and determines the local mesh region affected by each node. This region is centered on the node and defined by the change in potential function value or spatial distance as the radius. The system performs a refinement reconstruction operation on the local region of each identified branch and terminal node. Within this local region, the system increases the sampling density, reduces the original voxel size by a factor, and re-establishes the local voxel mesh. Based on a local subset of the training point set used in steps S22 and S23 and the corresponding potential function, the system recalculates the potential function within the local region and calls the traveling cube algorithm to generate a higher-resolution set of triangular facets. The system stitches the newly generated local triangular facet set with the original denoised triangular mesh at the boundary, ensuring geometric continuity through vertex matching and edge splitting operations. After completing the refinement reconstruction of all topologically complex regions, the system outputs a topology-enhanced triangular mesh. This mesh has higher geometric resolution in morphologically abrupt regions such as ore body branches and pinch-outs, and can more accurately represent topological connectivity relationships.

[0091] S43: Perform potential-oriented hole filling processing on the topology-enhanced triangular mesh, detect all hole boundaries in the mesh, construct a sequence of boundary vertices for each hole boundary, calculate the projection point of each boundary vertex in the gradient direction of the potential function, connect the projection point with the boundary vertex to form a new triangular patch through a constrained triangulation algorithm, iterate filling until all hole boundaries are completely closed, and generate a manifold triangular mesh in which each edge is shared by at most two triangular patches.

[0092] Specifically, the system acquires a topology-enhanced triangular mesh and a complete potential function. The system traverses all edges in the mesh, counting the number of times each edge is referenced by a triangular facet. If an edge is referenced by only one triangular facet, the system determines that edge is a boundary edge. The system connects the boundary edges end-to-end, constructing a closed sequence of boundary loops, each corresponding to a hole. For each detected hole boundary loop, the system extracts the vertex sequence constituting the loop and sorts it according to the loop's direction. The system calculates the average potential function value of the hole boundary loop, which is typically close to zero but may deviate slightly due to mesh discretization errors. The system then performs a process on each vertex of the boundary loop... Calculate the gradient of the potential function at the current vertex. The gradient direction is perpendicular to the isosurface of the potential function, i.e., pointing towards the normal direction of the ore body boundary. The system sets the iterative projection step size factor, which is determined based on the local average side length of the cavity boundary loop. The system performs iterations on each boundary vertex. Move one step along the gradient direction to obtain the projection point. .

[0093] To constrain the projection points to the isosurface where the potential function is zero, the system employs Newton's iteration method for correction. The core calculation formula for the correction is:

[0094]

[0095] in, The vertex position after the k-th iteration. This represents the potential function value at that point. Let be the gradient vector of the potential function at that point. This formula projects the vertex directly onto a surface where the potential function is zero. The direction of movement in each iteration is either opposite to or the same as the gradient direction, depending on the sign of the potential function value. The distance moved is determined by the ratio of the current potential function value to the gradient magnitude. The system repeats the iterative process until the absolute value of the potential function is less than a set tolerance. The system collects all projected points into a new vertex set. The system uses the original hole boundary vertices and projected points as a point set, and employs a constrained triangulation algorithm to generate triangular patches within the hole region. The constraints are that the original boundary edges must be preserved, and the lines connecting projected points and the original boundary vertices cannot intersect.

[0096] The system adds newly generated triangular facets to the mesh. Since the addition of new facets may introduce new boundary edges, the system repeats the above detection and filling process, projecting and triangulating the newly added boundary loops again until no boundary edges remain in the mesh. Finally, the system checks for non-manifold edges in the mesh, i.e., edges shared by three or more triangular facets. If non-manifold edges exist, the system splits the vertices along the edge, dividing the shared edge into two groups, each using different vertices after the split, thus decomposing the non-manifold edge into multiple manifold edges. The system outputs a manifold triangular mesh, in which each edge is shared by at most two triangular facets, with no holes and no dangling edges. The system uses this projection formula to correct the position of each hole boundary vertex, ensuring that the newly generated triangular facets are geometrically continuous with the original mesh and conform to the isosurface shape defined by the potential field.

[0097] In one embodiment, S5 of the automated ore body morphology modeling method based on three-dimensional topological reconstruction provided by the present invention specifically includes the following steps:

[0098] S51: Using a manifold triangular mesh as the spatial boundary constraint, generate a grid of block center points arranged according to a set sampling resolution inside the grid. For each block center point, use the ray method to detect whether the block center point is located inside the manifold triangular mesh. Only retain the block center points located inside the grid to generate the set of effective block locations in the ore body domain.

[0099] Specifically, the system acquires a manifold triangular mesh, which defines the closed boundary surface of the ore body. The system sets a sampling resolution parameter, which is related to the average borehole spacing and the desired level of detail in the block model, to control the distribution density of block center points in three spatial directions. Based on the lengths of the bounding box of the manifold triangular mesh along the X, Y, and Z axes and the set sampling resolution, the system calculates the sampling step size in each direction. Starting from the smallest point of the bounding box, the system generates a series of regularly arranged points in the three directions using the calculated sampling step size, forming a three-dimensional point matrix. Each point corresponds to the center point of a potential block. The system performs an inclusion test on each generated block center point to determine whether the point is located inside the closed body enclosed by the manifold triangular mesh.

[0100] Preferably, the inclusion test can employ a ray casting method, generating a ray from the center point of the block in any direction, such as the positive X-axis, and calculating the number of intersections between this ray and all triangular faces of the manifold triangular mesh. For each triangular face, the system solves for the intersection points between the ray and the plane containing the triangle, and then determines whether the intersection point lies inside the triangle. This step of determining the intersection points between the ray and the triangular face requires solving for the intersection parameter t and the centroid coordinates. The calculation formula is:

[0101]

[0102] in, The coordinates of the ray origin, i.e., the center point of the block. Let be the ray direction vector. , , Let be the coordinates of the three vertices of the triangular facet. For ray parameters, and These are the components of the centroid coordinate system. Solving the above system of linear equations yields... , , ,like and If the ray intersects the triangular facet, then the system performs this calculation on all triangular facets and counts the number of intersections that satisfy the condition.

[0103] The system records all valid intersection points. If the number of intersection points is odd, the block center point is determined to be inside the grid; if the number is even, it is determined to be outside. In its implementation, the system uses random ray directions or multiple rays emitted along the coordinate axes for voting to avoid ambiguity caused by rays passing exactly through vertices or edges. The system retains block center points determined to be inside the grid and records their 3D coordinates, while discarding block center points determined to be outside. After detecting all potential block center points, the system collects all retained internal points into a set of valid block locations. Each element in this set corresponds to the center location of a block within the ore body domain.

[0104] S52: Spatial variation analysis and Kriging interpolation are performed on the grade data within the internal and boundary point sets. Based on all borehole grade data, experimental variation function values ​​under different hysteresis distances are calculated. The weighted least squares method is used to fit a spherical model to the experimental variation function point cloud to obtain the nugget value, sill value, and range parameter. The fitted variation parameters are used to establish and solve the ordinary Kriging equation system for each block center point in the effective block location set. The optimal linear unbiased estimation of the grade at each block center point is performed using the borehole grade data as input to generate the block grade model.

[0105] Specifically, the system acquires the grade values ​​and three-dimensional coordinates of each sampling point in the internal and boundary point sets. The system performs spatial variation analysis to characterize the variation of grade values ​​with spatial distance. The system defines a series of lag intervals. For each lag interval, the system searches for all pairs of sampling points within the allowable range of that lag interval, calculates the squared difference in grade values ​​for each pair, and uses half of the average of all these differences as the experimental variation function value. The calculation formula is:

[0106]

[0107] in, The lag distance, This represents the number of all point pairs that are within the allowable hysteresis distance. and These represent the grade values ​​of the two sampling points in the point pair. The system calculates the directional experimental variability function along the strike, dip, and perpendicular directions of the ore body to reflect the anisotropy of mineralization, resulting in a scatter set consisting of the hysteresis distance and its corresponding experimental variability function value.

[0108] The system selects a spherical model as the theoretical variogram model. This model increases throughout the range, reaches a sill value at the range, and then remains constant. The mathematical expression for the spherical model is:

[0109]

[0110] in, Value of a nugget. The value is the partial sill value, and the sill value is... and The sum of For variable range, the system uses weighted least squares to fit the experimental variation point cloud, assigning higher weights to each point pair with a larger number of small lag intervals to obtain model parameters. The fitted parameters include nugget value, sill value, and variable range.

[0111] The system obtains a set of valid block locations, where each element corresponds to the center point of a block for which the grade needs to be estimated. The system performs ordinary kriging interpolation on the center point of each block. It searches the sample point set for the few closest sampling points to \mathbf{x}_0. If the number of sampling points is less than the minimum requirement, the search radius is expanded or an ellipsoidal search region is used. Based on the fitted variogram model, the system calculates the variogram values ​​between sampling points and between sampling points and the point to be estimated, constructing the ordinary kriging equations:

[0112]

[0113] in, Let be the weight coefficient for the j-th sampling point. For Lagrange multipliers, To fit the obtained variation function model, and These are the coordinates of the sampling point. The coordinates of the center point of the block to be estimated are: Let be the number of sampling points in the neighborhood. The system solves this system of equations to obtain the weight coefficient for each sampling point. The grade value of each sampling point is multiplied by its corresponding weight coefficient and then summed to obtain the grade estimate for the center point of the block. The formula for calculating the grade estimate is as follows: The system simultaneously calculates the estimated variance as a measure of estimation uncertainty. The system associates and stores the coordinates of all block center points with the corresponding grade estimates to form a block grade model.

[0114] S53: Perform resource-weighted accumulation processing on the block grade model. Traverse each effective block in the block grade model, and perform block-by-block weighted accumulation of the estimated grade value of each effective block with the corresponding block volume and ore density. At the same time, independently accumulate the volume of each block. Summarize the accumulation results according to ore quantity and metal quantity respectively to generate a resource reserve assessment report of the ore body. The manifold triangular mesh is used as the geometric skeleton of the three-dimensional solid model and directly output as an STL or OBJ format file to output the three-dimensional solid model of the ore body.

[0115] Specifically, the system acquires a set of valid block locations and a block grade model. Each valid block corresponds to a block center point. The system calculates the volume of each block based on the sampling resolution parameter; the volume value is the product of the sampling step sizes in three directions, and all blocks have the same volume value. The system sets the ore density parameter, which is determined based on measured data from the mining area or empirical formulas. The system iterates through each valid block in the block grade model. For the i-th block, the system reads its estimated grade value. Block volume and ore density parameters The core formula for calculating resource quantity is:

[0116]

[0117] in, The total number of valid blocks. This represents the total ore quantity. The system calculates the ore quantity of this block, which is equal to the product of volume and density. The system adds the current ore quantity of this block to the total ore quantity accumulator. The system also calculates the metal content of this block, which is equal to the ore quantity multiplied by the estimated grade. The corresponding calculation formula is:

[0118]

[0119] in, To determine the total metal content, the system adds the metal content of the current block to the total metal content accumulator. Simultaneously, the system independently accumulates the total volume of all blocks for subsequent average grade calculation. After traversing all blocks, the system obtains the total ore quantity, total metal content, and total volume. The system then calculates the average grade of the ore body using the following formula:

[0120]

[0121] in, The average grade of the ore body is calculated by averaging the estimated grades of all blocks when their volumes are equal. The system uses a weighted average formula to ensure generality. The system generates a resource reserve assessment report, which includes the total ore quantity, total metal quantity, average grade, and statistical information categorized by grade range. Categorical statistics divide the blocks into different grade levels by setting multiple grade thresholds, summarizing the ore quantity and metal quantity within each level. The system outputs manifold triangular meshes in standard 3D model file formats. The STL format stores the mesh as a series of triangular facets with vertex coordinates and normal vectors, without color or texture information. The OBJ format supports detailed descriptions of vertices, facets, and normal vectors, offering wider compatibility. Both formats are output for use by subsequent geological modeling or mine design software. The system uses the manifold triangular mesh as the geometric skeleton of the ore body's 3D solid model, accurately representing the spatial morphological boundaries of the ore body.

[0122] It should be understood that although the steps in the flowcharts of the embodiments described above are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages of other steps.

[0123] Based on the same inventive concept, this application also provides an automated ore body morphology modeling device based on three-dimensional topology reconstruction for implementing the aforementioned automated ore body morphology modeling method based on three-dimensional topology reconstruction. The solution provided by this device is similar to the solution described in the above method. Therefore, the specific limitations of one or more embodiments of the automated ore body morphology modeling device based on three-dimensional topology reconstruction provided below can be found in the limitations of the automated ore body morphology modeling method based on three-dimensional topology reconstruction described above, and will not be repeated here.

[0124] Preferably, such as Figure 2 As shown, the present invention provides an automated ore body morphology modeling system 600 based on three-dimensional topology reconstruction, which is configured with the following modules:

[0125] The borehole data classification module 610 is used to perform spatial processing on the acquired discrete borehole data, and classify the sampling points into internal points of the ore body, boundary crossing points and surrounding rock points according to the boundary grade threshold, generating internal point sets, boundary point sets and surrounding rock point sets;

[0126] The implicit potential field construction module 620 is used to construct a radial basis function implicit potential field by applying positive constraints to points in the internal point set, zero constraints to points in the boundary point set, and negative constraints to points in the surrounding rock point set, and to solve for the radial basis function coefficients and drift term coefficients of the radial basis function implicit potential field to generate the potential function.

[0127] The orebody mesh initial construction module 630 is used to call the moving cube algorithm to perform voxel-level tracking of the potential function at the isosurface level. For each cube cell that intersects with the isosurface in the modeling space, the connection method of the triangular facets inside the cube cell is determined by a lookup table based on the positive and negative combinations of the potential function values ​​at the eight vertices of the cube cell, thereby generating an initial triangular mesh that reflects the orebody boundary.

[0128] The mesh topology optimization module 640 is used to perform three-dimensional topology reconstruction on the initial triangular mesh. It sequentially removes noise fragments through connected component analysis, performs local densification and reconstruction of branches and pinch-out regions through skeleton extraction, fills the hole boundaries through iterative projection guided by potential function gradient, and after multiple local potential field reconstructions and iterative reconnection of mesh patches, the non-manifold mesh is transformed into a manifold triangular mesh in which each edge is shared by at most two triangular patches, thus generating a manifold triangular mesh.

[0129] The orebody reserve modeling module 650 is used to build block models within a spatially constrained domain using a manifold triangular mesh. Based on the grade data within the internal point set and boundary point set, spatial variation analysis is performed to obtain variation function parameters. Spatial interpolation methods are called to estimate the grade of the center point of each block, generating a block grade model. Based on the block grade model, all valid blocks are traversed to perform a weighted summation of volume and grade, generating a three-dimensional solid model of the orebody and a resource reserve assessment report.

[0130] In one embodiment, this application also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the above-described automated modeling method for ore body morphology based on three-dimensional topological reconstruction.

[0131] In one embodiment, this application also provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the above-described automated ore body morphology modeling method based on three-dimensional topological reconstruction.

[0132] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of this application. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of those different embodiments or examples.

[0133] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to in the description of the method embodiments. The device embodiments described above are merely illustrative. The components described as separate parts may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this disclosure according to actual needs. Those skilled in the art can understand and implement this without creative effort.

[0134] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various variations or substitutions within the technical scope disclosed in this application, and these should all be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. An automated modeling method for ore body morphology based on three-dimensional topological reconstruction, characterized in that, Includes the following steps: S1: Spatialize the acquired discrete borehole data and classify the sampling points into internal points of the ore body, boundary crossing points, and surrounding rock points according to the boundary grade threshold, generating internal point sets, boundary point sets, and surrounding rock point sets; S2: Apply positive constraints to the points in the internal point set, zero constraints to the points in the boundary point set, and negative constraints to the points in the surrounding rock point set to construct a radial basis function implicit potential field, and solve for the radial basis function coefficients and drift term coefficients of the radial basis function implicit potential field to generate the potential function. S3: Call the moving cube algorithm to perform voxel-level tracing of the potential function at the isosurface level. For each cube cell that intersects with the isosurface in the modeling space, determine the connection method of the triangular facets inside the cube cell by looking up the positive and negative combinations of the potential function values ​​at the eight vertices of the cube cell, and generate an initial triangular mesh that reflects the boundary of the ore body. S4: Perform three-dimensional topological reconstruction on the initial triangular mesh, successively remove noise fragments through connected component analysis, perform local densification and reconstruction of branches and pinch-out regions through skeleton extraction, fill the hole boundaries through iterative projection guided by potential function gradient, and after multiple local potential field reconstructions and iterative reconnection of mesh patches, transform the non-manifold mesh into a manifold triangular mesh in which each edge is shared by at most two triangular patches, thus generating a manifold triangular mesh; S5: Using the manifold triangular mesh as the spatial constraint domain, a block model is established within the spatial constraint domain. Spatial variation analysis is performed based on the grade data within the internal point set and the boundary point set to obtain the variation function parameters. The spatial interpolation method is called to estimate the grade of the center point of each block, generating a block grade model. Based on the block grade model, all valid blocks are traversed to perform a weighted summation of volume and grade, generating a three-dimensional solid model of the ore body and a resource reserve assessment report.

2. The method according to claim 1, characterized in that, S1 includes: S11: Spatialize the acquired discrete borehole data, bind the spatial coordinates of each sampling point with the grade value to generate a three-dimensional spatial point set, and divide the sampling points into points inside the ore body and points in the surrounding rock based on the preset boundary grade threshold to generate preliminary classification results. S12: Perform boundary crossing detection processing on the preliminary classification results, compare the relative relationship between the grade values ​​of adjacent sampling points and the boundary grade threshold segment by segment along the depth direction of each borehole, identify the boundary points where the grade crosses the boundary grade threshold from high to low and the boundary points where the grade crosses the boundary grade threshold from low to high, and generate a set of boundary positions to be interpolated. S13: Perform linear interpolation on each boundary position in the set of boundary positions to be interpolated. Take the sampling points on both sides of the boundary position as high-grade points and low-grade points, respectively. Calculate the straight line equation between the high-grade points and the low-grade points and substitute the boundary grade threshold to obtain the scaling factor. Use the scaling factor to perform a weighted average of the spatial coordinates of the high-grade points and the spatial coordinates of the low-grade points to obtain the precise three-dimensional coordinates of the boundary crossing points, generate the boundary point set, and classify all internal points of the ore body into the internal point set and all surrounding rock points into the surrounding rock point set.

3. The method according to claim 1, characterized in that, S2 includes: S21: Configure constraint conditions for the internal point set, the boundary point set, and the surrounding rock point set. Based on the geometric constraint principle that the potential function value inside the ore body is positive, the potential function value at the ore body boundary is zero, and the potential function value outside the ore body is negative, configure a positive constraint value for each point in the internal point set, a zero constraint value for each point in the boundary point set, and a negative constraint value for each point in the surrounding rock point set, thereby generating a training data point set with constraint values. S22: The training data point set with attached constraint values ​​is processed by radial basis function implicit potential field construction. The compact support radial basis function is used as the interpolation kernel function. A radial basis function is assigned to each training point. All radial basis functions are weighted and superimposed to form the main part of the potential field. It is combined with the linear drift function to form a composite implicit potential field model and generate a potential function expression containing unknown coefficients. S23: Solve the potential function expression using a system of linear equations. Establish a system of linear algebraic equations based on the condition that the potential function value at each training point is equal to the corresponding constraint value. Rearrange the coefficient matrix and right-hand vector of the system of equations into a symmetric saddle point system. Solve the saddle point system using the conjugate gradient method to obtain all radial basis function coefficients and drift term coefficients. Substitute the solved radial basis function coefficients and drift term coefficients back into the potential function expression to generate a potential function with a complete analytical expression.

4. The method according to claim 3, characterized in that, The expression for the composite implicit potential field model is as follows: in, This is the functional expression for the composite implicit potential model. Represents a three-dimensional spatial coordinate vector. To tightly support the Wendland kernel function, , Let J be the adaptive weight coefficient for the j-th training point. Let be the coefficients to be solved for the j-th radial basis function. For linear drift basis vectors, The vector of drift term coefficients. Let be the Euclidean distance between the j-th training point and its nearest neighbor borehole sampling point. The scaling parameter is set based on the variance of the distribution of all training point locations. The shape parameter is used to adjust the tilt of the sigmoid curve.

5. The method according to claim 1, characterized in that, S3 includes: S31: The potential function is modeled by spatial meshing. The bounding box of the modeling space is determined according to the spatial distribution range of the borehole data. The inside of the bounding box is divided into uniform cubic unit meshes according to the set mesh resolution. The potential function values ​​at the eight vertices of each cubic unit are calculated to generate a voxel mesh field with potential function values. S32: Perform isosurface intersection detection processing on the voxel grid field, traverse each cube cell, check the sign distribution of the potential function values ​​of the eight vertices of the cube cell, if all eight vertices have the same sign, the cube cell does not intersect with the ore body boundary and is skipped directly, if there are opposite signs, it is determined that the cube cell intersects with the isosurface, and a set of cube cells to be generated for isosurfaces is generated. S33: Perform isosurface tracing processing on each cube unit in the cube unit set. Based on the positive and negative distribution pattern of the potential function values ​​at the eight vertices of the cube unit, obtain the corresponding triangular facet connection topology from the predefined connection method lookup table. Construct one or more triangular facests inside the cube unit according to the found connection method. Determine the vertex positions of the triangular facests by linear interpolation of the potential function values ​​and output them uniformly as a continuous mesh structure to generate the initial triangular mesh.

6. The method according to claim 1, characterized in that, S4 includes: S41: Perform connectivity component analysis and noise filtering on the initial triangular mesh, construct a graph connectivity model based on the edge-sharing relationship of triangular facets, divide the mesh into several independent connected subgraphs through breadth-first search, estimate the three-dimensional volume value of each connected subgraph by calling the signed tetrahedral volume method, and delete the connected subgraphs with volume values ​​lower than the preset noise fragment volume threshold from the mesh as a whole to generate a denoised triangular mesh. S42: Based on the skeleton extraction technology, the topologically complex regions of the denoised triangular mesh are locally densified. The topological skeleton of the mesh is represented by a Reeb graph using the gradient field of the potential function and the connectivity of the level set. All branch nodes and the topological structure of the branches in the topological skeleton are extracted. The branch compound and pinch-out reproduction regions in the ore body morphology are identified. The triangular patches in the branch compound and pinch-out reproduction regions are densified and reconstructed to generate a topologically enhanced triangular mesh. S43: Perform potential field-guided hole filling processing on the topology-enhanced triangular mesh, detect all hole boundaries in the mesh, construct a boundary vertex sequence for each hole boundary, calculate the projection point of each boundary vertex in the gradient direction of the potential function, connect the projection point with the boundary vertex to form a new triangular patch through a constrained triangulation algorithm, iterate filling until all hole boundaries are completely closed, and generate a manifold triangular mesh in which each edge is shared by at most two triangular patches.

7. The method according to any one of claims 1-6, characterized in that, S5 includes: S51: Using the manifold triangular mesh as a spatial boundary constraint, generate a block center point mesh arranged according to a set sampling resolution inside the mesh. Use the ray method to detect whether the block center point is located inside the manifold triangular mesh for each block center point. Only retain the block center points located inside the mesh to generate a set of effective block positions in the ore body domain. S52: Perform spatial variation analysis and Kriging interpolation on the grade data within the internal point set and the boundary point set. Calculate the experimental variation function values ​​under different lag distances based on all borehole grade data. Use the weighted least squares method to fit a spherical model to the experimental variation function point cloud to obtain the nugget value, sill value, and range parameter. Use the fitted variation parameters to establish and solve the ordinary Kriging equation system for each block center point in the effective block location set. Use the borehole grade data as input to perform optimal linear unbiased estimation of the grade at each block center point, generating a block grade model. S53: Perform resource-weighted accumulation processing on the block grade model. Traverse each effective block in the block grade model, and perform block-by-block weighted accumulation of the estimated grade value of each effective block with the corresponding block volume and ore density. At the same time, independently accumulate the volume of each block. Summarize the accumulation results according to ore quantity and metal quantity to generate a resource reserve assessment report of the ore body. The manifold triangular mesh is used as the geometric skeleton of the three-dimensional solid model and directly output as an STL or OBJ format file to output the three-dimensional solid model of the ore body.

8. An automated ore body morphology modeling system based on three-dimensional topological reconstruction, characterized in that, The system includes: The borehole data classification module is used to spatialize the acquired discrete borehole data. Based on the boundary grade threshold, the sampling points are classified into points inside the ore body, boundary crossing points, and surrounding rock points, generating internal point sets, boundary point sets, and surrounding rock point sets. An implicit potential field construction module is used to construct a radial basis function implicit potential field by applying positive constraints to points in the internal point set, zero constraints to points in the boundary point set, and negative constraints to points in the surrounding rock point set, and to solve for the radial basis function coefficients and drift term coefficients of the radial basis function implicit potential field to generate a potential function. The orebody mesh initialization module is used to call the moving cube algorithm to perform voxel-level tracing of the potential function at the isosurface level. For each cube cell that intersects with the isosurface in the modeling space, the connection mode of the triangular facets inside the cube cell is determined by a lookup table based on the positive and negative combinations of the potential function values ​​at the eight vertices of the cube cell, thereby generating an initial triangular mesh that reflects the orebody boundary. The mesh topology optimization module is used to perform three-dimensional topology reconstruction on the initial triangular mesh. It sequentially removes noise fragments through connected component analysis, performs local densification and reconstruction of branches and pinch-out regions through skeleton extraction, fills the hole boundaries through iterative projection guided by potential function gradient, and after multiple local potential field reconstructions and iterative reconnection of mesh patches, the non-manifold mesh is transformed into a manifold triangular mesh in which each edge is shared by at most two triangular patches, thus generating a manifold triangular mesh. The ore body reserve modeling module is used to establish a block model within the spatial constraint domain using the manifold triangular mesh as the spatial constraint domain. Based on the grade data in the internal point set and the boundary point set, spatial variation analysis is performed to obtain the variation function parameters. The spatial interpolation method is called to estimate the grade of the center point of each block, generating a block grade model. Based on the block grade model, all valid blocks are traversed to perform a weighted summation of volume and grade, generating a three-dimensional solid model of the ore body and a resource reserve assessment report.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 7.

Images