地矿工程可视化仿真中品位与储量的计算实现

地矿工程可视化仿真中品位与储量的计算实现1李翠平，李仲学，郝秀强，郝晋会（北京科技大学金属矿山高效开采与安全教育部重点实验室 北京（100083））Email: cpli@摘 要： 品位与储量计算是实现地矿工程可视化仿真系统的重要组成部分。

其实现首先需对原始采样数据进行预处理，以构造满足需要的规则体数据；之后基于不同插值方法下的规则体素信息，来计算任意指定范围下的矿石品位、储量等属性。

通过品位、储量计算过程的实现，也为矿床的三维可视化仿真工作提供了体素处理手段。

关键字： 地矿工程，空间插值，体素，品位、储量计算建立地矿工程的三维可视化仿真系统，实现矿床及其采矿工程的计算机三维及动态表示，对于优化开采设计方案、改善矿床开采设计精度、消除安全及灾害隐患、降低矿石开采成本、提高采矿生产经济效益与社会效益，具有重要的促进作用。

对于计算矿石品位与储量，则是地矿工程可视化仿真系统中不可缺少的重要内容[1,2]。

本文重在探讨仿真系统中品位与储量等矿石属性的计算与表达的实现工作。

1. 原始数据预处理地矿工程可视化仿真中要实现品位、储量计算以及矿体的真三维显示，需对原始的采样数据进行预处理，使之形成满足需要的体数据，作为实现矿体的真三维显示及属性（品位、储量）的计算与处理等工作的直接数据来源。

其预处理工作分成两个阶段，即由原始的采样数据转换为空间点数据和由初始的体数据转化为规则体数据两个阶段。

前者的数据处理一方面是通过原始数据的空间信息来取得各个钻孔及刻槽的空间坐标，进而得到各空间点的坐标值，另一方面是根据样品数据取得各个样品的品位、岩性、矿体号等属性值，进而得到各个空间点的属性值，从而确定各空间点的拓扑信息和属性信息的完整描述，即构成初始的体数据。

后者的数据处理是由初始体数据通过空间插值细化处理，得到规则的空间体数据形式，作为矿床可视化的基础数据来源[3]。

经过前一阶段的工作，将具有一定长度的地勘、生勘的钻孔采样数据及进路、沿脉、切割巷中的槽探数据处理成分布在空间中具有属性值的各离散点（这里的属性值就是品位值）。

而由地质勘探与生产勘探所获取的原始采样数据以及对其进行正则化处理所得到的离散的空间点数据，均具有明显的离散性、稀疏性及有限性。

而为实现地矿数据的体视化、地矿工程的真三维显示及赋存矿石平均品位、储量的计算，均需首先进行数据的空间属性插值处理，通过估算未知点属性得到想要的具有规则分布的数据。

因此插值计算的任务实际就是基于原始数据处理阶段得到的数据做进一步的数据处理，形成具有规则分布的数据，为矿石储量的计算以及为矿体三维显示提供基础。

其规则网格结构化体数据构型可表示为： ),,(),,(z k y j x i f z y x f ∆×∆×∆×= (1)1本课题得到了教育部高等学校博士学科点专项科研基金项目(20020008006)的资助。

- 1 -式中，i,j,k 分别为任意一个插值点(x,y,z)所在的行、列、层的序号；Δx,Δy,Δz 分别为网格点在三个方向上的间距，由此构成规则体素的三维尺寸；f(x,y,z)为体数据在网格点(i,j,k )处的属性值。

2. 空间插值方案设计2.1 泰森多边形法泰森（Thiessen）多边形法（又称Drichlet、Voronoi多边形法、最近邻点法）插值简单易行，在地学、资源、环境、气象等方面的研究中作为一种由点到面的插值方法得到了广泛的应用[4]。

其插值原理为：将整个数据平面按已知采样点的位置分割成若干个由泰森多边形表示的子区域，而每个泰森多边形的构成则是由相应的采样点与周围的所有邻域点间作垂直平分线，并将各垂直平分线依次连接组合而成，各泰森多边形内的每一点属性均由各多边形内的已知点确定。

由插值原理可见，泰森多边形法是一种基于二维平面的插值方法。

而将该插值方法引入到三维空间中，来解决三维空间地质体的属性插值（如品位插值），生成规则空间网格结构化的体数据，则需进行如下分析：由于地质勘探与生产勘探所获取的样品数据以及对其进行正则化处理所得到的离散的空间点数据，均具有沿取样方向数据相对其它方向分布集中特性，且地质勘探和生产勘探的取样方向常为大地坐标系的Z 轴方向。

故针对样品数据在Z 方向上分布相对集中甚至连续的特性，则在对原始采样数据进行正则化处理时，按照规则体素的Δz 间隔对原始的样品数据进行正则化，这样正则化后的空间离散点数据已是在Z 方向上按Δz 间隔分层的数据。

而此时对于Z 方向上每一分层的数据，其分布是散乱、稀疏的，其离散稀疏的程度则取决于勘探网度的大小，为此需经插值处理而使各分层的数据呈规则分布。

由于此时的插值只是在二维平面上进行，故采用泰森多边形法按规定的Δx、Δy 间隔对分层数据进行插值，沿Z 方向顺次处理各分层，直至所有分层均被插值处理。

这样，对于原本三维空间的插值问题，已转变为在若干二维平面上的插值问题，从而较轻松的将泰森多边形法插值引入到了三维空间中。

通过综合分析，利用泰森多边形法插值，待插值点的属性确定可由如下关系表达：⎪⎩⎪⎨⎧≥≤×−×=∑=0010011d d d d )O (Z )d d (n )A (Z n k k (2)式中,d 为待插值点与其控制点间的距离；d 0代表临界影响距离，其取值应根据所研究的具体问题来具体确定；n 为公共泰森多边形的个数，即公共控制点的个数；Z k (O)为第k 个控制点单独对未知点算出的属性值。

可见，无需真正构造泰森多边形，便可进行泰森多边形法的插值。

2.2 距离幂次反比法距离幂次反比法是由气象学家和地质工作者提出的，后来由于D. Shepard的工作，故又 - 2 -被称为Shepard法[5]。

其主要思想是认为未知属性被分布在周围的已知属性影响，这种影响与距离相关，距离越近影响越大，反之影响越小。

实践表明，该方法尤其对于构造不甚复杂的层状矿床及埋藏条件平稳的矿床，效果应用良好，故本文中设计了基于该方法的插值计算工作。

通过对距离幂次反比法的原理与应用情况的分析，该方法用于地矿工程需考虑两个问题：一是估值距离，即距未知点多远的已知点参加估值最为合理；二是幂指数，虽在应用中常取其值为2，不过应根据实际情况具体确定。

同时，由于该方法要求空间属性变量各向同性，若直接用于三维空间插值，其插值精度只有在空间属性变量各向同性时才具有可信度。

故采用与多边形插值法的同样考虑，把三维的插值计算转化成平面的二维插值，分别在每个层上进行插值计算，当每一层的插值计算都完成后，整个空间网格点上的属性即可得到。

在一个平面上对未知属性点A 实行距离幂次反比插值计算，需要在以A 点的坐标为圆心，以指定的插值距离为半径进行搜索，在搜索范围内进行对已知点的收集，采用下式计算A 点的属性Z *(A)，若在搜索范围内找不到已知属性点，则认为A 点的属性为零。

∑∑===N i Ni i i i d d g A Z 11*1)(µµ (3) 式中，N 是参与估值的已知属性点个数；d i 是第i 个已知属性点到未知属性点的距离；g i 是第i 个已知属性点的属性。

2.3 克立格法克立格（Kriging）法源于数学地质研究，是由法国学者G. Matheron等提出的、以南非采矿工程师D. G. Krige命名的一种地质统计学方法[6]。

该方法以区域化变量理论为基础、以变异函数为基本工具，能够用来研究分布于空间并呈现出一定结构性和随机性的数据。

克立格法处理的是区域化变量，区域化变量具有空间的连续性、变异性、结构性和相关性等特点。

故克立格法既能反映变量的空间结构特性，又能反映变量的随机分布特性。

同时还能充分考虑各采样点属性之间的相关性和变异性，轻松解决采样点之间的“丛聚效应”、各向异性等问题。

根据克立格法，研究对象空间某一未知点属性的估计值Z *(x 0)可以通过该未知点影响范围内n 个有效观察数据（即原始离散数据）的线性组合来得到，即∑==n i i i x Z x Z 10*)()(λ (4)式中，Z(x i )是指在观察点或实测点x i (=1,2,…,n)处的采样值，λi 是权系数，表示各采样值Z(x i )对估计值Z *(x 0)影响的大小，其求解要保证所求估计值无偏（即使所有待求点的估计值与其真实值之间偏差的期望值为0），并且估计方差最小。

对于无偏条件应满足各权系数之和等于1，即∑==n i i 11λ(5)- 3 -对于方差最小由下式给定∑=+=n i i i ex x 102),(ϕγλσ (6)其中j x x x x n i j j ii ∀=+∑=10),(),(γϕγλ (7)式中，Υ(x i ,x j )是Z 在已知点x i 和x j 之间的半方差；Υ(x j ,x 0)是Z 在已知点x j 和待求点x 0之间的半方差；Ф是极小化处理时的拉格朗日乘数。

这样，通过确定权系数λi 来求出Z *(x 0)，实现克立格法插值。

3. 储量算法设计通过上述三种插值方法的任何一种对原始数据进行预处理，均可得到体素大小相同的结构化体数据。

从而只要获得每一体素的具体属性，就得到了整个模型的属性。

对于各体素体积，只要确定了网格间距便为已知，而对于体素内任意一点的属性可通过八个角点的属性进行三元线性插值确定，对于各体素的平均品位可取体素中心点的品位。

而对于每一体素的体重值，则存在着与品位之间的函数关系：D=f(C)。

对于这一关系获得，严格上讲是通过对已知样品化验所得的品位和体重，来绘制品位体重函数图，根据此函数关系，来求出具体品位对应的体重值。

而在实际工作中，常依经验公式来推断，以满足现场生产的需要。

本文则采用后者，其关系为：u h uh i h h i C C G G C C G G −−•−−=)( (8)式中，G i 为第i 个体素的待求体重值；C i 为第i 个体素品位值；C h 、C u 分别为经验中高于和低于C i 品位的临界品位值；G h 、G u 分别为与C h 、C u 对应的临界体重值。

这样，求得了所有体素的体积、品位及对应体重值，便可计算出整个开采水平、任意矿块、甚至任意指定区域的矿石储量和金属量，其计算可利用下面的关系求得∑∑==i i i ii C V G Q VG T (9)式中，T 表示所求区域的矿石储量；Q 表示所求区域的矿石金属量；G i 表示第i 个体素的矿石体重；V i 表示第i 个体素的体积；C i 表示第i 个体素的品位。

4. 品位与储量计算的实现基于多边形法、距离幂次反比法和克立格法三种插值方法的结果数据，可以进行指定范围的矿石平均品位、储量及金属量的计算。

本文中提供了两种指定范围的方案，一是按指定的水平分段进行计算，一是用户自定义范围进行计算。