Layer-controlled lead-zinc deposit resource reserve evaluation method and system based on three-dimensional model

Abstract

The present application relates to the technical field of three-dimensional geological modeling, in particular to a strata-bound lead-zinc deposit resource reserve evaluation method and system based on a three-dimensional model, comprising the following steps: calculating a stratum grid normal and dip vector, constructing a unit set and giving a vector to generate an occurrence vector field, screening the same lithology sample to construct a list and establish a local coordinate system, analyzing the sample relative position vector mapping distance components to the local coordinate system and normalizing and synthesizing structural distance, solving weight coefficients based on the structural distance and weightedly estimating grade, and combining density and volume to calculate resource quality and accumulate and output. In the present application, the grade evaluation is constrained in the stratum extension direction by constructing the occurrence vector field and the local coordinate system, the spatial correlation is corrected by using anisotropic structural distance, the evaluation main shaft is adaptively rotated with the fold, the cross-layer diffusion phenomenon is eliminated, the natural occurrence state of the ore body is restored, and the reserve estimation accuracy and model fidelity are improved.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional geological modeling technology, and in particular to a method and system for assessing the resource reserves of stratabound lead-zinc deposits based on three-dimensional models. Background Technology

[0002] The field of 3D geological modeling technology encompasses the use of computer graphics principles and spatial data interpolation algorithms to digitally reconstruct and visualize underground geological structures, stratigraphic distribution, and ore body morphology in three dimensions. Traditional methods for assessing the resource reserves of stratabound lead-zinc deposits involve geological engineers first obtaining core and cuttings samples at different depths of the deposit through drilling, trenching, or pitting. Lead and zinc grades are then determined through chemical analysis. Subsequently, primarily using the geological block method or parallel section method, ore body boundaries are manually delineated on paper or electronic 2D geological profiles based on boundary grades and minimum industrial indicators. The cross-sectional area of ​​the ore body within each profile is measured using a planimeter or auxiliary design software. The volume of the ore body between adjacent profiles is then calculated using trapezoidal or truncated cone formulas. Finally, the resource reserves of the deposit are estimated by combining this with the average ore weight.

[0003] Existing technologies mainly rely on two-dimensional profiles to artificially delineate ore bodies and combine them with geometric formulas to estimate reserves. This approach severs the three-dimensional spatial connection of geological data and ignores the anisotropic characteristics of ore bodies in stratabound deposits, which are strictly distributed along the bedding of strata and undergo bending deformation with tectonic folds. This makes it difficult to accurately restore the true shape of ore bodies when dealing with structurally complex areas. Using simple geometric formulas to connect adjacent profiles cannot adapt to drastic changes in the occurrence of ore bodies, which can easily lead to blurred ore body boundaries and deviations in reserve calculations, making it difficult to meet the needs of refined resource assessment. Summary of the Invention

[0004] To achieve the above objectives, the present invention adopts the following technical solution: a method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model, comprising the following steps: S1: Calculate the normal vector and dip vector by calling the stratigraphic triangular mesh data, construct a set of elements, match associated patches for the elements in the set of elements, and assign the normal vector and dip vector of the associated patches to the elements to generate a three-dimensional attitude vector field; S2: Extract target units from the unit set, filter core sample data that match the lithology of the target units to construct a list of valid samples, and use the normal vector and the dip vector to construct a local coordinate system centered on the target unit; S3: Calculate the relative position vector of the effective sample list relative to the target unit, map the relative position vector to the local coordinate system to analyze the distance components in each direction, and use the range parameter to normalize the distance components in each direction to synthesize the constructed distance; S4: Substitute the structural distance into the spherical variogram model to calculate the semivariogram value, solve the weighting coefficient based on the semivariogram value, use the weighting coefficient to calculate the estimated grade value by weighting the grade values ​​of the effective sample list, combine the average rock density value, the estimated grade value and the target unit volume to calculate the metal resource quality and accumulate and output the total reserve value.

[0005] As a further aspect of the present invention, step S1 specifically comprises: S11: Traverse all triangular patches contained in the formation triangular mesh data, extract the spatial coordinate data of the vertices of the triangular patches, use the cross product algorithm to calculate the unit normal vector perpendicular to the triangular patch, and project the unit normal vector onto the horizontal plane to solve the unit dip vector indicating the dip direction of the formation. S12: Discretize the three-dimensional space based on the boundary range of the ore body to generate the unit set composed of regular hexahedral voxels, calculate the three-dimensional coordinates of the geometric center point of each voxel in the unit set, and establish a KD tree data structure for accelerating spatial indexing. S13: For each voxel in the unit set, use the KD tree data structure to search for the triangular facet with the closest Euclidean distance in the stratigraphic triangular mesh data as the associated facet, extract the unit normal vector and the unit dip vector of the associated facet and assign them to the corresponding voxel, and traverse the entire region to generate the three-dimensional attitude vector field.

[0006] As a further aspect of the present invention, step S2 specifically comprises: S21: Obtain the geological attribute code of the target unit, traverse the core sample records in the borehole database of the whole area, compare the lithology code of the core sample with the geological attribute code of the target unit, retain samples with consistent codes and spatial distance within the preset search radius, and establish the list of effective samples. S22: Obtain the normal vector and dip vector corresponding to the target unit in the three-dimensional attitude vector field, define the dip vector as the local U-axis, define the normal vector as the local W-axis, and generate the local V-axis by using the cross product operation of the normal vector and the dip vector, thereby constructing the local coordinate system that can characterize the local anisotropy of the stratabound ore body.

[0007] As a further aspect of the present invention, step S3 specifically comprises: S31: Obtain the spatial coordinates of each sample in the list of valid samples and the center coordinates of the target unit respectively. Obtain the global relative position vector of each sample relative to the target unit through vector subtraction. Project the global relative position vector onto the U-axis, V-axis and W-axis of the local coordinate system respectively. Analyze the distance components along the dip, along the direction and along the normal. S32: Obtain the preset variation function range ellipsoid parameters, extract the main range, secondary range and vertical range as normalization factors respectively, and use the normalization factors to perform dimensionless processing on the distance components in each direction according to the anisotropic distance synthesis formula and calculate the Euclidean norm to generate the constructed distance. S33: Store the calculated constructed distances one-to-one with the samples in the effective sample list, and establish a distance weight calculation dataset containing the mapping relationship between sample identifiers and constructed distances, so as to provide a constructed corrected distance parameter basis for subsequent spatial interpolation using the spherical variogram model.

[0008] As a further aspect of the present invention, step S4 specifically comprises: S41: Based on the constructed distance and the preset nugget value and sill value parameters, call the spherical variogram model to calculate the semivariance values ​​between each sample and the target unit and between each sample, and construct the semivariance matrix and right-hand covariance vector required for the Kriging linear equation system. S42: Introduce a Lagrange multiplier factor to satisfy the unbiased weight constraint, perform inverse matrix operation or Gaussian elimination on the Kriging equation system containing the semivariance matrix to obtain the weight coefficient corresponding to each sample in the effective sample list, and ensure that the algebraic sum of all the weight coefficients is 1. S43: Using the weighting coefficients, perform a linear weighted summation of the test grades of the corresponding core samples in the list of valid samples to obtain the estimated grade value of the target unit. Multiply the estimated grade value by the volume of the target unit and the average density of the rock to obtain the metal resource quality of the unit. Iterate through all units and sum them up to obtain the total reserves of the deposit.

[0009] As a further aspect of the present invention, the process of constructing the local coordinate system capable of characterizing the local anisotropic features of the stratabound ore body as described in S22 specifically includes: Obtain the normal vector and the yaw vector of the target unit, and normalize the two vectors to eliminate the difference in magnitude; The normalized directional vector is set as the primary extension direction axis, and the normalized normal vector is set as the thickness variation direction axis. By calculating the vector product of the primary extension direction axis and the thickness variation direction axis, a secondary extension direction axis orthogonal to the first two axes is generated. An orthogonal rotation matrix is ​​constructed using the primary extension direction axis, the secondary extension direction axis, and the thickness variation direction axis, and the center coordinates of the target unit are set as the origin to establish the local coordinate system.

[0010] As a further aspect of the present invention, the process of generating the constructed distance according to the anisotropic distance synthesis formula in step S32 is specifically performed according to the following formula: ; in, Represents the constructed distance, Represents the distance component along the dip. Represents the distance component along the direction. Represents the distance component along the normal direction. Represents the main variable range parameter. Represents the secondary range parameters. This represents the vertical range parameter.

[0011] As a further aspect of the present invention, the process of calculating the semivariogram value by calling the spherical variogram model in step S41 is specifically performed according to the following formula: ; in, This represents the semivariance value. Represents the constructed distance, Represents the nugget constant. Represents the height of the arch. This represents the normalized range threshold.

[0012] As a further aspect of the present invention, the process of calculating and estimating the grade value by weighting the grade values ​​of the effective sample list using the weighting coefficients in step S43 specifically includes: Extract the lead and zinc grade test values ​​of each sample in the list of valid samples, and multiply the test value by the weight coefficient of the corresponding sample obtained in S42 to obtain the weighted component. The weighted components of all samples are summed. If the summation result is negative, it is truncated and corrected using a preset lower limit value to generate the estimated grade value of the target unit. Determine whether the estimated grade value is greater than the preset boundary grade index. If it is greater, mark the unit as an industrial ore body unit and participate in the subsequent metal resource quality calculation. Otherwise, mark the unit as a surrounding rock unit and exclude it from the reserve calculation range.

[0013] A 3D model-based system for assessing the resource reserves of stratabound lead-zinc deposits is provided. This system is used to implement the aforementioned 3D model-based method for assessing the resource reserves of stratabound lead-zinc deposits. The system includes: The vector field generation module is used to call the stratigraphic triangular mesh data to calculate the normal vector and dip vector, construct a set of cells, match associated patches for the cells in the set of cells, and assign the normal vector and dip vector of the associated patches to the cells to generate the three-dimensional attitude vector field. The coordinate system construction module is used to extract target units from the unit set, filter core sample data that match the lithology of the target units to construct the effective sample list, and construct a local coordinate system centered on the target unit using the normal vector and the dip vector. The construction distance analysis module is used to calculate the relative position vector of the effective sample list relative to the target unit, map the relative position vector to the local coordinate system to analyze the distance components in each direction, and use the range parameter to normalize the distance components in each direction to synthesize the construction distance. The reserve estimation module is used to substitute the structural distance into the spherical variogram model to calculate the semivariogram value, solve the weighting coefficient based on the semivariogram value, use the weighting coefficient to weight the grade values ​​of the effective sample list to calculate the estimated grade value, combine the average rock density value, the estimated grade value and the target unit volume to calculate the metal resource quality and accumulate and output the total reserve value.

[0014] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, by constructing a three-dimensional attitude vector field representing the morphology of strata and establishing a local coordinate system that varies with the undulation of strata, the grade estimation calculation is forcibly constrained to the direction of strata extension. By using anisotropic structural distance to dynamically correct the spatial correlation between sample points, the main axis of the estimation search is made to adaptively rotate with the attitude of strata folds. This effectively eliminates the grade diffusion phenomenon caused by fixed direction at structural turning points. The natural occurrence state of the ore body with the curvature of strata is accurately restored in the digital model, improving the accuracy of resource reserve estimation in complex structural areas and the fidelity of geological models. Attached Figure Description

[0015] Figure 1 This is the main flowchart of the stratabound lead-zinc deposit resource reserve assessment method based on a three-dimensional model according to the present invention; Figure 2This is a flowchart of the three-dimensional attitude vector field generation process of the present invention; Figure 3 This is a flowchart illustrating the local coordinate system construction and sample screening process of this invention. Figure 4 A flowchart for distance synthesis calculation is provided for this invention. Figure 5 This is a flowchart illustrating the resource reserve estimation and classification process of this invention. Detailed Implementation

[0016] To make the objectives, technical solutions, and advantages of this invention clearer, the software-based technical solution is described in detail below with reference to system architecture diagrams and embodiments. It should be understood that the specific embodiments described herein are only for explaining the technical solutions of this invention and do not constitute a limitation on the scope of protection.

[0017] In the description of this invention, the system architecture relationships or data processing flows indicated by terms such as "layer," "module," "interface," "data flow," "client," and "server" are all defined based on the architecture diagram or flowchart corresponding to the embodiments. This way of describing is only used to clearly illustrate the logical relationships between the elements in the technical solution, and not to limit the physical deployment form. The term "multiple" includes two or more technical units, including but not limited to multiple data nodes, processing threads, service instances, or functional components and other scalable elements. The specific number is determined according to the actual business scenario and needs to be specifically specified.

[0018] Please see Figure 1 and Figure 2 This invention provides a technical solution: a method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model, comprising the following steps: S1: Calculate the normal and dip vectors by calling the stratigraphic triangular mesh data, construct a set of elements, match associated patches for the elements in the set of elements, and assign the normal and dip vectors of the associated patches to the elements to generate a three-dimensional attitude vector field; The specific steps of S1 are as follows: S11: Traverse all triangular patches contained in the stratigraphic triangular mesh data, extract the spatial coordinate data of the vertices of the triangular patches, use the cross product algorithm to calculate the unit normal vector perpendicular to the triangular patch, and project the unit normal vector onto the horizontal plane to solve the unit dip vector indicating the dip direction of the stratigraphy. S12: Discretize the three-dimensional space based on the boundary range of the ore body to generate a set of cells composed of regular hexahedral voxels. Calculate the three-dimensional coordinates of the geometric center point of each voxel in the set of cells and establish a KD tree data structure to accelerate spatial indexing. S13: For each voxel in the unit set, use the KD tree data structure to search for the triangular facet with the closest Euclidean distance in the stratigraphic triangular mesh data as the associated facet. Extract the unit normal vector and unit dip vector of the associated facet and assign them to the corresponding voxel. Traverse the entire region to generate a three-dimensional attitude vector field.

[0019] For processing stratigraphic triangular mesh data and constructing a 3D attitude vector field, the system first initializes the graphics workstation memory area and loads the stratigraphic interface triangular mesh file covering the target mining area. The system reads the file header information and establishes a list of topological structures containing triangular facets. For step S11, the system starts a computation thread to index and traverse each triangular facet in the list. Taking the triangular facet numbered T_4025 as an example, the system extracts its three vertices. The three-dimensional spatial coordinates are respectively , and .

[0020] System-defined vector with vector Performing vector subtraction yields and .

[0021] Then, the cross product function from the linear algebra library is called to calculate the normal vector. The calculation process involves performing determinant calculations: Quantity ; Quantity ; Quantity The non-unit normal vector is obtained. .

[0022] The system further calculates the magnitude of the vector. Then perform a normalized division operation to obtain the unit normal vector. This unit normal vector represents the vertical variation trend of a local stratum. Next, the system will... Projected onto the horizontal plane ,make When the component is 0, the projection vector is obtained. Again After normalization, the unit dip vector indicating the maximum dip direction of the strata is obtained. This process is repeated until all faces have been traversed, generating a face orientation attribute table. The aforementioned cross product function is a mathematical algorithm that operates on two vectors in a vector space to generate a third vector. The generated third vector is perpendicular to the plane formed by the first two vectors, its direction follows the right-hand rule, and its magnitude is equal to the area of ​​the parallelogram formed by the first two vectors.

[0023] In step S12, the system sets the discretized boundary in three-dimensional space based on the maximum boundary range delineated by the ore body exploration line profile: Axis range , Axis range , Axis range The system reads the preset voxel resolution parameters and sets them. The directional step length is 10 meters. The directional step length is 10 meters. The directional step size is 5 meters. The program executes a triple nested loop to divide the area into regular hexahedral meshes within the aforementioned boundaries. (Index number) Taking a voxel as an example, the system calculates the coordinates of its geometric center point using a formula. After the loop completes, the system generates a unit set object containing 800,000 voxels. To address the efficiency issue of matching massive amounts of data, the system invokes a spatial indexing algorithm, using the geometric center points of all triangular facets extracted in S11 as the dataset to construct a KD-tree data structure. During the construction process, the system calculates the dataset's... The variances of the three dimensions are used to select the dimension with the largest variance as the splitting surface. The median point is selected to split the dataset into a left subtree and a right subtree. This process is recursively executed until the number of faces contained in the leaf node is less than a preset threshold (e.g., 10), thereby establishing a binary tree structure for fast nearest neighbor search.

[0024] The aforementioned KD-tree data structure refers to a binary tree partitioning data structure used to organize a set of points in a k-dimensional space. It recursively divides the multidimensional space into multiple sub-regions using hyperplanes perpendicular to the coordinate axes, thereby achieving logarithmic time complexity optimization in range search and nearest neighbor search operations.

[0025] In step S13, the system initiates a parallel batch processing task, distributing voxels from the cell set to multiple computing cores. Regarding the aforementioned center point... For the target voxel, the system performs a nearest neighbor search in the KD-tree. The search algorithm starts from the root node and calculates... The distance from the dividing plane is used to determine the search path. During the backtracking process, the search branches are pruned using the current nearest distance, and finally the triangle with the closest Euclidean distance is located. Assume the nearest face found is the one described above. System read Associated unit normal vector With unit tendency vector The two depth vectors are directly copied to the attribute field of the target voxel. For voxels located in areas of drastic attitude change, such as faults or fold axes, if the distance to the nearest facet exceeds a preset confidence threshold (e.g., 50 meters), the system triggers an interpolation correction subroutine to search for the five nearest facets, synthesize the average normal and dip vectors using an inverse distance weighting method, and assign them to the voxel. After traversing all voxels in the entire region, a complete three-dimensional attitude vector field is formed in memory.

[0026] Please see Figure 1 and Figure 3 S2: Extract target units from the unit set, filter core sample data that match the lithology of the target units to construct a list of valid samples, and use normal vectors and dip vectors to construct a local coordinate system centered on the target units; The specific steps of S2 are as follows: S21: Obtain the geological attribute code of the target unit, traverse the core sample records in the borehole database of the whole area, compare the lithology code of the core sample with the geological attribute code of the target unit, retain the samples with the same code and the spatial distance within the preset search radius, and establish a list of valid samples. S22: Obtain the normal vector and dip vector corresponding to the target unit in the three-dimensional attitude vector field, define the dip vector as the local U-axis, define the normal vector as the local W-axis, and generate the local V-axis by using the cross product operation of the normal vector and the dip vector, thereby constructing a local coordinate system that can characterize the local anisotropic characteristics of the stratabound ore body; The process of constructing a local coordinate system, S22, capable of characterizing the local anisotropy of stratabound ore bodies, specifically includes: Obtain the normal vector and dip vector of the target element, and normalize the two vectors to eliminate the difference in magnitude; The normalized dip vector is set as the primary extension direction axis, and the normalized normal vector is set as the thickness variation direction axis. By calculating the vector product of the primary extension direction axis and the thickness variation direction axis, a secondary extension direction axis orthogonal to the first two axes is generated. An orthogonal rotation matrix is ​​constructed using the primary extension direction axis, the secondary extension direction axis, and the thickness variation direction axis, and the center coordinates of the target element are set as the origin to establish a local coordinate system.

[0027] After constructing the three-dimensional attitude vector field, the system proceeds to step S2, which focuses on lithological matching and the accurate construction of the local coordinate system. First, step S21 is executed, where the system locates the "target unit" to be estimated from the unit set and reads the geological attribute code assigned to that target unit during the geological model building stage. Let the current target unit be... The geological attribute code is "D3z" (representing the Upper Devonian Zaige Formation dolomite). The system then establishes a database connection and accesses the entire area's borehole database, which contains the analysis records of core samples. The system executes a query, setting the filter condition to "Lithology_Code=='D3z'" to ensure lithological homogeneity. After passing the initial lithological screening, the system extracts the center coordinates of the target unit. It also introduces a preset search radius parameter. This parameter is determined based on the type of mineral deposit exploration, setting the horizontal search radius. meters, vertical search radius Meters. The system calculates each lithological sample. and Spatial distance, if it satisfies and and The judgment logic then includes the sample in the "valid sample list". Table 1 shows some of the valid sample data after screening.

[0028] Table 1. List of valid samples - data snippet ; As shown in Table 1, the system successfully screened out sample sets that were spatially adjacent and had consistent lithology. The system then proceeded to step S22 and the specific construction process within S22. The system read the target unit... The unit normal vector obtained in S1 With unit tendency vector To construct a strictly orthogonal local coordinate system, the system performs the following orthogonalization steps: First, the unit normal vector... Defined as the W-axis (thickness variation direction axis) of the local coordinate system, i.e. Secondly, using the W-axis and the horizontal dip vector Calculate the orthogonal orientation direction. Because Located on the horizontal plane, the system calculates As a local V-axis (along the direction axis), and normalized. Calculate Axis components: ; ; .

[0029] Normalization : .

[0030] After unitization Finally, the cross product operation is performed. Calculate the true direction along the dip axis U to ensure the orthogonality of the coordinate system. Calculate Axis components: ; ; .

[0031] Thus, the system constructs an orthogonal rotation matrix. and the center coordinates of the target unit Set as the local coordinate origin.

[0032] The aforementioned local coordinate system refers to a Cartesian coordinate system constructed based on the natural occurrence characteristics of ore body units. One axis is perpendicular to the stratum bedding plane, and the other two axes are parallel to the dip and strike lines of the strata, respectively. It is used to convert Euclidean distance in global space into relative distances that reflect geological anisotropy.

[0033] Please see Figure 1 and Figure 4 S3: Calculate the relative position vector of the effective sample list relative to the target cell, map the relative position vector to the local coordinate system to analyze the distance components in each direction, and use the range parameter to normalize the distance components in each direction to synthesize the constructed distance. The specific steps for S3 are as follows: S31: Obtain the spatial coordinates of each sample in the list of valid samples and the center coordinates of the target unit respectively. Obtain the global relative position vector of each sample relative to the target unit through vector subtraction. Project the global relative position vector onto the U-axis, V-axis and W-axis of the local coordinate system respectively. Analyze the distance components along the dip, along the directional and along the normal directions. S32: Obtain the preset variation function range ellipsoid parameters, extract the main range, secondary range and vertical range as normalization factors respectively, and use the normalization factors to perform dimensionless processing on the distance components in each direction according to the anisotropic distance synthesis formula and calculate the Euclidean norm to generate the construction distance. The process of generating the construction distance based on the anisotropic distance synthesis formula in S32 is specifically performed according to the following formula: ; in, Represents the construction distance. Represents the distance component along the dip. Represents the distance component along the direction. Represents the distance component along the normal direction. Represents the main variable range parameter. Represents the secondary range parameters. Represents the vertical range parameter; S33: Store the calculated constructed distances one-to-one with the samples in the effective sample list, and establish a distance weight calculation dataset containing the mapping relationship between sample identifiers and constructed distances, so as to provide a constructed corrected distance parameter basis for subsequent spatial interpolation using the spherical variogram model.

[0034] Step S3 primarily performs anisotropic mapping and structural correction of spatial distances. The system first executes step S31, iterating through the list of valid samples in Table 1. Taking sample S_101 as an example, its global coordinates are... The center coordinates of the target unit are The system first performs vector subtraction to calculate the global relative position vector. Next, the system uses the local coordinate axis vector constructed in S2 to... Project onto the local axis. Calculate the distance component along the dip direction. .

[0035] Calculate the distance component along the direction .

[0036] Calculate the normal distance component .

[0037] The system takes the absolute value of the above results to obtain the transformed anisotropic components: m, m, m.

[0038] Then, in step S32, the system reads the preset range parameters of the variogram function. These parameters are obtained by fitting the semivariogram function of the exploration line profile and experimental samples in the mining area in the early stage, and can quantitatively characterize the range of mineralization continuity. In this embodiment, the main range parameter is set as follows: Meters (with the best mineralization continuity along the dip), secondary range parameters Meters (secondary along the direction), vertical range parameter Meters (where the direction of cross-layer change is fastest). The system generates the structural distance based on the anisotropic distance synthesis formula referenced in S32, as follows:

[0039] in, Representing tectonic distance, it is a dimensionless scalar value used to characterize relative spatial correlation in a geostatistical sense; Represents the distance component along the dip; Represents the distance component along the direction; Represents the distance component along the normal direction; Represents the main range parameter, used to normalize the tendency distance; This represents the secondary range parameter, used to normalize the travel distance; This represents the vertical range parameter, used to normalize the normal distance.

[0040] Substituting the calculated components and range parameters of sample S_101 into the formula, a practical calculation example is performed: First term squared: ; The square of the second term: ; The third term squared: ; Summation and square root: .

[0041] The calculation result This indicates that although sample S_101 is only about 47 meters away from the target unit in physical space, its relative distance after structural correction is significantly increased because its component in the normal direction (thickness direction) exceeds the vertical range, which means that its geological relevance is greatly reduced.

[0042] After the calculation is completed, step S33 is executed. The system creates a "distance weight calculation dataset" object in memory, stored using a key-value pair structure. The system executes the above steps S31 to S32 on all samples in Table 1. As shown in Table 2, the anisotropic component analysis and final tectonic distance results for some samples are presented. These data not only achieve the scalarization of spatial distances but also embed the geological law of stratabound deposits: "easy to extend along bedding planes, difficult to extend across bedding planes."

[0043] Table 2 Calculation results of sample construction distance ; Please see Figure 1 and Figure 5 S4: Substitute the structural distance into the spherical variogram model to calculate the semivariogram value, solve the weighting coefficient based on the semivariogram value, use the weighting coefficient to calculate the estimated grade value by weighting the grade values ​​of the effective sample list, combine the average rock density value, the estimated grade value and the target unit volume to calculate the metal resource quality and accumulate to output the total reserve value.

[0044] The specific steps for S4 are as follows: S41: Based on the construction distance and the preset nugget value and sill value parameters, call the spherical variogram model to calculate the semivariance values ​​between each sample and the target unit and between each sample, and construct the semivariance matrix and right-hand covariance vector required for the Kriging linear equation system. The process of S41 calling the spherical variogram model to calculate the semivariogram is performed according to the following formula: ; in, Represents the semivariance value. Represents the construction distance. Represents the nugget constant. Represents the height of the arch. Represents the normalized range threshold; S42: Introduce the Lagrange multiplier factor to satisfy the unbiased weight constraint, perform inverse matrix operation or Gaussian elimination on the Kriging equation system containing the semivariance matrix to obtain the weight coefficient corresponding to each sample in the effective sample list, and ensure that the algebraic sum of all weight coefficients is 1. S43: Use weighted coefficients to linearly weight and sum the test grades of the corresponding core samples in the effective sample list to obtain the estimated grade value of the target unit. Multiply the estimated grade value by the volume of the target unit and the average density of the rock to obtain the metal resource quality of the unit. Iterate through all units and sum them up to obtain the total reserves of the deposit. The process of S43 in calculating and estimating grade values ​​by weighting grade values ​​of the valid sample list using weighting coefficients specifically includes: Extract the lead and zinc grade test values ​​of each sample in the list of valid samples, and multiply the test value by the weight coefficient of the corresponding sample obtained in S42 to obtain the weighted component. The weighted components of all samples are summed. If the summation result is negative, it is truncated and corrected using a preset lower limit value to generate the estimated grade value of the target unit. Determine whether the estimated grade value is greater than the preset boundary grade index. If it is, mark the unit as an industrial ore body unit and include it in the subsequent metal resource quality calculation. Otherwise, mark the unit as a surrounding rock unit and exclude it from the reserve calculation range.

[0045] Step S4 is the core step in resource estimation, involving geostatistical interpolation and reserve accumulation. First, step S41 is executed, where the system uses the structural distance calculated in step S3. The semivariogram is calculated by combining the preset variogram model parameters. The preset parameters are derived from the fitting of the semivariogram of the mining area experiment: nugget value. Arch height Normalized range threshold Calling the formula for the spherical variogram model: ; in, Represents the semivariance value, used to measure the degree of variation between pairs of points in space; Represents the construction distance; This represents the nugget constant, reflecting random variability over extremely short distances; This represents the arch height, which is the difference between the total slab value and the nugget value. The normalized range threshold represents the distance limit at which the variogram reaches the sill value.

[0046] A calculation example is performed based on the sample data in Table 2: For sample S_101, the construction distance is... .because The range has been exceeded, so a value is assigned directly. .

[0047] This result indicates that S_101 has extremely low geostatistical correlation with the target unit, exhibiting only random variance. For sample S_102, the tectonic distance... .because Substitute into the formula in the first half to calculate: .

[0048] For sample S_103, the construction distance Calculated .

[0049] Proceeding to step S42, the system constructs a system of Kriging linear equations. The system uses the calculated semivariance values ​​of each sample and the target cell to form the right-hand vector. The semivariance values ​​of each sample are then used to construct the left-hand matrix. The system then adds Lagrange multiplier constraints to the rows and columns of the matrix. The Gaussian elimination method from the numerical analysis library is used to solve the system of equations. The weighting coefficients corresponding to each sample are obtained through the solution. In this embodiment, the weighting coefficients are as follows: , , After verification, the sum of the weight coefficients is 1.0, which satisfies the unbiased constraint condition.

[0050] The aforementioned spherical variogram model is a theoretical model widely used in geostatistics. Its characteristic is that the variogram value increases monotonically with increasing distance until it reaches the range and then remains at the sill value level. It can accurately describe the spatial continuity and random structure of geological variables.

[0051] Perform step S43 to estimate the grade using weighting coefficients. Extract the grade from Table 1: , , .

[0052] Calculate the weighted components: ; ; .

[0053] The estimated grade value is obtained by summing. .

[0054] The system introduces a preset lower grade limit (e.g., 0%) for truncation correction to ensure the result is non-negative. Subsequently, a boundary grade indicator (e.g., 0.7%) is introduced. [The system then makes a judgment.] The target unit was determined to be an industrial ore body unit. Finally, the metal resource quality of this unit was calculated. The target unit volume was also determined. .

[0055] Read the preset average rock density value Formula for calculating metal content: ton.

[0056] The system accumulates this value into the total reserves register of the deposit. Table 3 shows the final estimation results. By traversing all units in the entire area, the system finally outputs the total ore quantity and metal quantity data of the entire lead-zinc deposit, completing the entire process from geometric modeling to geostatistical estimation.

[0057] Table 3. Target Unit Valuation and Resource Quantity Calculation Results ; A 3D model-based system for assessing the resource reserves of stratabound lead-zinc deposits is provided. This system is used to execute the aforementioned 3D model-based method for assessing the resource reserves of stratabound lead-zinc deposits. The system includes: The vector field generation module is used to call the stratigraphic triangular mesh data to calculate the normal vector and dip vector, construct the cell set, match the associated patches for the cells in the cell set, and assign the normal vector and dip vector of the associated patches to the cells to generate a three-dimensional attitude vector field. The coordinate system construction module is used to extract target units from the unit set, filter core sample data that match the lithology of the target units to construct a list of valid samples, and use normal vectors and dip vectors to construct a local coordinate system centered on the target units. The construction distance analysis module is used to calculate the relative position vector of the effective sample list relative to the target cell, map the relative position vector to the local coordinate system to analyze the distance components in each direction, and use the range parameter to normalize the distance components in each direction to synthesize the construction distance. The reserve estimation module is used to substitute the structural distance into the spherical variogram model to calculate the semivariogram value, solve the weighting coefficient based on the semivariogram value, use the weighting coefficient to calculate the estimated grade value by weighting the grade values ​​of the effective sample list, and combine the average rock density value, the estimated grade value and the target unit volume to calculate the metal resource quality and accumulate and output the total reserve value.

[0058] The above embodiments illustrate preferred embodiments of the present invention. Any equivalent adjustments to the technical solution based on software engineering methods are within the scope of protection, including but not limited to: implementing algorithm logic using different programming languages, refactoring functional modules into services, adjusting data interaction protocols, and optimizing resource scheduling strategies. Any implementation scheme derived from reasonable modifications to the data processing flow, service call chain, or system architecture layer without departing from the core technology of the present invention should be considered within the scope of protection of the present invention.

Claims

1. A method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model, characterized in that, Includes the following steps: S1: Calculate the normal vector and dip vector by calling the stratigraphic triangular mesh data, construct a set of elements, match associated patches for the elements in the set of elements, and assign the normal vector and dip vector of the associated patches to the elements to generate a three-dimensional attitude vector field; S2: Extract target units from the unit set, filter core sample data that match the lithology of the target units to construct a valid sample list, and use the normal vector and the dip vector to construct a local coordinate system centered on the target unit; S3: Calculate the relative position vector of the effective sample list relative to the target unit, map the relative position vector to the local coordinate system to analyze the distance components in each direction, and use the range parameter to normalize the distance components in each direction to synthesize the constructed distance; S4: Substitute the structural distance into the spherical variogram model to calculate the semivariogram value, solve the weighting coefficient based on the semivariogram value, use the weighting coefficient to calculate the estimated grade value by weighting the grade values ​​of the effective sample list, combine the average rock density value, the estimated grade value and the target unit volume to calculate the metal resource quality and accumulate and output the total reserve value.

2. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 1, characterized in that, The specific steps of S1 are as follows: S11: Traverse all triangular patches contained in the formation triangular mesh data, extract the spatial coordinate data of the vertices of the triangular patches, use the cross product algorithm to calculate the unit normal vector perpendicular to the triangular patch, and project the unit normal vector onto the horizontal plane to solve the unit dip vector indicating the dip direction of the formation. S12: Discretize the three-dimensional space based on the boundary range of the ore body to generate the unit set composed of regular hexahedral voxels, calculate the three-dimensional coordinates of the geometric center point of each voxel in the unit set, and establish a KD tree data structure for accelerating spatial indexing. S13: For each voxel in the unit set, use the KD tree data structure to search for the triangular facet with the closest Euclidean distance in the stratigraphic triangular mesh data as the associated facet, extract the unit normal vector and the unit dip vector of the associated facet and assign them to the corresponding voxel, and traverse the entire region to generate the three-dimensional attitude vector field.

3. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 1, characterized in that, The specific steps of S2 are as follows: S21: Obtain the geological attribute code of the target unit, traverse the core sample records in the borehole database of the whole area, compare the lithology code of the core sample with the geological attribute code of the target unit, retain samples with consistent codes and spatial distance within the preset search radius, and establish the list of effective samples. S22: Obtain the normal vector and dip vector corresponding to the target unit in the three-dimensional attitude vector field, define the dip vector as the local U-axis, define the normal vector as the local W-axis, and generate the local V-axis by using the cross product operation of the normal vector and the dip vector, thereby constructing the local coordinate system that can characterize the local anisotropy of the stratabound ore body.

4. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 1, characterized in that, The specific steps of S3 are as follows: S31: Obtain the spatial coordinates of each sample in the list of valid samples and the center coordinates of the target unit respectively. Obtain the global relative position vector of each sample relative to the target unit through vector subtraction. Project the global relative position vector onto the U-axis, V-axis and W-axis of the local coordinate system respectively. Analyze the distance components along the dip, along the direction and along the normal. S32: Obtain the preset variation function range ellipsoid parameters, extract the main range, secondary range and vertical range as normalization factors respectively, and use the normalization factors to perform dimensionless processing on the distance components in each direction according to the anisotropic distance synthesis formula and calculate the Euclidean norm to generate the constructed distance. S33: Store the calculated constructed distances one-to-one with the samples in the effective sample list, and establish a distance weight calculation dataset containing the mapping relationship between sample identifiers and constructed distances, so as to provide a constructed corrected distance parameter basis for subsequent spatial interpolation using the spherical variogram model.

5. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 1, characterized in that, The specific steps of S4 are as follows: S41: Based on the constructed distance and the preset nugget value and sill value parameters, call the spherical variogram model to calculate the semivariance values ​​between each sample and the target unit and between each sample, and construct the semivariance matrix and right-hand covariance vector required for the Kriging linear equation system. S42: Introduce a Lagrange multiplier factor to satisfy the unbiased weight constraint, perform inverse matrix operation or Gaussian elimination on the Kriging equation system containing the semivariance matrix to obtain the weight coefficient corresponding to each sample in the effective sample list, and ensure that the algebraic sum of all the weight coefficients is 1. S43: Using the weighting coefficients, perform a linear weighted summation of the test grades of the corresponding core samples in the list of valid samples to obtain the estimated grade value of the target unit. Multiply the estimated grade value by the volume of the target unit and the average density of the rock to obtain the metal resource quality of the unit. Iterate through all units and sum them up to obtain the total reserves of the deposit.

6. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 3, characterized in that, The process of constructing the local coordinate system capable of characterizing the local anisotropic features of the stratabound ore body, as described in S22, specifically includes: Obtain the normal vector and the yaw vector of the target unit, and normalize the two vectors to eliminate the difference in magnitude; The normalized directional vector is set as the primary extension direction axis, and the normalized normal vector is set as the thickness variation direction axis. By calculating the vector product of the primary extension direction axis and the thickness variation direction axis, a secondary extension direction axis orthogonal to the first two axes is generated. An orthogonal rotation matrix is ​​constructed using the primary extension direction axis, the secondary extension direction axis, and the thickness variation direction axis, and the center coordinates of the target unit are set as the origin to establish the local coordinate system.

7. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 4, characterized in that, The process of generating the constructed distance according to the anisotropic distance synthesis formula described in S32 is specifically performed according to the following formula: ； in, Represents the constructed distance, Represents the distance component along the dip. Represents the distance component along the direction. Represents the distance component along the normal direction. Represents the main variable range parameter. Represents the secondary range parameters. This represents the vertical range parameter.

8. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 5, characterized in that, The process of calculating the semivariance value by calling the spherical variogram model as described in S41 is specifically performed according to the following formula: ； in, This represents the semivariance value. Represents the constructed distance, Represents the nugget constant. Represents the height of the arch. This represents the normalized range threshold.

9. The method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model according to claim 5, characterized in that, The process of calculating and estimating the grade value by weighting the grade values ​​of the valid sample list using the weighting coefficients as described in S43 specifically includes: Extract the lead and zinc grade test values ​​of each sample in the list of valid samples, and multiply the test value by the weight coefficient of the corresponding sample obtained in S42 to obtain the weighted component. The weighted components of all samples are summed. If the summation result is negative, it is truncated and corrected using a preset lower limit value to generate the estimated grade value of the target unit. Determine whether the estimated grade value is greater than the preset boundary grade index. If it is greater, mark the unit as an industrial ore body unit and participate in the subsequent metal resource quality calculation. Otherwise, mark the unit as a surrounding rock unit and exclude it from the reserve calculation range.

10. A resource reserve assessment system for stratabound lead-zinc deposits based on a three-dimensional model, characterized in that, The system is used to implement the method for assessing the resource reserves of stratabound lead-zinc deposits based on a three-dimensional model as described in any one of claims 1-9, and the system comprises: The vector field generation module is used to call the stratigraphic triangular mesh data to calculate the normal vector and dip vector, construct a set of cells, match associated patches for the cells in the set of cells, and assign the normal vector and dip vector of the associated patches to the cells to generate the three-dimensional attitude vector field. The coordinate system construction module is used to extract target units from the unit set, filter core sample data that match the lithology of the target units to construct the effective sample list, and construct a local coordinate system centered on the target unit using the normal vector and the dip vector. The construction distance analysis module is used to calculate the relative position vector of the effective sample list relative to the target unit, map the relative position vector to the local coordinate system to analyze the distance components in each direction, and use the range parameter to normalize the distance components in each direction to synthesize the construction distance. The reserve estimation module is used to substitute the structural distance into the spherical variogram model to calculate the semivariogram value, solve the weighting coefficient based on the semivariogram value, use the weighting coefficient to weight the grade values ​​of the effective sample list to calculate the estimated grade value, combine the average rock density value, the estimated grade value and the target unit volume to calculate the metal resource quality and accumulate and output the total reserve value.

Images