矿山离散钻孔样品变异函数模型计算与拟合.pdf
矿山离散钻孔样品变异函数模型计算与拟合 ① 荆永滨1, 王公忠1, 毕 林2 (1.河南工程学院 安全工程学院,河南 郑州 451191; 2.中南大学 资源与安全工程学院,湖南 长沙 410083) 摘 要 为建立用于矿产资源储量估算的变异函数模型,对矿山离散钻孔样品变异函数的计算与拟合方法进行了研究。 在定向实 验变异函数计算中,采用定向搜索棱柱代替搜索方向,用以搜索空间位置分布不规则的离散样品中满足滞后距条件的成对样品,给 出了对于离散样品的定向实验变异函数计算步骤及算法实现。 以实验变异函数图变程范围内的离散点及成对样品数目为主要依 据,对理论变异函数进行人工拟合,建立理论变异函数模型,最后通过交叉验证的方法进行最优理论模型的选择。 以某矿山铜元素 品位数据为例,利用组合样品计算矿体走向、倾向和厚度方向的实验变异函数,并分别采用 3 种理论模型进行拟合,根据交叉验证 结果计算平均误差、标准化平均误差、误差标准差,经过比较最终选择球状模型为最优理论模型。 结果表明,本文提出的变异函数 计算与拟合方法能准确建立矿山离散样品的变异函数模型,使矿石品位估值和资源储量计算结果更加可靠和准确。 关键词 地质统计学; 变异函数; 理论模型; 模型拟合 中图分类号 TD672; TP391.72文献标识码 Adoi10.3969/ j.issn.0253-6099.2017.04.005 文章编号 0253-6099(2017)04-0019-04 Calculating and Fitting Variogram Models for Dispersed Mine Drill Samples JING Yong-bin1, WANG Gong-zhong1, BI Lin2 (1. School of Safety Engineering, Henan Institute of Engineering, Zhengzhou 451191, Henan, China; 2.School of Resources and Safety Engineering, Central South University, Changsha 410083, Hunan, China) Abstract The calculating and fitting methods for variogram of dispersed drilling samples were studied for constructing a variogram model for geostatistical mineral resource/ reserve estimation. In the calculation of directional experimental variogram, directional search prism was proposed instead of search direction for searching those sample pairs meeting the lag distance requirement. Calculation steps and algorithm implementation of directional experimental variogram for dispersed drill samples were given. Based on the scattered points within the range and number of sample pairs in the graph of experimental variogram, the theoretical variogram was fitted manually to establish a theoretical variogram model. And an optimal theoretical model was selected through cross validation. With copper grade data from one mine taken as an example, the experimental variogram including orebody strike, dip direction and thickness of orebody, were calculated using a kind of composite sample. The obtained results were fitted with three theoretical models. Then, results from cross validation were taken for calculating mean error, standardized mean error and standard deviation of error. Finally, a spherical model was selected as the optimal theoretical model after comparison. It is shown that the proposed calculating and fitting method for variogram can exactly create a variogram model for dispersed drill samples in mines, leading to the grade estimation and reserve calculation with more reliability and accuracy. Key words geostatistic; variogram; theoretical model; fitting model 在矿石品位估值及储量计算中,对于矿山空间样 品数据相关性的研究,通常采用地质统计学的区域化 变量理论。 以区域化变量理论为基础,以变异函数为 工具,研究矿山空间样品数据的随机性和结构性。 以 变异函数为基础的各种地质统计学估值方法广泛应用 于固体矿产资源储量计算中,包括铜矿床、铁矿床、金 矿床、铅锌矿床、多金属矿床以及煤层等[1-6]。 变异函 数反映矿床变化程度的优劣影响到地质统计学方法矿 床品位和储量计算的准确性[7]。 Matheron 提出了经典 变异函数的计算方法,对于不同的滞后距,在滞后距容 差范围内搜索满足条件的成对样品,计算不同滞后距 对应的实验变异函数值。 针对样品数据在空间不规则 ①收稿日期 2017-02-17 基金项目 国家自然科学基金(41572317);河南省科技攻关项目(172102310643);河南工程学院博士基金(D2013021) 作者简介 荆永滨(1981-),男,河南郑州人,讲师,博士,主要研究方向为矿山三维建模及可视化、地质统计学。 第 37 卷第 4 期 2017 年 08 月 矿矿 冶冶 工工 程程 MINING AND METALLURGICAL ENGINEERING Vol.37 №4 August 2017 万方数据 分布的特点,为了保证计算的变异函数有足够的样品 点,研究人员提出对空间任一方向指定容差角,使该方 向的样品点搜索范围变成一个圆锥体。 对于理论变异 函数拟合的研究包括自动拟合和人工拟合 2 类。 自动 拟合方法主要包括加权最小二乘法、线性规划法和多 项式回归法等[8]。 自动拟合对于具有较好空间结构 性的样品能够得到理想的拟合结果,多用在环境、气象 等领域。 因此,本文对实验变异函数计算中满足滞后 距条件的成对样品,采用定向搜索棱柱的方法进行搜 索,并对相关搜索参数进行研究。 在成对样品搜索的 基础上,给出了实验变异函数计算的具体步骤和算法。 研究了人工拟合并通过交叉验证选择最优的理论变异 函数模型的方法。 1 变异函数 变异函数建模有 2 个步骤首先绘制实验变异函 数图,然后根据实验变异函数图选择合适的理论变异 函数模型对其进行拟合,得到理论变异函数。 实验变 异函数分析和确定理论变异函数的基本流程如图 1 所示。 Y N 样品数据是否足够 仅全向模型 可用 模型可用 模型不可用 模型不可用 模型可用 退出 计算全向变异函数 拟合理论变异函数模型 确定合适的变异函数计算方法 常规变异函数 指示变异函数 计算定向变异函数。 使用不同的搜素方向、角度 容差、带宽和滞后距 确定主轴、次轴和短轴方向 的实验变异函数模型 选项 1.获取更多数据 2.不使用地质统计学方法 选项 1.获取更多数据 2.使用各向异性模型 3.不使用地质统计学方法 图 1 变异函数建模流程 1.1 计算公式 对于有限的样品数据集 z(xi),i= 1,2,,n,变异 函数计算公式为 γ∗(h) = 1 2N(h)∑ N(h) i = 1 z(xi) - z(xi+ h)[] 2 (1) 式中 γ∗(h)为变异函数值;N(h)为由滞后距 h 分隔的 成对样品个数;z(xi)和 z(xi+h)分别为成对样品的底 部值和顶部值。 计算时,指定一组滞后距 h 的值,从而 得到相应的一组变异函数值。 1.2 变异函数图 变异函数图是根据一组滞后距 h 和相应 γ(h)值 绘制的二维图形,如图 2 所示,变异函数图的结构包括 块金 C0、变程 a、基台 C0+C(其中 C 为先验方差)。 h 0.20 0.15 0.10 0.05 0.00 200406080100120 块金 基台 变程 理论变异函数实验变异函数 h γ 图 2 变异函数图及其组成 实验变异函数图是离散的数据点,实验变异函数 不能直接用于地质统计学中的克立格和条件模拟等估 值方法,要根据实验变异函数图拟合出合适的理论变 异函数,获得相应类型理论变异函数参数,从而用于后 续的品位估值等计算。 理论变异函数有多种数学模 型,利用数学公式及相应的参数值来表示从实验变异 函数图中分析出的样品空间变异性,进一步用于地质 统计学估值算法。 2 离散样品变异函数计算 2.1 成对样品搜索 矿山地质勘探以及生产勘探的钻孔样品数据都是 在空间不规则分布的,根据式(1)计算变异函数时,对 于指定的滞后距,满足方向和距离的成对样品非常少, 甚至没有。 滞后距 h 的方向和距离应确定为一个范 围,从而有更多的成对样品数据参与变异函数的计算。 因此,在实际计算中,使用滞后距容差、容差角度、水平 和竖直带宽来实现[9]。 如图 3 所示,通过增加滞后距容差和方位角容差, 距离变为[4n-2,4n+2],方向变为[37.5,82.5],即 图中 2 条圆弧之间的区域。 再通过带宽减小距离较远 北 东 方 向 方 位 角 =60 滞后距 容差 2 m 带宽5 m 角 度容差22.5 滞后距4 m S0 S1 S2 S3 图 3 成对样品搜索二维图 02矿 冶 工 程第 37 卷 万方数据 样品的方向范围,最终满足条件的样品被限制在一个 棱锥范围内,对于图中的样品 S0,S1和 S2均满足要 求,S3则因带宽而不满足要求。 当容差角度增大到 90、带宽无限大时,计算出的为全向变异函数。 对于空间样品数据,计算某个方向的变异函数则 需增加滞后距容差、方位角容差、倾角容差、水平带宽 和竖直带宽,形成一个搜索棱柱,如图 4 所示。 方位角容差 倾角 容差 竖 直 带 宽 水 平 带 宽 图 4 空间离散样品定向搜索棱柱 2.2 变异函数值计算 设定变异函数滞后距的方向方位角和倾角分别为 α 和 β,滞后距距离为 d,滞后距数目为 n;方位角容差 和倾角容差分别为 α′和 β′,滞后距容差为 dt,水平带 宽和竖直带宽分别为 du和 dv。 循环取出全部成对样 品,对于每一对样品,其空间坐标记为 P1(x1,y1,z1) 和 P2(x2,y2,z2)。 样品之间的向量 h=P2 -P 1。 对向 量 h 的距离和方向进行判断,确定是否满足条件。 距离 h 为 h = Δx2+ Δy2+ Δz2(2) 式中 Δx=x2 -x 1,Δy=y2 -y 1,Δz=z2 -z 1。 判断 h 哪一 个滞后距范围,[d i-dt,d i+dt],i=1,,n。 判断方向的水平分量,即向量 h 方向的水平分量 与计算方向方位角的夹角小于方位角容差,且在水平 带宽范围内。 arccos Δxsinα + Δycosα Δx2+ Δy2 < αt (3) Δysinα - Δxcosα < u(4) 判断方向的垂直分量,即向量 h 方向的垂直分量 与计算方向倾角的夹角小于倾角容差,且在竖直带宽 范围内。 acos Δx2+ Δy2cosβ + Δzsinβ h < βt (5) Δzcosβ -Δx2 + Δy2sinβ < v (6) 对于向量 h 满足条件的成对样品,根据其底部值 和顶部值,计算变异函数的值。 循环全部样品后,计算 各滞后距满足条件的成对样品数目、滞后距平均值和 变异函数平均值,得到实验变异函数最终计算结果。 2.3 计算结果 数据来源为云南某铜矿,该矿山原始地质数据包 括地表钻探和生产勘探钻孔数据,共计 586 个钻孔。 铜元素特高品位临界值确定为 4.4%,去除特高品位样 品9 个。 样品组合长度取原始样品长度的平均值 1.2 m, 组合样品数量为 11 631。 计算组合样品铜元素品位全 向实验变异函数,结果见图 5。 h γ h 0.30 0.25 0.20 0.15 0.10 0.05 0.00 1020030405060708090 38028 57734 50440 55696 69334 95126 118780 128208 156908 183762 317606 394220 377402 382380 373298 371928 384872 方差 图 5 铜品位全向实验变异函数图 按地质要素法分别计算 Cu 元素在矿体走向、倾 向和厚度方向的实验变异函数,各方向参数见表 1。 各方向方位角容差和倾角容差均为 15,水平带宽和 竖直带宽均为 20 m,滞后距距离为 5 m,滞后距容差为 2.5 m,滞后距数目为 16。 计算结果见图 6。 表 1 实验变异函数计算方向 方向方位角/ ()倾角/ () 走向2850 倾向195 -29 厚度19561 h γ h 0.30 0.25 0.20 0.15 0.10 0.05 0.00 1020030405060708090 走向方向 倾向方向 厚度方向 ■ ■ ■ ■ ■ ■ ■■ ■ ■ ■ ■ ■ ■ ■ ■ ■ ■ 图 6 铜品位定向实验变异函数图 3 理论变异函数拟合 实验变异函数只是计算了一组离散滞后距对应的 变异函数估计值。 要通过变异函数对矿石品位等区域 化变量进行结构分析,或者是利用克立格方法进行属 12第 4 期荆永滨等 矿山离散钻孔样品变异函数模型计算与拟合 万方数据 性插值,必须根据这些离散的估计值拟合出理论变异 函数模型及模型相应的块金、基台和变程等参数。 理 论变异函数模型分为有基台和无基台 2 类,有基台的 模型主要包括球状模型、高斯模型和指数模型,无基台 模型主要包括纯块金效应模型、对数函数模型、幂函数 模型及空穴效应模型等。 在矿石品位的结构分析和品 位估值中,常用的有球状模型、高斯模型、指数模型和 幂函数模型[10]。 3.1 拟合方法 目前,在矿山工程领域,进行矿产矿化研究、矿石 品位估值时人工拟合理论变异函数仍然是主要做法。 依靠地质专家根据离散样品的变异函数图,选择合适 的理论模型,并改变模型参数,通过观察绘制出的变异 函数曲线的拟合程度人工判断理论模型拟合结果的好 坏。 人工拟合的过程实际上是一种获取专家知识库中 专业经验和对地质信息理解的综合判断方法。 理论变异函数人工拟合的步骤如下计算各滞后 距满足条件的成对样品数目、滞后距平均值和变异函 数平均值,得到实验变异函数最终计算结果。 1) 以实验变异函数计算结果中的滞后距平均值 和变异函数平均值绘制变异函数图,通常用折线连接 多个离散点,并在每个点的一侧标注成对样品数目。 2) 根据实验变异函数图的折线形态选择变异函 数理论模型的类型及初始参数,绘制理论变异函数 曲线。 3) 观察理论曲线与离散点的拟合程度,实际拟合 过程中对拟合程度的判断主要参考变程范围以内的离 散点,通常为前 3 个点。 重复步骤 2)得到多个较好的 拟合结果。 3.2 拟合结果 分别采用球状模型、高斯模型和指数模型对实验 变异函数进行拟合,其中走向方向的理论变异函数曲 线拟合结果如图 7 所示,理论变异函数模型参数拟合 结果见表 2。 h γ h 0.30 0.25 0.20 0.15 0.10 0.05 0.00 1020030405060708090 实验变异函数 球状模型 高斯模型 指数模型 166 327 428 354 644 2366 1072 797 2145 2175 20839 16509 3428 3480 2787 4620 8770 图 7 走向方向理论变异函数拟合 表 2 理论变异函数模型 类型方向 C0 Ca 走向45 球状模型倾向0.0330.1635 厚度18 走向35 高斯模型倾向0.0330.1630 厚度16 走向55 指数模型倾向0.0130.1835 厚度12 4 理论模型交叉验证 通过计算离散样品的实验变异函数并绘制变异函 数图,采用人工方式进行拟合,得到多个拟合程度较好 的理论变异函数模型。 此时,要采用交叉验证的方法 对多个模型进行比较,然后选择一个最好的理论模型。 交叉验证是通过地质统计学估值方法,利用变异函数 理论模型对样品进行属性估值,然后将估计值和真实 值进行比较分析。 需要注意的是,交叉验证只是一种 探索性方法,并不能确定性衡量变异函数模型的优劣, 只是通过它辅助选择变异函数模型。 根据估计值 Z∗(xi)和真实值 z(xi)计算统计值平均误差 ME 和平 均标准化误差 MSE[11] ME = 1 N∑ N i = 1 [z(xi) - Z∗(xi)](7) MSE = 1 N∑ N i = 1 z(xi) - Z∗(xi) σ∗(xi) (8) 式中 σ∗(xi)为估计值对应的克立格方差。 从统计值 上分析,ME 应为 0,与 ME 相比,MSE 考虑了克立格方 差,其值应接近于 0。 此外,还可以利用探索性数据分析方法对模型进 行选择,绘制诊断图形并计算相关统计特征值。 1) 样品估计值和真实值的散点图,散点应集中分 布在过原点且斜率为 1 的直线两侧,具有较少的离群 值,且相关系数较大。 2) 误差直方图,图形应以 0 为中心对称分布,且 标准偏差较小。 3) 估计值和误差的散点图,衡量估计值的无偏 性,散点应集中分布在零误差线上下。 4) 误差位置分布图,样品的正负误差应在空间均 匀分布。 对理论变异函数 3 个方向上的各向异性模型进行 套合,采用普通克立格法进行交叉验证,根据交叉验证 结果进行统计分析,结果见表 3。 由表 3 可知,球状模 型理论变异函数的交叉验证统计结果更优。 (下转第 27 页) 22矿 冶 工 程第 37 卷 万方数据 参考文献 [1] 王元汉,李卧东,李启光,等. 岩爆预测的模糊数学综合评判方法 [J]. 岩石力学与工程学报, 1998,17(5)493-501. [2] 何 正,李晓红,卢义玉. BP 神经网络模型在深埋隧道岩爆预测 中的应用[J]. 地下空间与工程学报, 2008,4(3)494-498. [3] 胡 炜,杨兴国,周宏伟,等. 基于灰色关联分析的岩爆预测方法 研究[J]. 人民长江, 2011,42(9)38-43. [4] 冯夏庭,赵洪波. 岩爆预测的支持向量机[J]. 东北大学学报, 2002,23(1)57-60. [5] 宫凤强,李夕兵. 岩爆发生和烈度分级预测的距离判别方法及应 用[J]. 岩石力学与工程学报, 2007,26(5)1012-1017. [6] 徐 飞,徐卫亚. 岩爆预测的粒子群优化投影寻踪模型[J]. 岩土 工程学报, 2010,32(5)719-725. [7] 葛启发,冯夏庭. 基于 AdaBoost 组合学习方法的岩爆分类预测研 究[J]. 岩土力学, 2008,29(4)943-948. [8] Breiman L. Random Forests[J]. Mach Learn, 2001,45(1)5-32. [9] Breiman L, Friedman J H, Olshen R A, et al. Classification and Regression Trees[M]. New YorkChapman and Hal, 1984. [10] 董师师,黄哲学. 随机森林理论浅析[J]. 集成技术, 2013,2(1) 1-7. [11] 李 亭,田 原,邬 伦,等. 基于随机森林方法的滑坡灾害危险 性区划[J]. 地理与地理信息科学, 2014,30(6)25-30. [12] 明均仁,肖 凯. 基于 R 语言的面向需水预测的随机森林方法 [J]. 统计与决策, 2012(9)81-83. [13] 方匡南,朱建平,谢邦昌. 基于随机森林方法的基金收益率方向 预测与交易策略研究[J]. 经济经纬, 2010(2)61-65. [14] 张 雷,王琳琳,张旭东,等. 随机森林算法基本思想及其在生态 学中的应用 以云南松分布模拟为例[J]. 生态学报, 2014 (3)650-659. [15] 温廷新,张 波,邵良杉. 煤与瓦斯突出预测的随机森林模型 [J]. 计算机工程与应用, 2014,50(10)233-237. [16] 周科平,雷 涛,胡建华. 深部金属矿山 RS-TOPSIS 岩爆预测模 型及其应用[J]. 岩石力学与工程学报, 2013,32(2)3706- 3711. [17] 郝 杰,侍克斌,王显丽,等. 基于模糊 C-均值算法粗糙集理论 的云模型在岩爆等级评价中的应用[J]. 岩土力学, 2016,37 (3)857-866. [18] Dong Longjun, Li Xibing, Peng Kang. Prediction of rockburst classi- fication using Random Forest[J]. Transactions of Noferrous Metals Society of China, 2013,23(2)472-477. [19] 贾义鹏,吕 庆,尚岳全. 基于粒子群算法和广义回归神经网络 的岩爆预测[J]. 岩石力学与工程学报, 2013,32(2)343-348. [20] Feng Xiating, Wang Lina. Rockburst prediction based on neural net- works[J]. Transactions of N F Soc, 1994,4(1)7-14. 引用本文 杨悦增,邓红卫,虞松涛. 基于随机森林模型的岩爆等级预 测研究[J]. 矿冶工程, 2017,37(4)23-27. �������������������������������������������������������������������������������������������������� (上接第 22 页) 表 3 交叉验证结果统计 类型 平均 误差 标准化 平均误差 误差 标准差 球状模型1.3910 -4 4.5110 -4 0.39 高斯模型-4.1510 -3 -1.2810- 2 0.41 指数模型2.2210 -4 8.0010 -4 0.39 5 结 论 1) 针对矿山钻孔样品空间分布不规则特点,在指 定方向的实验变异函数计算中,利用滞后距距离容差、 角度容差以及水平和竖直带宽形成定向搜索棱柱,用 以搜索离散样品中满足滞后距条件的成对样品。 给出 对于离散样品的实验变异函数计算步骤及算法实现。 2) 绘制实验变异函数图,在其基础上进行理论变 异函数的人工拟合,对样品通过克里格估值进行交叉 验证,选择最优理论模型。 3) 通过某矿山铜元素品位数据,计算实验变异函 数并采用 3 种理论模型进行拟合,对交叉验证结果进 行统计分析,确定球状模型为最优理论模型。 参考文献 [1] 揣媛媛,范继璋,肖克炎,等. 西岔金矿普通克里格法可视化储量 计算应用研究[J]. 石油天然气学报(江汉石油学院学报), 2006 (5)63-65. [2] 贾福聚,秦德先,黎应书,等. 变异函数在都龙锡多金属矿床的应 用[J]. 地质与勘探, 2008(2)77-81. [3] 罗周全,王中民,刘晓明,等. 基于地质统计学与 Surpac 的某铅锌 矿床储量计算[J]. 矿业研究与开发, 2010(2)4-6. [4] 王炯辉,李 毅,黄冬梅,等. 基于普通克里格法的泥河铁矿床资 源储量估算研究[J]. 地质与勘探, 2013(6)1108-1113. [5] 荆永滨,孙光中. 矿体块段模型品位估值快速搜索算法研究[J]. 矿冶工程, 2016(4)8-10. [6] 刘晓明,吕太含冰,陈建宏,等. 基于地质统计学的资源储量估算 参数优选研究[J]. 矿冶工程, 2015(2)23-28. [7] 侯景儒. 地质统计学在我国的应用及其发展[J]. 地质与勘探, 1991,27(4)36-38. [8] 徐武平,邱 峰,徐爱萍. 空间数据插值的自动化方法研究[J]. 武汉大学学报(信息科学版), 2016(4)498-502. [9] Deutsch C V, Journel A G. Gslib Geostatistical Software Library and User′s Guide[M]. New York, NY Oxford University Press, 1998. [10] 熊俊楠,马洪滨. 变异函数的自动拟合研究[J]. 测绘信息与工 程, 2008,33(1)27-29. [11] 徐爱萍,舒 红,译. 空间数据分析与 R 语言实践[M]. 北京清 华大学出版社, 2013. 引用本文 荆永滨,王公忠,毕 林. 矿山离散钻孔样品变异函数模型 计算与拟合[J]. 矿冶工程, 2017,37(4)19-22. 72第 4 期杨悦增等 基于随机森林模型的岩爆等级预测研究 万方数据