A kriging modeling method of regular optimization and topological collapse

Abstract

The application discloses a kind of regular optimization and topological collapse's Kriging modeling method, method includes: obtaining borehole sample data, generates sample point set and voxel grid;According to the binning of lag distance, generate experimental half variance set and point pair count set, construct regularization objective function based on trust weight and global variance constraint term, and output variation parameter set containing range by L-BFGS-B optimizer;With range as cKDTree truncation boundary extraction neighborhood sample set, according to matrix scale and available physical memory generation batch capacity, solve ordinary Kriging equation set in batches, obtain predicted grade voxel set;Again generate grade mask and render voxel set, encode as Trace tensor and static grid object after writing into WebGL buffer object.The application realizes the stable, efficient and lightweight display of ore body Kriging modeling by regularization optimization, range-driven computing scheduling and topological collapse rendering.

Description

Technical Field

[0001] This application relates to the field of digital mining and 3D geological modeling technology, and in particular to a kriging modeling method with regularization optimization and topological collapse. Background Technology

[0002] Three-dimensional orebody resource estimation typically employs ordinary kriging interpolation, which fits spatial autocorrelation through a variogram and solves for neighborhood sample weights on a voxel grid.

[0003] Existing methods often separate the fitting of variogram functions, neighborhood search, matrix solving and front-end rendering. Although they can complete modeling of conventional scale, they still have underlying defects that are easily overlooked when there are sparse boreholes, dense local sample clusters and millions of voxels displayed on the same screen.

[0004] First, existing automatic fitting methods often only pursue the minimum residual, ignoring the differences in the number of point pairs within bins with different lag distances. A small number of distant outliers in the semivariogram can create spurious gradients, causing the optimizer to output mathematically convergent but geophysically distorted nugget values, sill values, or ranges, and even abnormal results such as sill exceeding global variance. Second, existing massive voxel kriging calculations usually set the search radius and batch size empirically, without mapping the range to cKDTree pruning boundaries and memory scheduling factors. When the target voxel falls into a high-density sample cluster, the size of the covariance matrix will increase quadratically with the number of neighboring samples, leading to RAM overflow, thread crashes, and interpolation interruptions.

[0005] Furthermore, if the generated heterogeneous grade voxels continue to be rendered using SVG discrete objects or GPU instance matrices, it will cause rendering call congestion and UI freezing. Therefore, there is an urgent need for a 3D orebody kriging modeling method that integrates adaptive regularization optimization, variable-range driven neighborhood pruning, dynamic batch processing, and topology collapse rendering. Summary of the Invention

[0006] To address the aforementioned problems, embodiments of the present invention provide a Kriging modeling method for regularization optimization and topology collapse, the method comprising:

[0007] Obtain borehole sample data, perform coordinate cleaning, quality verification and spatial registration on the borehole sample data, generate sample point set, and generate voxel mesh based on the spatial bounding box of the sample point set;

[0008] Based on the sample point set, an experimental semivariance set and a point-pair count set are generated by binning according to the lag distance. The trust weight is determined according to the point-pair count set. The trust weight, global variance constraint term and spherical variability function are written into the regularization objective function. The L-BFGS-B optimizer with boundary mapping outputs a set of variability parameters containing nugget value, sill difference and range.

[0009] Using the range as the cutoff boundary of the cKDTree retrieval topology, the neighborhood sample set of the target voxel is extracted. The matrix size of the ordinary kriging equation system is determined based on the sample size of the neighborhood sample set. The batch processing capacity is generated based on the matrix size and available physical memory. The ordinary kriging equation system is assembled and solved according to the batch processing capacity to obtain the predicted grade voxel set.

[0010] Based on the predicted quality voxel set, a quality mask and a rendering voxel set are generated. The rendering voxel set is encoded into a Trace tensor containing NaN breakpoints and a static mesh object containing an absolute vertex pool and a face index pool. The static mesh object is then written into a WebGL vertex buffer object and an index buffer object.

[0011] Further, borehole coordinates, survey trajectory, and test grade are extracted from the borehole sample data; a spatial coordinate sequence is generated based on the borehole coordinates and the survey trajectory; the spatial coordinate sequence is matched with the test grade by segment to generate a grade coordinate table; the grade coordinate table is subjected to null value removal, coordinate unification, and grade validity screening to obtain the sample point set; the spatial bounding box is generated based on the spatial extrema of the sample point set, and the voxel mesh is generated within the spatial bounding box.

[0012] Furthermore, the method for generating the experimental semivariance set and the point pair count set includes: generating a lag distance sequence based on the Euclidean distance between sample points in the sample point set; establishing binning units according to the lag distance sequence; counting the number of sample point pairs in each binning unit and generating the point pair count set; generating a representative value of the experimental semivariance based on the squared grade difference of the sample point pairs in each binning unit, and generating the experimental semivariance set from the representative value of the experimental semivariance.

[0013] Furthermore, the regularization objective function includes a weighted residual term and a variance constraint term; the weighted residual term consists of the experimental semivariance set, the confidence weight, and the function value of the spherical variogram; the variance constraint term consists of the global variance constraint term, the nugget value, and the abutment difference value; when solving the regularization objective function, the L-BFGS-B optimizer with boundary mapping maps the candidate variogram parameter set to the nugget boundary set, the abutment boundary set, and the range boundary set.

[0014] Further, the method for generating the batch processing capacity includes: writing the sample point set into the cKDTree retrieval topology; querying the neighborhood sample set of the target voxel using the truncated boundary; generating a neighborhood size value based on the sample size of the neighborhood sample set; determining the matrix size of the ordinary Kriging equation system based on the neighborhood size value; reading the available physical memory of the computing node and generating a memory reserve value; and generating the batch processing capacity based on the matrix size, the memory reserve value, the security factor, and the matrix storage unit.

[0015] Further, the method for assembling and solving the ordinary kriging equations according to the batch processing capacity includes: dividing the voxel grid into voxel batches according to the batch processing capacity; generating, for the target voxel within the voxel batch, an inter-sample semivariance matrix, a target sample semivariance vector, and unbiased constraint rows and columns based on the neighborhood sample set and the variation parameter set; and assembling the inter-sample semivariance matrix, the target sample semivariance vector, the unbiased constraint rows and columns, and the Lagrange multipliers into the ordinary kriging equations.

[0016] Furthermore, the method for obtaining the predicted grade voxel set includes: solving the ordinary Kriging equations to obtain a weight vector; generating a grade vector based on the neighborhood sample set; performing a dot product operation between the weight vector and the grade vector to obtain the predicted grade value of the target voxel; writing the spatial coordinates of the target voxel and the predicted grade value into a voxel record, and generating the predicted grade voxel set from the voxel record.

[0017] Further, the method for generating a quality mask and a rendering voxel set based on the predicted quality voxel set includes: generating the quality mask based on the predicted quality values ​​and quality thresholds in the predicted quality voxel set; the quality mask includes a retention flag and a filtering flag; writing the voxels corresponding to the retention flags into a retention rendering queue; generating a rendering quota based on the rendering capacity, and writing the voxels corresponding to the filtering flags into a filtering rendering queue according to the rendering quota; merging the retention rendering queue and the filtering rendering queue to obtain the rendering voxel set.

[0018] Furthermore, encoding the rendering voxel set into a Trace tensor containing NaN breakpoints and a static mesh object containing an absolute vertex pool and a face index pool includes: extracting the corner sequence of the rendering voxel set according to the two-dimensional projection plane, writing the NaN breakpoints between the closed and starting ends of adjacent voxel projection polygons to form the Trace tensor; generating the absolute vertex pool based on the voxel centers and mesh stride of the rendering voxel set, generating the face index pool based on the vertex indices of the absolute vertex pool and the set wrap order, and encapsulating the absolute vertex pool and the face index pool into the static mesh object.

[0019] Furthermore, a kriging modeling method for regularization optimization and topology collapse further includes: after extracting the neighborhood sample set and before generating the inter-sample semivariance matrix, generating neighborhood configuration labels based on the sample spacing sequence and orientation occupancy sequence, and performing isomorphic sample merging or orientation compensation selection according to the neighborhood configuration labels, and updating the neighborhood size value, matrix size and batch processing capacity.

[0020] The technical effects and advantages of the Kriging modeling method for regularization optimization and topological collapse provided by this invention are as follows:

[0021] This invention achieves stable, efficient, and lightweight display of orebody kriging modeling through regularized optimization, variable-range driven computation scheduling, and topology collapse rendering. It generates trust weights by counting point pairs and combines global variance constraints and an L-BFGS-B optimizer with boundary mapping. This ensures that the solution processes for nugget values, sill differences, and variable ranges are simultaneously constrained by sample support and physical boundaries, reducing the interference of sparse binning outlier semivariogram points on the optimization of variation parameters. The variable range is used as the truncation boundary for cKDTree retrieval topology, and the neighborhood sample size is converted into the matrix size of ordinary kriging equations. Combined with available physical memory to generate batch processing capacity, this creates a continuous constraint relationship between neighborhood search, matrix solution, and memory scheduling, reducing the risk of computational interruptions caused by high-density sample clusters. Finally, it generates rendering voxel sets through grade masks and encodes them as Trace tensors with NaN breakpoints and static mesh objects with absolute vertex pools and face index pools, reducing discrete primitive submissions and voxel-by-voxel transformation data transmission, and improving the front-end display continuity of heterogeneous grade voxels. Attached Figure Description

[0022] Figure 1 This is a schematic diagram of a Kriging modeling method for regularization optimization and topology collapse in Example 1.

[0023] Figure 2 This is a schematic diagram illustrating the joint display of the rendered voxel set, formed by processing the predicted grade voxel set with grade masking in Example 1, on a two-dimensional orthographic projection plane and a three-dimensional view.

[0024] Figure 3 This is a schematic diagram of the spherical variogram optimization based on trust weights and variance constraints in Example 1;

[0025] Figure 4 This is a flowchart of the neighborhood configuration verification and matrix size update method in Example 2. Detailed Implementation

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

[0027] Example 1:

[0028] Please see Figure 1 As shown, embodiments of the present invention provide a kriging modeling method for regularization optimization and topology collapse, the method comprising:

[0029] S1. Obtain borehole sample data, perform coordinate cleaning, quality verification and spatial registration on the borehole sample data, generate sample point set, and generate voxel mesh based on the spatial bounding box of the sample point set.

[0030] S2. Based on the sample point set, generate an experimental semivariance set and a point-pair count set by binning according to the lag distance. Determine the trust weight according to the point-pair count set. Write the trust weight, global variance constraint term and spherical variability function into the regularization objective function. The L-BFGS-B optimizer with boundary mapping outputs a set of variability parameters containing nugget value, sill difference and range.

[0031] S3. Using the variable range as the truncation boundary of the cKDTree retrieval topology, extract the neighborhood sample set of the target voxels. Determine the matrix size of the ordinary kriging equation system based on the sample size of the neighborhood sample set. Generate the batch processing capacity based on the matrix size and available physical memory. Assemble and solve the ordinary kriging equation system according to the batch processing capacity to obtain the predicted grade voxel set.

[0032] S4. Generate a quality mask and a rendering voxel set based on the predicted quality voxel set, encode the rendering voxel set into a Trace tensor containing NaN breakpoints, and a static mesh object containing an absolute vertex pool and a face index pool, and write the static mesh object into a WebGL vertex buffer object and an index buffer object.

[0033] In one implementation, such as Figure 2As shown, after completing the 2D projection encoding and 3D mesh encoding, the rendering voxel set can form a joint display result of the 2D orthogonal projection surface and the 3D view. Specifically, based on the spatial coordinates and predicted grade values ​​recorded by each voxel in the rendering voxel set, corner point sequences of the corresponding projection surfaces are generated in the top view, front view, and side view, respectively. NaN breakpoints are written between adjacent voxel projection polygons, so that the discrete voxel projection polygons are displayed in the form of Trace tensors. At the same time, an absolute vertex pool is generated based on the voxel centers and mesh step size in the same rendering voxel set, and a face index pool is generated based on the vertex index of the absolute vertex pool and the set winding order. After encapsulating the absolute vertex pool and the face index pool into a static mesh object, the voxelized ore body model is displayed in the 3D view. Figure 2 The grade classification in the code is used to represent the correspondence between predicted grade values ​​and display colors. Figure 2 The two-dimensional orthographic projection plane and the three-dimensional view in the image both originate from the same set of rendering voxels.

[0034] In this embodiment, the received borehole sample data includes borehole coordinates for identifying the starting position of the borehole, inclinometer trajectory for describing the borehole's underground extension direction, and chemical grade for characterizing the mineralization properties of the sampling section. The borehole coordinates are used to determine at least the starting position of the borehole in a unified spatial coordinate system; the inclinometer trajectory includes at least dip angle and azimuth information corresponding to the borehole depth to determine the spatial orientation of the borehole axis as it changes with depth; and the chemical grade is the mineralization content data corresponding to the sampling section, used as attribute input for subsequent kriging modeling.

[0035] First, the borehole number is used as the associated field to match the borehole coordinates with the survey trajectory. For the same borehole, the survey trajectory is organized according to the order of survey depth. Based on the depth interval, dip angle and azimuth information between adjacent survey records, the incremental direction of the borehole axis in three-dimensional space is determined. Then, starting from the borehole coordinates, the coordinates are accumulated segment by segment along the incremental direction to generate the spatial coordinate sequence corresponding to the borehole. The spatial coordinate sequence refers to the three-dimensional position sequence formed along the borehole depth direction, which is used to convert the original borehole record, which is represented by depth, into a three-dimensional coordinate record that can participate in spatial calculation.

[0036] After generating the spatial coordinate sequence, the test grade is mapped to the spatial coordinate sequence according to the sampling segment. Specifically, the starting depth and ending depth corresponding to the test grade are read, and the representative spatial position of the sampling segment is determined in the spatial coordinate sequence. The representative spatial position can be the midpoint of the segment or the weighted position of the coordinates within the segment. The representative spatial position is bound to the corresponding test grade to form a grade coordinate table. The grade coordinate table includes at least the borehole number, spatial coordinates, segment identifier, and test grade, which is used as the basic data table for subsequent sample point sets.

[0037] Subsequently, the grade coordinate table was subjected to null value removal, coordinate unification, and grade validity screening. Null value removal included deleting records that lacked spatial coordinates, laboratory grades, or corresponding segment relationships. Coordinate unification included converting coordinate fields from different sources to the same spatial coordinate system while maintaining consistent coordinate units. Grade validity screening included removing records that did not conform to the mineralization attribute value rules and merging or selecting duplicate or overlapping segments. After the above processing, a sample point set was generated, and each sample point in the sample point set had three-dimensional spatial coordinates and a corresponding laboratory grade.

[0038] The three-dimensional spatial coordinates of the sample point set are further traversed to obtain the minimum and maximum boundary values ​​in each coordinate direction, and a spatial bounding box covering the sample point set is generated accordingly. The spatial bounding box is used to limit the spatial range of subsequent voxelization calculations. Within the spatial bounding box, the data is discretized according to a preset grid step size to generate a voxel grid. The voxel center in the voxel grid serves as the target voxel position for the subsequent ordinary Kriging equations. For example, when the borehole coordinates, inclinometer trajectory, and test section of a certain borehole have been matched, the midpoint depth of the test section can be converted into spatial coordinates, and the spatial coordinates and the corresponding grade value can be written into the grade coordinate table. However, this example does not limit the selection method of representing the spatial position.

[0039] In this embodiment, after obtaining the sample point set, the three-dimensional spatial coordinates and test grade of each sample point in the sample point set are read first. The three-dimensional spatial coordinates are used to calculate the spatial distance between sample points, and the test grade is used to calculate the attribute difference between sample points. For any pair of sample points, the Euclidean distance is calculated based on their three-dimensional spatial coordinates. The Euclidean distance represents the straight-line interval of the sample point pair in the three-dimensional ore body space. At the same time, the test grade of the sample point pair is read, and the square of the grade difference is calculated as the basis for subsequent semivariance calculation.

[0040] After calculating the distance between sample point pairs, a lag distance sequence is generated based on the Euclidean distance range of the sample point pairs. The lag distance sequence is a distance interval benchmark arranged in the direction of increasing distance, used to characterize the statistical position of spatial correlation as distance changes. Subsequently, binning units are established according to the lag distance sequence, with each binning unit corresponding to a lag distance range. For any sample point pair, if its Euclidean distance falls within the corresponding lag distance range, the sample point pair is assigned to the corresponding binning unit, thereby enabling sample point pairs with close spatial distances to participate in the calculation of the same statistical unit.

[0041] Within each binning unit, the number of sample point pairs assigned to that binning unit is counted to obtain the point pair count value corresponding to that binning unit. The point pair count values ​​of each binning unit are arranged in order of lag distance to form a point pair count set. The point pair count set is used to characterize the data support level of each lag distance position in the experimental semivariogram set, and can be used to determine the confidence weights later.

[0042] Furthermore, within each binning unit, the squared grade difference is calculated for each sample point pair belonging to that binning unit, and the experimental semivariance representative value of that binning unit is obtained based on the squared grade difference. The experimental semivariance representative value is used to characterize the average degree of difference in mineralization grade with spatial distance within the lag distance range. The experimental semivariance representative values ​​of each binning unit are arranged in order of lag distance to form an experimental semivariance set. For example, if there are several sample point pairs within the same binning unit, the squared grade differences of each sample point pair can be summed and normalized according to the number of point pairs to obtain the experimental semivariance representative value of that binning unit. This example does not limit the specific statistical form of the experimental semivariance representative value.

[0043] In one implementation, the first Representative experimental semivariance of each binning unit It can be represented as:

[0044] ;

[0045] in, For the first The center lag distance of each compartment; This represents the total number of data point pairs that fall within the lag interval. It is a three-dimensional spatial coordinate vector; For spatial points The mineralization grade at the location; Distance The mineralization grade at the two spatial coordinate points.

[0046] In this embodiment, after forming the experimental semivariance set and the point-pair count set, the trust weight is first determined based on the point-pair count set. The trust weight is used to represent the degree of data support of the experimental semivariance representative value in each bin. When the number of sample point pairs in a bin is small, the influence of the experimental semivariance representative value corresponding to that bin in the subsequent objective function is correspondingly limited. When the number of sample point pairs in a bin is sufficient, the experimental semivariance representative value corresponding to that bin has a corresponding weight in the subsequent objective function. Thus, the point-pair count set obtained by lag binning is no longer just used as a statistical result, but is still used to construct the regularization objective function.

[0047] The spherical variogram is used to calculate the theoretical semivariogram values ​​for each lag distance position based on the candidate variogram parameter set. The candidate variogram parameter set includes nugget value, sill difference, and range. The nugget value is used to characterize grade variation or measurement disturbance at extremely short distances, the sill difference is used to characterize the increment of semivariogram with increasing distance, and the range is used to characterize the distance scale before spatial correlation reaches a steady state. After each update of the candidate variogram parameter set, the function value of the spherical variogram is calculated based on the candidate variogram parameter set and the lag distance sequence, and the residual is calculated by comparing the function value with the corresponding experimental semivariogram representative value in the experimental semivariogram set.

[0048] In this embodiment, the spherical variogram is used to calculate the theoretical semivariance values ​​for each lag distance position based on the nugget value, sill difference, and range. Subsequently, the experimental semivariance set, confidence weights, global variance constraint term, and spherical variogram are jointly written into the regularization objective function to form the basis for optimizing the candidate variogram parameter set. The theoretical spherical model equation... And define a loss function that aims to minimize the sum of squared residuals. .

[0049] ;

[0050] When solving the regularization objective function, the L-BFGS-B optimizer with boundary constraints is invoked, and the nugget boundary set, sill boundary set, and range boundary set are set respectively. When a candidate value in the candidate mutation parameter set exceeds the corresponding boundary set, the candidate value is mapped back into the corresponding boundary set.

[0051] Gold Value (Indicating measurement error or extremely minute variation): Boundary constraints are defined as ;

[0052] abutment difference (Used to calibrate the population variance magnitude): Boundary constraints are defined as follows ;

[0053] Variable range (Representing the radius of maximum influence of spatial correlation): Boundary constraints are defined as follows .

[0054] The optimizer updates the candidate mutation parameter set based on the gradient information generated during the iteration process, and outputs a mutation parameter set that satisfies the convergence condition of the regularization objective function within the nugget boundary set, sill boundary set, and range boundary set. The range parameters in the mutation parameter set are then used as the truncation boundary for the subsequent cKDTree retrieval topology.

[0055] like Figure 3As shown, in one embodiment, the experimental semivariance set can be represented as discrete semivariance points varying with lag distance, and the function value of the spherical variogram can be represented as a fitted curve corresponding to the discrete semivariance points. Figure 3 The horizontal axis represents the lag distance, and the vertical axis represents the semivariance. The intercept of the curve at the initial position characterizes the nugget value; the increment magnitude corresponding to the curve stabilizing as the lag distance increases characterizes the sill difference; and the distance scale corresponding to the curve before reaching a steady state characterizes the range. Figure 3 It is not only used to represent ordinary variogram curves, but also to illustrate that in this embodiment, the experimental semivariogram set, trust weights, global variance constraint terms, and boundary mappings jointly participate in the optimization process of the candidate variogram parameter set, so that the output range can continue to serve as the truncation boundary of the subsequent cKDTree retrieval topology.

[0056] The regularization objective function includes a weighted residual term and a variance constraint term. The weighted residual term consists of the experimental semivariance set, the confidence weight, and the function value of the spherical variogram. It is used to compare the representative experimental semivariance value of each bin with the theoretical semivariance value in a weighted manner. The variance constraint term consists of a global variance constraint term, nugget value, and sill difference. The global variance constraint term is obtained from the analysis grade statistics of the sample point set and is used to constrain the variance scale after combining the nugget value and the sill difference. Through this construction, the regularization objective function simultaneously includes the experimental semivariance fitting relationship and the overall grade variation constraint relationship.

[0057] To achieve efficient fitting of the theoretical model to the experimental scatter points, a loss function is constructed with the objective of minimizing the sum of squared residuals. And set boundary constraints. :

[0058] ;

[0059] in, This represents the total global variance of the sample obtained from statistics in memory. The penalty coefficient is... This is the adaptive trust weight coefficient.

[0060] When solving the regularization objective function, the L-BFGS-B optimizer with boundary mapping uses the candidate mutation parameter set as the variable to be optimized, and sets the nugget value, stilt difference, and range boundary set respectively. After obtaining a new candidate mutation parameter set in each iteration, it first determines whether the candidate value falls into the corresponding boundary set. For candidate values ​​that exceed the corresponding boundary set, they are mapped back into the corresponding boundary set and then participate in the next objective function calculation. After iterative solution, the output is the mutation parameter set that satisfies the boundary set and makes the regularization objective function converge. This mutation parameter set is then used for subsequent processing of the cKDTree retrieval topology truncation boundary with the range.

[0061] In this embodiment, after obtaining the set of variation parameters, the sample point set is written into the cKDTree retrieval topology. The cKDTree retrieval topology is a tree-shaped spatial retrieval structure built with the three-dimensional spatial coordinates of the sample points as the index key. Each sample point retains its spatial coordinates and corresponding test grade in this structure. When writing, the spatial coordinates in the sample point set are first extracted to form a searchable coordinate sequence. Then, the coordinate sequence is associated with the sample point identifier and test grade, so that the results obtained from subsequent neighborhood queries can be traced back to the corresponding sample point.

[0062] When querying a target voxel in a voxel grid, the range in the mutation parameter set is used as the cutoff boundary. The cutoff boundary represents the maximum spatial correlation distance participating in the ordinary kriging calculation. For any target voxel, a radius query is performed in the cKDTree retrieval topology based on its voxel center coordinates. Only sample points located within the cutoff boundary are extracted as the neighborhood sample set. Thus, the source of the neighborhood sample set is no longer determined by a fixed empirical radius, but by the range limit obtained by the optimization of the preceding mutation function.

[0063] After obtaining the neighborhood sample set, the number of sample points is counted to generate a neighborhood size value. The neighborhood size value is used to determine the matrix size of the ordinary Kriging equation system. When the neighborhood size value increases, the size of the equation system formed by the semivariance matrix between samples, the unbiased constraint rows and columns, and the Lagrange multipliers also increases. Thus, the neighborhood size value is converted into the basis for the storage size required for subsequent matrix solving, rather than just being used as the number of query results.

[0064] Then, the physical memory currently available for this computation on the computing node is read, and a memory margin value is generated. The memory margin value is used to represent the storage space that can be allocated to solving the Kriging matrix while maintaining the stability of the current operating environment. Combined with the safety factor and matrix storage unit, the estimated storage occupancy of the equation system corresponding to a single target voxel is determined, and the batch processing capacity is calculated accordingly. The matrix storage unit refers to the basic storage unit of matrix elements in memory, which, together with the matrix size, determines the storage requirements of a single equation system. The safety factor is used to reserve operating redundancy in the memory margin value.

[0065] The batch processing capacity obtained from the above processing is used to limit the number of target voxels participating in the assembly and solution of the ordinary Kriging equation system at one time. For example, when the borehole samples in a certain region are relatively dense, the neighborhood size value corresponding to the target voxel is larger, the matrix size increases accordingly, and the batch processing capacity decreases accordingly; when the neighborhood size value corresponding to the target voxel is smaller, the batch processing capacity can be increased accordingly. This example is only used to illustrate the linkage between batch processing capacity and neighborhood size value and memory margin value, and does not limit the specific values.

[0066] In this embodiment, after obtaining the batch processing capacity, the voxel grid is first divided according to the batch processing capacity to form voxel batches that participate in the calculation in sequence. A voxel batch refers to the set of target voxels that participate in the processing in the same round of assembly and solution of the ordinary Kriging equations. Its capacity is determined by the aforementioned matrix size and memory margin. For each voxel batch, the voxel center coordinates of each target voxel in the batch are read sequentially, and the aforementioned neighborhood sample set is called as the known sample input of the target voxel.

[0067] For any target voxel within a voxel batch, first read the coordinates of the sample points and the assay grade in its neighborhood sample set, and determine the spherical variogram function based on the nugget value, sill difference, and range in the variogram parameter set. Then, based on the spatial distance between any two sample points in the neighborhood sample set, substitute them into the spherical variogram function to obtain the corresponding inter-sample semivariogram value. According to the arrangement order of the sample points in the neighborhood sample set, write each inter-sample semivariogram value into the corresponding position in the matrix to form the inter-sample semivariogram matrix. The inter-sample semivariogram matrix is ​​used to represent the spatial variation relationship between the neighborhood samples participating in the valuation.

[0068] Furthermore, based on the spatial distance between the center coordinates of the target voxel and the coordinates of each sample point in the neighborhood sample set, the spherical variogram function is substituted to obtain the sample target semivariogram value; according to the arrangement order of the sample points in the neighborhood sample set, the sample target semivariogram value is written into the corresponding position of the column vector to form the sample target semivariogram vector, which is used to represent the spatial variation relationship between the target voxel and each neighborhood sample.

[0069] After forming the inter-sample semivariance matrix and the target semivariance vector, an unbiased constraint column is added. The unbiased constraint column is used to constrain the summation relationship of the ordinary kriging weights, so that the obtained weights satisfy the unbiased condition of ordinary kriging. Corresponding to the unbiased constraint column, Lagrange multipliers are introduced into the equation system. The Lagrange multipliers are added as additional unknowns in the constraint solution, and participate in the solution together with the weight unknowns. Thus, the inter-sample semivariance matrix, the target semivariance vector, the unbiased constraint column, and the Lagrange multipliers are assembled into the ordinary kriging equation system.

[0070] In one example, if the neighborhood sample set of the target voxel contains several sample points, the semivariance matrix between samples is constructed according to the pairwise semivariance values ​​between these sample points, and the target semivariance vector is constructed according to the semivariance values ​​from these sample points to the target voxel. This example is only used to illustrate the source of the matrix and vector, and does not limit the number of neighborhood sample points. After solving the equation system for one voxel batch, the next voxel batch is read according to the batch processing capacity to continue processing.

[0071] In this embodiment, after the ordinary kriging equation system is assembled, the ordinary kriging equation system is solved to obtain the weight vector corresponding to the target voxel. Each weight value in the weight vector corresponds one-to-one with the sample point in the neighborhood sample set, and is used to represent the contribution ratio of the test grade of the corresponding sample point in the prediction of the target voxel. Since the ordinary kriging equation system already contains unbiased constraint rows and columns and Lagrange multipliers, the weight vector obtained by solving inherits the spatial variation relationship expressed by the aforementioned inter-sample semivariance matrix and sample-target semivariance vector.

[0072] After generating the weight vector, the test grade of each sample point is read according to the arrangement order of the sample points in the neighborhood sample set when assembling the ordinary Kriging equation system, forming a grade vector. The element order of the grade vector is consistent with the element order of the weight vector, so that each weight value corresponds to the test grade of its source sample point. If the neighborhood sample set contains sample point identifiers in the previous query, the grade coordinate table is backtracked according to the sample point identifiers, the corresponding test grade is read and written into the grade vector.

[0073] Subsequently, the weight vector and the grade vector are multiplied by a dot product to obtain the predicted grade value of the target voxel. The dot product operation multiplies the test grade of each sample point with its corresponding weight value and combines the resulting products to form the grade estimate at the location of the target voxel. This predicted grade value and the spatial coordinates of the target voxel together constitute a voxel record. The voxel record includes at least the target voxel identifier, the spatial coordinates of the target voxel, and the predicted grade value. When it is necessary to retain the calculation source, the neighboring sample set identifier or the voxel batch identifier can also be recorded.

[0074] After performing the above processing on each target voxel in a voxel batch, a corresponding voxel record set is formed. After each voxel batch is processed, the voxel record sets are summarized according to the spatial positional relationship in the voxel grid to generate a predicted grade voxel set. For example, if the neighborhood sample set corresponding to a target voxel contains several sample points, the predicted grade value of the target voxel is obtained by merging the test grade of the several sample points with the corresponding weight value. This example is only used to illustrate the correspondence between the weight vector and the grade vector and does not limit the number of neighborhood sample points.

[0075] In this embodiment, after generating the predicted grade voxel set, the predicted grade value and spatial coordinates in each voxel record are read first. The predicted grade voxel set is the set of voxel records obtained by the previous ordinary kriging calculation. Each voxel record includes at least voxel spatial coordinates and predicted grade value. The grade threshold is a grade boundary condition predetermined according to the mineral evaluation rules, modeling task or display requirements, used to distinguish voxels that need to maintain spatial morphology first and voxels that can participate in screening.

[0076] Based on the predicted grade value and grade threshold, a grade mask is generated for each voxel in the predicted grade voxel set. The grade mask is a marker sequence corresponding to the voxel record, including a retention marker and a filtering marker. For voxels whose predicted grade value meets the grade threshold condition, a retention marker is written; for voxels whose predicted grade value does not meet the grade threshold condition, a filtering marker is written. Through this marker sequence, the grade attribute of the voxel is converted into the selection criterion for subsequent rendering queues.

[0077] After the quality mask is formed, queue writing is performed according to the correspondence between the quality mask and the predicted quality voxel set. Voxels corresponding to the reserved markers are written to the reserved rendering queue. The reserved rendering queue is used to store voxel records that need to participate in subsequent rendering encoding in their entirety. For voxels corresponding to the filter markers, they are not directly written to the rendering result. Instead, a rendering quota is generated based on the rendering capacity. The rendering capacity represents the range of voxels that the current rendering end can receive, and the rendering quota represents the capacity that can be allocated to the voxels with the filter markers after deducting the reserved rendering queue.

[0078] Subsequently, based on the rendering quota, voxels to be rendered are selected from the voxels corresponding to the filter markers and written into the filter rendering queue. The selection can be performed according to the spatial order of voxels in the voxel grid or according to the sorting result of the voxel records, so as to maintain the correspondence between the filter rendering queue and the spatial distribution of voxels. After the queue writing is completed, the reserved rendering queue and the filter rendering queue are merged to obtain the rendering voxel set. The rendering voxel set continues to be used for subsequent encoding processing of Trace tensors containing NaN breakpoints and static mesh objects.

[0079] For example: when the rendering capacity is insufficient to accommodate all predicted quality voxels, voxels that meet the quality threshold conditions enter the reserved rendering queue, and the remaining voxels enter the filtered rendering queue according to the remaining rendering quota; this example is only used to illustrate the queue allocation relationship between the reserved marker and the filtered marker, and does not limit the specific threshold or quota value.

[0080] In this embodiment, after obtaining the rendering voxel set, the voxel center, mesh step size and predicted quality value of each voxel record in the rendering voxel set are read first. The rendering voxel set is obtained by merging the reserved rendering queue and the filtered rendering queue. Its voxel records have completed quality mask filtering and serve as the common input for subsequent two-dimensional projection coding and three-dimensional mesh coding.

[0081] For 2D projection encoding, the rendering voxel set is read according to the preset 2D projection plane. The 2D projection plane can correspond to any orthogonal direction of the voxel mesh and is used to map 3D voxels into planar polygons. For each voxel, the corner coordinates of the voxel on the 2D projection plane are calculated according to its voxel center and mesh step size, and written into the corner sequence according to the polygon closure order. The corner sequence is the coordinate arrangement result of each voxel projection polygon. In order to avoid adjacent but unconnected projection polygons being continuously connected in the same trajectory, NaN breakpoints are written between the closed end of one voxel projection polygon and the starting end of the next voxel projection polygon. NaN breakpoints are non-numerical identifiers used to separate adjacent projection polygons during trajectory parsing. The corner sequence and NaN breakpoints together form the Trace tensor. The Trace tensor is used to carry discrete projection polygons with continuous coordinate sequences.

[0082] For 3D mesh encoding, the corner coordinates of each voxel in 3D space are calculated based on the voxel center and mesh step size of each voxel in the rendering voxel set. The corner coordinates are then added to the absolute vertex pool in the order they are written. The absolute vertex pool is a set of vertices represented by a unified spatial coordinate system, where the vertex coordinates do not depend on individual voxel transformation matrices. Subsequently, based on the vertex index in the absolute vertex pool and the set wrap order, a corresponding face index is generated for each voxel. The set wrap order is used to determine the connection order of the face vertices. The face index pool is used to record the vertex indexes referenced by each face. Thus, the geometric position of the voxel is saved by the absolute vertex pool, and the voxel surface connection relationship is saved by the face index pool.

[0083] After completing the above processing, the absolute vertex pool and face index pool are encapsulated into a static mesh object. The static mesh object contains vertex data and index data, but does not contain a per-voxel independent transformation matrix. Then, the vertex data in the static mesh object is written to a WebGL vertex buffer object, and the index data is written to a WebGL index buffer object, so that subsequent rendering calls can draw a 3D voxel model based on the spatial coordinates in the vertex buffer object and the face connectivity in the index buffer object. For example, a cube voxel can be represented by several corner coordinates and several face indices. This example is only used to illustrate the correspondence between the absolute vertex pool and the face index pool, and does not limit the specific number of voxels or the mesh step size.

[0084] Example 2:

[0085] like Figure 4As shown, this embodiment further improves upon the design of Embodiment 1. The difference lies in the fact that, in the actual operation of Embodiment 1, it was found that the neighborhood sample set extracted by cKDTree topology retrieval exhibits situations where sample points are concentrated along a single borehole trajectory, sample point spacing is too close, or sample direction distribution is biased in local dense borehole areas or spatial boundary areas. This leads to enhanced row and column correlation in the semivariance matrix generated based on this neighborhood sample set. Although the matrix size of the ordinary Kriging equation system is limited by batch processing capacity, its numerical stability may still decrease, failing to maintain a stable weight solution relationship between locally high-density sample clusters and boundary sparse voxels. Based on this, a Kriging modeling method based on regularization optimization and topology collapse further includes:

[0086] After extracting the neighborhood sample set of the target voxel using the variable range as the cutoff boundary, and before generating the inter-sample semivariance matrix, a neighborhood configuration check is performed on the neighborhood sample set to generate a neighborhood configuration label. The neighborhood configuration check includes: generating a sample spacing sequence based on the Euclidean distance between sample points in the neighborhood sample set, generating an orientation occupancy sequence based on the spatial direction from the target voxel to each sample point, and determining whether there is a sample close clustering or orientation concentration state in the neighborhood sample set based on the sample spacing sequence and orientation occupancy sequence. Sample close clustering refers to the presence of sample points in the neighborhood sample set that are spatially too close and form similar row and column contributions in the inter-sample semivariance matrix; orientation concentration state refers to the sample points in the neighborhood sample set being mainly distributed on the same spatial lateral direction or the same borehole extension zone of the target voxel.

[0087] When the neighborhood configuration marker indicates that the neighborhood sample set is in a state of close sample clustering, isomorphic sample merging is performed on the neighborhood sample set. The isomorphic sample merging process includes: classifying sample points that are spatially close and have continuous sources in the test area into isomorphic sample groups; generating group representative coordinates based on the spatial coordinates of sample points within the isomorphic sample group; generating group representative grades based on the test grades of sample points within the isomorphic sample group; and replacing the isomorphic sample group with the group representative coordinates and group representative grades to participate in subsequent matrix assembly. The group representative coordinates and group representative grades are used to maintain the representative relationship of the isomorphic sample group in terms of spatial location and mineralization properties.

[0088] When the neighborhood configuration marker indicates that the neighborhood sample set is in a state of orientation concentration, without breaking the truncation boundary corresponding to the range, orientation compensation selection is performed on the neighborhood sample set according to the spatial position of the target voxel. Orientation compensation selection includes: partitioning the spatial orientation around the target voxel, identifying orientation partitions that lack sample point support, and selecting sample points located in the orientation partition from the candidate sample points within the truncation boundary to add to the neighborhood sample set; if there are no candidate sample points in the orientation partition, the original neighborhood sample set is retained and an orientation missing marker is written.

[0089] After isomorphic sample merging or direction compensation selection, the neighborhood size value is regenerated based on the updated neighborhood sample set, and the matrix size of the ordinary kriging system is determined by the regenerated neighborhood size value. If the updated matrix size changes the batch processing capacity, the voxel batch is re-divided according to the updated batch processing capacity. Then, the semivariance matrix between samples, the semivariance vector of the sample target, and the unbiased constraint rows and columns are generated based on the updated neighborhood sample set and the mutation parameter set, and the assembly and solution of the ordinary kriging system are continued.

[0090] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

[0091] The above description is merely a preferred embodiment of the present application, but the scope of protection of the present application is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present application, based on the technical solution and concept of the present application, should be covered within the scope of protection of the present application.

Claims

1. A Kriging modeling method for regularization optimization and topological collapse, characterized in that, The methods include: Obtain borehole sample data, perform coordinate cleaning, quality verification and spatial registration on the borehole sample data, generate sample point set, and generate voxel mesh based on the spatial bounding box of the sample point set; Based on the sample point set, an experimental semivariance set and a point-pair count set are generated by binning according to lag distance. The trust weight is determined according to the point-pair count set. The trust weight is used to represent the degree of data support for the representative value of the experimental semivariance in each lag distance bin. The global variance constraint term is obtained based on the laboratory grade statistics of the sample point set. The global variance constraint term is used to constrain the variance scale after combining the nugget value and the abutment difference. The trust weight, the global variance constraint term, and the spherical variogram function are written into the regularization objective function, and the L-BFGS-B optimizer with boundary mapping outputs a set of variogram parameters containing the nugget value, the abutment difference, and the range. Using the range as the cutoff boundary of the cKDTree retrieval topology, the neighborhood sample set of the target voxel is extracted. The matrix size of the ordinary kriging equation system is determined based on the sample size of the neighborhood sample set. The batch processing capacity is generated based on the matrix size and available physical memory. The ordinary kriging equation system is assembled and solved according to the batch processing capacity to obtain the predicted grade voxel set. A quality mask and a rendering voxel set are generated based on the predicted quality voxel set. Topology collapse is performed on the rendering voxel set, which includes two-dimensional projection encoding and three-dimensional mesh encoding. The two-dimensional projection encoding includes extracting the corner sequence of the rendering voxel set according to the two-dimensional projection plane, writing NaN breakpoints between the closed and starting ends of adjacent voxel projection polygons, and forming a Trace tensor from the corner sequence and the NaN breakpoints. The NaN breakpoints are used to separate adjacent voxel projection polygons during trajectory resolution. The three-dimensional mesh encoding includes generating an absolute vertex pool based on the voxel centers and mesh stride of the rendering voxel set, generating a face index pool based on the vertex indices of the absolute vertex pool and the set wrap order, encapsulating the absolute vertex pool and the face index pool into a static mesh object, and writing the static mesh object into a WebGL vertex buffer object and an index buffer object.

2. The Kriging modeling method for regularization and topological collapse according to claim 1, characterized in that, The borehole coordinates, inclination trajectory, and test grade are extracted from the borehole sample data; a spatial coordinate sequence is generated based on the borehole coordinates and the inclination trajectory; the spatial coordinate sequence is matched with the test grade by segment to generate a grade coordinate table; the grade coordinate table is subjected to null value removal, coordinate unification, and grade validity screening to obtain the sample point set; The spatial bounding box is generated based on the spatial extrema of the sample point set, and the voxel mesh is generated within the spatial bounding box.

3. The Kriging modeling method for regularization and topological collapse according to claim 1, characterized in that, The method for generating the experimental semivariance set and the point pair count set includes: generating a lag distance sequence based on the Euclidean distance between sample points in the sample point set; establishing binning units according to the lag distance sequence; counting the number of sample point pairs in each binning unit and generating the point pair count set; generating a representative value of the experimental semivariance based on the squared grade difference of the sample point pairs in each binning unit, and generating the experimental semivariance set from the representative value of the experimental semivariance.

4. The Kriging modeling method for regularization and topological collapse according to claim 3, characterized in that, The regularization objective function includes a weighted residual term and a variance constraint term; the weighted residual term consists of the experimental semivariance set, the confidence weight, and the function value of the spherical variogram; the variance constraint term consists of the global variance constraint term, the nugget value, and the abutment difference value; when solving the regularization objective function, the L-BFGS-B optimizer with boundary mapping maps the candidate variogram parameter set to the nugget boundary set, the abutment boundary set, and the range boundary set.

5. The Kriging modeling method for regularization and topological collapse according to claim 1, characterized in that, The batch processing capacity generation method includes: writing the sample point set into the cKDTree retrieval topology; querying the neighborhood sample set of the target voxel using the truncated boundary; generating a neighborhood size value based on the sample size of the neighborhood sample set; determining the matrix size of the ordinary Kriging equation system based on the neighborhood size value; reading the available physical memory of the computing node and generating a memory reserve value; and generating the batch processing capacity based on the matrix size, the memory reserve value, the security factor, and the matrix storage unit.

6. The Kriging modeling method for regularization and topological collapse according to claim 5, characterized in that, The method for assembling and solving the ordinary kriging equations according to the batch processing capacity includes: dividing the voxel grid into voxel batches according to the batch processing capacity; generating, for the target voxel within the voxel batch, a sample semivariance matrix, a sample target semivariance vector, and unbiased constraint rows and columns based on the neighborhood sample set and the variation parameter set; and assembling the sample semivariance matrix, the sample target semivariance vector, the unbiased constraint rows and columns, and the Lagrange multipliers into the ordinary kriging equations.

7. The Kriging modeling method for regularization and topological collapse according to claim 6, characterized in that, The method for obtaining the predicted grade voxel set includes: solving the ordinary Kriging equations to obtain a weight vector; generating a grade vector based on the neighborhood sample set; performing a dot product operation between the weight vector and the grade vector to obtain the predicted grade value of the target voxel; writing the spatial coordinates of the target voxel and the predicted grade value into a voxel record, and generating the predicted grade voxel set from the voxel record.

8. The Kriging modeling method for regularization and topological collapse according to claim 1, characterized in that, The method for generating a quality mask and a rendering voxel set based on the predicted quality voxel set includes: generating the quality mask based on the predicted quality values ​​and quality thresholds in the predicted quality voxel set; the quality mask includes a retention flag and a filtering flag; writing the voxels corresponding to the retention flags into a retention rendering queue; generating a rendering quota based on the rendering capacity, and writing the voxels corresponding to the filtering flags into a filtering rendering queue according to the rendering quota; merging the retention rendering queue and the filtering rendering queue to obtain the rendering voxel set.

9. The Kriging modeling method for regularization and topological collapse according to claim 6, characterized in that, Also includes: After extracting the neighborhood sample set and before generating the inter-sample semivariance matrix, a neighborhood configuration label is generated based on the sample spacing sequence and orientation occupancy sequence. Isomorphic sample merging or orientation compensation selection is then performed according to the neighborhood configuration label, and the neighborhood size value, matrix size, and batch processing capacity are updated.

Images