A geological statistics-based deep gold mine resource quantity estimation method

CN122883252APending Publication Date: 2026-10-09SHANDONG GEOLOGICAL EXPLORATION INST OF SINOCHEM GEOLOGY & MINING ADMINISTRATION
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202611382170.X
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2026-09-08
Publication Date
2026-10-09

AI Technical Summary

Technical Problem

然而,受地质构造控制,矿化分布在地下三维空间内常表现出明显的各向异性和空间相关性,这种空间结构的变异在均质单元假设下无法被有效表征和利用,从而造成了估算结果与真实空间分布规律之间的偏差

Benefits of technology

1.本方法将变异函数模型沿主要控矿构造的走向方向和倾向方向分别构建,使成矿元素品位在三维空间中不同方向上的变异结构得以独立表征,避免了均质单元假设下对方向性差异的忽略,克里金插值依据各向异性的变异函数模型分配权重系数,插值计算能够与矿化的实际空间相关性保持一致,品位估计值在走向方向和倾向方向上均能反映成矿元素分布的延展特征,插值结果与矿化空间分布规律之间的系统性偏差得到抑制,同时插值不确定度的计算同时纳入留一法交叉验证残差与克里金插值方差两重信息,交叉验证残差提供了各局部区域内插值估计与实际值偏离程度的直接度量,克里金插值方差提供了基于数据空间构形的估计稳定性评价,两重信息融合后形成的插值不确定度对品位估计值可信度的表征更为完整,单一指标可能遗漏的低置信度区域能够被有效识别,降低了不可靠估计值参与后续圈定所带来的不确定性传递。

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN122883252A_ABST
    Figure CN122883252A_ABST
Patent Text Reader

Abstract

The application discloses a deep gold mine resource quantity estimation method based on geostatistics, and particularly relates to the technical field of statistical interpolation processing of geological data, and is used to solve the problem that the existing data processing method does not consider the anisotropic variation structure of mineralization properties in three-dimensional space controlled by geological structure, and thus systematic deviation exists between the grade interpolation result and the actual spatial distribution characteristics of mineralization; the three-dimensional geological model data of a target mining area is acquired, a variation function model is constructed along the strike direction and the tendency direction of the main ore-controlling structure, the grade estimation value is obtained by performing Kriging interpolation on the cubic unit by using the variation function model, the interpolation uncertainty is calculated based on the leave-one-out cross-validation residual and the Kriging interpolation variance, the spatial distribution range of the ore body is circled by using the grade estimation value and the interpolation uncertainty as double constraints, and the resource quantity is calculated according to the circled range.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of statistical interpolation processing technology for geological data, and more specifically, to a method for estimating deep gold resources based on geological statistics. Background Technology

[0002] In quantitative mineral resource assessment, dividing geological bodies into three-dimensional spatial units and estimating resource quantity based on unit attributes is a standard data processing procedure. Current technologies often rely on discretizing the underground space into independent homogeneous cubic units, assigning geological attribute values ​​to each unit, and then using statistical methods such as information content calculation or volume weighting to obtain the estimation results. However, controlled by geological structures, mineralization distribution in the underground three-dimensional space often exhibits significant anisotropy and spatial correlation. This variation in spatial structure cannot be effectively characterized and utilized under the assumption of homogeneous units, resulting in a deviation between the estimated results and the actual spatial distribution patterns. The aforementioned data processing procedure fails to incorporate the spatial variation structure between units into the interpolation estimation framework, leading to systematic errors in the inference of grade attribute distribution.

[0003] Existing data processing methods divide the three-dimensional space into homogeneous cubic units and assign and statistically analyze them independently, without considering the differences in the degree of variation of mineralization attributes along different directions of the ore-controlling structure. In other words, they ignore the anisotropic variation structure of mineralization attributes in three-dimensional space controlled by geological structures. Since the strength of variation in direction cannot be distinguished, the weight allocation cannot be adjusted according to the actual spatial correlation during interpolation calculation, resulting in a systematic deviation between the grade interpolation results and the actual spatial distribution characteristics of mineralization. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a method for estimating deep gold resources based on geostatistics to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: A geostatistical method for estimating deep gold resources includes the following steps: S1: Obtain the three-dimensional geological model data of the target mining area; S2: Based on three-dimensional geological model data, search and construct variation function models that characterize the spatial variation structure of ore-forming element grades along the strike and dip directions of the main ore-controlling structures. S3: Using the variogram model, kriging interpolation is performed on the ore-forming element grades of the cubic units divided in the three-dimensional space of the target mining area to obtain grade estimates. S4: Calculate the interpolation uncertainty of each cube element based on the leave-one-out cross-validation residuals of known sample points and the variance of Kriging interpolation; S5: Delineate the spatial distribution range of the ore body based on the grade estimate and interpolation uncertainty; S6: Calculate the resource quantity of deep ore bodies based on the delineated spatial distribution range of the ore bodies.

[0006] Furthermore, obtain three-dimensional geological model data of the target mining area, including: Borehole logging data, exploration line profile data, and structural geological maps were extracted from geological data of the study area. Based on borehole logging data and exploration line profile data, determine the spatial location of stratigraphic boundaries, rock mass boundaries, and known ore body outlines; Extract the strike, dip, and occurrence information of the main ore-controlling fault structures from structural geological maps; Spatial correlation is formed by linking the spatial locations of stratigraphic boundaries, rock mass boundaries, and known ore body outlines with the strike, dip, and occurrence information of major ore-controlling fault structures to create three-dimensional geological model data.

[0007] Furthermore, based on three-dimensional geological model data, variogram models characterizing the spatial variation structure of ore-forming element grades are searched and constructed along the strike and dip directions of the main ore-controlling structures, including: Extracting the ore-forming element grade values ​​of known ore bodies from three-dimensional geological model data; Based on the strike and occurrence information of the main ore-controlling fault structures, the data point pairs corresponding to the grade values ​​are searched along the strike direction with lag distance, the experimental variability function of the strike direction is calculated and fitted into the theoretical variability function of the strike direction. Based on the dip and occurrence information of the main ore-controlling fault structures, the data point pairs corresponding to the grade values ​​are searched along the dip direction with lag distance, the experimental variability function of the dip direction is calculated and fitted into the theoretical variability function of the dip direction. The variogram of the direction theory and the variogram of the tendency theory are combined to form a variogram model.

[0008] Furthermore, fitting the theoretical variogram of the directional direction includes: performing least-squares fitting on the discrete points of the experimental variogram of the directional direction using a spherical model to obtain the theoretical variogram of the directional direction; fitting the theoretical variogram of the tendency direction includes: performing least-squares fitting on the discrete points of the experimental variogram of the tendency direction using a spherical model to obtain the theoretical variogram of the tendency direction.

[0009] Furthermore, using a variogram model, kriging interpolation is performed on the ore-forming element grades of the cubic units divided in the three-dimensional space of the target mining area to obtain grade estimates, including: Extract the ranges of the directional theory variogram and the tendency theory variogram from the variogram model; The ratio of the range of the variogram of the direction theory to the range of the variogram of the tendency theory is used as the anisotropy ratio. The side lengths of the initial cube elements are scaled according to the anisotropy ratio so that the ratio of the side length in the orientation direction to the side length in the dip direction is equal to the anisotropy ratio. The three-dimensional space of the target mining area is divided using scaled cubic units; For each cubic cell, the grade of ore-forming elements is estimated by kriging interpolation using a variogram model.

[0010] Furthermore, based on the leave-one-out cross-validation residuals of known sample points and the variance of Kriging interpolation, the interpolation uncertainty of each cubic element is calculated, including: Leave-one-out cross-validation is performed on the known sample points, and the known sample points are removed one by one. Kriging interpolation is performed on the removal positions using the variogram model to obtain the interpolated estimates of the removal positions. The difference between the interpolated estimate at the removal location and the actual grade value of the removed known sample point is calculated as the residual. Centered on each cubic unit, the residuals falling into the search ellipsoid determined by the variogram model are collected, and the root mean square value of the falling residuals is calculated as a local residual index. Extracting the Kriging interpolation variance generated during the Kriging interpolation process; The interpolation uncertainty is obtained by multiplying the local residual index by the Kriging interpolation variance.

[0011] Furthermore, the search ellipsoid determined by the variogram model satisfies the following: the major axis radius of the search ellipsoid takes the range of the theoretical variogram in the strike direction; the minor axis radius of the search ellipsoid takes the range of the theoretical variogram in the dip direction; and the major axis azimuth of the search ellipsoid is consistent with the strike and occurrence information of the main ore-controlling fault structures.

[0012] Furthermore, based on the grade estimate and interpolation uncertainty, the spatial distribution range of the ore body is delineated, including: The grade estimate is compared with the set grade threshold, and cubic cells with a grade estimate greater than the set grade threshold are selected. The interpolation uncertainty is compared with a set uncertainty threshold, and cubic cells with interpolation uncertainty less than the set uncertainty threshold are selected. The spatial locations of cubic units that simultaneously satisfy the condition that the grade estimate is greater than the set grade threshold and the interpolation uncertainty is less than the set uncertainty threshold are merged to form the spatial distribution range of the ore body.

[0013] Furthermore, the uncertainty threshold is determined based on the interpolation uncertainty corresponding to the known ore body: the interpolation uncertainty corresponding to each cubic unit within the spatial distribution of the known ore body is extracted; the maximum value of the extracted interpolation uncertainty is used as the uncertainty threshold.

[0014] Furthermore, the resource quantity of the deep ore body is calculated based on the delineated spatial distribution range of the ore body, including: Grade estimates for each cubic unit within the spatial distribution range of the ore body are extracted; Get the volume of each cube cell; The resource quantity of the deep ore body is obtained by summing the product of the grade estimates of all cubic units within the spatial distribution range of the ore body and the unit volume.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. This method constructs variogram models along the strike and dip directions of the main ore-controlling structures, enabling independent characterization of the variation structure of ore-forming element grades in different directions in three-dimensional space. This avoids neglecting directional differences under the homogeneous unit assumption. Kriging interpolation assigns weight coefficients based on the anisotropic variogram model, ensuring that the interpolation calculation remains consistent with the actual spatial correlation of mineralization. The grade estimates reflect the extension characteristics of ore-forming element distribution in both the strike and dip directions, suppressing systematic deviations between the interpolation results and the spatial distribution patterns of mineralization. Furthermore, the calculation of interpolation uncertainty incorporates both leave-one-out cross-validation residuals and Kriging interpolation variance. The cross-validation residuals provide a direct measure of the deviation between the interpolation estimate and the actual value in each local region, while the Kriging interpolation variance provides an evaluation of the estimation stability based on the spatial configuration of the data. The interpolation uncertainty formed by fusing these two pieces of information provides a more complete characterization of the reliability of the grade estimate. Low-confidence areas that might be missed by a single indicator can be effectively identified, reducing the uncertainty transmission caused by unreliable estimates participating in subsequent delineation.

[0016] 2. The delineation of the spatial distribution range of ore bodies is subject to dual constraints: grade estimates and interpolation uncertainty. Grade constraints ensure that the delineated range meets the content requirements of ore-forming elements, while interpolation uncertainty constraints ensure that areas with insufficient estimation reliability are excluded. These two types of constraints screen cubic units from two dimensions: mineralization intensity and estimation reliability. The basis for determining the ore body boundary has been expanded from a single grade threshold to a joint judgment of grade and reliability. The distribution location of the delineation results in three-dimensional space is consistent with the interpolation uncertainty level of known ore bodies. The set of cubic units on which resource quantity calculation is based eliminates parts with poor estimation quality, thereby improving the reliability and robustness of the resource quantity calculation results. Attached Figure Description

[0017] Figure 1 This is a flowchart of a method for estimating deep gold resources based on geostatistics according to the present invention. Detailed Implementation

[0018] 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 of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0019] Example: Figure 1 This invention presents a method for estimating deep gold resources based on geostatistics, which includes the following steps: S1: Obtain the three-dimensional geological model data of the target mining area; S2: Based on three-dimensional geological model data, search and construct variation function models that characterize the spatial variation structure of ore-forming element grades along the strike and dip directions of the main ore-controlling structures. S3: Using the variogram model, kriging interpolation is performed on the ore-forming element grades of the cubic units divided in the three-dimensional space of the target mining area to obtain grade estimates. S4: Calculate the interpolation uncertainty of each cube element based on the leave-one-out cross-validation residuals of known sample points and the variance of Kriging interpolation; S5: Delineate the spatial distribution range of the ore body based on the grade estimate and interpolation uncertainty; S6: Calculate the resource quantity of deep ore bodies based on the delineated spatial distribution range of the ore bodies.

[0020] In the specific implementation of S1, borehole logging data, exploration line profile data, and structural geological maps were extracted from the geological data of the study area. The geological data of the study area consisted of digitized data accumulated during previous geological exploration and mining production, including borehole logging results from previous explorations, exploration line profile measurement results, and structural geological mapping results of the mining area and surrounding areas. Borehole logging data was obtained by reading the borehole information tables stored in the geological database of the study area. Each borehole logging data included the three-dimensional coordinates of the borehole opening, borehole trajectory inclination data, borehole opening azimuth and dip angle, and a lithological stratification description record along the borehole axis. The lithological stratification description record recorded the start and end depths and lithological names of each lithology layer. For example, the lithology from a depth of 250 meters to 310 meters was named sericite-quartzized granitic breccia. Exploration line profile data were obtained by reading the exploration line profile information table stored in the geological database of the study area. Each exploration line profile includes the three-dimensional coordinates of the profile endpoints, the profile azimuth, the geometric shape of each geological boundary on the profile, and the corresponding depth value. The structural geological maps are geological maps of the study area that have been digitized after field measurements and indoor interpretation. They include the surface outcrops of the main ore-controlling fault structures and the fault attitude markings measured at each measuring point.

[0021] Based on borehole logging data and exploration line profile data, the spatial locations of stratigraphic boundaries, rock mass boundaries, and known orebody outlines are determined. For each borehole logging data, lithological names in the lithological stratification description record are read sequentially along the borehole trajectory according to depth. Depth points where lithological names change are marked as stratigraphic boundary points or rock mass boundary points. Simultaneously, depth intervals identified in the lithological names as mineralized altered rocks or industrial ore bodies are marked as known orebody outline points. The three-dimensional coordinates corresponding to each stratigraphic boundary point, rock mass boundary point, and known orebody outline point are recorded. The three-dimensional coordinates are calculated by combining the three-dimensional coordinates of the borehole opening with borehole trajectory inclination data and the corresponding depth. In the exploration line profile data, the geometric morphology of each geological boundary on the profile is read. The depth and location information of the same geological interface on different exploration line profiles are correlated and connected to determine the distribution morphology of each geological interface in three-dimensional space, including the spatial locations of stratigraphic boundaries, rock mass boundaries, and known orebody outlines.

[0022] The strike and dip orientation information of major ore-controlling faults were extracted from the structural geological map. The map marked the surface outcrops of the major ore-controlling faults and the measured fault attitudes at each measuring point. Fault attitudes were expressed using three sets of angle values: strike azimuth, dip azimuth, and dip angle. The strike azimuth and dip azimuth data were read point by point along the surface outcrops of the major ore-controlling faults. The strike azimuth represents the angle of extension of the major ore-controlling fault relative to true north on the horizontal plane; the dip azimuth represents the pointing angle of the dip direction of the major ore-controlling fault relative to true north; and the dip angle represents the angle between the dip surface of the major ore-controlling fault and the horizontal plane. The strike azimuth data read from each measuring point along the surface outcrop line are statistically averaged. The statistical averaging method is to normalize the strike azimuth of each measuring point to the range of 0 degrees to 180 degrees and then take the arithmetic mean. This yields the strike and attitude information representing the overall distribution characteristics of the main ore-controlling fault structures. The dip azimuth data read from each measuring point along the surface outcrop line are also statistically averaged. The statistical averaging method is to normalize the dip azimuth of each measuring point to the range of 0 degrees to 360 degrees and then take the arithmetic mean. This yields the dip and attitude information representing the overall distribution characteristics of the main ore-controlling fault structures.

[0023] The spatial relationships between stratigraphic boundaries, rock mass boundaries, and known ore body outlines with the strike and dip information of major ore-controlling fault structures are spatially correlated to form three-dimensional geological model data. This spatial correlation process is performed within a three-dimensional coordinate framework. The spatial locations of stratigraphic boundaries, rock mass boundaries, and known ore body outlines are used as geometric elements of the geological bodies, while the strike and dip information of major ore-controlling fault structures is used as structural constraint elements. These geometric elements and structural constraint elements are registered in the same three-dimensional space, ensuring geometric consistency between the spatial distribution of stratigraphic boundaries and rock mass boundaries and the attitude of major ore-controlling fault structures. The three-dimensional geological model data includes the three-dimensional spatial distribution of stratigraphic boundaries, rock mass boundaries, and known ore body outlines, as well as the strike and dip information of major ore-controlling fault structures. All information is organized in three-dimensional coordinate data format, collectively constituting the three-dimensional geological model data of the target mining area.

[0024] In the specific implementation of S2, the ore-forming element grade values ​​of known ore bodies are extracted from the 3D geological model data. The spatial distribution of the known ore body outline in the 3D geological model data defines the 3D spatial range of the known ore body. Within the spatial distribution range of the known ore body outline, the corresponding lithological stratification description records from the borehole logging data are extracted. The sample analysis data, obtained from core samples collected at fixed sampling intervals within the known ore body outline during borehole drilling, are the ore-forming element grade values ​​obtained through laboratory analysis. The ore-forming element grade value represents the content of ore-forming elements per unit weight of ore; the ore-forming element is gold, and the unit of ore-forming element grade value is grams per ton. The ore-forming element grade value of each sample is correlated with the sampling position coordinates in 3D space. The sampling position coordinates are calculated from the 3D coordinates of the corresponding borehole opening, borehole trajectory inclination data, and sampling depth, forming a data point set of ore-forming element grade values.

[0025] Based on the strike and occurrence information of the main ore-controlling fault structures, data point pairs corresponding to the ore-forming element grades are searched along the strike direction using lag distances. The experimental variability function of the strike direction is calculated and fitted to the theoretical variability function of the strike direction. The strike and occurrence information of the main ore-controlling fault structures is recorded in the three-dimensional geological model data in the form of strike azimuths. The strike azimuths define the spatial directional reference for searching data point pairs corresponding to the ore-forming element grades along the strike direction. The initial value of the lag distance is determined based on the minimum distance between data points in the ore-forming element grade data point set. For example, if the minimum distance between data points in the ore-forming element grade data point set is 20 meters, then the initial value of the lag distance is set to 20 meters, and the tolerance range of the lag distance is set to half of the initial value of the lag distance, i.e., 10 meters. Along the strike direction specified by the strike azimuth, using the lag distance as the search step size, search for all data point pairs corresponding to the ore-forming element grade values ​​that satisfy the spatial distance condition within the lag distance tolerance range. For each data point pair corresponding to the searched ore-forming element grade values, calculate the square of the difference between the two ore-forming element grade values ​​in each data point pair. The specific expression for calculating the variogram value of the strike direction experimental variogram at the lag distance can be expressed as follows: γ(h)=1 / (2N(h))×Σ[Z(xi)-Z(xi+h)] 2Where γ(h) represents the value of the strike direction experimental variogram at a lag distance h, h represents the lag distance, N(h) represents the number of data point pairs corresponding to the searched ore-forming element grade values, Z(xi) and Z(xi+h) represent the two ore-forming element grade values ​​in the i-th data point pair, i is the data point number. The lag distance is increased, for example, sequentially increasing to 40 meters, 60 meters, and 80 meters, or continuing to increase the lag distance in increments of 20 meters. The process of searching for data point pairs corresponding to ore-forming element grade values ​​and calculating the value of the strike direction experimental variogram at the lag distance is repeated to obtain a series of discrete points of the strike direction experimental variogram at different lag distances. For the discrete points of the strike direction experimental variogram, a spherical model is used for least-squares fitting. The spherical model is defined by three parameters: sill value, range, and nugget constant. The specific expression of the spherical model is: γ(h) = C0 + C[1.5(h / a) - 0.5(h / a)]. 3 When h ≤ a; γ(h) = C0 + C, when h > a; where γ(h) represents the model value of the spherical model at the lag distance h, h represents the lag distance, C0 represents the nugget constant, C represents the sill value, and a represents the range. The least squares fitting process involves adjusting the values ​​of C0, C, and a to minimize the sum of squared deviations between the model value of the spherical model at the lag distance corresponding to the discrete point of the experimental variogram of the directional direction and the variogram value of the experimental variogram of the directional direction at the discrete point of the experimental variogram of the directional direction. After fitting, the specific values ​​of C0, C, and a are obtained, thereby determining the theoretical variogram of the directional direction.

[0026] Based on the dip and occurrence information of the main ore-controlling fault structures, data point pairs corresponding to the ore-forming element grades are searched along the dip direction using lag distances. The experimental variability function of the dip direction is calculated and fitted to the theoretical variability function of the dip direction. The dip and occurrence information of the main ore-controlling fault structures is recorded in the three-dimensional geological model data in the form of dip azimuth and dip angle. The dip azimuth and dip angle together define the spatial directional reference for searching data point pairs corresponding to the ore-forming element grades along the dip direction. The dip direction is the direction along the dip angle in the vertical plane passing through the dip azimuth. The initial value of the lag distance is determined based on the minimum distance between data points in the ore-forming element grade data point set. For example, if the minimum distance between data points in the ore-forming element grade data point set is 20 meters, then the initial value of the lag distance is set to 20 meters, and the tolerance range of the lag distance is set to half of the initial value of the lag distance, i.e., 10 meters. Along the dip direction specified by the dip azimuth and dip angle, with the lag distance as the search step size, search for all data point pairs corresponding to the ore-forming element grade values ​​that meet the spatial distance conditions within the lag distance tolerance range.

[0027] Calculate the square of the difference between the grades of the two ore-forming elements in each data point pair. The variogram value of the experimental variogram in the dip direction at the lag distance is calculated as: γ(h) = 1 / (2N(h)) × Σ[Z(xi) - Z(xi+h)] 2 Where γ(h) represents the variogram value of the dip direction experimental variogram at a lag distance h, h represents the lag distance, N(h) represents the number of data point pairs corresponding to the searched ore-forming element grade values, and Z(xi) and Z(xi+h) represent the two ore-forming element grade values ​​in the i-th data point pair, respectively. The lag distance is increased, for example, sequentially increasing it to 40 meters, 60 meters, and 80 meters, or continuing to increase the lag distance in increments of 20 meters. The process of searching for data point pairs corresponding to ore-forming element grade values ​​and calculating the variogram value of the dip direction experimental variogram at the lag distance is repeated to obtain a series of discrete points of the dip direction experimental variogram at different lag distances. For the discrete points of the dip direction experimental variogram, a spherical model is used for least-squares fitting. The spherical model is defined by three parameters: sill value, range, and nugget constant. The expression of the spherical model is: γ(h) = C0 + C[1.5(h / a) - 0.5(h / a)]. 3 When h ≤ a; γ(h) = C0 + C, when h > a; where γ(h) represents the model value of the spherical model at the lag distance h, h represents the lag distance, C0 represents the nugget constant, C represents the sill value, and a represents the range. The least squares fitting process involves adjusting the values ​​of C0, C, and a to minimize the sum of squared deviations between the model value of the spherical model at the lag distance corresponding to the discrete point of the experimental variogram in the trend direction and the variogram value at the discrete point of the experimental variogram in the trend direction. After fitting, the specific values ​​of C0, C, and a are obtained, thereby determining the theoretical variogram in the trend direction.

[0028] It should be noted that in geostatistics, experimental variogram values ​​and fitted values ​​from theoretical variogram models are represented by the same symbols. The forms of the experimental variogram calculation expressions for the strike direction and the dip direction are the same, but the data point pairs they rely on are different. The strike direction experimental variogram calculation expression uses data point pairs corresponding to the ore-forming element grade values ​​obtained by searching along the strike direction, while the dip direction experimental variogram calculation expression uses data point pairs corresponding to the ore-forming element grade values ​​obtained by searching along the dip direction. Therefore, the variogram values ​​calculated by the two expressions correspond to different spatial directions. Similarly, the expressions for the spherical model are the same in the strike and dip directions, but the discrete points of the experimental variogram used for fitting are different. The fitting of the theoretical variogram for the strike direction uses the discrete points of the experimental variogram for the strike direction, while the fitting of the theoretical variogram for the dip direction uses the discrete points of the experimental variogram for the dip direction. The fitted sill values, range, and nugget constants also correspond to their respective tectonic directions. The symbols of the same form in the expressions refer to the experimental variogram values ​​or spherical model values ​​in their respective directions in the context of the strike direction and the dip direction, respectively. Those skilled in the art can clearly distinguish their specific meanings based on the context.

[0029] The strike-direction theoretical variability function and the dip-direction theoretical variability function are combined to form a variability model. The variability model includes two sets of parameters: the strike-direction theoretical variability function, its range, and its nugget constant. The dip-direction theoretical variability function also includes these parameters. The strike-direction theoretical variability function characterizes the spatial variation structure of ore-forming element grades along the strike direction of the main ore-controlling fault structure, i.e., the correlation distance and degree of variation of ore-forming element grades along the strike direction. The dip-direction theoretical variability function characterizes the spatial variation structure of ore-forming element grades along the dip direction of the main ore-controlling fault structure, i.e., the correlation distance and degree of variation of ore-forming element grades along the dip direction. The variability model comprehensively describes the anisotropic characteristics of the spatial variation structure of ore-forming element grades in both the strike and dip directions.

[0030] In the specific implementation of S3, the ranges of the strike direction theoretical variogram and the dip direction theoretical variogram are extracted from the variogram model. The variogram model includes parameters for both the strike direction and dip direction theoretical variograms. The strike direction theoretical variogram parameters include the range of the strike direction theoretical variogram, and the dip direction theoretical variogram parameters include the range of the dip direction theoretical variogram. The range represents the maximum distance of spatial correlation between ore-forming element grade values; that is, beyond the range range, there is no longer spatial correlation between ore-forming element grade values. The ranges of the strike direction and dip direction theoretical variograms are directly obtained from the parameters of the variogram model.

[0031] The anisotropy ratio is defined as the ratio of the range of the theoretical variogram in the strike direction to the range of the theoretical variogram in the dip direction. The anisotropy ratio is calculated by dividing the range of the theoretical variogram in the strike direction by the range of the theoretical variogram in the dip direction; the resulting ratio must be greater than zero. For example, if the range of the theoretical variogram in the strike direction is 120 meters and the range of the theoretical variogram in the dip direction is 60 meters, then the anisotropy ratio is 2. The anisotropy ratio reflects the difference in the extent of the spatial variation structure of ore-forming element grades along the strike and dip directions. An anisotropy ratio greater than 1 indicates that the continuity along the strike direction is greater than the continuity along the dip direction, while an anisotropy ratio equal to 1 indicates that the continuity along the strike and dip directions is the same.

[0032] The initial cube unit's side length is scaled according to the anisotropy ratio so that the ratio of the side length in the strike direction to the side length in the dip direction equals the anisotropy ratio. The side length of the initial cube unit is determined based on the average spacing of data points in the ore-forming element grade data set. The average spacing is calculated by sorting all pairwise spatial distances between data points in the ore-forming element grade data set and taking the median of the sorted distances as the average spacing. For example, if the average spacing is 20 meters, then the side length of the initial cube unit in all three directions is set to 20 meters. During scaling, the side length in the dip direction remains unchanged from the initial cube unit's side length, while the side length in the strike direction is adjusted to the initial cube unit's side length multiplied by the anisotropy ratio, so that the ratio of the side length in the strike direction to the side length in the dip direction equals the anisotropy ratio. For example, when the initial cube unit's side length is 20 meters and the anisotropy ratio is 2, after scaling, the side length in the strike direction is 40 meters, the side length in the dip direction is 20 meters, and the ratio of the side length in the strike direction to the side length in the dip direction is 2. Once the side lengths in the strike and dip directions are determined, the side lengths in the third direction perpendicular to the plane containing the strike and dip directions remain consistent with the side lengths of the initial cubic unit. Through scaling, the geometry of the cubic unit is adapted to the principal axis direction of the spatial variation structure of ore-forming element grades.

[0033] The target mining area's three-dimensional space is divided using scaled cubic units. The division is based on the smallest bounding cuboid of the target mining area's three-dimensional space. This smallest bounding cuboid encompasses all geological geometric elements and structural constraints within the target mining area's three-dimensional space. Starting from a corner of the smallest bounding cuboid, the division proceeds along the strike direction with a step size equal to the scaled side length in that direction; along the dip direction with a step size equal to the scaled side length in that direction; and along the planes perpendicular to the strike and dip directions with a step size equal to the side length of the initial cubic unit, forming a regularly arranged array of cubic units. Each cubic unit has a unique three-dimensional spatial location identifier, represented by its index ordinal number in three directions or by the three-dimensional coordinates of its center point. For cubic units located at the boundary of the target mining area's three-dimensional space, if part of the cubic unit's volume exceeds the target mining area's three-dimensional space, the cubic unit is retained but marked as a boundary unit, and subsequent calculations only apply to the portion of the cubic unit within the target mining area's three-dimensional space.

[0034] For each cubic cell, Kriging interpolation of the ore-forming element grade is performed using a variogram model to obtain an estimated grade value. The specific process of Kriging interpolation is as follows: For each cubic cell, the center point of the cubic cell is taken as the point to be estimated. The shape, orientation, and size of the search ellipsoid are determined according to the variogram model. The major axis direction of the search ellipsoid is consistent with the strike and attitude information of the main ore-controlling fault structure, and the minor axis direction of the search ellipsoid is consistent with the dip and attitude information of the main ore-controlling fault structure. The major axis radius of the search ellipsoid is taken from the range of the theoretical variogram of the strike direction, and the minor axis radius of the search ellipsoid is taken from the range of the theoretical variogram of the dip direction. The radii of the search ellipsoid in the plane directions perpendicular to the strike direction and the dip direction are taken as the side length of the initial cubic cell. Data points containing the grade values ​​of ore-forming elements within the search ellipsoid are collected as known sample points for Kriging interpolation. The number of known sample points is no less than a preset minimum number of sample points. This preset minimum number of sample points is set according to the basic requirements for the stability of Kriging interpolation in geostatistics, for example, a minimum of 4 sample points. If the number of known sample points within the search ellipsoid is insufficient, the radius of the search ellipsoid is increased by a set ratio, for example, 1.2 times, and this process is repeated until the minimum number of sample points is met. The covariance between the known sample points and the covariance between the known sample points and the point to be estimated are calculated using a variogram model. The covariance is obtained by converting the theoretical variogram value determined by the variogram model. The conversion relationship is that the covariance equals the sill value of the variogram model minus the variogram value of the variogram model at the distance between the known sample points, or the covariance equals the sill value of the variogram model minus the variogram value of the variogram model at the distance between the known sample points and the point to be estimated. A covariance matrix for the known sample points is constructed based on the covariance between the known sample points, and a covariance vector for the point to be estimated is constructed based on the covariance between the known sample points and the point to be estimated. A system of Kriging equations is established using the ordinary Kriging method. This system includes unbiased constraints, meaning it minimizes the estimation error variance of the point to be estimated under the condition that the sum of the weighting coefficients equals 1. The weighting coefficients for each known sample point are obtained by solving the Kriging equations. The grade values ​​of the ore-forming elements at each known sample point are multiplied by their corresponding weighting coefficients and then summed to obtain the grade estimate of the point to be estimated, i.e., the grade estimate of the cubic unit. During the Kriging interpolation process, the Kriging interpolation variance is simultaneously obtained. This variance represents the minimum estimation error variance after solving the Kriging equations and is used in subsequent steps to calculate the interpolation uncertainty. The Kriging interpolation process is repeated for all cubic units to obtain the grade estimate of each cubic unit in the three-dimensional space of the target mining area.

[0035] In implementing S4, leave-one-out cross-validation is performed on known sample points. Known sample points are removed one by one, and kriging interpolation is performed on the removed locations using a variogram model to obtain interpolated estimates for those locations. Known sample points are sampling points with ore-forming element grades from the data point set. Each known sample point contains the sampling location coordinates and the ore-forming element grade value. The leave-one-out cross-validation process is as follows: one known sample point is removed from all known sample points at a time, and the sampling location coordinates of the removed known sample point are used as the location to be estimated. Kriging interpolation is then performed on the removed locations using the remaining known sample points. The method for Kriging interpolation at the removal location is the same as the method for Kriging interpolation of ore-forming element grades at cubic units. That is, a search ellipsoid is constructed based on the variogram model with the removal location as the center. The search ellipsoid determined by the variogram model satisfies the following: the major axis radius of the search ellipsoid takes the range of the theoretical variogram in the strike direction, the minor axis radius of the search ellipsoid takes the range of the theoretical variogram in the dip direction, and the azimuth of the major axis of the search ellipsoid is consistent with the strike and occurrence information of the main ore-controlling fault structures. The remaining known sample points within the search ellipsoid are collected, and the covariance between the remaining known sample points and the covariance between the remaining known sample points and the removal location are calculated using the variogram model. The covariance is obtained by converting the theoretical variogram value in the strike direction or the theoretical variogram value in the dip direction determined by the variogram model. The covariance is equal to the sill value of the variogram model minus the variogram value of the variogram model at the corresponding distance. Construct the covariance matrix of the remaining known sample points and the covariance vector of the removal positions. Solve the Kriging equations under the unbiased constraint that the sum of the weight coefficients equals 1 to obtain the weight coefficients of each remaining known sample point. Multiply the ore-forming element grade values ​​of each remaining known sample point by their corresponding weight coefficients and sum them to obtain the interpolated estimate of the removal position. Repeat the leave-one-out cross-validation process for all known sample points, and an interpolated estimate of the removal position will be obtained for each known sample point.

[0036] The difference between the interpolated estimate at the removal location and the actual grade value of the removed known sample point is calculated as the residual. For each removed known sample point, the actual grade value is the original ore-forming element grade value of that known sample point. The difference between the interpolated estimate at the removal location and the actual grade value of the known sample point is the residual. The sign of the residual indicates whether the interpolated estimate is overestimated or underestimated relative to the actual grade value of the known sample point, and the absolute value of the residual indicates the magnitude of the interpolation error. A residual set is obtained for all removed known sample points, and the number of residuals in the residual set is equal to the total number of known sample points. Each residual is associated with the sampling location coordinates of the corresponding known sample point.

[0037] Centered on each cubic unit, residuals falling into the search ellipsoid determined by the variogram model are collected, and the root mean square value of the fallen residuals is calculated as a local residual index. For each cubic unit, the center point of the cubic unit is taken as the center of the search ellipsoid, and the search ellipsoid determined by the variogram model is constructed. The search ellipsoid determined by the variogram model satisfies the following: the major axis radius of the search ellipsoid takes the range of the theoretical variogram function in the strike direction, the minor axis radius of the search ellipsoid takes the range of the theoretical variogram function in the dip direction, and the azimuth of the major axis of the search ellipsoid is consistent with the strike and occurrence information of the main ore-controlling fault structures. Among all residuals associated with the sampling position coordinates, those whose sampling position coordinates fall within the search ellipsoid are collected. The method for determining whether a sampling position coordinate falls within the search ellipsoid is to calculate the normalized distance between the sampling position coordinates and the center point of the cube element along each axis of the search ellipsoid. The normalized distance is the projected distance between the sampling position coordinates and the center point of the cube element along each axis of the search ellipsoid divided by the radius of the corresponding axis. A normalized distance less than or equal to 1 indicates that the sampling position coordinates fall within the search ellipsoid. The collected residuals falling within the search ellipsoid form a local residual set. The number of residuals in the local residual set is denoted as the local residual count. The root mean square (RMS) value is calculated for all residuals in the local residual set. The RMS value is calculated by summing the squares of each residual in the local residual set, dividing by the local residual count, and then taking the square root. The resulting value is the local residual index of that cube element. If the local residual count is zero, meaning no residual falls within the search ellipsoid, the local residual index of that cube element is set to zero. The local residual index reflects the overall level of cross-validation residuals within a local spatial region centered on a cube cell. The larger the local residual index, the greater the deviation between the interpolated estimate and the actual grade value of the known sample point within that local spatial region.

[0038] Extracting the Kriging interpolation variance generated during the Kriging interpolation process. When performing Kriging interpolation of the ore-forming element grade for each cubic cell, solving the Kriging equations simultaneously generates the Kriging interpolation variance. The Kriging interpolation variance is the minimum estimation error variance of the point to be estimated, and its value is determined by the solution to the Kriging equations. For each cubic cell, the Kriging interpolation variance is directly extracted from the Kriging interpolation results. The dimension of the Kriging interpolation variance is the square of the unit of the ore-forming element grade value, i.e., grams per ton squared. The Kriging interpolation variance reflects the degree of uncertainty in the interpolation estimate evaluated solely based on the spatial configuration and variogram model of known sample points during the Kriging interpolation process.

[0039] The interpolation uncertainty is obtained by multiplying the local residual index by the Kriging interpolation variance. For each cube cell, the interpolation uncertainty is calculated by multiplying the local residual index of that cube cell by the Kriging interpolation variance of that cube cell; the product is the interpolation uncertainty. Both the local residual index and the Kriging interpolation variance are non-negative, therefore the interpolation uncertainty is also non-negative. The interpolation uncertainty combines the overall level of cross-validation residuals within the local spatial region reflected by the local residual index and the estimation uncertainty based on the spatial configuration reflected by the Kriging interpolation variance, forming the interpolation uncertainty used for subsequent confidence screening of cube cells.

[0040] In the specific implementation of S5, the grade estimate is compared with the set grade threshold, and cubic units with a grade estimate greater than the set grade threshold are selected. The set grade threshold is the lower limit value of the ore-forming element grade used to define the boundary of the ore body. The set grade threshold is obtained by comprehensively analyzing the ore's feed grade and economically recoverable grade based on the industrial indicators of similar gold deposits in the study area, combined with the beneficiation data accumulated during the mining process. For example, if the feed grade of similar gold deposits in the study area is 1.0 g / t, then the set grade threshold is set to 1.0 g / t. For the grade estimate of each cubic unit obtained by kriging interpolation of the ore-forming element grade, the grade estimate is compared with the set grade threshold one by one. If the grade estimate is greater than the set grade threshold, the cubic unit is marked as a qualified grade cubic unit; if the grade estimate is less than or equal to the set grade threshold, the cubic unit is marked as an unqualified grade cubic unit.

[0041] The interpolation uncertainty is compared with a set uncertainty threshold, and cubic cells with interpolation uncertainties less than the set uncertainty threshold are selected. The set uncertainty threshold is determined based on the interpolation uncertainty corresponding to the known ore body. The determination process involves extracting the interpolation uncertainty corresponding to each cubic cell within the spatial distribution of the known ore body. The spatial distribution of the known ore body is the spatial range defined by the spatial distribution of the known ore body outline in the 3D geological model data. Within the spatial distribution range of the known ore body, the position of each cubic cell is determined. The interpolation uncertainty value calculated for each cubic cell in the step of calculating the interpolation uncertainty of each cubic cell is read. The maximum value among all the read interpolation uncertainty values ​​is taken as the set uncertainty threshold. For example, if the read interpolation uncertainty values ​​are 0.08, 0.12, 0.06, and 0.15, then the set uncertainty threshold is 0.15. For each cube element in the three-dimensional space of the target mining area, the interpolation uncertainty is calculated and compared with the set uncertainty threshold one by one. If the interpolation uncertainty is less than the set uncertainty threshold, the cube element is marked as a qualified cube element. If the interpolation uncertainty is greater than or equal to the set uncertainty threshold, the cube element is marked as an unqualified cube element.

[0042] The spatial locations of cube units that simultaneously satisfy the condition that the grade estimate is greater than a set grade threshold and the interpolation uncertainty is less than a set uncertainty threshold are merged to form the spatial distribution range of the ore body. For each cube unit in the three-dimensional space of the target mining area, it is determined whether the cube unit belongs to both the grade qualified cube unit and the confidence qualified cube unit. That is, it is determined whether the cube unit simultaneously satisfies the two conditions that the grade estimate is greater than a set grade threshold and the interpolation uncertainty is less than a set uncertainty threshold. Cube units that simultaneously satisfy both conditions are selected, and the spatial locations of these cube units are merged. The spatial location merging method is to connect adjacent cube units to form a continuous ore body space. Adjacent cube units refer to cube units that are adjacent in position after translation along the strike direction, or adjacent in position after translation along the dip direction, or adjacent in position after translation along the plane perpendicular to the strike and dip directions. The continuous orebody space formed after the merger is the orebody spatial distribution range. The orebody spatial distribution range contains one or more orebody spatial volumes that are not connected to each other. Each orebody spatial volume corresponds to an independent orebody. The orebody spatial distribution range defines the three-dimensional spatial location, shape and scale of the deep orebody.

[0043] In the specific implementation of S6, the grade estimate of each cubic unit is extracted from the spatial distribution range of the ore body. The spatial distribution range of the ore body is a set of continuous ore body spatial volumes formed by merging the spatial locations of cubic units that simultaneously satisfy the condition that the grade estimate is greater than a set grade threshold and the interpolation uncertainty is less than a set uncertainty threshold. The spatial distribution range of the ore body contains one or more ore body spatial volumes that are not connected to each other, and each ore body spatial volume is composed of several adjacent cubic units connected together. For each cubic unit within the spatial distribution range of the ore body, the grade estimate is the ore-forming element grade estimate calculated for that cubic unit when performing kriging interpolation of the ore-forming element grade of the cubic units divided in the three-dimensional space of the target mining area using a variogram model. The unit of the grade estimate is grams per ton. The method for extracting the grade estimate is to read the three-dimensional spatial location identifier of each cubic unit within the spatial distribution range of the ore body, look up the corresponding grade estimate value in the kriging interpolation result data based on the three-dimensional spatial location identifier, and associate the found grade estimate value with that cubic unit.

[0044] Obtain the unit volume of each cube element. The unit volume is the spatial volume occupied by the scaled cube element in the three-dimensional space of the target mining area. The unit volume is calculated by multiplying the side lengths of the cube element in three directions. The side lengths of the scaled cube element in the three directions are the side length in the strike direction, the side length in the dip direction, and the side length in the plane perpendicular to the strike and dip directions. The side length in the strike direction is the original cube element's side length multiplied by the anisotropy ratio. The side lengths in the dip direction and the plane perpendicular to the strike and dip directions are the original cube element's side lengths. The original cube element's side length and anisotropy ratio were determined in the step of scaling the original cube element's side length according to the anisotropy ratio. The unit volume of the unit volume is cubic meters. For example, if the initial cubic unit has a side length of 20 meters and an anisotropy ratio of 2, then the side length in the strike direction is 40 meters, the side length in the dip direction is 20 meters, and the side length in the plane perpendicular to the strike and dip directions is 20 meters. The unit volume is 40 meters × 20 meters × 20 meters = 16,000 cubic meters. The calculation method for the unit volume of each cubic unit within the spatial distribution range of the ore body is the same, that is, the unit volume is determined based on the side length of the scaled-down cubic unit.

[0045] The resource quantity of the deep ore body is obtained by summing the product of the grade estimates and the unit volumes of all cubic units within the spatial distribution range of the ore body. For each cubic unit within the spatial distribution range of the ore body, the product of the grade estimate and the unit volume is calculated. The grade estimate is in grams per ton, and the unit volume is in cubic meters. The product of the grade estimate and the unit volume is then multiplied by the ore density. The ore density is obtained by measuring the density of core samples collected from known ore bodies and taking the arithmetic mean. The ore density is in tons per cubic meter. The product of these is then used to obtain the amount of ore-forming element metal contained in the cubic unit, in grams. For example, if the ore density is 2.7 tons per cubic meter, the amount of ore-forming element metal is calculated by multiplying the grade estimate by the unit volume by 2.7. The total metal content of the deep ore body is obtained by summing the metal content of ore-forming elements in all cubic units within the spatial distribution range of the ore body. The unit of the total metal content of the deep ore body is grams. Divide the total metal content of the deep ore body by 1000 to convert it to the resource quantity of the deep ore body in kilograms, or divide the total metal content of the deep ore body by 1,000,000 to convert it to the resource quantity of the deep ore body in tons. The resource quantity of the deep ore body is the total mass of ore-forming elements contained within the spatial distribution range of the ore body.

[0046] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0047] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0048] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and inventive constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0049] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0050] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0051] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for estimating deep gold resources based on geostatistics, characterized in that, Includes the following steps: S1: Obtain the three-dimensional geological model data of the target mining area; S2: Based on three-dimensional geological model data, search and construct variation function models that characterize the spatial variation structure of ore-forming element grades along the strike and dip directions of the main ore-controlling structures. S3: Using the variogram model, kriging interpolation is performed on the ore-forming element grades of the cubic units divided in the three-dimensional space of the target mining area to obtain grade estimates. S4: Calculate the interpolation uncertainty of each cube element based on the leave-one-out cross-validation residuals of known sample points and the variance of Kriging interpolation; S5: Delineate the spatial distribution range of the ore body based on the grade estimate and interpolation uncertainty; S6: Calculate the resource quantity of deep ore bodies based on the delineated spatial distribution range of the ore bodies.

2. The method for estimating deep gold resources based on geostatistics according to claim 1, characterized in that, Obtain 3D geological model data of the target mining area, including: Borehole logging data, exploration line profile data, and structural geological maps were extracted from geological data of the study area. Based on borehole logging data and exploration line profile data, determine the spatial location of stratigraphic boundaries, rock mass boundaries, and known ore body outlines; Extract the strike, dip, and occurrence information of the main ore-controlling fault structures from structural geological maps; Spatial correlation is formed between the spatial location of stratigraphic boundaries, rock mass boundaries, and known ore body outlines and the strike, dip, and occurrence information of major ore-controlling fault structures to create three-dimensional geological model data.

3. The method for estimating deep gold resources based on geostatistics according to claim 1, characterized in that, Based on three-dimensional geological model data, variogram models characterizing the spatial variation structure of ore-forming element grades were searched and constructed along the strike and dip directions of the main ore-controlling structures, including: Extracting the ore-forming element grade values ​​of known ore bodies from three-dimensional geological model data; Based on the strike and occurrence information of the main ore-controlling fault structures, the data point pairs corresponding to the grade values ​​are searched along the strike direction with lag distance, the experimental variability function of the strike direction is calculated and fitted into the theoretical variability function of the strike direction. Based on the dip and occurrence information of the main ore-controlling fault structures, the data point pairs corresponding to the grade values ​​are searched along the dip direction with lag distance, the experimental variability function of the dip direction is calculated and fitted into the theoretical variability function of the dip direction. The variogram of the direction theory and the variogram of the tendency theory are combined to form a variogram model.

4. The method for estimating deep gold resources based on geostatistics according to claim 3, characterized in that, The fitting of the theoretical variogram of the direction of travel includes: using a spherical model to perform least squares fitting on the discrete points of the experimental variogram of the direction of travel to obtain the theoretical variogram of the direction of travel; the fitting of the theoretical variogram of the tendency direction includes: using a spherical model to perform least squares fitting on the discrete points of the experimental variogram of the tendency direction to obtain the theoretical variogram of the tendency direction.

5. The method for estimating deep gold resources based on geostatistics according to claim 1, characterized in that, Using a variogram model, kriging interpolation is performed on the ore-forming element grades of cubic units in the three-dimensional spatial division of the target mining area to obtain grade estimates, including: Extract the ranges of the directional theory variogram and the tendency theory variogram from the variogram model; The ratio of the range of the variogram of the direction theory to the range of the variogram of the tendency theory is used as the anisotropy ratio. The side lengths of the initial cube elements are scaled according to the anisotropy ratio so that the ratio of the side length in the orientation direction to the side length in the dip direction is equal to the anisotropy ratio. The three-dimensional space of the target mining area is divided using scaled cubic units; For each cubic cell, the grade of ore-forming elements is estimated by kriging interpolation using a variogram model.

6. The method for estimating deep gold resources based on geostatistics according to claim 1, characterized in that, Based on the leave-one-out cross-validation residuals of known sample points and the variance of Kriging interpolation, the interpolation uncertainty of each cube element is calculated, including: Leave-one-out cross-validation is performed on the known sample points, and the known sample points are removed one by one. Kriging interpolation is performed on the removal positions using the variogram model to obtain the interpolated estimates of the removal positions. The difference between the interpolated estimate at the removal location and the actual grade value of the removed known sample point is calculated as the residual. Centered on each cubic unit, the residuals falling into the search ellipsoid determined by the variogram model are collected, and the root mean square value of the falling residuals is calculated as a local residual index. Extracting the Kriging interpolation variance generated during the Kriging interpolation process; The interpolation uncertainty is obtained by multiplying the local residual index by the Kriging interpolation variance.

7. The method for estimating deep gold resources based on geostatistics according to claim 6, characterized in that, The search ellipsoid determined by the variogram model satisfies the following: the major axis radius of the search ellipsoid is taken as the range of the theoretical variogram in the strike direction; the minor axis radius of the search ellipsoid is taken as the range of the theoretical variogram in the dip direction; and the major axis azimuth of the search ellipsoid is consistent with the strike and occurrence information of the main ore-controlling fault structures.

8. The method for estimating deep gold resources based on geostatistics according to claim 1, characterized in that, Based on the grade estimate and interpolation uncertainty, the spatial distribution range of the ore body is delineated, including: The grade estimate is compared with the set grade threshold, and cubic cells with a grade estimate greater than the set grade threshold are selected. The interpolation uncertainty is compared with a set uncertainty threshold, and cubic cells with interpolation uncertainty less than the set uncertainty threshold are selected. The spatial locations of cubic units that simultaneously satisfy the condition that the grade estimate is greater than the set grade threshold and the interpolation uncertainty is less than the set uncertainty threshold are merged to form the spatial distribution range of the ore body.

9. The method for estimating deep gold resources based on geostatistics according to claim 8, characterized in that, The uncertainty threshold is determined based on the interpolation uncertainty corresponding to the known ore body: extract the interpolation uncertainty corresponding to each cubic unit within the spatial distribution of the known ore body; and use the maximum value of the extracted interpolation uncertainty as the set uncertainty threshold.

10. The method for estimating deep gold resources based on geostatistics according to claim 1, characterized in that, The resource quantity of deep ore bodies is calculated based on the delineated spatial distribution range of the ore bodies, including: Grade estimates for each cubic unit within the spatial distribution range of the ore body are extracted; Get the volume of each cube cell; The resource quantity of the deep ore body is obtained by summing the product of the grade estimates of all cubic units within the spatial distribution range of the ore body and the unit volume.