(一)实验变差函数的计算

我们知道实验变差函数的计算公式为

地质统计学(空间信息统计学)基本理论与方法应用

计算时,是根据数据的空间构形进行的。不同的数据空间构形按以下3种不同的方法计算。

(1)列线等间距的数据空间构形的计算

这种数据空间构形依据研究对象的实际情况,可分作一维、二维和三维空间进行计算。这种构形计算比较简单,只需按某一方向的点距等于基本滞后距,逐点对取值按公式 便可计算实验变差函数,如下图。

设沿α方向和β方向的基本滞后距分别为a和b,实验变差函数分别为 其中k1和k2等于1,2,…,n。

地质统计学(空间信息统计学)基本理论与方法应用

(2)列线非等间距数据构形的计算

对于这一构形主要是合理地解决好数据对问题,通常的作法是,将非等间距数据点按区间[hk-ε(h),hk+ε(h)]组合成距离组,然后用[kh±e(h)]内的所有数据对来计算γ 。譬如对于点xi,与xi的距离在[kh-ε(h)k,h+ε(h)]范围内的数据点(钻孔点),却可以被认作是xi距离为kh的点,ε(h)称作距离误差,对于ε(h)的选择,一般小于样品间距(h)的一半,

[在列线非等间距的构形中,因数据点的间距(h)不等,通常的作法是以统计多数样品的间距为准,取其 若统计多数样品的间距在100m左右,则确定h=100m,ε(h)≤50m]

这里需补充的是,在计算 为了使两钻孔点的距离接近kh,通常选用:

地质统计学(空间信息统计学)基本理论与方法应用

其中hj为任一钻孔距x点的实际距离。

此处的hj≥kh+ε(h)或hj≤k-hε(h)。

(3)非列线不等间距数据构形的计算

对于这种数据构形,其数据点对的确定,通常是按一定的角度和距离误差来组合数据对。

具体步骤如下:首先将数据组合成角度组。设dα为x方向的角度允许误差限,这样,α±αd内的数据都可认为是沿x方向的数据(dα一般为两相邻方向的夹角的 最大不能超过其夹角的一半第二步将数据点按距离组[hk±ε(h)]进行组合,在角度范围(α±αd)及距离范围[hk±ε(h)]内的数据点都可认为是xi点在α方向上相距为hk的数据点(如下图)。

地质统计学(空间信息统计学)基本理论与方法应用

主要是你的数据问题，从最后一张变异函数图可以看出来，原数据最少有三个套合结构，整体看不出变程和基台值，不符合球状模型。看上去又不像指数模型，surfur中不知道可不可以选一般线性模型，如果可以选择，那么可以尝试用一般线性模型模拟。

从原数据的图形中可以明显看出呈线性，所以建议用一般线性模型做进一步尝试。

另一方面，因为我不知道你的原数据来源，是否具有物理特性，如果就是一组纯数字，建议把步长（滞后距）分别调整到10，15，20，相应的最大滞后距也调整到200左右，再看看图形是否符合球状模型特征？滞后个数就设置20几个就差不多。如果你得数据是有特定物理意义的，那么上述尝试就不用了。

以丁家山矿区的组合样品数据为原始数据，分别以铅品位和锌品位为区域化变量，计 算了三维空间中各个方位的实验半变异函数，绘制主要平面方向上的实验半变异函数等值 线图，构造三维变异椭球体，揭示矿体标志空间变化的结构性特征。

对三维空间中各个方向均计算了实验半变异函数。在任意一个方向面上，计算了24 个方向的实验半变异函数，每一个方向按照10m的滞后距计算实验半变异函数值γ（h）。 在每一个方向面上，根据实验半变异函数的计算结果，可绘制实验半变异函数值γ（h）的 等值线图。

在图59中，59a和59b是矿体走向平面上铅和锌品位的实验半变异函数等值线图，从图中可知实验半变异函数等值线的延伸方向为NE－SW（约NE39°），表明沿NE－SW方 向，铅和锌品位在空间上变化相对较慢，具有较大的空间自相关距离；而在与该方向垂直 的NW－SE方向上，铅和锌品位在空间上变化相对较快，具有较小的空间自相关距离。由 此可以初步确定，三维变异椭球体的长轴方向为NE39°、倾角0°。

为了进一步确定变异椭球体的次长轴、短轴的方向，对坐标轴进行旋转，计算并绘制 出与变异椭球体长轴方向相垂直的平面（矿体横剖面）方向上的变异函数等值线图（图 510），选择NW309°、倾角0°为椭球体次长轴方向，SE 129°、倾角90°为椭球体短轴方向。 图511所示为椭球体长轴与短轴所在平面（NE39°∠0°，矿体纵剖面）方向上的实验半 变异函数等值线图，从该图可知，等值线的分布形态与长轴、短轴的方向延伸基本一致。

根据上述分析，反映铅品位和锌品位空间变异方向的椭球体的三个轴的方向为：(1)长轴，NE39° ∠0°（或SW219° ∠0°）； (2)次长轴（次短轴 ），NW309° ∠0°（或 SE 129° ∠0°）；(3)短轴，NW309° ∠-90°（或SE 129° ∠90°）

对变异椭球体三个轴方向的实验半变异函数，采用球状模型进行拟合（孙洪泉，1990），可得到该三个方向的理论半变异函数模型（图512与图513）。

图512和图513分别是铅和锌的实验及理论变异函数的曲线图，其中虚线是实验变 异函数曲线，实线是与之具有相同颜色的实验变异函数的拟合曲线。图512中，X方向的 拟合曲线和Y方向的拟合曲线相重合，所以图中只显示了X方向的拟合曲线。

根据上述实验半变异函数的理论模型（球状模型）拟合结果，得到了铅品位和锌品位 空间变异椭球体的参数，如表55和图514、图515所示。表55中的方向、变程分别为图514 和图515所示椭球体的三个轴的延伸方向、长度。

表55 球状变异函数模型主要参数

铅品位和锌品位在空间上的变化具有明显的结构性特征。上述半变异函数模型和椭球 体参数表明，在不同方向上的空间相关性或空间连续性有明显的差别，沿NE39° ∠ 0°方向，铅品位和锌品位的变化速度最慢，连续性较好，分别在达到约241m、约181m的距离范围 内，铅品位和锌品位均具有一定程度的空间相关性；而沿SE129°∠90°方向，铅品位和锌 品位的变化速度最快，连续性较差，分别在达到约100m、约90m的距离后，铅品位和锌 品位的空间相关性消失。

图59 铅（a）和锌（b）品位在矿体走向平面上的实验变异函数等值线图（图面方向：上NE39°，左NW309°，下SW219°，右SE129°；图面倾角：0°）

图510 铅（a）和锌（b）品位在矿体横剖面上的实验变异函数等值线图（图面方向：左NW309°，右SE129°；图面倾角：90°）

图511 铅（a）和锌（b）品位在矿体纵剖面上的实验变异函数等值线图（图面方向：左SW219°，右NE39°；图面倾角：90°）

图512 铅品位的实验及理论变异函数的曲线图

图513 锌品位的实验及理论变异函数的曲线图

图514 铅品位空间变异椭球体（X轴—东，Y轴—北，Z轴—上）

图515 锌品位空间变异椭球体（X轴—东，Y 轴—北，Z轴—上）

将铅品位和锌品位的空间变异椭球体参数与主要铅矿体和锌矿体的形态产状对比可 知，铅品位和锌品位的变异函数模型与变异椭球体较好地反映了铅矿体和锌矿体的空间分 布形态、产状和品位变化。其一，椭球体的长轴、次长轴、短轴的方向分别与矿区最大的 矿体Ⅲ1及Ⅲ2的走向、倾向及倾角、厚度方向基本一致（Ⅲ1及Ⅲ2矿体向南西倾伏，呈 北东高，南西低，南东翼向南东倾，北西翼向北西倾，倾角0°～25° ）；其二，椭球体的形 态及空间方位与主要矿体Ⅲ1及Ⅲ2的空间展布形态和延伸方向相似；其三，变异函数模型 表达的不同方位的空间相关距离和连续性差异，与主要矿体的铅品位和锌品位的空间分布 及变化相似（矿体沿走向连续性好、沿厚度方向连续性差）。

对于关兜矿区的矿体空间变化结构分析，由于该区已有的组合样品数量过少，无法有 效取得具有统计意义的结构分析结果，因此，在此未对其进行地质统计学意义上的分析。 其后续的资源量估算也只采用了常用的距离平方反比法来进行。

341 桂林市区岩溶地下水水位观测网分布特征及存在问题

3411 桂林市区岩溶地下水水位观测网分布特征

图311 桂林市区岩溶地下水水位观测网分布示意图

（根据1990年9月水位资料）

至1990年，在区内不同地下水类型、环境地质亚区和环境地质问题突出段先后布置了60个水位观测孔，其中钻孔35个，天然水点14个，大口井11个；孔隙水11处，岩溶地下水49处；桂林市区岩溶地下水水位观测网分布状况详见图311。从分布图上可看出，观测孔主要集中分布在市区中心，研究区周边区域很少或没有分布观测孔。

3412 桂林市区岩溶地下水水位观测网存在的主要问题

桂林市区岩溶地下水水位观测网始建于1981年，当时很多观测孔来源于地质勘察孔或其他目的的施工钻孔，存在一些不合理、不完善的地方。其存在的主要问题如下：

1）观测网提供的信息数据具有冗余性。根据水位资料信息和分布图分析可知，一些相邻的观测孔的水位及水位动态变化情况基本相同。如观测孔GⅢ46、GⅢ47（图312）。

图312 1990年观测孔GⅢ46、GⅢ47地下水水位动态曲线图

2）监测点分布不合理。在排泄区布置的观测孔太多，而在一些径流区布置的观测孔过少。如在朝阳北侧、研究区西南部为地下水径流区，水位变化较大，水力坡度变化也大，布置的观测孔都太少。

3）地下水监测中水位监测与水质监测不同步，长期以来重水位而轻水质，不能满足全面评价地下水资源质量的要求。

342 桂林市区岩溶地下水水位观测网优化研究

3421 研究对象与目标

桂林市区岩溶地下水水分布范围广，水质好，是目前主要的开采层。但在强岩溶发育地区，岩溶地下水的过度开采往往会造成地面沉降、塌陷等地质灾害。因此，必须建立合理的地下水观测网，对岩溶地下水位进行实时观测，才能真正做到地下水的合理开采，防止漏斗面积的扩散及保护地下水水质。以岩溶地下水水位为研究对象，对桂林市区岩溶地下水水位观测网进行优化设计。

选用地下水位估计误差标准差为代用目标函数，用整个观测网平均估计误差标准差作为评价观测网整体精度的标准，优化的目标是在节省一定费用情况下，使优化后平均估计误差标准差与原观测网的平均估计误差标准差相差不大。

3422 地下水观测网点布设原则

观测孔总的布置原则是：对于面积较大的监测区域，以顺沿地下水流向为主与垂直地下水流向为辅相结合布设监测网；对于面积较小的监测区域，根据地下水的补给、径流、排泄条件来布置控制性监测点。

1）国家级区域地下水动态监测点，应在水文地质单元和含水层层序划分的基础上，依据地质环境背景和水文地质条件进行布置。主要布置在：

a岩溶水具有供水意义的地区，以及已经产生或可能产生岩溶塌陷的地区。

b已经或将要形成区域环境地质问题的地区。

2）省级地市区地下水动态监测点，要在国家级城市区地下水动态监测网的基础上布置。布置时应考虑以下几个方面：

a在城市供水水源地的补给、径流、排泄区，污染源附近和水源地保护区，均应布设监测点。

b在水源地应在平行和垂直于地下水流向布设两条监测线，以监测地下水位下降漏斗的形成和发展趋势。

c在查明水源地之间的相互影响或附近矿区排水对水源地的影响时，应于连接两个开采区的地带布设监测点。

d为建立城市地下水均衡计算模型或地下水管理模型的需要，可在边界处及计算分区内布设控制性监测点。

3423 桂林市区岩溶地下水位观测网优化设计研究

（1）岩溶水文地质半定量分析

考虑到克立格法人为增减观测孔的不准确性，从影响观测网优化设计布设的主要因素（岩溶地下水防护条件、岩溶地下水富水性、岩溶发育特征以及岩溶地下水开采条件等）角度分析提出了岩溶水文地质半定量分析法，根据其生成的半定量分区图与水位估计误差方差等值线图合成一个综合图形，通过综合图形对观测孔进行增减，不但能更准确地布置观测孔的位置，而且能提高优化效率。其方法介绍如下。

1）赋值。按各因素对桂林市区岩溶地下水位观测网布设的影响效果，分为大、较大、中、小四级标准。根据桂林市区实际情况，选取的主要影响因素及其赋值标准如表33所示。

表33 各影响因素赋值标准

2）权重的确定。权重反映各影响因素的相对重要程度，它既是决策者的主观评价，又是指标本质的物理属性的客观反映，是主客观综合度量的结果。权重主要决定于两个方面：第一，指标本身在决策中的作用和指标价值的可靠程度；第二，决策者对该指标的重视程度。目前，对指标的权重的确定存在着不少主观随意性，严重影响着评价结果的客观性。建立了最优化模型，对主客观赋权进行线性加权法求取权重，其表达式为

岩溶地区地下水与环境的特殊性研究

求解上述模型，即得综合集成权系数：

岩溶地区地下水与环境的特殊性研究

式中： 为专家所评各因素的权重；hk为各专家所评结果的重要程度， hk=1，q为专家人数；α，β分别表示主、客观赋权方法的相对重要程度，α+β=1；m为评价方案的个数；αij为第i个评价方案的第j个因素值； 为各评价方案中第j个因素值的最优值，即α中最优值，一般取其平均值；ij 为综合集成后的权重向量。

根据多名深入了解研究区水文地质条件的水利专家所建议各影响因素的初始权重，求得各影响因素综合权重值如表34所示。

表34 各影响因素权重赋值表

3）综合评价。基于MapGIS软件具有任意图形边界的多图层矢量数据叠加的功能，对岩溶地下水防护条件分区图、岩溶发育强度分区图等影响桂林市区岩溶地下水水位观测网优化设计的各因素的分区图进行图层叠加，在叠加后所形成的小分区内采用加权平均模型计算各小分区的综合评价值，其数学模型表达式为

岩溶地区地下水与环境的特殊性研究

式中：R为各区的综合评价值； 为各区因素的权重；Yj为各区内因素的赋值；n为各区内所有因素的数量。

根据计算的结果综合评价值基于MapGIS软件生成岩溶水文地质半定量分析图。其形成过程如下：

a各影响因素综合权重值的确定。根据多名水利专家所建议各影响因素的权重值，求得各影响因素综合权重值，在岩溶地下水富水性、岩溶地下水防护条件、岩溶发育强度等各主要影响因素的分区图中给各分区赋予因素值及综合权重值。

b基于MapGIS图层叠加功能，把各分区图叠加成一个综合图层。

c导出综合图的属性值（各分区图的因素值及权重），利用综合评价模型，计算出综合图中的各小分区的综合评价值（R），再导入到叠加图层中。计算结果在0～7之间，将其分成7个区间，每个区间赋予一个颜色并在综合图中表示出来，从而得到岩溶水文地质半定量分析图，如图313所示。

图313 岩溶水文地质半定量分析示意图

4）评价结果及分析。从岩溶水文地质半定量分析图上可看出，原观测网大多数观测孔分布在评价值高的区域，评价值在3以上的区域就分布有41个观测孔（表35），评价值越高，布孔率也就越大。

表35 原观测网各区间的分布面积及布孔率

（2）桂林市区岩溶地下水位观测网优化设计定量分析

1）岩溶水文地质条件分析。在岩溶地区，地下水含水介质的强非均质性给地下水观测网的布局和优化方法的选择上带来一定的困难。经过前面对研究区水文地质条件分析总结，克立格法适用于桂林市区岩溶地下水观测网的优化设计，其依据如下：

a观测孔集中分布在峰林平原或漓江两岸，其孔深都在100m以内，而在其深度范围内岩溶发育强烈，溶洞之间通过裂隙、微裂隙相互连通，从而使岩溶地下水保持密切的水力联系，构成了一个统一的岩溶含水层。

b研究区地下水系统含水介质以岩溶裂隙、溶洞为主，局部发育管道，具非均质性，但不强；地下水运动方式以散流为主，水的流态以层流为主，呈各向异性，但也不强。在研究区东区，广泛发育的溶洞、裂隙构成“地下溶洞-裂隙网络”，使岩溶含水层表现出类似于多孔含水介质的特征，可近似简化为相对均质各向异性含水层。

综上所述，研究区地下水系统具有非均质性及各向异性，但均不强，在局部范围内甚至表现出均质性，因此，克立格法适用于桂林市区岩溶地下水观测网优化设计的研究。

2）克立格法在桂林市区岩溶地下水观测网优化设计的应用。

A桂林市区各岩溶地下水水位观测孔坐标及水位值。桂林市区岩溶地下水水位观测网各观测孔1990年9月水位均值如表36所示。

表36 1990年9月桂林市区岩溶水位观测网各观测孔水位均值

B实验变差函数的计算。运用实验变差函数公式计算各方向实验变差函数值。经过分析，沿N—S、S—W、NE—SW、NW—SE4个方向分别计算的实验变差函数值相差不大，这也说明在观测孔分布范围内地下水系统含水介质非均质性不强。

由于观测孔分布的不规则性，可采用平均距离法计算实验变差函数。首先计算出每个观测点与其他观测点的距离hij，把它们分成8类，基于Matlab编程计算，算出的每类hij的个数、平均值，然后计算的平均实验变差函数 （h）如表37所示。

表37 原观测网每类hij的平均值及平均实验变差函数值

根据每一类的hij平均值及对应的平均实验变差函数值 （h），作出实验变差函数曲线，即图314中的折线。

图314 变差函数曲线拟合图

C实验变差函数的拟合。运用加权线性规划法对变差函数进行自动拟合，选用的球状模型为

岩溶地区地下水与环境的特殊性研究

将球状模型变差函数的拟合问题转化为多元线性回归问题。

根据加权线性规划法，运用Matlab编程计算得到b0=186、b1=061、b2=002；然后计算得到：C0=186，a=317，C=129，验前标准差C（0）=C0+C=315。

表38 加权线性规划法X1j、X1j、ωj值

续表

注：A为放大系数，在拟合时可根据人机互动方式设定。

根据表38计算结果，可得最佳拟合的理论变差函数为

岩溶地区地下水与环境的特殊性研究

D估计误差标准差及水位估值的计算。变差函数模型确定后，即可运用克立格方程组计算各观测孔估计误差标准差值和水位估值。

首先应用普通克立格方程组计算估计误差标准差值及水位估值。在计算中，由于计算时克立格系数出现了负数值，因此运用改进的克立格模型进行计算。从图313可知，当各观测孔两点间距离大于变程317km时，变差函数γ（h）几乎不变。这说明当两点间距离大于变程时，观测网内其他观测孔对估计点的估值有影响，但影响不大，可忽略不计。为了计算简单且符合精度要求，取距离小于变程范围内的所有观测孔参加克立格权系数和估计误差方差的计算。下面详细介绍观测孔GⅢ1、GⅢ63的克立格权系数和估计误差方差计算过程。

aGⅢ1。经过计算可知，有GⅢ3、GⅢ32、GⅢ35与GⅢ38观测孔与GⅢ1距离小于317km，因此只有它们4个观测孔参与计算，它们之间变差函数值及它们与GⅢ1观测孔间的变差函数值如表39所示。

表39 变差函数γi19、γij和权系数λi及拉格朗日乘子μ

利用表中数据，计算GⅢ1估计误差方差方程：

岩溶地区地下水与环境的特殊性研究

GⅢ1的水位估计值为

岩溶地区地下水与环境的特殊性研究

水位估计值与实际水位值相差003m。

bGⅢ63。依据变差函数计算公式，计算得到与GⅢ63距离小于变程所有观测孔之间的变差函数值和这些观测孔及GⅢ63的变差函数值的结果如表310所示。

表310 变差函数γi63、γij和权系数λi及拉格朗日乘子μ

根据表中变差函数值，计算GⅢ63水位估计误差方差为48852，把权系数代入水位估计公式求的其水位估计值为14495m，与实际值相差090m。

原岩溶地下水水位观测网其他各观测孔计算过程与GⅢ19、GⅢ63观测孔相同，各观测孔水位估计误差标准差值及水位估值如表311所示。

从表311中可看出，研究区地下水观测网平均估计误差标准差为39820。因此，给定估计误差标准差的临界值为39820。根据估计误差标准差等值线图，在水位估计误差标准差较大的区域，适当增加观测孔；反之，在水位估计误差标准差较小的区域则减少观测孔。然后再拟定几个优化方案，并计算各方案下平均估计误差标准差值和水位估值，经过精度比较和费用分析，选择其中最佳方案。

表311 原观测网各观测井估计误差标准差值

注：Δh为估计水位值与实际偏差表中，由于观测孔GⅢ21、GⅢ71周围317km内没观测孔，不能参加统计计算。

此外，由计算得到的各观测孔水位估值可知，在观测网布设密集区，水位估值较为准确，而在一些周边或水位变化较大的区域，由于分布观测孔较少，使得估计值不准确，甚至偏差很大。如观测孔GⅢ41，其变程范围内只分布有观测孔GⅢ45，故估计值即为GⅢ45的水位值，与实际相差1115m。在观测孔GⅢ13变程范围内分布的GⅢ11、GⅢ14、GⅢ23、GⅢ25等16个观测孔观测水位都比观测孔GⅢ13要大，其中观测孔GⅢ57水位值最小，仍比观测孔GⅢ13大072m，由于采用的克立格空间插值估计法是一种线性加权法，其算得的估值自然不会比参与计算的观测孔中的最小水位值小，即其估计水位值与真实观测水位值之间的偏差应在072m以上。因此，应在岩溶地下水位变化较大的观测孔周围布设观测孔以更全面地提取研究区水文地质信息量。

根据计算得到的原观测网岩溶地下水水位估计误差标准差值及水位估值，运用克立格空间插值技术，基于MapGIS软件生成估计误差标准差等值线图等值线图，如图315、图316所示。结果显示，在MapGIS中运用克立格空间插值法生成的估计误差标准差等值线能真实的反映实际情况。奇峰镇、研究区东北角处观测孔与周围观测孔距离太远，不能计算出估计误差标准差值，故缺少等值线的分布，需要增加观测孔以更全面反映整个桂林市区岩溶地下水水位动态；在岩溶所、湘桂铁路附近，观测孔分布较多，估计误差标准差值也小，可以减少一定量的观测孔。

图315 原岩溶地下水观测网水位估计误差标准差等值线示意图

从应用估计水位所作的地下水流场图与原地下水流场图可看出，在朝阳处流场发生了很大的变化，由于布置的观测孔少，GⅢ41 与GⅢ45 观测孔只有两者之间相互估计对方的水位，造成水位估计误差过大，从而使得估计的地下水流场不能真实反映实际情况，这也说明了原观测网布置得不合理。

E拟定被选方案。由于桂林市区岩溶地下水位观测网布置还相对较稀疏，因此在适当减少少量观测孔的条件下，以调整为主。将原观测网水位估计误差标准差等值线图和岩溶水文地质半定量分析图合成一个综合图，在水位估计误差标准差较大的区域，经过水文地质条件分析，在其区域内评价值较大的分区里增加观测孔；估计误差标准差较小的区域内，则在该区域内评价值较小的分区里减少观测孔，根据这个原则，拟定出两个优化方案；然后运用改进克立格法计算各方案的估计误差标准差值和水位估值；最后根据费用分析及精度比较，取其中较优者。

图316 原岩溶地下水观测网水位估计等值线示意图

（根据1990年9月水位资料）

方案1：

a分布在湘桂铁路处各观测孔估计误差标准差值都较小，其中最大的估计误差标准差为37941，因此根据岩溶水文地质评价分区图在该处评价值1～2减少观测孔GⅢ47、GⅢ48，保留与GⅢ47提供相同水位信息的观测孔GⅢ46；在2～3区间内减少观测孔GⅢ42，由于该处地下水位等值线分布较密集，水力坡度大，是地下水径流区，故还保留GⅢ24、GⅢ42两个观测孔。长海机械厂及其附近仅分布观测孔GⅢ6，其估计误差标准差较大为55694，在评价值为6～7的区间内增加观测孔Z1；在岩溶研究所及市区中心漓江两岸虽为地下水排泄区，但布置的观测孔很较多，估计误差标准差值很小，故可按一定的距离间隔减少分布在评价值大的区域内的观测孔GⅢ53、GⅢ54、GⅢ55、GⅢ59、GⅢ60、GⅢ68。

b拓木镇为地下水径流区，在评价值高的区域内增设观测孔Z2；瓦窑为人口密集区，发生的岩溶塌陷灾害的密度较大，而且估计误差标准差值也大，因此在评价值为5～6的区域新增观测孔Z3；朝阳北侧为地下水排泄区，该区域估计的水位流场已不能真实的反映实际情况，估计误差标准差值也大，故在评价值为5～6的区间内增加观测孔Z4。

c观测孔GⅢ20、GⅢ21布置的距离周围观测孔距离较远，因此将其稍作挪动，使其能参与估计误差标准差值和水位估值的计算；GⅢ29、GⅢ71观测水位与其他观测孔相差较大，为了更全面反映水文地质信息，在两者中间塘家湾处增加观测孔Z5。

此外由于观测孔GⅢ17与GⅢ65之间、GⅢ31与GⅢ64之间距离很近，水位值及水位动态变化情况也大致相同，因此减少分布在评价值较低的观测孔GⅢ17、GⅢ31。本方案在原观测网基础上减少观测孔11个，挪移位置2个，新增观测孔5个，优化后共43个观测孔，方案1中桂林市区岩溶地下水水位观测网分布情况见图317。

图317 方案1岩溶地下水观测网水位估计方差标准差等值线示意图

调整优化后，对各观测孔进行估计方差标准差及水位估计的计算。由于球状模型描述的是区域化变量结构的变化，而是在原观测网基础上进行的局部调整，整个观测网空间结构变化不大，因此调整优化后的观测网仍用原观测网拟合的球状模型。根据调整后的观测网各观测孔的坐标及水位值，运用Matlab编制的改进克立格程序模型计算的估计误差标准差值及水位估值如表313所示，并基于MapGIS软件生成的估计误差标准差等值线、估计水位等值线图如图317和图318所示。

图318 方案1岩溶地下水观测网水位估计等值线示意图

（根据1990年9月水位资料）

方案2：

在方案1中，GⅢ41、GⅢ65、GⅢ71、GⅢ20估计误差标准差值仍较大，因此方案2在方案1的基础上，在这四个观测孔附近综合评价值较高区域内增加Z6、Z7、Z8、Z9观测孔。计算结果如表312所示，根据计算结果所生成的估计误差标准差及估计水位等值线图如图319，图320所示。

F备选方案分析与确定。根据原观测网、方案1和方案2的计算结果和估计误差标准差、估计水位等值线图可知，虽然优化后减少的观测孔数目不多，但平均估计误差标准差有所增大，而估计水位值与实际水位值的平均偏差却有所减少。最主要原因是原观测网大多观测孔集中布设在市区中心地带，周边分布的观测孔很少，有些甚至不能进行统计计算，这就使得原观测网平均估计误差标准差值较小。而方案1和方案2是针对整个研究区进行优化调整的，在观测孔集中区减少观测孔，在周边及水位变动较大的地区增加观测孔以全面反映岩溶地下水水位信息，使得优化后所有的观测孔都参与了统计计算，观测网的观测范围也变大了，而观测孔数目又有所减少，所以估计误差标准差也相应地增加；由于在水位变动较大的区域增加设了观测孔，使得这些观测孔估计水位与实际水位偏差变小，不会出现水位估计有较大的偏差情况，所以优化后观测网平均水位估计偏差值比原观测网有所减少。

图319 方案2岩溶地下水观测网水位估计误差标准差等值线示意

优化后的观测网在顺沿地下水流和垂直地下水流方向上都布设了观测孔，在地下水补给、径流、排泄区都设有观测孔，在漓江河谷排泄区减少提供冗余信息的观测孔，在水力坡度较大区域增加了观测孔，因此在空间布局上优化方案较原观测网合理。从流场图上分析，原观测网估计水位流场已与原观测网水位流场变化较大，而优化方案估计的水位流场更符合实际。

图320 方案2岩溶地下水观测网水位估计误差标准差等值线示意图

（根据1990年9月水位资料）

表312 方案1和方案2计算结果

续表

原观测网、方案1和方案2的精度及费用分析比较如表313所示。

表313 方案1与方案2的比较

注：新建一个井约需8万元，每个井运行管理费及修理费每年按工程投资的5%计，各方案的设计年限为10年。

由表可知，方案1虽然可以节约24万元，但平均估计误差标准差相对较大，比原观测网增加158%，说明优化后整个观测网的布置很稀疏，达不到精度要求；方案2 是在方案1的基础上增加了4个观测孔，虽然费用只节省8万元，但观测精度提高了很多。根据桂林市区现有观测网布置密度相对较稀疏的实际情况，为了保证观测精度，优选方案2。

一、矿体的空间形态及其地质概念模型

阿舍勒铜锌矿床范围内的地质特征在第三章中已叙述，矿区内的地层及主体构造线呈近南北向展布，以F16断层破碎带为界，东西两侧分别为Ⅱ号和I号矿化蚀变带，一号铜锌矿床分布于Ⅰ号矿化蚀变带内。在构造上，一号铜锌矿床分布于矿区4号次级同斜倒转向斜的回转端及两翼，呈近南北向延伸，与构造线一致。赋矿地层岩石为中泥盆统阿舍勒组第二岩性段中亚段中上部的蚀变火山碎屑岩：主要是硅化的石英角斑质凝灰岩、沉凝灰岩，含角砾凝灰岩。工业矿体集中分布于矿区北17勘探线至南16勘探线之间（图8-1），其产出严格受地层和向斜构造的控制，与中酸性火山（沉积）岩有关，以火山喷气—沉积成因为主并经历了后期变形变质热液迭加改造形成的块状硫化物铜锌矿床。经勘探工作，已探明铜储量达大型规模，共生锌达中型规模，共生硫储量达大型规模，此外该矿床还伴生有Au、Ag、Pb、Ga、Cd、Se等多种有益组分，可综合利用。

图8-1 阿舍勒铜锌矿床勘查工程布置

（据新疆地矿局第四地质大队资料编制）

Q—第四系； —齐也组； —阿舍勒组第二岩性段上亚段； 阿舍勒组第二岩性段中亚段； —阿舍勒组第二岩性段下亚段；λ φ—钠长石英斑岩；Sq—次生石英岩；1—细碧岩；2—铁帽；3—黄铁矿化绢云母化强硅化蚀变带；4—地质界线和不整合界线；5—岩性段和岩性亚段界线；6—断裂；7—勘探线及孔位；8—Ⅰ号矿体地表水平投影范围

一号铜锌矿床共由三个矿体组成，编号分别为Ⅰ、Ⅱ、Ⅲ号矿体。矿体数量少但单体规模大，是本矿床的特点。Ⅰ号矿体为本矿床的主矿体，形态严格受地层、褶皱构造控制；Ⅱ、Ⅲ号矿体形态则受与I号矿体同生成矿期的同生裂隙的控制。矿体在空间展布上有侧伏现象，侧伏方向为北北东，侧伏角为45—65。

I号矿体分布于18—17勘探线间，总体南北向展布，为半隐伏—隐伏矿体，构造上位于4号向斜回转端内侧及两翼。矿体上覆地层为阿舍勒组第二岩性段上亚段（D2as2-3）组第二岩性段中亚段（D2as2-2）的蚀变火山—沉积碎屑岩。主矿体呈层状、似层状或透镜状，产于细碧岩与石英角斑质火山碎屑岩之间的接触界面上，在水平断面上形态呈似月牙状，横断面上呈鱼钩状，与上下地层整合接触，同步褶皱，形态严格受向斜构造的控制，呈向北倾伏、向南扬起、东翼向西倒转的紧闭向斜形态（图8-2）。

Ⅱ号矿体产于阿舍勒组第二岩性段中亚段（D2as2-2）中上部位，位于I号矿体之下，矿体近南北向展布，地表与Ⅰ号矿体平行分布，深部与其倒转翼及回转端部位斜交，逐渐向东偏离Ⅰ号矿体，并呈逐渐尖灭的趋势。

Ⅲ号矿体赋存在阿舍勒组第二岩性段上亚段（D2as2-3）的玄武熔岩中，隐伏于1—5线附近，呈一脉状展布。

为了从三维空间的角度真正的体现矿区的地质形态，需要建立矿区地质三维概念模型。模型的建立与矿区大量的基础地质工作分不开，特别是在矿区投入的大量钻探、槽探及硐探等山地工程，不仅为矿床储量的探明提供基本条件，而且也有效地控制了矿体的空间赋存状态，并详细查明了一号铜锌矿体的成矿地质条件、控矿因素、矿床规模及矿体的空间赋存规律。研究区位于16勘探线至13勘探线之间，区内施工有60个钻孔及部分探槽，有效地控制了矿体的空间形态及展布，就为建立矿区的地质三维概念模型提供了详实的基础数据。

1矿体三维概念模型

矿区内的矿体共三个，其中Ⅲ号矿体只在2—5勘探线间分布，呈脉状。Ⅱ号矿体，分布于16—1勘探线之间，由于受同生裂隙的控制，其在空间的三维模型呈一板状矿体，并逐渐往深部尖灭。在4—8勘探线间，Ⅱ号矿体近南北向展布，8—16线之间，矿体走向转向北西。

图8-2 5号勘探线矿体分布

（据新疆地矿局第四地质大队资料编制）

1—细碧岩；2—绢云母化凝灰岩；3—多金属重晶石矿石；4—多金属矿石；5—铜锌黄铁矿矿石；6 含铜黄铁矿矿石；7—块状黄铁矿矿石；8—浸染状黄铜黄铁矿矿石；9—浸染状黄铁矿矿石；10—角砾凝灰岩；11—绢云母化硅化凝灰岩；12—硅化角砾凝灰岩；13—大理岩；14—次生石英岩；15—硅化凝灰岩；16—角斑岩；17—断层及推测断层；18—钻孔及编号

Ⅰ号矿体为主矿体。在空间上，总体形态为一倒转的向形，严格受地层及4号向斜构造的控制（图8-2）。褶皱变形后的形态为一东翼向西倒转，向北倾伏的紧闭向斜形态。以4勘探线为界，矿体产状南缓北陡。4勘探线以南，向斜枢纽向南扬起，先是矿体的倒转翼部分出露地表，在10勘探线以南矿体正常翼也相继出露地表，在18勘探线附近矿体的回转端也完全出露地表。4勘探线往北，向斜枢纽向北倾伏，而且倾角较陡，呈隐伏状，且隐伏深度逐渐加大。矿体在其回转断部位受强烈的褶皱变形作用致使厚度明显膨大，形成典型的顶厚褶皱形态。Ⅰ号矿体两翼明显不对称，东翼矿体埋藏相对较浅，但倾斜延伸较大，呈稳定的厚层状。矿体的连续性较好，沿矿体倾斜方向厚度逐渐增大，变化较有规律。西翼矿体较东翼矿体平均埋深大，矿体沿倾向、走向延伸不稳定，变化大，呈复合分支状，从16线剖面线至13剖面线，矿体形态由紧闭向形变为吊钩状。垂直剖面方向，矿体大致呈向北倾伏（图8-3）。

2地质三维概念模型

本次研究的地质三维概念模型的范围主要局限于16勘探线至13勘探线之间。本区出露的地层有齐也组（D2q1）粗火山碎屑岩和中酸性火山熔岩、阿舍勒组第二岩性段上亚段（D2as2-3）、中亚段（D2as2-2）、下亚段（D2as2-1）中酸性火山碎屑岩为主体的火山—沉积岩夹碳酸盐岩建造（λφD2）、强硅化的细碧岩（sq）以及阿舍勒组旋回钠长石英斑岩（图8-4）。

齐也组只出露在9—13勘探线之间，为一半圆锥状，平行不整合于阿舍勒组第二岩性段上亚段，在研究区内，出露范围小。

阿舍勒组第二岩性段上亚段，分布于12—13勘探线之间，为向斜的核部，向南扬起，向北倾伏，轴向从13勘探线往12勘探线由南北向转至北北西向。地质体由南向北逐渐增大，也逐渐向下延伸，厚度增大。

图8-3 阿舍勒铜锌矿区Ⅰ号矿体联合水平断面

（据新疆地矿局第四地质大队资料编制）

1—铜硫矿石；2—硫矿石；3—铜锌矿石

阿舍勒组第二岩性段中亚段地层，在平面上走向为近南北向，为4号向斜的翼部，分布于细碧岩两侧。在剖面上，西翼为正常的地层，倾角相对较缓，而东翼为倒转翼，往南扬起，向北倾伏，在向斜的转折端增厚。在垂直剖面方向，呈现向北倾伏的趋势，空间上连续性较好。

图8-4 阿舍勒铜锌矿区地质三维概念模型

（据新疆地矿局第四地质大队资料编制）

D3q1—齐也组；D2a s2-1—阿舍勒组第二岩性段上亚段；D2as2-2—阿舍勒组第二岩性段中亚段；D2as2-2—阿舍勒组第二岩性段下亚段；λφD2—钠长石英斑岩；Sq—细碧岩；Ⅰ、Ⅱ—矿体代号；16—勘探线编号

阿舍勒组第二岩性段下亚段地层。从平面上东翼分布范围较西翼大，基本呈长条状分布，受断裂构造的影响，但大致与地层的倾斜方向一致，构成4号向斜翼部的最外层。西翼受断裂控制，出露范围较小。从剖面上看，东翼在空间上较连续，呈倒转，而西翼为正常翼。在垂直剖面方向，地层向北倾伏，向南扬起，并有向倾伏方向增大的趋势。

由于强硅化的细碧岩为矿体的主要赋存空间，故将其单列为一层。其在空间上也是由南往北倾伏，向南扬起，紧贴着上亚段，西翼由8线起减薄，但下部变大；东翼由4线起逐渐减薄，到9线后紧贴上亚段，往下部逐渐增大。空间上基本呈一倒转向斜。Ⅰ号、Ⅱ号矿体主要赋存在其中，连续性较好。

阿舍勒组钠长石英斑岩（λφD2）。分布于8—13勘探线之间，主要侵入中亚段地层。在平面上受风化剥蚀和地形的影响，与上亚段地层的形态相似，由南往北逐渐增大。在空间上，岩体呈葫芦状，两头大，中间小。剖面上，由8线起往北，岩脉逐渐增大。从垂直剖面方向，岩体向北倾伏，向南扬起，连续性较好。

总之，由于受4号紧闭倒转向斜的影响，地质体在空间上形成紧闭倒转褶皱，呈现向南扬起，向北倾伏的趋势。

二、阿舍勒铜锌矿床三维立体模型的建立

1空间数据库的建立

空间数据库是三维模型建立的基础。空间数据库是指在地球表面某一范围内与空间地理相关的，反映某一主题信息的数据集合。与传统的数据库既有相同点，又有其自身的特点。作为数据库，它必须满足如系统的灵活性、数据输入与更新的标准化、系统的有效性检验、数据的安全性、最小冗余等原则，而其特点表现为数据包含三维空间的信息。三维模型系统功能的好坏，直接与空间数据库的合理组织有关的，由于空间数据都是密切相关的，因此有必要将它们合理地组织起来。

矿区数据资料收集的完整性及准确性关系到空间数据库的建立，也影响到三维模型是否符合实际情况，因此充分收集原始资料是首要的任务。为了建立空间三维模型，MRES软件系统需要使用矿区的山地工程资料，包括：钻孔、槽探、平硐、沿脉、穿脉等，为了形象、方便及简单的描述，以下均将山地工程简称为钻孔。所需原始资料主要包括三类：一是钻孔的空间总体位置信息，即钻孔的测量数据，包括钻孔在三维空间的起点坐标（X,Y,Z）以及钻孔的长度；二是钻孔在空间的位置变化信息，即钻孔在空间的倾斜方向和倾角，这两个关于钻孔空间位置信息的资料描述了钻孔在空间的形态；三是对钻孔的 *** 作及有关的地质描述，即采样信息：包括采样位置、样品代号、分析结果、样品长度，地质描述则是在充分了解钻孔揭露的地质内容，划分钻孔揭露的地质界线，包括岩性代号及地质代号。表8-1为钻孔的空间总体位置信息，表8-2为钻孔在空间的位置变化信息，表8-3为采样信息及地质描述。

表8-1 钻孔的空间总体位置信息

表8-2 钻孔在空间的位置变化信息

续表

表8-3 采样信息及地质描述

续表

续表

续表

续表

续表

续表

续表

续表

续表

当收集了原始资料后，需要进行空间数据库的建立。首先是建立数据表，这是数据交换的前提。数据表包括多项记录，而每一项记录又包括多个信息段。第二步是添加记录，图8-5为记录界面示意图，每一个记录本身包含一些信息，包括记录的名字及其描述、记录的内容、记录数据输入的次数等。表8-4为记录内容表。

表8-4 记录内容

图8-5 记录界面示意

其中DCOLLAR记录是钻孔的空间总体位置信息，DSAMPT记录是钻孔空间的位置变化信息，DSURVEY记录为采样信息及地质描述。

第三步是添加信息段（字段）。信息段的界面如图8-6，每一个信息段包括以下信息：

图8-6 信息段界面示意

（1）信息段的名称及描述；

（2）数据类型有：字符型、整型、浮点型、双精度型；

（3）字符信息包括：起始位置、终点位置、字符宽度；

（4）小数点位数；

（5）属性特征（每个信息段有一个别名，代表信息段的属性）；本次研究使用的记录及信息段具体内容见表8-5。

当原始资料及数据表创建完成后，即可建立数据库文件（如test），然后导入相关的数据信息，则空间数据库的建立工作完成了。

表8-5 信息字段内容

2钻孔三维可视模型

如上所述，钻孔是山地工程的一种简称，它包含有孔口三维坐标（X,Y,Z）以及钻孔的长度、钻孔的三维空间变化信息（空间延伸方向和倾角）。运用MRES软件系统，根据这些三维空间信息，可以创建钻孔的三维模型，图8-7为铜元素钻孔的三维模型。要显示钻孔，就要为钻孔在屏幕确立一定的显示空间，即建立矿区在屏幕显示的三维空间（X,Y,Z的区间）框架，主要依据矿区的公里网坐标：东西方向（X）范围为29450500—29452000，南北方向（Y）范围为535000—5352000，高程范围（Z）为100—1000。则钻孔的三维模型显示在屏幕的中心，同样，下文要描述的矿体、地质三维模型均将显示在这个范围内。

为了显示钻孔及元素空间分布，需要首先建立显示图例。矿区主要进行Cu、Zn、S、Ag等四种元素的研究，建立了Cu、Zn、S、Ag四种元素的图例。图例的建立，主要依据矿区元素的分析成果、边界品位、工业品位等信息划分五个区间，如：Cu元素分－999—00,00—03,03—05,0 5—10,10—999;Zn元素分－999—00,00—05,05—10,10—30,30—999;S元素分－999—00,00—80,80—120,120—250,250—999,Ag元素分－999—00,00—10,10—400,400—800,800—1500,150—9999999。每一个区间给以一定的颜色表示。这样就可显示样品在空间分布情况，也可以显示地质描述特征。而通过颜色的不同，大致显示了元素在空间富集的形态。

图8-7 铜元素钻孔的三维模型

3矿体三维模型

（1）矿体三维模型的建立：

使用MRES单工程自动圈定功能，输入边界品位和最低可采厚度，应用单工程原理进行矿体边界的划分。同时在原采样信息表中增加一个字段，标明该段是否是矿体。

（2）剖面编辑：

采用MRES剖面编辑CAD技术进行剖面矿体地层的圈定，通过交互方法输入地质体边界线，并通过放大、缩小和移动数据点等功能进行精细编辑，在剖面编辑形成的编辑都有拓扑关系属性。在圈定过程中，不同品位的矿体可以用不同的颜色表达。（图8-8）

图8-8 矿体剖面圈定界面

（3）曲面连接和体生成：

每条剖面上的矿体圈定完成后（图8-9），就要进行剖面间矿体的空间连接，即将每个剖面上的矿体形态连接成一个整体的三维空间模型。当矿体由于形态特别复杂无法自动的拟合成完整的空间三维模型时，就需要人为进行交互 *** 作，即给系统一个连接三维模型的趋势，来引导系统完成模型的建立（图8-10）。本次研究首先建立了铜矿体的三维模型。

图8-9 矿体剖面圈定结果

图8-10 矿体剖面连接示意

（4）铜矿体三维模型简介：

由于金属矿体呈不规则形态，因此，矿体的圈定主要依据钻孔中样品元素分析结果并结合地质规律。本次共圈定了2个铜矿体，Cul为Ⅰ号主矿体的主要部分，Cu2为Ⅰ号主矿体在深部的分支脉体，并分布于2勘探线至13勘探线之间，Ⅰ号主矿体在空间上呈一倒转向形，东翼长而厚，西翼短，回转端部位矿体厚度明显增大。矿体总体向北倾伏，向南扬起（图8-11,8-12）。

图8-11 矿体三维显示图

图8-12 矿体联合剖面圈

在三维空间上，Cu 1矿体在16线—12线呈两翼大致相等，而到8线发生突变，变成东翼长西翼短。Cu l矿体向北恻伏，在矿体南部分布标高在500米上，而在北部矿体在500～50米之间。Cu2矿体主要分布于16勘探线—8勘探线之间，呈一板状，紧贴在Cul矿体的东翼，并向下逐渐尖灭。

4地质体三维模型

在充分研究前人工作的基础上，对矿区的地质情况进行研究，分析地质体的空间展布，开始地质体的三维建模。本次研究选择的是16～13勘探线之间的地质体，主要包括：齐也组（D3q1）、阿舍勒组第二岩性段上亚段（D2as2-3）、中亚段（D2as2-2）、下亚段（D2as2-1）以及矿体。借鉴矿体三维模型建立的方法，同时又考虑到不同地质体间的连续性界面，利用“层”的方法，即不同的地质体为不同的层，层与层之间的面是共同的面，下层地质体的上顶面是其上一层的下底面。从剖面入手，在剖面上圈定地质体的界线；然后，同一界线连接成一过渡层，两个不同的过渡层，即一层上顶面和一层下底面通过建立另一过渡层——墙；最后，利用这三个过渡层建立地质体的三维空间模型。逐个的建立地质体，完成矿区三维空间模型的建立。对于地质体在空间的显示，考虑到地质体不像矿体有化学分析结果，因此采用给地质体变量赋予一定的缺省值的办法来达到显示地质体的目的。

本次所建立的地质三维空间模型能较好地反映该地区的地质情况，但与该区的地质体三维空间形态仍存在一些区别，一是该地区的地质情况非常复杂，要完全模拟还存在一些不小的困难，这与对地质的认识程度和简化地质体的分层有关系；二是钻孔主要控制的是矿体的空间形态，对地质体的控制程度不够。在进行实验时为了能在计算机上描述，进行了一些简化工作。

从地质体三维模型整体上看（图8-13），由于矿体及矿区地层的空间展布严格受构造（特别是4号向斜构造）的控制，地质体三维模型在空间中主要形成向斜的空间形态，无论从平面还是从三维模型的两侧剖面（16、13勘探线剖面）来看，由向斜核部到翼部依次为：齐也组（D3q1）、阿舍勒组第二岩性段上亚段（D2as2-3）、矿体、中亚段（D2as2-2）、下亚段（D2as2-1）。

（1）三维地质体显示：

图8-13 地质体三维立体

图8-14 DEM数据与遥感及三维数据的复合显示

①Ⅰ号矿体主要赋存于中泥盆统阿舍勒组第二岩性段中亚段（D2as2-2）、上亚段（D2as2-3）、强硅化细碧岩中并靠近中亚段一侧，矿体呈半隐伏、隐伏状。矿体形态严格地层、向斜构造的控制，褶皱变形后为东翼向西倒转，西翼正常，剖面上呈鱼钩状，深部出现分支现象。②地层严格受褶皱构造的控制，形成倒转的紧闭向斜形态，阿舍勒组第二岩性段是构成褶皱的岩性段。③阿舍勒组旋回钠长石英斑岩侵入向斜倒转翼的阿舍勒组第二岩性段中亚段地层中，并靠近上亚段。

（2)DEM数据与遥感及三维数据的复合显示。利用DEM数据及配准好的遥感数据可以和三维矿体地质体构成真三维图像，方便从地上到地下全三维显示。

三、储量计算

1概述

（1）工业指标的确定：

阿舍勒铜锌矿区储量计算的工业指标主要参照《新疆哈吧河县阿舍勒铜矿区一号铜锌矿床勘探地质报告》中采用的工业指标，见表8-6。

表8-6 矿床工业指标

（2）储量计算方法的选择：

软件系统提供了进行储量计算的方法，主要利用地质统计学的方法进行品位估算和储量的计算。在本次实验中主要采用地质统计学的方法计算阿舍勒铜锌矿区的储量，主要计算了铜元素的储量。

2储量计算过程

（1）概述：

在软件系统中进行储量计算时，主要分以下几个步骤：

①样品组合：进行样品的归一化，创建组合样品信息数据库；②创建块段系统（BLOCK）：为每一个矿种创建其显示空间范围，所建的空间框架能包围矿体，并用超级块段及子块段来拟合矿体及其外边界；③变异函数：进行实验半变异函数的计算及理论变异函数的拟合，寻找块金值、基台值、变程、搜索椭球体的参数；④品位估算：依据理论变异函数拟合提供的结果，利用软件提供的计算方法，估算矿体的品位；⑤储量计算：根据品位估算的结果、矿体的体积，依据品位条件，计算矿体的矿石量、金属量以及可进行资源的评估。

（2）分析数据简单统计：

对原始铜数据可以进行简单的分析统计，得到铜平均值162、最大值148、最小值001。铜元素分析结果的累积概率曲线如图8-15。

图8-15 铜数据概率分布曲线

矿体铜元素分析数据PP图表明矿体铜元素含量分布较复杂，不服从标准正态分布（图8-16）。

（3）变异函数：

对区域化变量进行克立格估值需要研究它们的变异函数。MRES软件中提供了以下变异函数以供选择：半变异函数、对数半变异函数、偏相关半变异函数、相关半变异函数、指示半变异函数、正交半变异函数、协变异函数、协同变异函数。本次研究主要采用半变异函数进行计算，其表达式为：

地质勘探三维可视化技术及系统开发

式中：N(h）为滞后距为h时，参加实验变异函数计算的样品个数；h为滞后距。

图8-16 铜数据概率PP分布曲线

选择矿体在平面上及垂直方向进行实验半变异函数的计算。考虑到工程中样品不是严格按照网格采取的，为了充分利用所采集样品，在计算某一个方向a上的变异函数时，给出该方向上的角度容许误差限da，距离容许误差限dh。当计算某一个方向a上的实验变异函数时，以方向a为基准线，在角度范围a±da及距离范围h±dh圆锥形所夹区域内的所有样品都参加实验变异函数的计算。

实验变异函数计算结束时，得到了元素的实验变异函数散点图。根据这些散点图或是实验变异函数曲线图（通过连接散点得到）可以进行理论半变异函数的拟合。MRES提供了以下的理论模型：球状模型、指数模型、高斯模型、德文金斯模型、线形模型、周期模型等。依据散点图的散点分布情况选择拟合的理论模型。拟合过程完成后，将得到理论变异函数的曲线以及相关的信息：块金值、基台值、变程、搜索椭球体的轴参数和角度参数。

（4）铜矿体的储量计算：

采用沿水平方向和垂直方向的划分原则，在水平方向上：以南北为Y轴，东西方向为X轴，垂直方向上以高程方向为Z轴。自南：2945118539，至北：2945148779；自西：535071353，至东535144853;标高为300～91975；整个计算体积为302m×735m×916m。

①块段系统（BLOCK）原点大地坐标为：

x0=2945118539M;

yo=535071353M;

z0=300M；

块段的大小为：X方向为100M,Y方向为100M,Z方向为10OM。

块段的数目为：X方向为31,Y方向为74,Z方向为92。

块段的数目为211 048。但扣除矿体外的块段，实际参加品位估算的块段的数目为25 100个。

图8-17 铜元素变差函数计算界面

②变异函数计算时所选用的参数如下：

滞后距：20M

距离容许误差限：2M

分别计算了水平方向0°、45°、90°和垂直方向90°和315°方向0°、45°、90。和45°方向45°、90°的变异函数值，其实验变异函数曲线图8-18～图8-26（由于垂直方向的变异函数曲线大致相同只放一张）。数据值见附表8-7。

图8-18 铜数据实验变差函数曲线

图8-19 方向角45°，倾角0°变差函数拟合

图8-20 方向角90°，倾角0°变差函数拟合

图8-21 方向角135°，倾角0°变差函数拟合

图8-22 方向角0°，倾角90°变差函数拟合

图8-23 方向角45°，倾角45°变差函数拟合

图8-24 方向角45°，倾角90°变差函数拟合

图8-25 方向角315°，倾角0°变差函数拟合

图8-26 方向角315°，倾角45°变差函数拟合

表8-7 不同方向理论实验函数拟合结果

理论变异模型的拟合，在水平方向采用球状模型，在垂直方向用球状模型，拟合的参数见表8-8。

表8-8 理论实验函数拟合结果

③品位估算：铜矿的品位估算分别应用了普通克立格和距离倒数法，其参数主要利用了理论变异函数拟合的结果。在参与品位估算的样品数选择最少为一个样品，最多为10个样品。估算的品位存放与变量Cu-ivd 中。

④储量计算：分别用普通克立格和距离倒数法估算的品位进行铜矿的储量计算。其计算公式为：矿石量TONNAGE=VOLUMESG，金属量METAVALUE=TONNAGECu-ivd。铜矿的体重（SG）为36(t/m3）。

二种方法进行铜储量计算，其结果见表8-9。与勘探报告结果的对比见表8-10。

表8-9 铜矿不同方法储量计算结果比较

表8-10 储量计算结果对比

⑤储量管理。

3任意盘区储量统计

在储量计算过程中记录了盘区号，这样可以方便统计各盘区的储量。

4色块图分布

MRES提供储量盘区块段可视化功能，可以在三维空间观看储量分布，也可以进行切片分析（图8-27）。

图8-27 1号矿体储量估算Y方向联合剖面

5等值线分布

将切片结果存成GRD文件，然后调用等值线进行储量二维显示（图8-28）。

图8-28 1号矿体500米标高储量等值线分布