第五节 品位、矿量计算的垂直断面法
垂直断面法是传统的手工计算矿量的常用方法，其一般步骤为： 第一步：沿勘探线做垂直剖面，将勘探线上的钻孔及其取样品位标在剖面图上（图1-11）。
图 1-11 标有取样品位的剖面图
第二步：根据给定的边界品位进行矿体圈定。简单地讲，矿体圈定的过
程就是将相邻钻孔上高于边界品位的样品点相连的过程。当一条矿体被一个钻孔穿越，而在相邻的钻孔消失时，一般将矿体延伸到两钻孔的中点，或是根据矿体的自然尖灭趋势在两钻孔之间实行自然尖灭。在矿体圈定过程中，要充分考虑矿床的地质构造（如断层和岩性）和成矿规律。图1-12是当边界品位等于25％时根据图1-11中的取样品位圈定的矿体示意图。
第三步：矿体圈定完成后，可用求积仪求得每个断面上的矿石面积，然
后就可以进行矿量计算。
(a)当一条矿体在两个相邻断面上的面积（S1和S2）相差不到40
％时，两断面之间的矿体体积用下式计算：
V=
s1+s2
2
L
(1-31)
式中，L为断面间距。
(b)当两个相邻断面上的面积相差大于40％时，采用下式计算：
V=
s1+s2+1⋅s2
3
L
(1-32)
(c)当矿体在二断面间是楔形尖灭时，计算公式为：
V=
sL 2
计算出两断面间矿石块段体积后，矿石块段的矿量为：
T=V⋅γ
(1-33)
式中，γ为本块段矿石体重。然后将所有块段的矿量相加，即得矿床的总矿量。
图 1-12 边界品位为25％时矿体圈定示意图
第四步：计算矿体的平均品位： (a)对穿越矿体的每一钻孔的样品进行“矿段样品组合”，求出
组合样品的品位。
(b)求出每一组合样品的影响面积。该面积是以钻孔为中线向两
侧各外推二分之一钻孔间距得到的矿体面积。
(c)对组合样品品位以其影响面积为权值进行加权平均计算，求
出矿体在断面上的平均品位。
(d)一条矿体的总平均品位是该条矿体在各断面上的平均品位以
断面所代表的矿量为权值的加权平均值。
上述矿体圈定与矿量、品位计算过程，可通过编制程序由计算机完成，但矿体圈定的完全计算机化具有较高的难度。垂直断面法是一种传统的手工方法，在我国仍广泛采用，但在发达国家由于计算机的广泛应用已基本成为过去。
第六节 矿量、品位计算的水平断面法
在露天矿山，矿石的开采是分台阶进行的，因此用于矿量、品位计算的一个水平断面即为一个台阶。常用的水平断面法有多边形法和三角形法。 600N
300N
300E
600E
图 1-13 台阶水平面上钻孔取样分布示意图
一、多边形法(polygonal method)
在多边形法中，水平断面上的每一钻孔取样（即进行台阶样品组合后的组合样品），位于一个多边形的中心。多边形的形成由以下步骤完成：
第一步：把穿越水平面的钻孔根据钻孔坐标，绘于水平面上，并将本平
面的组合样品品位标注在图上（图1-13）。
第二步：根据经验和地质统计学分析，确定影响半径R，R的确定将在后
面的地质统计学部分介绍。
第三步：以每一样品为中心，确定其相邻样品，一般情况下相邻样品是
落在半径为2R的圈内的样品。以样品S-41为例，其相邻样品如图1-14a所示。
第四步：用直线将中心样品和相邻样品连接起来（图1-14a）。
第五步：在中心样品与每一相邻样品的连线中点作垂直于连线的直线
（称为二分线），这些二分线相交围成的多边形即为所求的多边形。当两条二分线近于平行时，在二者相交之前，将与另外一条二分线相交，这时，取二者中离中心样品最近者作为多边形的边。以样品S-41为中心的多边形如图1-14b所示。
( b ) ( a )
图 1-14 中心样品四周均有样品时多边形的形成过程
如果在某些区域，钻孔间距大于2R，就以中心样品为中心，以R为半径做一八边形。位于边缘上的样品，只在其一侧有相邻样品而在另一侧没有相邻样品。这时，在没有相邻样品的一侧以中心样品为中心，画一半径为R的圆。然后在0o,90o与±45o方向上分别作圆的切线，这些切线与有相邻样品一侧的二分线相交就形成了所求的多边形。图1-15是以边缘样品S－14为例的多边形形成过程。
( a )
( b )
图 1-15 边缘样品的多边形的形成过程
以每一样品为中心形成多边形后，在每一多边形内，品位被看作是常数且多边形的品位等于其中心样品的品位。第i个多边形的重量Ti由下式计算：
Ti=SiγH
(1-34)
式中，Si为第i个多边形的面积，γ为体重；H为台阶高度。 如果多边形的品位xi大于边界品位，该多边形即为矿石多边形，所有矿石多边形的集合就形成了该水平面上的矿体。基于图1-13中的样品分布，应用多边形法进行品位计算和矿体圈定的结果如图1-16所示。这里假设边界品位为0.6％。 水平断面上的矿石总量T为矿石多边形的重量之和，矿石的平均品位为矿石多边形品位的面积加权平均值，即：
T=∑Ti
i=1
n
(1-35)
=∑xisi
i=1
n
∑si
i=1
n
(1-36)
式中，xi为矿石多边形i的品位，n为矿石多边形个数。 在某些矿床中，往往不同区域的成矿作用与矿石种类不同，或是地质构造使矿体发生错动。在这种情况下，使用多边形法时，应注意多边形不跨越区域界限和构造线。
图1-16 多边形法品位估算与矿体圈定结果
图1-16 多边形品位估算与矿体圈定结果
二、三角形法
在多边形法中，中心样品的品位被外推到整个多边形，即一个多边形的平均品位只用了一个已知样品作估计。一般情况下，对一给定体积的品位进行估算时，利用的已知数据越多，估计误差越小。从而就产生了三角形法。 在三角形法中，每一样品是三角形的一个顶点。图1-17是将图1-13中的样品相连形成的三角形。很显然，在一定条件下顶点连法不同，所得到的三角形也不同。一般性的连接原则是使每个三角形的边尽可能短，面积尽可能小。 每一三角形的品位xi是位于其顶点上的三个样品品位的平均值，即：
xi=
xi1+xi2+xi3
3
(1-37)
品位高于边界品位的三角形是矿石三角形，水平面上所有矿石三角形的集合形成该水平面上的矿体。一个水平面上的总矿量为矿石三角形的重量之和，矿石的平均品位为矿石三角形品位的重量或面积加权平均值。 600N
300N
300E
图 1-17 三角形法品位估算与矿体圈定结果
三角形法的优点在于它利用三个样品的品位来估计一个三角形的平均品位。从理论上讲，利用三角形法求得的品位、矿量较多边形法误差小。在应用三角形法时，也要注意三角形不超越区域界限和地质构造线。 第七节 三维块状模型
三维块状模型是将矿床划分为许多单元块形成的离散模型（图1-18）。单元块一般是尺寸相等的长方体。随着计算机在矿山的普及应用和计算机的容量与速度的不断提高，三维块状模型在国际上得到越来越广泛的应用。三维块状模型不仅被广泛用于品位、矿量计算，也被用于
图1-18 三维块段模型示意图
露天矿最终开采境界优化和开采计划优化。实际上，许多优化方法是由于三维块状模型的引入而产生的，而且当今矿山优化界正在进行的各种优化方法的研究，绝大多数以三维块状模型为基础。在我国，虽然三维块状模型的概念引入较早，但由于观念的陈旧与国家有关储量计算和矿山设计规范的束缚等原因，三维块状模型除在矿山设计院有少量应用外，至今没有发挥应有的作用。 三维块状模型中单元块的高度等于露天矿台阶高度，单元块在水平方向一般取正方形，其边长视具体情况而定。有人错误地认为，单元块越小，品位、矿量计算结果越精确。但是，从地质统计学理论可知，在已知数据（即样品数与样品在空间的分布）一定的条件下，单元块越小，对其品位的估计误差越大。另外，当单元块小到一定程度时，相邻的几个单元块的品位估计值会非常接近，与它们的平均品位相差无几，这种现象称为平滑作用，故用一个大块代替几个小块，品位与矿量的计算结果变化很小。而这样做可以降低对计算机容量的要求，加快计算速度。一般的经验规则是，单元块在水平方向上的边长不应小于钻孔平均间距的1/4或1/5。对于100米的钻孔间距，单元块的边长一般取30米左右。 将矿床分为单元块后，需要应用某种方法对每一小块的平均品位进行估计。常用的方法有三，即最近样品法、距离N次方反比法和地质统计学法（即克里金法）。三者均基于样品加权平均的概念，即对落在以单元块为中心的影响范围内的样品品位进行加权平均求得单元块的品位。三种方法的根本区别在于所用权值不同。本节介绍前两种方法，第三种方法将在下节给予较详细的论述。 一、最近样品法
所谓最近样品法就是将距离某一单元块最近的样品品位作为该单元块的品位估计值。前面介绍的多边形法其实是不规则单元块情况下的最近样品法，最近样品法的一般步骤为：
第一步：以被估计的单元块的中心为圆心，做半径为影响半径R的圆。 第二步：计算落入影响范围内的每一样品与单元块中心点的距离。
第三步：选取离单元块中心最近的样品，其品位即为被估单元块的品位。 当没有样品落入影响范围时，被估单元块的品位是未知的。一般情况下，未知单元块的品位取0值，当废石处理。如果有理由认为未知单元块所处的区域可能有矿石，未知单元块的出现表明该区域的数据量太少，要想确定其品位与矿量，需要增加钻孔。 在三维状态下，上述的影响范围，由二维状态下的圆变为球，需要对落入球内的所有样品进行考查，找出离单元块中心最近的样品。 求得矿床中所有单元块的品位以后，品位大于边界品位的单元块的集合组成矿体。矿石量和矿石平均品位可由矿石单元块重量的简单累加
和品位的平均求得。基于图1-13所示的样品分布和品位值，用最近样品法求得该台阶上单元块的品位如图1-19所示。这里采用的影响半径(R)为75m，边界品位(gc)为0.6％。将图1-19与图1-16作比较可见，用多边形法与最近样品法所得的矿体形态相近。
600E
图1-19 最近样品法品位估算与矿体圈定结果
二、距离N次方反比法(Inverse Distance Method)
在多边形法和最近样品法中，只有一个样品参与单元块品位的估值，如果落入影响范围的样品都参与单元块的品位估值，估值结果会更为精确。然而，由于各样品距单元块中心的距离不同，其品位对单元体的影响程度也不同。显然，距离单元块越近的样品，其品位对单元体品位的影响也就越大。因而在计算中，离单元体近的样品的权值应比离单元体远的样品的权值大。距离N次方反比法就是基于这一思想产生的。在此法中，一个样品的权值等于样品到单元块中心距离的N次方的倒数（1/dN）。 参照图1-20，距离N次方反比法的一般步骤如下：
第一步：以被估单元块中心为圆心、以影响半径R为半径做圆，确定影
响范围（在三维状态下，圆变为球）。
第二步：计算落入影响范围内每一样品与被估单元块中心的距离。 第三步：利用下式计算单元块的品位Xb：
xi
x=∑N b
i=1di
n
∑dN
i=1
i
n
1
(1-38)
式中，xi为落入影响范围的第i个样品的品位；di为第i个样品到单元块中心的距离。
图1-20 距离N次方反比法示意
在实际应用中，有时采用所谓的角度排除，即当一个样品与被估单元块中心的连线与另一个样品与被估单元块中心的连线之间的夹角小于某一给定值α时，距单元块较远的样品将不参与单元块的估值运算（如图1-20中的G3与G5）。α值一般在15度左右。如果没有样品落入影响范围之内，单元块的品位为零。公式(1-38)中的指数N对于不同的矿床取值不同。假设有两个矿床，第一个矿床的品位变化程度较第二个矿床的品位变化程度大，即第二个矿床的品位较第一个矿床连续性好。那么在离单元体同等距离的条件下，第一个矿床中样品对单元块品位的影响应比第二个矿床小。因此，在估算某一单元块的品位时第一个矿床中样品的权值在同等距离条件下应比第二个矿床中样品的权值小。也就是说在品位变化小的矿床，N取值较小；在品位变化大的矿床，N取值较大。在铁、镁等品位变化较小的矿床中，N一般取2；在贵重金属（如黄金）矿床中，
N的取值一般大于2，有时高达4或5。如果有区域异性存在，不同区域中品位的变化不同，则需要在不同区域取不同的N值。同时，一个区域的样品一般不参与另一区域的单元块品位的估值运算。以图1-20为例，若n=2，则被估单元块的品位为0.628。 600N
300N
300E 600E
图1-21 距离平方反比法品位计算与矿体圈定结果
基于图1-13所示的样品分布和品位值，应用上述方法，当N＝2，R=75m，gc=0.6％时，品位计算和矿体圈定结果如图1-21所示。将图1-16，1-19和1-21所示的计算结果作比较，应用距离平方反比法求得的矿体形态与前两种方法（多边形法与最近样品法）得出的矿体形态之间的差别较为明显。
第八节 地质统计学法（Geostatistical Method）
地质统计学是60年代初期出现的一个新兴应用数学分支，其基本思想是由南非的Danie Krige在金矿的品位估算实践中提出来的，后来由法国的Georges Matheron经过数学加工，形成一套完整的理论体系。在过去的三十多年中，地质统计学不仅在理论上得到发展与完善，而且在
实践中得到日益广泛的应用。如今，地质统计学在国际上除被用于矿床的品位估算外，也被用于其他领域中研究与位置有关的参数变化规律和参数估计。如农业中农作物的收成、环保中污染物的分布等等。本节将从矿床的品位估算的角度，简要介绍地质统计学的基本概念、原理和方法。
一、区域化变量、协变异函数与半变异函数
应用传统统计学（“传统”二字是相对于地质统计学而言的）可以对矿床的取样数据进行各种分析，并估计矿床的平均品位及其置信区间。在给定边界品位时，传统统计学也可用于初步估算矿石量和矿石平均品位。然而，传统统计学的分析计算均基于一个假设，即样品是从一个未知的样品空间随机选取的，而且是相互独立的。根据这一假设，样品在矿床中的空间位置是无关紧要的，从相隔上千米的矿床两端获取的两个样品与从相隔几米的两点获取的两个样品从理论上讲是没有区别的，它们都是一个样本空间的两个随机取样而已。 但是在实践中，相互独立性是几乎不存在的，钻孔的位置（即样品的选取）在绝大多数情况下也不是随机的。当两个样品在空间的距离很小时，样品间会存在较强的相似性，而当距离很大时，相似性就会减弱或不存在。也就是说，样品之间存在着某种联系，这种联系的强弱是与样品的相对位置有关的。这样就引出了“区域化变量”的概念。
（一）区域化变量及协变异函数： 如果以空间一点为中心获取一样品，样品的特征值X(z)是该点的空间位置z的函数，那么变量X即为一区域化变量。 显然，矿床的品位是一个区域化变量，而控制这一区域化变量之变化规律的是地质构造和矿化作用。区域化变量的概念是整个地质统计学理论体系的核心，用于描述区域化变量变化规律的基本函数是协变异函数和半变异函数。 设有两个随机变量X1与X2，如果X1与X2之间存在某种相关性，那么从传统统计学可知，这种相关关系由X1与X2的协方差σ(x1,x2)表示：
σ(x1,x2)=E[(x1-E(x1))(x2-E(x2))] (1-39)
22
让σx和分别表示X1和X2的方差，则： σx12
2
σx=E[(x1-E(x1))2]
1
(1-40) (1-41)
2σx=E[(x2-E(x2))2]
2
式中，E[*]表示随机变量[*]的数学期望。 X1与X2之间的相关系数为：
ρx⋅x=
1
2
σ(x1,x2)σx⋅σx
1
2
(1-42)
当X1和X2互相独立时，即二者之间不存在任何相关性时，协方差与相关系数均为零。当X1和X2“完全相关”时，相关系数为1.0（或-1.0）。 如果X1和X2不是一般的随机变量，而是区域化变量X在矿体Ω中的取值，即： X1代表 X(z)：区域化变量X在矿体Ω中z点的取值； X2代表X(z+h)：区域化变量X在矿体Ω中距z点h处的取值。
那么，由式(1-39)可计算X(z)与X(z+h)在矿体Ω中的协方差：
σ(x(z),x(z+h))=σ(h)
](x(z+h)-E[x(z+h))]]=E[(x(z)-E[x(z))
(1-43)
σ(h)称为区域化变量X在Ω中的协变异函数(Covariogram)。
22
让σ1和σ2分别表示X(z)与X(z+h) 在矿体Ω中的方差,则：
2
σ1=E[(x(z)-E[x(z)])2] 2σ2=E[(x(z+h)-E[x(z+h)])2]
(1-44) (1-45)
那么X(z)与X(z+h)之间的相关系数为：
σ(h)ρ(h)=
σ1σ2
(1-46)
ρ(h)称为区域化变量X在Ω中的相关函数(Correlogram)。 对于任何矿床，都可能计算出其协变异函数σ(h)，但在利用σ(h)对
矿床中单元体的品位进行估值时，需满足二阶稳定性条件。
[二阶稳定性条件](Second order stationarity conditions):
● X(z)的数学期望与空间位置z无关，即对任意位置z0 ：
E[x(z0)]= μ
(1-47)
● 协变异函数与空间位置无关只与距离h有关，即对于任何位置
z0 ：
E[(x(z)-μ)(x(z
+h)
-μ)]=σ(h)
(1-48)
22
当(1-48)成立时，X(z)与X(z+h)的方差相等，即：σ1式(1-46)=σ2=σ2，变为
σ(h)ρ(h)=2
σ
(1-49)
（二）半变异函数 用于描述区域化变量变化规律的另一个更具实用性的函数是半变异函数( Semivariogram)。半变异函数的定义为：
γ(h)=E[(x(z)-x(z+h))2]
12
(1-50)
如果满足二阶稳定性条件，半变异函数和协变异函数之间存在以下关系：
γ(h)=σ2-σ(h)
(1-51)
证明：
γ(h)=
1
E[(x(z)-x(z+h))2]21
=E[((x(z)-μ)-(x(z+h)-μ))2]2
11
=E[(x(z)-μ)2]+E[(x(z+h)-μ)2]-E[(x(z)-μ)(x(z+h)-μ)]2211
=σ12+σ22-σ(h)22=σ2-σ(h)
图1-22是关系式(1-51)的示意图。
图 当h=0时，点z和z+h变为一点，区域化变量X的取值X(z)与X(z+h)应变为同一取值。从以上各式可以看出σ(0)=σ2，γ(0)=0。实际上，在
同一位置获得两个完全相同的样品是几乎不可能的。如果我们从紧挨着的两点(h=0)取两个样品，由于取样过程中的误差和微观矿化作用的变化，两个样品不会完全相同；即使是把同一个样品化验两次，由于化验过程中的误差，化验结果也难以完全相同。因此，半变异函数在原点附近实际上不等于零，这种现象称为块金效应。块金效应的大小用块金值N表示：
N=limitγ(h)=σ2-limit[σ(h)]
h→0
h→0
(1-52)
应用半变异函数进行参数估计时，需满足内蕴假设。
[内蕴假设](Intrinsic Hypothesis)
● 区域化变量X的增量的数学期望与位置无关，即对于区域Ω内的
任意位置z0 :
E[x(z)-x(z
+h)
]=m(h)
(1-53)
● 半变异函数与位置无关，即对于区域Ω内的任意位置z0 :
1
E[(x(z0)-x(z02
+h)
)2]=γ(h)
(1-54)
内蕴假设的内涵是：区域化变量的增量，在给定区域Ω内的所有位置上具有相同的概率分布。内蕴假设要求的条件要比二阶稳定性条件宽松得多，当满足后者时，前者自然得到满足。 二、实验半变异函数及其计算
象普通随机变量的概率分布特征值一样，半变异函数对任一给定矿床Ω是未知的，需要通过取样值对之进行估计。 设从矿床Ω中获得一组样品，相距h的样品对数为n(h)，那么半变异函数γ(h)可以用下式估计：
γ(h)=
1
n(h)
2n(h)i=1
∑[xz
(
i
)
-x(zi
+h)
]2
(1-55)
式中，X(zi)是在zi处的样品值，X(zi+h)是在与zi相距h处的样品值。由(1-55)计算的半变异函数称为实验半变异函数。下面举例说明实验半变异函数的计算过程。
[例1-2]在一条直线上取得10个样品，其位置如图1-23所示，试计算实验半变异函数。
５ ５ ７
１２ １１ ８ ７ ２ ３ ３
０ １ ５ １０ １２
图1-23一维取样分布
表1-3 基于图1-23中数据的半变异函数计算结果
表1-4 h=3时γ(h)的计算过程
γ
图 1-24 实验半变异函数
解：样品是一个离散集，因此我们只能对几个离散h值计算γ(h)。应用
公式(1-55)计算结果列于表1-3中。以h=3为例，计算过程列于表1-4中。表1-3中的计算结果绘于图1-24。 上例中样品落于一直线上，是一个在一维空间计算实验半变异函数的问题。在二维或三维空间，半变异函数是具有方向性的，即在不同的方向上，半变异函数可能不一样。下面是一个二维空间下求半变异函数的算例。
[例1-3]如图1-25所示，在矿床的某一台阶取样31个，样品位于间距为1的规则网格点上，各样品的品位如图中的数字所示。试求在4个方向上的实验半变异函数。
方向
2 1
图 1-25 二维取样分布
解：在任一方向上计算过程与例1-2相同。只是在一给定方向上选取间距为h的样品对时，只能在该方向上选取。在方向1和2上的实验半变异函数计算结果列在表1-5中，在方向3和4上的计算结果列于表1-6中。
表1-5 例1-3中在方向1和2上的实验半变异函数计算结果
表1-6 例1-3中在方向3和4上的实验半变异函数计算结果
若将平面上所有方向上相距为h的样品对用于计算γ(h)，得到的实验半变异函数称为该平面上的平均实验半变异函数，例1-3中的平均实验半变异函数计算结果列于表1-7。
表1-7例1-3中平均实验半变异函数计算结果
所有4个方向上的实验半变异函数与平均半变异函数计算结果绘于图1-26。
γ
h
图 1-26 不同方向上的实验半变异函数
在实践中，样品在平面上的分布可能很不规则，不可能所有样品都位于规则的网格点上，样品间的距离也不会是一个基数的整数倍，而且往往需要计算任意方向的实验半变异函数。因此，恰好落在某一给定方向的方向线上和间距恰好等于某一给定h的样品对很少（或几乎不存在）。所以如图1-27所示，在计算实验半变异函数时，我们需要确定一个最大方向角偏差∆α和距离偏差∆h，如果一对样品X(zi)和X(zj)所在的位置所组成的向量Zi→
Zj的方向落于α－∆α和α＋∆α之间，那么就可以认为
X(zi)和X(zj)是在方向α上的一个样品对；如果样品X(zi)和X(zj)之间的距离落于h -∆h和h +∆h之间，就可认为这两个样品是相距h的一个样品对。2∆α称为窗口(window)。在实际计算中，往往以2∆h作为h的增量，以∆h作为最小h值（即偏移量Offset）。例如，当2∆h＝10米时，h取5米，15米，25米„„。
图1-27 二维半变异函数的实用计算方法
在三维空间，图1-27中的扇形变为图1-28中的锥体，空间的某一方
向由方位角φ与倾角Ψ表示。另外，在三维空间，一个样品不是一个二维点，而是具有一定长度的三维体，所以在计算半变异函数前，需要将样品进行组合处理，形成等长度的组合样品。在实验半变异函数的实际计算中，首先要对所有样品对进行矢量运算，找出落于方向与间距最大偏差范围内的样品对，然后对这些样品对应用公式(1-55)进行计算，获得半变异函数γ(h)曲线上的一个点。需要说明的是，距离h是所有这些样品对的平均距离。
三、半变异函数的数学模型
实验半变异函数由一组离散点组成，在实际应用时很不方便。因此常常将实验半变异函数拟合为一个可以用数学解析式表达的数学模型。常见的半变异函数的数学模型有以下几种：
图1-28 三维半变异函数的实用计算方法
（一）球状模型(Spherical Model) 实验半变异函数在大多数情况下可以拟合成球状模型。因此，球状模型是应用最广的一种半变异函数模型，其数学表达式为：
⎧3hh3
⎪C(-3)h≤aγ(h)=⎨2a2a
h>a⎪⎩C
(1-56)
2
2
式中，C称为槛值或台基值(sill)，一般情况下可以认为C＝σ(σ为
样品的方差），α称为变程(range)。 图1-29是球状模型的图示。从图中可以看出，γ(h)随h的增加而增加，当h达到变程时，γ(h)达到槛值C；之后γ(h)便保持常值C。这种特征的物理意义是：当样品之间的距离小于变程时，样品是相互关联的，关联程度随间距的增加而减小，或者说，变异程度随间距的增加而增大；当间距达到一定值时，样品之间的关联性消失，变为完全随机，这时γ(h)即为样品的方差。因此，变程实际上代表样品的影响范围。
（二）随机模型(Random Model) 当区域化变量X的取值是完全随机的，即样品之间的协方差σ(h)对
于所有h都等于0时，半变异函数是一常量：
γ(h)=C
(1-57)
这一模型称为随机模型，其图示为一水平直线（图1-30）。随机模型表明样品之间互不相关。随机模型有时也被称为纯块金效应模型(Pure nugget effect model)。
（三）指数模型(Exponential Model) 指数模型的数学表达式为：
γ(h)=C(1-e-h)
(1-58)
指数模型的特征与球状模型相似（图1-31），变异速率较小。式(1-58)中的α是原点处的切线达到C时的h值。
（四）高斯模型(Gaussian Model) 高斯模型的数学表达式为：
γ(h)=C(1-e-h
2
a2
)
(1-59)
如图1-32所示，高斯模型在原点的切线为水平线，表明γ(h)在短距离内变异很小。
（五）线性模型(Linear Model) 线性模型的数学表达式为一线性方程，即：
γ(h) ＝( P2/2 )h
(1-60 )
式中，p2为一常量，且
p2=E[(x(zi)-x(zi))2]
+1
(1-61)
如图1-33所示，线性模型没有槛值，γ(h)随h无限增加。
(六)对数模型(Logarithmic Model) 对数模型的表达式为：
γ(h)=3αlogeh
(1-62)
式中，α为常量。当h取对数坐标时，对数模型为一条直线（图1-34）。对数模型没有槛值。当h
γ(
hγ(
hCC
图1-29 球状模型示意图
图1-30 随机模型
γ(
hγ(hC
C
图1-31 指数模型
图1-32 高斯模型
γ(
hγ(h)
图1-33 线性模型
图1-34 对数模型
除对数模型和随机模型外，均有γ(0)=0。但由于取样、化验误差和矿化作用在短距离内（小于最小取样间距）的变化，在绝大多数情况下半变异函数在原点不等于零，即存在块金效应。因此，在实践中应用最广的模型是具有块金效应的球状模型，其数学表达式为：
3⎧3hh⎪N+C(-3)h
γ(h)=⎨2a2a
h≥a⎪⎩N+C
(1-63)
式中，N为块金效应；N+C为槛值。
在某些情况下，区域化变量的结构特性较复杂，难以用单一结构的数学模型描述。这时，往往采用几个结构的数学组合来描述较复杂的结构。带有块金效应的球模型（式1-63）实质上就是由两个结构组成的：一个是纯块金效应结构（或随机结构），另一个是球形结构。由多个半变异函数组成的结构称为嵌套结构(nested structures)。实践中较常见的嵌套结构由块金效应与两个球模型组成，即：
套结构的示意图：
(1-64) γ(h)=N+γ1(h)+γ2(h)
式中，γ1(h)和γ2(h)为具有不同α和C值的球模型。图1-35是这一嵌
γ(
hN+C1+C2N+C2N+C1N 2 1
12h
图1-35 球模型的嵌套结构示意图
四、半变异函数的拟合
实践中半变异函数是根据有限数目的地质取样建立的，而通过取样我们只能得到由一些离散点组成的实验半变异函数。因此，需要对实验半变异函数进行加工获得实验半变异函数的数学模型。将实验半变异函数加工成数学模型的过程称为半变异函数的拟合。这里只讲球模型的拟合。 图1-36是从一组样品得到的实验半变异函数。虽然数据点的分布不很规则，但仍可看出γ(h)随h首先增加，然后趋于稳定的特点。因此，其
数学模型应为具有块金效应的球模型。如果能确定块金效应N，槛值N+C和变程a，拟合也就完成了。 Y(h) 图1-36 实验半变异函数的拟合 首先确定槛值。从数据点的分布很难看出γ(h)稳定在何值，但从理论
2
上讲，可以认为槛值等于样品的方差σ。因此，在实际拟合时，往往取σ2为槛值。这里σ2＝0.135，故N+C＝0.135。 其次确定块金效应。根据槛值以下靠近原点的数据点变化趋势，作一条斜线，斜线与纵轴的截距即为块金效应N，从图中可以看出N≈0.02。这样C＝0.135-0.02=0.115。 最后确定变程。根据球模型的数学表达式可知，γ(h)在h=0处的切线
2
斜率为C/ ( a)。所以上面求块金效应时所作的斜线与等于槛值的水平
线的交点之横坐标为a。从图中可以看出，a约为100米，所以变程约为150米。 利用实际数据进行半变异函数的拟合通常是个十分复杂的过程，需要对地质特征有较好的了解和拟合经验。当取样间距较大时，变程以内的数据点很少，很难确定半变异函数在该范围内的变化趋势，而恰恰这部分曲线是半变异函数最重要的组成部分。在这种情况下，常常求助于“沿钻孔实验半变异函数”(down-hole variogram)，即沿钻孔方向建立的实验半变异函数。因为沿钻孔取样间距小，沿钻孔半变异函数可以捕捉短距离内的结构特征，帮助确定半变异函数的块金效应和变化趋势。但必须注意，当存在各向异性时，沿钻孔半变异函数只代表区域化变量沿钻孔方向的变化特征，并不能完全代表其他方向上半变异函数在短距离的变化特征。
第五节 品位、矿量计算的垂直断面法
垂直断面法是传统的手工计算矿量的常用方法，其一般步骤为： 第一步：沿勘探线做垂直剖面，将勘探线上的钻孔及其取样品位标在剖面图上（图1-11）。
图 1-11 标有取样品位的剖面图
第二步：根据给定的边界品位进行矿体圈定。简单地讲，矿体圈定的过
程就是将相邻钻孔上高于边界品位的样品点相连的过程。当一条矿体被一个钻孔穿越，而在相邻的钻孔消失时，一般将矿体延伸到两钻孔的中点，或是根据矿体的自然尖灭趋势在两钻孔之间实行自然尖灭。在矿体圈定过程中，要充分考虑矿床的地质构造（如断层和岩性）和成矿规律。图1-12是当边界品位等于25％时根据图1-11中的取样品位圈定的矿体示意图。
第三步：矿体圈定完成后，可用求积仪求得每个断面上的矿石面积，然
后就可以进行矿量计算。
(a)当一条矿体在两个相邻断面上的面积（S1和S2）相差不到40
％时，两断面之间的矿体体积用下式计算：
V=
s1+s2
2
L
(1-31)
式中，L为断面间距。
(b)当两个相邻断面上的面积相差大于40％时，采用下式计算：
V=
s1+s2+1⋅s2
3
L
(1-32)
(c)当矿体在二断面间是楔形尖灭时，计算公式为：
V=
sL 2
计算出两断面间矿石块段体积后，矿石块段的矿量为：
T=V⋅γ
(1-33)
式中，γ为本块段矿石体重。然后将所有块段的矿量相加，即得矿床的总矿量。
图 1-12 边界品位为25％时矿体圈定示意图
第四步：计算矿体的平均品位： (a)对穿越矿体的每一钻孔的样品进行“矿段样品组合”，求出
组合样品的品位。
(b)求出每一组合样品的影响面积。该面积是以钻孔为中线向两
侧各外推二分之一钻孔间距得到的矿体面积。
(c)对组合样品品位以其影响面积为权值进行加权平均计算，求
出矿体在断面上的平均品位。
(d)一条矿体的总平均品位是该条矿体在各断面上的平均品位以
断面所代表的矿量为权值的加权平均值。
上述矿体圈定与矿量、品位计算过程，可通过编制程序由计算机完成，但矿体圈定的完全计算机化具有较高的难度。垂直断面法是一种传统的手工方法，在我国仍广泛采用，但在发达国家由于计算机的广泛应用已基本成为过去。
第六节 矿量、品位计算的水平断面法
在露天矿山，矿石的开采是分台阶进行的，因此用于矿量、品位计算的一个水平断面即为一个台阶。常用的水平断面法有多边形法和三角形法。 600N
300N
300E
600E
图 1-13 台阶水平面上钻孔取样分布示意图
一、多边形法(polygonal method)
在多边形法中，水平断面上的每一钻孔取样（即进行台阶样品组合后的组合样品），位于一个多边形的中心。多边形的形成由以下步骤完成：
第一步：把穿越水平面的钻孔根据钻孔坐标，绘于水平面上，并将本平
面的组合样品品位标注在图上（图1-13）。
第二步：根据经验和地质统计学分析，确定影响半径R，R的确定将在后
面的地质统计学部分介绍。
第三步：以每一样品为中心，确定其相邻样品，一般情况下相邻样品是
落在半径为2R的圈内的样品。以样品S-41为例，其相邻样品如图1-14a所示。
第四步：用直线将中心样品和相邻样品连接起来（图1-14a）。
第五步：在中心样品与每一相邻样品的连线中点作垂直于连线的直线
（称为二分线），这些二分线相交围成的多边形即为所求的多边形。当两条二分线近于平行时，在二者相交之前，将与另外一条二分线相交，这时，取二者中离中心样品最近者作为多边形的边。以样品S-41为中心的多边形如图1-14b所示。
( b ) ( a )
图 1-14 中心样品四周均有样品时多边形的形成过程
如果在某些区域，钻孔间距大于2R，就以中心样品为中心，以R为半径做一八边形。位于边缘上的样品，只在其一侧有相邻样品而在另一侧没有相邻样品。这时，在没有相邻样品的一侧以中心样品为中心，画一半径为R的圆。然后在0o,90o与±45o方向上分别作圆的切线，这些切线与有相邻样品一侧的二分线相交就形成了所求的多边形。图1-15是以边缘样品S－14为例的多边形形成过程。
( a )
( b )
图 1-15 边缘样品的多边形的形成过程
以每一样品为中心形成多边形后，在每一多边形内，品位被看作是常数且多边形的品位等于其中心样品的品位。第i个多边形的重量Ti由下式计算：
Ti=SiγH
(1-34)
式中，Si为第i个多边形的面积，γ为体重；H为台阶高度。 如果多边形的品位xi大于边界品位，该多边形即为矿石多边形，所有矿石多边形的集合就形成了该水平面上的矿体。基于图1-13中的样品分布，应用多边形法进行品位计算和矿体圈定的结果如图1-16所示。这里假设边界品位为0.6％。 水平断面上的矿石总量T为矿石多边形的重量之和，矿石的平均品位为矿石多边形品位的面积加权平均值，即：
T=∑Ti
i=1
n
(1-35)
=∑xisi
i=1
n
∑si
i=1
n
(1-36)
式中，xi为矿石多边形i的品位，n为矿石多边形个数。 在某些矿床中，往往不同区域的成矿作用与矿石种类不同，或是地质构造使矿体发生错动。在这种情况下，使用多边形法时，应注意多边形不跨越区域界限和构造线。
图1-16 多边形法品位估算与矿体圈定结果
图1-16 多边形品位估算与矿体圈定结果
二、三角形法
在多边形法中，中心样品的品位被外推到整个多边形，即一个多边形的平均品位只用了一个已知样品作估计。一般情况下，对一给定体积的品位进行估算时，利用的已知数据越多，估计误差越小。从而就产生了三角形法。 在三角形法中，每一样品是三角形的一个顶点。图1-17是将图1-13中的样品相连形成的三角形。很显然，在一定条件下顶点连法不同，所得到的三角形也不同。一般性的连接原则是使每个三角形的边尽可能短，面积尽可能小。 每一三角形的品位xi是位于其顶点上的三个样品品位的平均值，即：
xi=
xi1+xi2+xi3
3
(1-37)
品位高于边界品位的三角形是矿石三角形，水平面上所有矿石三角形的集合形成该水平面上的矿体。一个水平面上的总矿量为矿石三角形的重量之和，矿石的平均品位为矿石三角形品位的重量或面积加权平均值。 600N
300N
300E
图 1-17 三角形法品位估算与矿体圈定结果
三角形法的优点在于它利用三个样品的品位来估计一个三角形的平均品位。从理论上讲，利用三角形法求得的品位、矿量较多边形法误差小。在应用三角形法时，也要注意三角形不超越区域界限和地质构造线。 第七节 三维块状模型
三维块状模型是将矿床划分为许多单元块形成的离散模型（图1-18）。单元块一般是尺寸相等的长方体。随着计算机在矿山的普及应用和计算机的容量与速度的不断提高，三维块状模型在国际上得到越来越广泛的应用。三维块状模型不仅被广泛用于品位、矿量计算，也被用于
图1-18 三维块段模型示意图
露天矿最终开采境界优化和开采计划优化。实际上，许多优化方法是由于三维块状模型的引入而产生的，而且当今矿山优化界正在进行的各种优化方法的研究，绝大多数以三维块状模型为基础。在我国，虽然三维块状模型的概念引入较早，但由于观念的陈旧与国家有关储量计算和矿山设计规范的束缚等原因，三维块状模型除在矿山设计院有少量应用外，至今没有发挥应有的作用。 三维块状模型中单元块的高度等于露天矿台阶高度，单元块在水平方向一般取正方形，其边长视具体情况而定。有人错误地认为，单元块越小，品位、矿量计算结果越精确。但是，从地质统计学理论可知，在已知数据（即样品数与样品在空间的分布）一定的条件下，单元块越小，对其品位的估计误差越大。另外，当单元块小到一定程度时，相邻的几个单元块的品位估计值会非常接近，与它们的平均品位相差无几，这种现象称为平滑作用，故用一个大块代替几个小块，品位与矿量的计算结果变化很小。而这样做可以降低对计算机容量的要求，加快计算速度。一般的经验规则是，单元块在水平方向上的边长不应小于钻孔平均间距的1/4或1/5。对于100米的钻孔间距，单元块的边长一般取30米左右。 将矿床分为单元块后，需要应用某种方法对每一小块的平均品位进行估计。常用的方法有三，即最近样品法、距离N次方反比法和地质统计学法（即克里金法）。三者均基于样品加权平均的概念，即对落在以单元块为中心的影响范围内的样品品位进行加权平均求得单元块的品位。三种方法的根本区别在于所用权值不同。本节介绍前两种方法，第三种方法将在下节给予较详细的论述。 一、最近样品法
所谓最近样品法就是将距离某一单元块最近的样品品位作为该单元块的品位估计值。前面介绍的多边形法其实是不规则单元块情况下的最近样品法，最近样品法的一般步骤为：
第一步：以被估计的单元块的中心为圆心，做半径为影响半径R的圆。 第二步：计算落入影响范围内的每一样品与单元块中心点的距离。
第三步：选取离单元块中心最近的样品，其品位即为被估单元块的品位。 当没有样品落入影响范围时，被估单元块的品位是未知的。一般情况下，未知单元块的品位取0值，当废石处理。如果有理由认为未知单元块所处的区域可能有矿石，未知单元块的出现表明该区域的数据量太少，要想确定其品位与矿量，需要增加钻孔。 在三维状态下，上述的影响范围，由二维状态下的圆变为球，需要对落入球内的所有样品进行考查，找出离单元块中心最近的样品。 求得矿床中所有单元块的品位以后，品位大于边界品位的单元块的集合组成矿体。矿石量和矿石平均品位可由矿石单元块重量的简单累加
和品位的平均求得。基于图1-13所示的样品分布和品位值，用最近样品法求得该台阶上单元块的品位如图1-19所示。这里采用的影响半径(R)为75m，边界品位(gc)为0.6％。将图1-19与图1-16作比较可见，用多边形法与最近样品法所得的矿体形态相近。
600E
图1-19 最近样品法品位估算与矿体圈定结果
二、距离N次方反比法(Inverse Distance Method)
在多边形法和最近样品法中，只有一个样品参与单元块品位的估值，如果落入影响范围的样品都参与单元块的品位估值，估值结果会更为精确。然而，由于各样品距单元块中心的距离不同，其品位对单元体的影响程度也不同。显然，距离单元块越近的样品，其品位对单元体品位的影响也就越大。因而在计算中，离单元体近的样品的权值应比离单元体远的样品的权值大。距离N次方反比法就是基于这一思想产生的。在此法中，一个样品的权值等于样品到单元块中心距离的N次方的倒数（1/dN）。 参照图1-20，距离N次方反比法的一般步骤如下：
第一步：以被估单元块中心为圆心、以影响半径R为半径做圆，确定影
响范围（在三维状态下，圆变为球）。
第二步：计算落入影响范围内每一样品与被估单元块中心的距离。 第三步：利用下式计算单元块的品位Xb：
xi
x=∑N b
i=1di
n
∑dN
i=1
i
n
1
(1-38)
式中，xi为落入影响范围的第i个样品的品位；di为第i个样品到单元块中心的距离。
图1-20 距离N次方反比法示意
在实际应用中，有时采用所谓的角度排除，即当一个样品与被估单元块中心的连线与另一个样品与被估单元块中心的连线之间的夹角小于某一给定值α时，距单元块较远的样品将不参与单元块的估值运算（如图1-20中的G3与G5）。α值一般在15度左右。如果没有样品落入影响范围之内，单元块的品位为零。公式(1-38)中的指数N对于不同的矿床取值不同。假设有两个矿床，第一个矿床的品位变化程度较第二个矿床的品位变化程度大，即第二个矿床的品位较第一个矿床连续性好。那么在离单元体同等距离的条件下，第一个矿床中样品对单元块品位的影响应比第二个矿床小。因此，在估算某一单元块的品位时第一个矿床中样品的权值在同等距离条件下应比第二个矿床中样品的权值小。也就是说在品位变化小的矿床，N取值较小；在品位变化大的矿床，N取值较大。在铁、镁等品位变化较小的矿床中，N一般取2；在贵重金属（如黄金）矿床中，
N的取值一般大于2，有时高达4或5。如果有区域异性存在，不同区域中品位的变化不同，则需要在不同区域取不同的N值。同时，一个区域的样品一般不参与另一区域的单元块品位的估值运算。以图1-20为例，若n=2，则被估单元块的品位为0.628。 600N
300N
300E 600E
图1-21 距离平方反比法品位计算与矿体圈定结果
基于图1-13所示的样品分布和品位值，应用上述方法，当N＝2，R=75m，gc=0.6％时，品位计算和矿体圈定结果如图1-21所示。将图1-16，1-19和1-21所示的计算结果作比较，应用距离平方反比法求得的矿体形态与前两种方法（多边形法与最近样品法）得出的矿体形态之间的差别较为明显。
第八节 地质统计学法（Geostatistical Method）
地质统计学是60年代初期出现的一个新兴应用数学分支，其基本思想是由南非的Danie Krige在金矿的品位估算实践中提出来的，后来由法国的Georges Matheron经过数学加工，形成一套完整的理论体系。在过去的三十多年中，地质统计学不仅在理论上得到发展与完善，而且在
实践中得到日益广泛的应用。如今，地质统计学在国际上除被用于矿床的品位估算外，也被用于其他领域中研究与位置有关的参数变化规律和参数估计。如农业中农作物的收成、环保中污染物的分布等等。本节将从矿床的品位估算的角度，简要介绍地质统计学的基本概念、原理和方法。
一、区域化变量、协变异函数与半变异函数
应用传统统计学（“传统”二字是相对于地质统计学而言的）可以对矿床的取样数据进行各种分析，并估计矿床的平均品位及其置信区间。在给定边界品位时，传统统计学也可用于初步估算矿石量和矿石平均品位。然而，传统统计学的分析计算均基于一个假设，即样品是从一个未知的样品空间随机选取的，而且是相互独立的。根据这一假设，样品在矿床中的空间位置是无关紧要的，从相隔上千米的矿床两端获取的两个样品与从相隔几米的两点获取的两个样品从理论上讲是没有区别的，它们都是一个样本空间的两个随机取样而已。 但是在实践中，相互独立性是几乎不存在的，钻孔的位置（即样品的选取）在绝大多数情况下也不是随机的。当两个样品在空间的距离很小时，样品间会存在较强的相似性，而当距离很大时，相似性就会减弱或不存在。也就是说，样品之间存在着某种联系，这种联系的强弱是与样品的相对位置有关的。这样就引出了“区域化变量”的概念。
（一）区域化变量及协变异函数： 如果以空间一点为中心获取一样品，样品的特征值X(z)是该点的空间位置z的函数，那么变量X即为一区域化变量。 显然，矿床的品位是一个区域化变量，而控制这一区域化变量之变化规律的是地质构造和矿化作用。区域化变量的概念是整个地质统计学理论体系的核心，用于描述区域化变量变化规律的基本函数是协变异函数和半变异函数。 设有两个随机变量X1与X2，如果X1与X2之间存在某种相关性，那么从传统统计学可知，这种相关关系由X1与X2的协方差σ(x1,x2)表示：
σ(x1,x2)=E[(x1-E(x1))(x2-E(x2))] (1-39)
22
让σx和分别表示X1和X2的方差，则： σx12
2
σx=E[(x1-E(x1))2]
1
(1-40) (1-41)
2σx=E[(x2-E(x2))2]
2
式中，E[*]表示随机变量[*]的数学期望。 X1与X2之间的相关系数为：
ρx⋅x=
1
2
σ(x1,x2)σx⋅σx
1
2
(1-42)
当X1和X2互相独立时，即二者之间不存在任何相关性时，协方差与相关系数均为零。当X1和X2“完全相关”时，相关系数为1.0（或-1.0）。 如果X1和X2不是一般的随机变量，而是区域化变量X在矿体Ω中的取值，即： X1代表 X(z)：区域化变量X在矿体Ω中z点的取值； X2代表X(z+h)：区域化变量X在矿体Ω中距z点h处的取值。
那么，由式(1-39)可计算X(z)与X(z+h)在矿体Ω中的协方差：
σ(x(z),x(z+h))=σ(h)
](x(z+h)-E[x(z+h))]]=E[(x(z)-E[x(z))
(1-43)
σ(h)称为区域化变量X在Ω中的协变异函数(Covariogram)。
22
让σ1和σ2分别表示X(z)与X(z+h) 在矿体Ω中的方差,则：
2
σ1=E[(x(z)-E[x(z)])2] 2σ2=E[(x(z+h)-E[x(z+h)])2]
(1-44) (1-45)
那么X(z)与X(z+h)之间的相关系数为：
σ(h)ρ(h)=
σ1σ2
(1-46)
ρ(h)称为区域化变量X在Ω中的相关函数(Correlogram)。 对于任何矿床，都可能计算出其协变异函数σ(h)，但在利用σ(h)对
矿床中单元体的品位进行估值时，需满足二阶稳定性条件。
[二阶稳定性条件](Second order stationarity conditions):
● X(z)的数学期望与空间位置z无关，即对任意位置z0 ：
E[x(z0)]= μ
(1-47)
● 协变异函数与空间位置无关只与距离h有关，即对于任何位置
z0 ：
E[(x(z)-μ)(x(z
+h)
-μ)]=σ(h)
(1-48)
22
当(1-48)成立时，X(z)与X(z+h)的方差相等，即：σ1式(1-46)=σ2=σ2，变为
σ(h)ρ(h)=2
σ
(1-49)
（二）半变异函数 用于描述区域化变量变化规律的另一个更具实用性的函数是半变异函数( Semivariogram)。半变异函数的定义为：
γ(h)=E[(x(z)-x(z+h))2]
12
(1-50)
如果满足二阶稳定性条件，半变异函数和协变异函数之间存在以下关系：
γ(h)=σ2-σ(h)
(1-51)
证明：
γ(h)=
1
E[(x(z)-x(z+h))2]21
=E[((x(z)-μ)-(x(z+h)-μ))2]2
11
=E[(x(z)-μ)2]+E[(x(z+h)-μ)2]-E[(x(z)-μ)(x(z+h)-μ)]2211
=σ12+σ22-σ(h)22=σ2-σ(h)
图1-22是关系式(1-51)的示意图。
图 当h=0时，点z和z+h变为一点，区域化变量X的取值X(z)与X(z+h)应变为同一取值。从以上各式可以看出σ(0)=σ2，γ(0)=0。实际上，在
同一位置获得两个完全相同的样品是几乎不可能的。如果我们从紧挨着的两点(h=0)取两个样品，由于取样过程中的误差和微观矿化作用的变化，两个样品不会完全相同；即使是把同一个样品化验两次，由于化验过程中的误差，化验结果也难以完全相同。因此，半变异函数在原点附近实际上不等于零，这种现象称为块金效应。块金效应的大小用块金值N表示：
N=limitγ(h)=σ2-limit[σ(h)]
h→0
h→0
(1-52)
应用半变异函数进行参数估计时，需满足内蕴假设。
[内蕴假设](Intrinsic Hypothesis)
● 区域化变量X的增量的数学期望与位置无关，即对于区域Ω内的
任意位置z0 :
E[x(z)-x(z
+h)
]=m(h)
(1-53)
● 半变异函数与位置无关，即对于区域Ω内的任意位置z0 :
1
E[(x(z0)-x(z02
+h)
)2]=γ(h)
(1-54)
内蕴假设的内涵是：区域化变量的增量，在给定区域Ω内的所有位置上具有相同的概率分布。内蕴假设要求的条件要比二阶稳定性条件宽松得多，当满足后者时，前者自然得到满足。 二、实验半变异函数及其计算
象普通随机变量的概率分布特征值一样，半变异函数对任一给定矿床Ω是未知的，需要通过取样值对之进行估计。 设从矿床Ω中获得一组样品，相距h的样品对数为n(h)，那么半变异函数γ(h)可以用下式估计：
γ(h)=
1
n(h)
2n(h)i=1
∑[xz
(
i
)
-x(zi
+h)
]2
(1-55)
式中，X(zi)是在zi处的样品值，X(zi+h)是在与zi相距h处的样品值。由(1-55)计算的半变异函数称为实验半变异函数。下面举例说明实验半变异函数的计算过程。
[例1-2]在一条直线上取得10个样品，其位置如图1-23所示，试计算实验半变异函数。
５ ５ ７
１２ １１ ８ ７ ２ ３ ３
０ １ ５ １０ １２
图1-23一维取样分布
表1-3 基于图1-23中数据的半变异函数计算结果
表1-4 h=3时γ(h)的计算过程
γ
图 1-24 实验半变异函数
解：样品是一个离散集，因此我们只能对几个离散h值计算γ(h)。应用
公式(1-55)计算结果列于表1-3中。以h=3为例，计算过程列于表1-4中。表1-3中的计算结果绘于图1-24。 上例中样品落于一直线上，是一个在一维空间计算实验半变异函数的问题。在二维或三维空间，半变异函数是具有方向性的，即在不同的方向上，半变异函数可能不一样。下面是一个二维空间下求半变异函数的算例。
[例1-3]如图1-25所示，在矿床的某一台阶取样31个，样品位于间距为1的规则网格点上，各样品的品位如图中的数字所示。试求在4个方向上的实验半变异函数。
方向
2 1
图 1-25 二维取样分布
解：在任一方向上计算过程与例1-2相同。只是在一给定方向上选取间距为h的样品对时，只能在该方向上选取。在方向1和2上的实验半变异函数计算结果列在表1-5中，在方向3和4上的计算结果列于表1-6中。
表1-5 例1-3中在方向1和2上的实验半变异函数计算结果
表1-6 例1-3中在方向3和4上的实验半变异函数计算结果
若将平面上所有方向上相距为h的样品对用于计算γ(h)，得到的实验半变异函数称为该平面上的平均实验半变异函数，例1-3中的平均实验半变异函数计算结果列于表1-7。
表1-7例1-3中平均实验半变异函数计算结果
所有4个方向上的实验半变异函数与平均半变异函数计算结果绘于图1-26。
γ
h
图 1-26 不同方向上的实验半变异函数
在实践中，样品在平面上的分布可能很不规则，不可能所有样品都位于规则的网格点上，样品间的距离也不会是一个基数的整数倍，而且往往需要计算任意方向的实验半变异函数。因此，恰好落在某一给定方向的方向线上和间距恰好等于某一给定h的样品对很少（或几乎不存在）。所以如图1-27所示，在计算实验半变异函数时，我们需要确定一个最大方向角偏差∆α和距离偏差∆h，如果一对样品X(zi)和X(zj)所在的位置所组成的向量Zi→
Zj的方向落于α－∆α和α＋∆α之间，那么就可以认为
X(zi)和X(zj)是在方向α上的一个样品对；如果样品X(zi)和X(zj)之间的距离落于h -∆h和h +∆h之间，就可认为这两个样品是相距h的一个样品对。2∆α称为窗口(window)。在实际计算中，往往以2∆h作为h的增量，以∆h作为最小h值（即偏移量Offset）。例如，当2∆h＝10米时，h取5米，15米，25米„„。
图1-27 二维半变异函数的实用计算方法
在三维空间，图1-27中的扇形变为图1-28中的锥体，空间的某一方
向由方位角φ与倾角Ψ表示。另外，在三维空间，一个样品不是一个二维点，而是具有一定长度的三维体，所以在计算半变异函数前，需要将样品进行组合处理，形成等长度的组合样品。在实验半变异函数的实际计算中，首先要对所有样品对进行矢量运算，找出落于方向与间距最大偏差范围内的样品对，然后对这些样品对应用公式(1-55)进行计算，获得半变异函数γ(h)曲线上的一个点。需要说明的是，距离h是所有这些样品对的平均距离。
三、半变异函数的数学模型
实验半变异函数由一组离散点组成，在实际应用时很不方便。因此常常将实验半变异函数拟合为一个可以用数学解析式表达的数学模型。常见的半变异函数的数学模型有以下几种：
图1-28 三维半变异函数的实用计算方法
（一）球状模型(Spherical Model) 实验半变异函数在大多数情况下可以拟合成球状模型。因此，球状模型是应用最广的一种半变异函数模型，其数学表达式为：
⎧3hh3
⎪C(-3)h≤aγ(h)=⎨2a2a
h>a⎪⎩C
(1-56)
2
2
式中，C称为槛值或台基值(sill)，一般情况下可以认为C＝σ(σ为
样品的方差），α称为变程(range)。 图1-29是球状模型的图示。从图中可以看出，γ(h)随h的增加而增加，当h达到变程时，γ(h)达到槛值C；之后γ(h)便保持常值C。这种特征的物理意义是：当样品之间的距离小于变程时，样品是相互关联的，关联程度随间距的增加而减小，或者说，变异程度随间距的增加而增大；当间距达到一定值时，样品之间的关联性消失，变为完全随机，这时γ(h)即为样品的方差。因此，变程实际上代表样品的影响范围。
（二）随机模型(Random Model) 当区域化变量X的取值是完全随机的，即样品之间的协方差σ(h)对
于所有h都等于0时，半变异函数是一常量：
γ(h)=C
(1-57)
这一模型称为随机模型，其图示为一水平直线（图1-30）。随机模型表明样品之间互不相关。随机模型有时也被称为纯块金效应模型(Pure nugget effect model)。
（三）指数模型(Exponential Model) 指数模型的数学表达式为：
γ(h)=C(1-e-h)
(1-58)
指数模型的特征与球状模型相似（图1-31），变异速率较小。式(1-58)中的α是原点处的切线达到C时的h值。
（四）高斯模型(Gaussian Model) 高斯模型的数学表达式为：
γ(h)=C(1-e-h
2
a2
)
(1-59)
如图1-32所示，高斯模型在原点的切线为水平线，表明γ(h)在短距离内变异很小。
（五）线性模型(Linear Model) 线性模型的数学表达式为一线性方程，即：
γ(h) ＝( P2/2 )h
(1-60 )
式中，p2为一常量，且
p2=E[(x(zi)-x(zi))2]
+1
(1-61)
如图1-33所示，线性模型没有槛值，γ(h)随h无限增加。
(六)对数模型(Logarithmic Model) 对数模型的表达式为：
γ(h)=3αlogeh
(1-62)
式中，α为常量。当h取对数坐标时，对数模型为一条直线（图1-34）。对数模型没有槛值。当h
γ(
hγ(
hCC
图1-29 球状模型示意图
图1-30 随机模型
γ(
hγ(hC
C
图1-31 指数模型
图1-32 高斯模型
γ(
hγ(h)
图1-33 线性模型
图1-34 对数模型
除对数模型和随机模型外，均有γ(0)=0。但由于取样、化验误差和矿化作用在短距离内（小于最小取样间距）的变化，在绝大多数情况下半变异函数在原点不等于零，即存在块金效应。因此，在实践中应用最广的模型是具有块金效应的球状模型，其数学表达式为：
3⎧3hh⎪N+C(-3)h
γ(h)=⎨2a2a
h≥a⎪⎩N+C
(1-63)
式中，N为块金效应；N+C为槛值。
在某些情况下，区域化变量的结构特性较复杂，难以用单一结构的数学模型描述。这时，往往采用几个结构的数学组合来描述较复杂的结构。带有块金效应的球模型（式1-63）实质上就是由两个结构组成的：一个是纯块金效应结构（或随机结构），另一个是球形结构。由多个半变异函数组成的结构称为嵌套结构(nested structures)。实践中较常见的嵌套结构由块金效应与两个球模型组成，即：
套结构的示意图：
(1-64) γ(h)=N+γ1(h)+γ2(h)
式中，γ1(h)和γ2(h)为具有不同α和C值的球模型。图1-35是这一嵌
γ(
hN+C1+C2N+C2N+C1N 2 1
12h
图1-35 球模型的嵌套结构示意图
四、半变异函数的拟合
实践中半变异函数是根据有限数目的地质取样建立的，而通过取样我们只能得到由一些离散点组成的实验半变异函数。因此，需要对实验半变异函数进行加工获得实验半变异函数的数学模型。将实验半变异函数加工成数学模型的过程称为半变异函数的拟合。这里只讲球模型的拟合。 图1-36是从一组样品得到的实验半变异函数。虽然数据点的分布不很规则，但仍可看出γ(h)随h首先增加，然后趋于稳定的特点。因此，其
数学模型应为具有块金效应的球模型。如果能确定块金效应N，槛值N+C和变程a，拟合也就完成了。 Y(h) 图1-36 实验半变异函数的拟合 首先确定槛值。从数据点的分布很难看出γ(h)稳定在何值，但从理论
2
上讲，可以认为槛值等于样品的方差σ。因此，在实际拟合时，往往取σ2为槛值。这里σ2＝0.135，故N+C＝0.135。 其次确定块金效应。根据槛值以下靠近原点的数据点变化趋势，作一条斜线，斜线与纵轴的截距即为块金效应N，从图中可以看出N≈0.02。这样C＝0.135-0.02=0.115。 最后确定变程。根据球模型的数学表达式可知，γ(h)在h=0处的切线
2
斜率为C/ ( a)。所以上面求块金效应时所作的斜线与等于槛值的水平
线的交点之横坐标为a。从图中可以看出，a约为100米，所以变程约为150米。 利用实际数据进行半变异函数的拟合通常是个十分复杂的过程，需要对地质特征有较好的了解和拟合经验。当取样间距较大时，变程以内的数据点很少，很难确定半变异函数在该范围内的变化趋势，而恰恰这部分曲线是半变异函数最重要的组成部分。在这种情况下，常常求助于“沿钻孔实验半变异函数”(down-hole variogram)，即沿钻孔方向建立的实验半变异函数。因为沿钻孔取样间距小，沿钻孔半变异函数可以捕捉短距离内的结构特征，帮助确定半变异函数的块金效应和变化趋势。但必须注意，当存在各向异性时，沿钻孔半变异函数只代表区域化变量沿钻孔方向的变化特征，并不能完全代表其他方向上半变异函数在短距离的变化特征。