矿化结构分析的变异椭球体构造及应用

毛先成，曹芳，邹艳红，赵莹

(中南大学 有色金属成矿预测教育部重点实验室，地球科学与信息物理学院，湖南 长沙，410083)

摘要：针对目前矿化结构分析中变异椭球体的构造仅考虑地质体走向、倾向、垂向上矿化结构的问题，提出一种以椭球体模型为基础，以变异函数为工具，综合考虑矿体的形态产状和品位的三维空间变异特征来构造变异椭球体的新方法。该方法通过计算铅、锌品位在空间中各个方位上的实验变异函数，绘制矿体走向倾斜面上的实验变异函数等值线图，研究矿体走向倾斜面上的铅、锌品位分布，确定变异椭球体的长轴方向，进而绘制与长轴相垂直平面上的实验变异函数等值线图，确定变异椭球体的次长轴和短轴方向；然后，分别用理论模型拟合这3个方向上的实验变异函数，得到其变程作为三轴的长度。将该方法应用于青海虎头崖矿床，构造矿床铅、锌品位的变异椭球体，并将其作为搜索椭球体来限制品位插值过程中已知样品点的搜索范围。研究结果表明：采用该方法构造的变异椭球体，不仅考虑了矿体形态对矿化空间分布的影响，而且有效地反映了铅、锌品位在空间中的分布及变异情况，为隐伏矿体定量预测提供了依据。

关键词：地质统计学；变异函数；变异椭球体; 倾伏向；倾伏角；旋转矩阵

中图分类号：P624             文献标志码：A         文章编号：1672-7207(2014)09-3076-09

Construction and application of variation ellipsoid for analysis of mineralization structure

MAO Xiancheng, CAO Fang, ZOU Yanhong, ZHAO Ying

(Key Laboratory of Metallogenic Prediction of Nonferrous Metals，Ministry of Education,

School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

Abstract: As the construction of variation ellipsoid only considered the mineralization structure in the trend, inclination, and vertical direction of geologic bodies in analysis of mineralization structure, a new method based on ellipsoid model was presented to construct variation ellipsoid by considering the shape and occurrence of ore-bodies and 3D spatial variation characteristics of ore grade comprehensively with the tool of variogram. Experimental variograms of grade variables in all directions in space were calculated and the contour map of experimental variogram on the inclined surface along the trend of ore-bodies was plotted, so grade distribution on the inclined surface along the trend of ore-bodies could be investigated and the long axis direction of variation ellipsoid was defined. Then the contour map of experimental variogram on the inclined surface perpendicular to the long axis was plotted, and so the secondary long axis direction, the short axis direction of variation ellipsoid could also be defined. Theoretical models were used to fit experimental variograms in the directions of the three axes, and their variation ranges were selected as lengths of the three axes. Finally, this method were applied in Hutouya ore deposit, Qinghai Province, and variation ellipsoids of Pb and Zn grades were constructed which were also taken as search ellipsoid to limit the search extent of known sample points in process of grade spatial interpolation. The results show that the variation ellipsoids constructed with this method not only consider the impact of the ore body morphology on mineralization spatial distribution, but also reflect spatial distribution and variability of grade variable effectively, which provides a basis for quantitative prediction of concealed ore bodies.

Key words: geostatistics; variogram; variation ellipsoid; pitch direction; pitch angle; rotation matrix

矿化是指地质体中供开采利用的物质含量超过地球化学背景值、有相对集中浓集分布且达到一定丰度的地质现象。研究矿体矿化特征、分布和变化对成矿信息定量提取及隐伏矿体定量预测有着极其重要的意义。随着20个世纪60年代初地质统计学的提出，利用地质统计学方法进行矿化结构分析便逐渐成为热点，国内外学者对此进行了很多研究，包括理论研究^[1-4]、方法技术研究^[5-8]、三维可视化研究^[9-11]等，并成功地应用于各个矿区^[12-16]。然而，这些方法更多注重地质统计学理论及各种方法技术的研究，在构造变异椭球体的过程中只考虑了矿体走向、倾向和垂向方向上地质的空间结构性，未全面考虑矿化在三维空间中任意方向的空间变异特征，有可能不能真实反映矿化的空间结构特征。随着计算机技术的发展，一批大型三维矿业软件被成功研发^[17-20]，使同时计算三维空间中各个方向上的实验变异函数变得很容易，也使得在三维空间中全面研究矿化结构成为可能。为了能够揭示矿化的三维空间结构特征，为成矿信息定量提取和隐伏矿体立体定量预测提供矿化结构模型支撑，本文作者提出一种以椭球体模型为基础，以变异函数为工具，综合考虑矿体的形态产状和品位的三维空间变异特征来构造变异椭球体的新方法。该方法首先通过地质分析或三维矿业软件建模获得矿体的形态和产状(走向、倾向、垂向)；然后，通过计算品位在空间中各个方位上的实验变异函数，绘制出矿体走向倾斜面上的实验变异函数等值线图，研究矿体走向倾斜面上的品位分布，确定变异椭球体的长轴方向；绘制出与长轴相垂直平面上的铅、锌品位的实验变异函数等值线图，确定变异椭球体的次长轴和短轴方向；最后，分别用理论模型拟合这3个方向上的品位的实验变异函数，得到其变程作为三轴的长度，即构造出变异椭球体模型，并将所得出的变异椭球体模型构造方法应用于青海虎头崖铅、锌矿床的矿化结构分析。

1  矿化结构分析的变异椭球体模型

由于矿化现象在空间中既具有随机性变化又具有结构性变化^[21-22]，本文以实验变异函数为工具进行矿化空间分析，获得反映矿体结构化信息的椭球体来反映矿体矿化的形态产状和空间展布规律。椭球体的各个轴在空间中的展布长度各不相同，并且从最长轴到最短轴是逐渐变化的，其数学模型如图1所示，其表面的数学表达式为

              (1)

式中：x，y和z为坐标变量；a，b和c为椭球体3个主要半轴的长度。

图1  椭球体数学模型

Fig. 1  Mathematical model of ellipsoid

变异椭球体是指定性描述区域变量在空间中各个方向上变异情况的一种椭球体模型，反映了区域变量的空间变异程度和各个方向上的自相关距离，模拟区域变量的分布情况。变异椭球体的轴尺寸表示区域变量沿此轴方向上的自相关距离，其中沿长轴方向区域变量变化最慢，具有较大的空间自相关距离，而沿短轴方向区域变量变化最快，空间自相关距离最小。

按照式(1)运用椭球体表面的数学模型来求取变异椭球体时，椭球体必须与坐标系处于一种相对水平的状态，即椭球体的3个主要轴分别与坐标系的X，Y和Z轴重合。而空间中反映区域变量变异情况的变异椭球体不一定如此，如图2所示，故需对坐标系进行旋转，使得坐标系的三轴与变异椭球体的3个主要轴重合。

矿体品位在不同的方向上具有不同的变异性^[23]，其空间分布也与矿体产状密切相关。在矿化空间分析过程中，在充分分析矿体产状的基础上，以变异函数为工具，构造品位的变异椭球体，并进行可视化建模，既模拟了各个方向上品位的变异情况反映铅、锌品位变异性质和变异程度，又形象而生动地展示了铅、锌矿体品位在空间中的分布情况。

图2  空间中的变异椭球体

Fig. 2  Variation ellipsoid in space

品位的变异椭球体反映了铅、锌品位在空间中各个方向上的自相关性，因此，在铅、锌品位空间插值过程中，可将变异椭球体作为搜索椭球体来限制空间中已知样品点的搜索范围，保证与待估点相关的已知样品点都位于该椭球体内，并在此椭球体内采用空间插值方法(如克立格法、反距离权重法等)对未知点品位进行估值，进而得到资源量估算结果。通过与矿山开采获得的铅、锌矿体平均品位和资源量进行对比可检验变异椭球体构造的合理性。

2  变异椭球体构造方法

矿化空间分布常常与矿体的形态存在着密切的空间关联关系^[24]，因此，在构造铅、锌品位变异椭球体时，应先充分考虑矿体的产状。一般沿矿体走向方向铅、锌品位变化较慢，具有较大的自相关性，沿矿体厚度方向铅、锌品位变化快，自相关性小，据此可以确定变异椭球体长轴、次长轴和短轴的大致方位。本文以矿体产状为基础，通过计算实验变异函数，在矿体走向倾斜面上研究铅、锌品位的分布来确定品位变异椭球体的长轴方向，在与长轴相垂直倾斜面上研究铅、锌品位分布来确定变异椭球体的次长轴和短轴方向，从而构造出铅、锌品位变异椭球体，反映空间中铅、锌矿体品位的分布和变异情况。

2.1  实验变异函数及其等值线图

变异函数是区域化变量增量的方差^[25]。实验变异函数是指在实际工作中由有限的实测样品值(品位、厚度等)构成的变异函数，能反映区域化变量的空间相关性和随机性，实验变异函数的计算公式为

       (2)

式中：h为滞后距；N(h)指滞后距为h时参加计算的样品个数；Z(x_i)为空间上x_i点的品位。通过计算某一方向上的实验变异函数，并用理论模型进行拟合，得出其空间变异的块金值、基台值和变程。变程反映了该方向上区域化变量自相关的最大距离，在距离范围内，区域化变量具有一定的自相关性，超过此距离时，自相关性消失。

根据实验变异函数的计算结果，将任意平面上每个计算方向上实验变异函数值相同的点连接成线，且相邻等值线的变异函数值之差相同，得到该平面上的实验变异函数等值线图。为了保证实验变异函数等值线图的可靠性，必须均匀地计算平面上各个方向上的实验变异函数，例如，要绘制水平面上的实验变异函数等值线图，可以0°为起始方位角，以15°水平移动增量，依次计算方位为0°，15°，30°，45°，60°，…，330°和345°共24个方向上的实验变异函数。

本文利用Datamine软件，对基于右手系的数学坐标系进行旋转，将XOY平面旋转到空间中的任一平面上，并根据实验变异函数计算结果绘制该平面上的实验变异函数等值线图。为了使生成的实验变异函数等值线图更符合地质人员的观察习惯，先将坐标轴绕Z轴旋转，再绕X轴旋转，并设置合理的旋转角，使得当观察者从倾斜面的倾向看倾斜面时，X轴指向观察者右方的斜面走向，Y轴指向观察者朝向的方向，即倾斜面倾向的反方向，且Y轴与水平面的夹角等于倾角，Z轴垂直倾斜面朝上。假设其旋转角分别为μ和ν，采用左手规则的约定来确定旋转角的正负号^[26]，即沿着坐标轴正向朝原点看时，顺时针旋转该坐标轴为正角，逆时针旋转为负角，则其旋转矩阵为

    (3)

实验变异函数等值线图上等值线的分布能反映区域化变量在某一平面不同方向上的变化速度：沿等值线延伸方向的区域化变量变化慢，该方向上的区域化变量具有较大的空间自相关性；与等值线延伸方向相垂直的方向区域化变量变化快，具有较小的空间自相关性。本文选择品位作为区域化变量，绘制不同平面上品位的实验变异函数等直线图，通过判读等值线分布了解矿体品位的空间分布及变异情况。

2.2  变异轴空间方位表达与计算

由式(1)、图1和图2可知：变异椭球体3个互相垂直主要轴(称为变异轴)的空间方位和长度唯一确定了变异椭球体的空间形态，同时在构造变异椭球体的过程中，可以通过间接构造变异轴的空间方位、长度来构造变异椭球体，由此可见，研究变异轴的空间方位表示方法是有必要的。变异轴类似于线状构造，其产状表达方法有以下2种。

1) 侧伏向、侧伏角和长度。侧伏向是指某线段所在的倾斜面走向线所指的该线段向下倾斜的那一端方位；侧伏角是指当线段在某一倾斜平面内时，此线与该平面走向线间所夹之锐角。如图3所示，OP为线段OQ在与水平面ABCD上的投影线，KM为倾斜面KMNJ的走向，则λ和为线段OQ在倾斜面KMNJ上的侧伏角，MK为OQ在面KMNJ上的侧伏向，其长度是该方位上变程的2倍。

2) 倾伏向、倾伏角和长度。倾伏向是指某一倾斜线段在水平面上的投影线所指的该线段向下倾斜的方位；倾伏角是指某一线段与其水平投影线间所夹的锐角。如图3所示，β即为OQ的倾伏角，PO为OQ的倾伏向，其长度是该方位上变程的2倍。

侧伏向和侧伏角是相对于线段所在的倾斜面而言的。使用侧伏向、侧伏角和长度表示变异轴的形态时，一定要指出参照平面，也就是说侧伏向、侧伏角和长度不能独立地表示变异轴的空间方位，使用倾伏向和倾伏角的表示方法显得简洁、明了，能独立表示变异轴的空间方位。本文通过判读实验变异函数等值线图只能得出变异轴在斜面上的侧伏角和侧伏向，故还必须推导公式求出变异轴的倾伏角和倾伏向来独立表示变异轴的空间方位。

线段OQ的倾伏向、倾伏角求解图如图3所示。要求得线段OQ倾伏向PO，必须先求出线段OP与倾斜面走向KM之间的夹角ε，根据ε和KM的方位及PO与KM的位置关系，可以快速得出PO的方位。已知θ为倾斜面KMNJ的倾角，则β和ε的推导公式如下：

；

；

；

；

故     (4)

      (5)

图3  线段OQ的倾伏向、倾伏角求解

Fig. 3  Solving of pitch direction, pitch angle of line segment OQ

2.3  变异椭球体的构造

变异椭球不仅反映了矿体品位的空间分布及变异情况，而且能作为搜索椭球体来限制品位插值过程中已知点的搜索范围，所以，变异椭球体的构造为分析矿化空间分布规律和资源量估算过程中的重要步骤。然而，一般研究不重视变异椭球体的构造过程和方法，仅仅依赖矿体产状来构造变异椭球体，无法模拟自然界中复杂的品位分布及变异情况。本文提出一种综合矿体产状和矿体品位分布的变异椭球体构造方法，具体步骤如下。

Step 1：分析研究区矿体产状，确定矿体的走向、倾向和倾角。

Step 2：建立综合地质数据库，并导入三维矿业软件建立矿床的勘探工程模型，确定样品的空间位置。

Step 3：为了保证统计参数是无偏估计量，参与统计的样品数据均应具有相同的支架或承载，即样品数据应该具有统一的样品长度^[27]，故需对原始样品进行组合。按照最少拆分原始取样长度的原则^[28]选取组合样长度对矿床的原始样品进行组合，得到长度相等的组合样数据，并对组合样品数据进行统计分析，了解品位的空间分布状况和矿体分布规律，剔除特高品位数据。

Step 4：运用式(3)所示的旋转矩阵R₁，将XOY平面旋转到矿体走向倾斜面上，以处理好的组合样品数据为原始数据，选择品位变量作为区域化变量。设置合适的滞后距、起始方位角、起始倾角、移动增量，计算矿体走向倾斜面上各个方向上的实验变异函数，绘制矿体走向倾斜面上的实验变异函数等值线图。

Step 5：对矿体走向倾斜面上的实验变异函数等值线分布进行判读，选择此面上等值线延伸方向作为变异椭球体的长轴方向。当等值线延伸方向不是很精确时，可以在一个范围内按一定的增量确定多个方向。在确定方向的同时，可以量出这些方向上的线段在该倾斜面上的侧伏角和侧伏向。

Step 6：运用式(4)和(5)求出矿体倾斜面上多个待定长轴方向上线段的倾伏角和倾伏向，计算这些方向上的实验变异函数，并分别用理论模型进行拟合，得到品位变量在这些方向上的理论变异函数模型。变异椭球体的长轴方向为其理论变异函数模型具有最大变程的那个方向，长半轴长度为其变程。

Step 7：用式(3)所示的旋转矩阵对坐标轴进行旋转变换，将XOY平面旋转到与长轴相垂直的平面上，计算该平面上各个方向上的实验变异函数，根据计算结果绘制该平面上的实验变异函数等值线图。重复Step 5和Step 6，确定变异椭球体的次长轴方向和长度，则该平面上与其相垂直的方向为变异椭球体的短轴方向，该方向上的理论变异函数模型的变程为短半轴的长度。

Step 8：将坐标轴旋转至与变异椭球体的3个主要轴重合，采用式(1)所示的椭球体数学模型生成变异椭球体。本研究先将坐标轴绕Z轴旋转使X轴移至变异椭球体长轴的水平面投影线上，再绕Y轴旋转使X轴移至变异椭球体长轴上，最后绕X轴旋转使Y轴、Z轴分别移至变异椭球次长轴、短轴上，使得X轴、Y轴、Z轴分别与长轴、次长轴和短轴重合。假设其旋转角分别为u，v和w，采用左手规则的约定来确定旋转角的正负号^[26]，则其旋转矩阵为

            (6)

3  应用实例

虎头崖铅锌多金属矿床由15个主要矿体组成近东西走向的铅锌矿带，矿体呈层状、透镜状产出。整个矿带大致可以分为2部分：(1) 67线以东包括矿体Ⅰ-17，Ⅰ-31，Ⅰ-32，Ⅰ-33，Ⅰ-34，Ⅰ-35和Ⅰ-36，矿体走向东西，倾向南，倾角为70°~85°；(2) 117线以西包括矿体Ⅰ-19，Ⅰ-20，Ⅰ-21，Ⅰ-1，Ⅰ-2，Ⅰ-3，Ⅰ-4和Ⅰ-5，矿体走向变为70°~250°，倾向南西，倾角为75°~85°。

本文以虎头崖铅锌矿床为例，对其进行矿化空间分析。采用所介绍的变异椭球体模型构造方法，利用三维矿业软件Datamine计算实验变异函数和绘制变异函数等值线图，实现对矿床铅、锌品位变异椭球体的构造，分析虎头崖矿床矿体铅、锌品位的空间分布结构及各个方位的自相关性，并将变异椭球体作为品位空间插值过程中的搜索椭球体来限制已知样品点的搜索范围，得到矿区的品位和资源量的估算结果。

3.1  数据预处理

虎头崖矿床的矿体分2个矿段分布于67线以东和117线以西，且距离较远，故分开对这2个矿段进行研究。对虎头崖矿床的相关原始资料进行数字化，建立矿床67线以东和117线以西的综合地质数据库，并分别导入到Datamine中建立其三维勘探工程模型，确定取样的空间位置。

对矿床的原始样品进行预处理，按样品平均值的8倍剔除特高品位样品，按边界品位剔除低品位样品，然后分别选择原始样品平均长度1 m作为样品组合长度对矿床67线以东和117线以西的样品进行组合，形成长度相等的组合样品。67线以东的原始样品经组合后得到铅、锌组合样2 660个，117线以西的原始样品经组合后得到长度为1 m的铅、锌组合样2 053个。分别对其组合样品铅锌品位进行统计分析，了解铅锌品位的空间分布特征，统计结果见表1。

3.2  矿床变异椭球体的构造

先研究矿床67线以东铅品位的变异椭球体构造。67线以东的矿体总体上走向东西，倾向南，倾角为75°。以67线以东的铅品位的组合样品数据为原始数据，选择铅品位作为区域化变量，每个方向按照10 m步长、15°容差角、100 m容差限来计算实验变异函数值。

表1  虎头崖矿床组合样品Pb和Zn品位统计表

Table 1  Statistical table of Pb and Zn grade of composite sample in Hutouya deposit

通过旋转矩阵R₁ (其中，u= 0°，v=-75°)，将XOY平面旋转到矿体走向倾斜面上，根据实验变异函数的计算结果绘制出该面上的实验变异函数等值线图，如图4所示。从图4可以看出：等值线沿北西—南东方向延伸，初步确定为变异椭球体长轴方向。沿此平面上侧伏角为50°，侧伏向为东的方向依次顺时针旋转5°至此平面上侧伏角为75°，侧伏向为东；依次逆时针旋转5°至此平面上侧伏角为30°，侧伏向为东。计算这些待定方向上的实验变异函数模型，并用球状模型分别进行拟合，得到这些方向上的块金值、基台值和变程。经过实验发现，此平面上侧伏角为60°，侧伏向为东的方向上的理论变异函数模型具有最大变程，选为变异椭球体的长轴方向，表明沿该方向铅、锌品位在空间上变化相对较慢，具有较大的空间自相关距离。

根据式(4)和式(5)算得长轴的倾伏角为；长轴在水平面的投影线与倾斜面走向之间的夹角为。故长轴的倾伏向为90°+24°=114°，则长轴的方向为(SE114°，∠56.8°)。

图4  虎头崖矿床67线以东Pb品位在走向倾斜面上的实验变异函数等值线图

Fig. 4  Contour of experimental variogram of Pb grade on inclined surface along trend in the east of exploration line 67 in Hutouya Mining Deposit

为了进一步确定变异椭球体的次长轴和短轴的方向，对坐标轴进行旋转，通过旋转矩阵R₁ (其中，u= 114°，v=-33.2°)，将XOY平面旋转到与长轴方向相垂直的平面上，根据实验变异函数的计算结果绘制出该平面上的实验变异函数等值线图，如图5所示。从图5可以看出：等值线的延伸方向为北东—南西向，侧伏角λ为53°，侧伏向为SW204°，该方向为变异椭球体的次长轴方向。

图5  虎头崖矿床67线以东Pb品位在长轴方向相垂直面的实验变异函数等值线图

Fig. 5  Contour of experimental variogram of Pb grade on inclined surface perpendicular to long axis in the east of exploration line 67 in Hutouya mining deposit

根据式(4)和(5)算得次长轴的倾伏角为；次长轴在水平面的投影线与该倾斜面走向之间的夹角为。故次长轴的倾伏向为204°+48°=252°，则次长轴的方向为(SW252°，∠25.9°)。该方向上的实验变异函数模型和理论变异函数模型如图6所示，根据拟合结果得到其变程为60 m，为次长半轴的长度。

图5中与次长轴相垂直的方向为短轴方向，其侧伏角为90°-53°=37°，侧伏向为NE24°，根据式(4)和(5)可算得短轴的倾伏角为 ；短轴在水平面的投影线与该斜面走向间的夹角为 ，故短轴的方位角为24.0°-32.2°=-8.2°，短轴方向为(NW351.8°，∠19.2°)。该方向上的实验变异函数模型和理论变异函数模型如图6所示，根据拟合结果得到其变程为25 m，为短半轴的长度。

图6  虎头崖矿床67线以东铅品位长轴、次长轴、短轴方向上的实验及理论变异函数曲线图

Fig. 6  Experimental and theoretical variogram graph of Pb grade in long axis,secondary long axis,short axis in the east of exploration line 67 in Hutouya Deposit

虎头崖矿床67线以东铅品位变异轴方向上的理论变异函数模型具体参数如表2所示。对其进行交叉验证，得到其平均差为0.009 3，标准差为1.150 3，如表3所示。从表3可知：品位的实际值均值与估计值均值相当接近，表明其估值是无偏的，理论变异函数模型可用于虎头崖矿床67线以东铅品位估算。

根据上述分析，反映铅、锌品位空间变异方向的椭球体的3个轴的方向和半径为：1) 长轴，(SE114°，∠56.8°)，85 m；2) 次长轴，(SW252°，∠25.9°), 60 m；3) 短轴，(NW351.8°，∠19.2°)，25 m。

利用式(6)所示的旋转矩阵R₂ (其中，u=24°，v=-56.8°，w=127.01°)，分别将X轴、Y轴和Z轴旋转至变异椭球体长轴、次长轴和短轴方向，采用椭球体数学模型，利用Datamine软件实现各向异性套合，生成变异椭球体如图7所示。

虎头崖矿床67线以东的铅、锌品位在空间上的变化具有明显的结构性特征。上述变异函数模型和椭球体参数表明：在不同方向上的空间相关性或空间连续性有明显的差别，沿(SE114°，∠56.8°)方向变化速度最慢，连续性较好，在达到约85 m的距离范围内，铅、锌品位均具有一定程度的空间相关性；沿(NW351.8°，∠19.2°)方向，铅、锌品位的变化速度最快，连续性较差，在达到约25 m的距离后，铅、锌品位的空间相关性消失。

图7  虎头崖矿床67线以东Pb品位空间变异椭球体

Fig.7  Spatial variation ellipsoid of Pb grade in the east of exploration line 67 in Hutouya Deposit

采用同样的方法构造虎头崖矿床67线以东的锌品位和117线以西的铅、锌品位变异椭球体，得到其变异轴方向上的理论变异函数模型如表2所示，其交叉验证结果如表3所示。

3.3  基于变异椭球体模型的资源量估算

以矿体块体模型为基础，采用对数正态克立格法对品位进行估值，并将得到的矿床铅、锌品位变异椭球体作为搜索椭球体来限制空间插值过程中已知样品点的搜索范围，得到矿体中每个块体单元铅、锌的平均品位及金属量。若变异椭球体构造合理，则与待估块体单元相关的已知样品点会都位于该搜索椭球范围内并参与估值，得到的品位和资源量的估算结果会更准确，更接近实际；反之，则估算结果与实际不相符。利用这一原理，可间接地评价变异椭球体构造的正确性。

表2  虎头崖矿床Pb和Zn品位变异轴方向上的理论变异函数模型

Table 2  Parameters of spatial variation ellipsoid of Pb, Zn grade in Hutouya deposit

表3  虎头崖矿床Pb和Zn品位理论变异函数模型交叉验证结果

Table3  Cross-validation results of theoretical semi-variogram model of Pb and Zn grade in Hutouya deposit

将本文得到的铅、锌品位和资源量估算结果与矿山实际开采获得的平均品位和资源量进行对比，发现本文的计算结果和矿山开采得到的结果基本一致，由此可知本文所构造的变异椭球体具有合理性。

4  结论

1) 利用旋转矩阵-，将XOY旋转到矿体走向倾斜面和与长轴相垂直的平面上，计算平面上各个方向上的实验变异函数，并绘制出符合地质人员观察习惯的实验变异函数等值线图。等值线的分布反映品位变量在该平面上各个方向上的变异情况。

2) 构造变异椭球体的关键是确定变异轴的空间方位和长度。通过判读实验变异函数等值线图来确定变异椭球体变异轴的空间方位，按照本文推导得到的变异轴倾伏向、倾伏角的计算公式，计算获得用倾伏向、倾伏角表示的变异轴的空间方位，并利用球状模型拟合变异轴方向上的实验变异函数，将球状模型的变程参数作为变异半轴的长度。

3) 将本文提出的变异椭球体构造方法应用到虎头崖矿床中，得到了研究区铅、锌品位的空间变异椭球体，反映了虎头崖矿床铅、锌品位的空间分布及各个方位上的变异情况，并将变异椭球体作为搜索椭球体应用到品位空间插值过程中，得到矿体的品位和资源量估算结果。通过与矿山开采得到的结果进行对比，发现本文构造的变异椭球体具有合理性，反映了品位变量在各个方位上的自相关性。

参考文献：

[1] Emery X. Variograms of order omega: A tool to validate a bivariate distribution model[J]. Mathematical Geology, 2005, 37(2): 163-181.

[2] 肖斌, 赵鹏大, 侯景儒. 地质统计学新进展[J]. 地球科学进展, 2000, 15(3): 293-296.

XIAO Bin, ZHAO Pengda, HOU Jingru. New development of geostatistics[J]. Advance in Earth Sciences, 2000, 15(3): 293-296.

[3] 吴胜和, 李文克. 多点地质统计学: 理论、应用与展望[J]. 古地理学报, 2005, 7(1): 137-144.

WU Shenghe, LI Wenke. Multiple-point geostatistics: Theory, application and perspective[J]. Journal of Palaeogeography, 2005, 7(1): 137-144.

[4] 余先川, 侯景儒, 姚力, 等. 时空域非参数和多元信息的地质统计学研究[J]. 自然科学进展, 2003, 13(11): 1217-1220.

YU Xianchuan, HOU Jingru, YAO Li, et al. Research of non-parametric and multi-informational geostatistics in space-time[J]. Advance in Nature Sciences, 2003, 13(11): 1217-1220.

[5] Yunsel T Y. A practical application of geostatistical methods to quality and mineral reserve modelling of cement raw materials[J]. Journal of the Southern African Institute of Mining and Metallurgy, 2012, 112(3): 239-249.

[6] Armony M. Field parametric geostatistics[J]. International Journal of Surface Mining, Reclamation and Environment, 2001, 15(2): 100-122.

[7] Sarkar B C, Roy I. A geostatistical approach to resource evaluation of Kalta iron ore deposit, Sundargarh district, Orissa[J]. Journal of the Geological Society of India, 2005, 65(5): 553-561.

[8] Modis K, Stavrou S, Terezopoulos N, et al. Geostatistics versus inverse distance squares in ore reserves estimation comparative case study in copper ore body in Cyprus[J]. Mining Technology, 2008, 117(1): 48-52.

[9] 胡平昭, 侯景儒, 朱礼, 等. 多元空间信息统计分析计算机系统[J]. 遥感信息, 1998(2): 2-5.

HU Pingzhao, HOU Jingru, ZHU Li, et al. Computer system of statistical analysis of multivariate spatial information[J]. Remote Information, 1998(2): 2-5.

[10] 肖斌, 侯景儒, 肖克炎, 等. 空间信息统计学可视化软件系统[J]. 地质论评, 2000, 46: 176-181.

XIAO Bin, HOU Jingru, XIAO Keyan, et al. Visually statistical software system of spatial information[J]. Geological Review, 2000, 46: 176-181.

[11] CHEN Jianping, TAN Xiaojin. Mining spatial association rules with geostatistics[C]// Zhang J X, Goodchild M F. Proceedings of the 8th international symposium on spatial accuracy assessment in natural resources and environmental sciences. Liverpool, England: World ACAD Union-Word ACAD Press, 2008: 369-376.

[12] Daya A A, Moarefvand P, Hassani H. Geostatistical modeling and reserve estimation of central part of north anomaly of choghart iron ore deposit through ordinary kriging method[C]// Baby S, Sandhu P S, Hang Y. Proceedings of 2010 International Conference on Environmental Science and Development. Liverpool, England: World ACAD Union-Word ACAD Press, 2010: 407-411.

[13] Roy D, Butt S D, Frempong P K. Geostatistical resource estimation for the Poura narrow-vein gold deposit[J]. Cim Bulletin, 2004, 97(1077): 47-51.

[14] Sarkar B C, Nair A M. A geostatistical modelling approach to gold mineralisation at Hutti mine, Raichur District, Karnataka[J]. Journal of the Geological Society of India, 2002, 60(6): 639-648.

[15] Fisher T R, Dagdelen K, Turner A K. Modelling 3D grade distributions on the Tarkwa paleoplacer gold deposit, Ghana, Africa[C]// Leuangthong O, Deutsch V C. Quantitative Geology and Geostatistics. Netherlands: Springer, 2005: 439-448.

[16] 余海军, 李文昌, 曾普胜, 等. 地质统计学在羊拉铜矿储量计算中的应用[J]. 地球学报, 2009, 30(5): 684-690.

YU Haijun, LI Wenchang, ZENG Pusheng, et al. The application of geostatistics to ore reserve calculation of the Yangla copper deposit[J]. Acta Geoscientica Sinica, 2009, 30(5): 684-690.

[17] 杨建宇, 秦德先, 康泽宁, 等. 北衙金矿三维模型的建立及研究[J]. 有色金属, 2006, 58(2): 17-20.

YANG Jianyu, QIN Dexian, KANG Zening, et al. Establishment and research of 3D model in Beiya gold mine[J]. Nonferrous Metals, 2006, 58(2): 17-20.

[18] 胡建明. 3DMine矿业软件在地勘工作中的应用[J]. 矿产勘查, 2010, 1(1): 78-80.

HU Jianming. Application of 3DMine software on geological exploration[J]. Mineral Exploration, 2010, 1(1): 78-80.

[19] 厉金龙. 基于Surpac 的三维地质建模及可视化研究[D]. 重庆: 重庆大学土木工程学院, 2011: 2-9.

LI Jinlong. Study of 3D geological modeling and visualization based on Surpac[D]. Chongqing: Chongqing University. College of Civil Engineering, 2011: 2-9.

[20] Johnson T. Structural block modeling at BHP iron ore’s Mt Whaleback Mine[J]. Australasian Institute of Mining and Metallurgy, 2001, 23: 169-186.

[21] 李岭. 地质统计学及其在某铜矿山储量计算中的应用[J]. 矿业研究与开发, 2004, 24(5): 21-23.

LI Ling. Geostatistics and the corresponding application of reserves estimation in some copper mine[J]. Mining Research and Development, 2004, 24(5): 21-23.

[22] 孙英君, 王劲峰, 柏延臣. 地质统计学方法进展研究[J]. 地球科学进展, 2004, 19(2): 268-274.

SUN Yingjun, WANG Jinfeng, BAI Yanchen. Study on progress of methods in geostatistics[J]. Advance in Earth Sciences, 2004, 19(2): 268-274.

[23] 张朔, 鲁学军, 赵英俊, 等. 空间变异函数等效应椭圆套合方法及其应用[J]. 地球信息科学学报, 2009, 11(3): 342-348.

ZHANG Shuo, LU Xuejun, ZHAO Yingjun, et al. Equivalent elliptic nested structure and its application in variogram nesting[J]. Journal of Geo-information Science, 2009, 11(3): 342-348.

[24] 毛先成, 唐艳华, 邓浩. 地质体的三维形态分析方法与应用[J]. 中南大学学报(自然科学版), 2012, 43(2): 588-595.

MAO Xiancheng, TANG Yanhua, DENG Hao. Three dimensional morphological analysis method for geologic bodies and its application[J]. Journal of Central South University (Science and Technology), 2012, 43(2): 588-595.

[25] 邹陈. 变异函数对矿体变化性的指示及其在矿床勘探网度优化中的应用[D]. 长沙: 中南大学地学与环境工程学院, 2010: 8-10.

ZOU Chen. Revealing variation of ore body by variogram and its application in optimizing exploration grid[D]. Changsha: Central South University. College of Geoscience and Environmental Engineering, 2010: 8-10.

[26] 张全伙, 张剑达. 计算机图形学[M]. 北京: 机械工业出版社, 2004: 9-25.

ZHANG Quanhuo, ZHANG Jianda. Computer graphics[M]. Beijing: China Machine Press, 2004: 9-25.

[27] 万昌林, 朱利平, 高祥. 应用地质统计学评估福建某铜矿资源[J]. 金属矿山, 2002(6): 140-143.

WAN Changlin, ZHU Liping, GAO Xiang. Applying geological statistics to evaluate a certain copper ore resource in Fujian[J]. Metal Mine, 2002(6): 140-143.

[28] 刘海英. 多种克里格方法在固矿储量估算中的应用研究[D]. 武汉: 中国地质大学信息工程学院, 2010: 8-10.

LIU Haiying. A study on application of multiple Kriging methods in solid mineral reserve estimation[D]. Wuhan: China University of Geosciences. School of Information Engineering, 2010: 8-10.

(编辑  陈灿华)

收稿日期：2013-09-19；修回日期：2013-11-23

基金项目：国家自然科学基金资助项目(41172297，41102204)；国家“十二五”科技支撑计划项目(2011BAB04B10)

通信作者：毛先成(1963-)，男，湖南澧县人，教授，博士生导师，从事隐伏矿床预测、地学信息技术研究；电话：0731-88877571；E-mail: xcmao@126.com