西藏甲玛铜多金属矿储量估算的新方法

距离为 ｄｉ，则式（２） 的最小化问题中的权值有如下
约束：
如果 ｄｉ ＜ｄｊ，则 λｋ，ｉ ＞λｋ，ｊ（ ｉ，ｊ ＝ １，…，ｍ） ，并且
Ｋｍ
������������λｋ，ｉ ＝ １。
（３）
ｋ＝１ ｉ＝１
２ 钻孔条带状效应分析
西藏甲玛铜多金属矿最主要的数据来源是沿着 钻孔获取到的数据，其中包括地面钻孔数据、钻孔化 验数据、小体重数据等。 原始数据为截止到 ２０１２ 年 甲玛矿区勘探工程得到的数据，包括 ５８ 个钻孔，共
但是，西藏甲玛铜多金属矿所处区域内的地层 比较复杂，地层曲面变化较大，甚至可能存在地层曲 面相交和错断的情况。 因此，运用协同克里格法计 算储量时，会 出 现 钻 孔 数 据 采 样 不 均、 呈 条 带 状 分 布［６ １１］ 的情况，即有限的条带两端数据权值偏大，导 致插值结果不合理。 在这种情况下，笔者采用拟牛 顿法［１２］ 对协同克里格插值得到的权值进行校正，在 估值方差基本不变的前提下使协同克里格估值的权 值更为合理，并且以西藏甲玛铜多金属矿的储量估 算为实例对该方法进行了验证。
５期
王琦标等：西藏甲玛铜多金属矿储量估算的新方法
·１００３·
３ 拟牛顿法权值校正
对权值系数进行约束需要满足式（３） 中的最小 化问题约束条件。 笔者采用拟牛顿法［１２］ ，其流程如 图 ５ 所示。
控制误差时，ｘｋ，１ ＝ ｘｋ，ｐ＋１ 即为最优解，Ｈｐ＋１ 可通过如下
计算公式求得：
Ｈｐ＋１
＝ Ｈｐ
－
Ｈｐ
ｙｐ
ｙ
Ｔ ｐ
Ｈｐ
ｙ
Ｔ ｐ
Ｈ
ｐ
ｙ
ｐ
＋
Ｓｐ ＳＴｐ 。 ｙＴｐ Ｓｐ
（４）
其中 ｙｐ ＝ ∇ｆ（ ｘｋ，ｐ＋１ ） －∇ｆ（ ｘｋ，ｐ ） ，Ｓｐ ＝ ｘｋ，ｐ＋１ －ｘｋ，ｐ ，计算
时的第一步迭代取 Ｈ１ ＝ Ｉ（单位矩阵）。
权值校正之后的协同克里格算法流程如图 ６
所示。
图 ６ 权值校正协同克里格插值算法流程
图 １ 地表及钻孔点三维模型
图 ２ 数据被断层边界截断 图 ３ 数据被搜索邻域截断
图 ４ 协同克里格引起的条带状效应
图 ４ 显示的这些权值体现了采用协同克里格估 值时的条带状效应：对于单个钻孔数据的估值，待估 点离钻孔数据最远的两端，其样品点所获得的权值 异常偏大，而离待估点较近位置的权值则普遍较小。 根据克里格估值方法中的屏蔽效应［１３ １８］ （ 即若块金 效应为 ０ 或很小时，则在待估区内的样本权系数最 大，稍远一点则其权系数显著减少），图 ４ 显示的权 值是不符合实际情况的，特别是当数据的分布在断 层或者搜索范围截出不平稳的数据时，就会出现端 点数据处的权值明显大于或小于中间数据的权值的 现象，所以需要对其进行校正。
西藏墨竹工卡县甲玛矽卡岩—角岩型铜多金属 矿床是目前我国同类矿床矿石品质最好、成矿元素 多、探明的资源量已经达到超大型规模的矿床［１］ 。 对于该地区的储量 估 算 采 用 过 传 统 断 面 法［２］ 、 ＳＤ 法［３］ 和地质统计学方法［４］ ，其中协同克里格法（ 地 质统计学方法中的一种） 最能反映该地区多种金属 之间的协同区域化［５］ 现象。
４ 处理结果
４．１ 数据预处理 采用权值校正协同克里格法估计矿床储量，首
先要选取区域化变量。 区域化变量必须具备随机性 和结构性的双重性质，即在受限的地质范围内既要 具有一定的相关性又要具有随机变化的统计特征。 由于甲玛矿床是 Ｃｕ、Ｍｏ、Ａｕ、Ａｇ、Ｐｂ、Ｚｎ 等的复合多 金属矿床，并且有大量的钻孔化验数据，所以选取金 属元素的含量作为区域化变量。 本次采用 Ｃｕ 含量 数据来进行说明。
表 １ 甲玛矿床 Ｃｕ 特异值处理结果
参数
均值 ／ ％ 标准差 ／ ％ 变化系数
特高品位代替前 ０．３２ １．０１４８ ２．１
特高品位代替后 ０．２１９３ ０．３７４５ １．７
甲玛矿区化学样的样长多为 １ ｍ 或 ２ ｍ，样品 长度各不相同。 在利用权值校正协同克里格方法进 行资源量估算时，要求所有参与计算的样品具有相 同的样长，即需要进行样品组合。 理论上，组合样品 的长度最好与实际采矿台阶的高度一致。 甲玛矿区 实际采矿台阶高度为 １０ ｍ，但考虑到该矿床为矽卡 岩型多金属矿床，矿体“ 分枝、复合” 现象比较明显， 而且地下开采部分夹石剔除厚度为 ４ ｍ，露天采矿 部分夹石剔除厚度为 ８ ｍ，因此，取 ６ ｍ 长度作为组 合样长，一方面可把夹石剔除，另一方面可把小范围 内（ 大于 ６ ｍ） 的变异也反映出来。 组合样中 Ｃｕ 特 征值的统计结果见表 ２。
０ １５ ５０ ０．５ ０．５ １５８
Ａｇ 次轴
０ － １０ ２０ ８０ ０．５ ０．５ １３２
短轴 ０ ８０ ２０ １０ ０．５ ０．５ ３２
图中所示，两种元素在各个方向上的变差函数
比较相似，因 此 推 断， 所 有 变 差 函 数 都 是 各 向 同 性
收稿日期： ２０１４⁃０６⁃０３ 基金项目： 国家自然科学基金重大科研仪器设备研制专项（４１２２７８０２） ；国家杰出青年基金（４１０２５０１５）
·１００２·
物 探 与 化 探
３９ 卷
计 １５ １８７ 件化学样。 矿床地表模型及钻孔点的三 维模型如图 １ 所示。
如果某一个钻孔遇到地层的错断情况，钻孔数 据就会出现边界截断，如图 ２ 所示。 另外如图 ３ 所 示，在进行协同克里格插值时首先会确定一个搜索 范围，该范围在三维空间内通常为一个椭球体的边 界截 断。 这 两 种 现 象 称 之 为 钻 孔 的 有 限 条 带 状 数据。
在矿产储量和品位计算之前，先对钻孔数据进 行预处理，包括特异值处理与样品组合处理。 首先 根据 Ｃｕ 的含量数据，在正态概率纸上作累积频率 曲 线（ 图７） 。从图上可以看出，Ｃｕ含量累积频率曲 线总体上为一直线，表明数据基本符合正态分布。
·１００４·
物 探 与 化 探
３９ 卷
图 ７ 甲玛矿床 Ｃｕ 含量累积频率曲线
表 ３ 组合样本中各元素与 Ｃｕ 的相关系数统计
元素 Ａｕ Ａｇ Ｍｏ Ｐｂ Ｚｎ
相关系数（ ｒ） ０．６３３８ ０．８８ ０．６５７９ ０．６３８２ ０．６１６６
进行权值校正协同克里格法插值的第一步是得 到 Ｃｕ 含量和 Ａｇ 含量的变差函数，以及两者的交叉 变差函数。 把所有钻孔点两两组合成点对，计算其 距离 值， 再 以 球 状 模 型 拟 合 理 论 实 验 变 差 函 数。 Ｃｕ、Ａｇ 三个变差函数轴的方向定位及其计算结果见 表 ４，变差函数见图 ８。
图 ５ 拟牛顿法校正权值流程
图中拟牛顿法的计算过程为：为了使距离最近 的已知点获得最大的权值，使最远的已知点获得最 小的权值，首先需要对距离 ｄ 按升序排列，然后用协 同克里格计算权值，在第 ｋ 个区域化变量的权值下， 给定初始迭代值 ｘｋ，０。 为求最优解 ｘｋ，１，通过下降方 向 Ｓｐ＋１ ＝ －Ｈｐ＋１ ∇ｆ（ ｘｋ，ｐ＋１ ） 得到迭代公式 ｘｋ，ｐ＋１ ＝ ｘｋ，ｐ ＋ Ｓｐ，其中 ｐ 为迭代次数，ｆ 函数为最小化方差计算公 式。 如式（２） 所示，当∇ｆ（ ｘｋ，ｐ＋１ ） 的范数达到给定的
表 ２ 甲玛矿床 Ｃｕ 样品组合处理的统计结果
元素 组合样数 最小值 ／ ％ 最大值 ／ ％ 平均值 ／ ％ 标准差 ／ ％
Ｃｕ ３５７０
０．００１４
９
０．８３０２
１．６１１６
４．２ 协同克里格法储量估算 ４．２．１ 变差函数及结构分析
对组合样本进行统计分析，从表 ３ 中可以看出， 组合样品中 Ｃｕ 含量数据与 Ａｇ 含量数据相关程度 最大，所以笔者以 Ａｇ 作为辅助变量对 Ｃｕ 的品位进 行估值。
在直线上部存在少量离散点，且偏离直线方向，这种 离散点才是真正的特异值。
本次估算采用频率曲线法进行特高品位代替， 即将曲线发生偏离直线方向处的含量值取为上截品 位，将所有高于上截品位的含量数据全部用上截品 位值替代。 笔者取 １．５％ 为上截品位，矿区 Ｃｕ 含量 高于 １．５％的样品共有 ７８２ 个，替代率为 ５．１５％。 如 果不对这 ７８２ 个样品含量进行替代，将会出现权值 普遍偏大的情况，导致插值结果不合理。 替代前、后 的统计分析见表 １。
估点 Ｚ（ｕ０）处的协同克里格的估计量为
Ｋｍ
������������ Ｚ∗（ｕ０） ＝
λｋ，ｉ Ｚ（ ｕｋ，ｉ ） ，
（１）
ｋ＝１ ｉ＝１
其中 λｋ，ｉ 是第 ｋ 个区域化变量中的第 ｉ 个样本点权
值，由最小化方差求出
σ２ ｕ０
ห้องสมุดไป่ตู้
＝
Ｅ｛ Ｚ（ ｕ０ ）
－
Ｚ∗（ ｕ０） ｝ ２
→ ｍｉｎ。
（２）
设样品数据点 ｕｉ（ ｉ ＝ １，２，…，ｍ） 与待估点 ｕ０ 之间的
西藏甲玛铜多金属矿储量估算的新方法
王琦标１ ，蒋鑫１ ，庹先国１，２ ，毛小波１ ，李怀良２
（１．成都理工大学 地质灾害防治与地质环境保护国家重点实验室，四川 成都 ６１００５９； ２．西南科技 大学 核废物与环境安全国防重点学科实验室，四川 绵阳 ６２１０１０）
摘 要： 西藏甲玛铜多金属矿的钻孔分析数据多呈条带状分布，采用协同克里格法对这种条带状数据进行插值时， 易出现端点数据权值偏大的情况。 笔者介绍了协同克里格算法的相关理论，以此为基础，对端点数据的权值进行 校正，使越远的数据点权值越小，并且给出了校正算法的程序代码。 以 ２０１２ 年甲玛矿区勘探工程的 Ｃｕ 含量数据 为例，分别进行了协同克里格法插值和权值校正的协同克里格法插值，实验结果表明，权值校正的协同克里格法能 取得更为良好的效果，对于西藏甲玛铜多金属矿的地质属性、储量估算等空间数据建模具有重要意义。 关键词： 储量估算；协同克里格法；条带状分布；权值校正；西藏甲玛铜多金属矿 中图分类号： Ｐ６３２ 文献标识码： Ａ 文章编号： １０００－８９１８（２０１５）０５－１００１－０７
１ 协同克里格算法原理
在用协同克里格法［１３ １８］ 进行储量估算时，对于 由 Ｋ 个变量构成的协同区域化变量 Ｚｋ（ ｘ），ｋ ＝ １，２，
…，Ｋ，有相邻 ｍ 个位于位置 ｕｉ ＝ （ ｉ ＝ １，２，…，ｍ） 处 的数据点，将中心点的品位数据设为 Ｚ（ ｕｋ，ｉ ） ，ｉ ＝ １， …，ｍ；ｋ ＝ １，…，Ｋ；Ｋ 为协同区域化变量个数。 在待
以 Ｃｕ 和 Ａｇ 协同区域化为例，对其中一列有限 的条带状数据进行协同克里格［１３ １８］ 插值计算。 估 值计算时 Ｃｕ 和 Ａｇ 的变差函数以及两种元素的交 叉变差函数均采用球状模型，通过计算可以得到块 金效应［１３ １８］（ 表征变差函数原点处的连续性） 为 ０ 的协同克里格权值系数，其结果如图 ４ 所示。
第 ３９ 卷第 ５ 期 ２０１５ 年 １０ 月
物 探 与 化 探
ＧＥＯＰＨＹＳＩＣＡＬ ＆ ＧＥＯＣＨＥＭＩＣＡＬ ＥＸＰＬＯＲＡＴＩＯＮ
Ｖｏｌ．３９，Ｎｏ．５ Ｏｃｔ．，２０１５
ｄｏｉ： １０．１１７２０ ／ ｗｔｙｈｔ．２０１５．５．２０ 王琦标，蒋鑫，庹先国，等．西藏甲玛铜多金属矿储量估算的新方法 ［ Ｊ］ ． 物探与化探，２０１５，３９（ ５） ：１００１ － １００７． ｈｔｔｐ： ／ ／ ｄｏｉ． ｏｒｇ ／ １０． １１７２０ ／ ｗｔｙｈｔ． ２０１５．５．２０ ＷＡＮＧ Ｑ Ｂ，Ｊｉａｎｇ Ｘ，Ｔｕｏ Ｘ Ｇ，ｅｔ ａｌ．Ａ ｎｅｗ ｒｅｓｏｕｒｃｅ ｅｓｔｉｍａｔｉｏｎ ｍｅｔｈｏｄ ｆｏｒ ｔｈｅ Ｊｉａｍａ ｐｏｌｙｍｅｔａｌｌｉｃ ｃｏｐｐｅｒ ｄｅｐｏｓｉｔ ｉｎ Ｔｉｂｅｔ［ Ｊ］ ．Ｇｅｏｐｈｙｓｉｃａｌ ａｎｄ Ｇｅｏｃｈｅｍｉ⁃ ｃａｌ Ｅｘｐｌｏｒａｔｉｏｎ，２０１５，３９（５） ：１００１－１００７．ｈｔｔｐ： ／ ／ ｄｏｉ．ｏｒｇ ／ １０．１１７２０ ／ ｗｔｙｈｔ．２０１５．５．２０
表 ４ Ｃｕ 和 Ａｇ 变差函数参数
参数
主轴
方位角 ／ （ °） ２７０
倾角 ／ （°） ０
容差 ／ （ °） １５
滞后距 ／ ｍ ５０
基台
０．４５
拱高
０．５５
变程 ／ ｍ １５０
Ｃｕ 次轴
０ － １０ ２０ ８０ ０．４５ ０．５５ １３５
短轴 ０ ８０ ２０ １０
０．４５ ０．５５ ２８
主轴 ２７０