micromine 基本原理与方法.pdf

第四章 基本原理与方法 第四章 基本原理与方法 4.1 矿床三维模型构建方法 4.1 矿床三维模型构建方法 运用计算机技术建立矿床三维模型的研究工作从六十年代为解决浸染状矿 床建模问题而采用三维块段模型以来，至今已经历了近四十年的发展。建模方法 也由早期简单的方块模型，发展到如今的实体模型。下面就三维矿床模型建模方 法分别进行简要的介绍。 4.1.1 线框模型 4.1.1 线框模型 矿体的地质形态复杂多变，很难用规则的几何体来描述。它需要一种灵活、 简便、 快速的方法来建立矿体的不规则几何模型。 目前， 比较知名的采矿 CAD 系 统均是采用表面模型来描述矿体的几何模型。 这种表面模型通常是由一系列的三 角面围成的表面。如 MICL 的 MICROMINE 的线框模型、Maptek 的 Vulcan 的模型 等均是表面模型。在不同的系统中表面模型的名称不同，但实质都是表面模型。 由于这种表面模型在未渲染前看似由线框构成， 因而在采矿 CAD 系统中多称为线 框模型。不过，这种表面模型可以进行体积估算、表面渲染、切制剖面、快速三 维显示等操作，比计算机图形学中的表面模型有所扩展。能满足矿山设计、生产 中地质制图的基本要求，也是建立矿体三维实体模型的基础。 线框模型的构建主要是采用了 TIN 技术（不规则三角网模型）中的 Voronoi 图与 Delaunay 三角形算法。TIN 是一种表示数字高程模型的方法，它既减少规 则网格方法带来的数据冗余，同时在计算效率（如坡度）方面又优于纯粹基于等 高线的方法。 TIN 模型根据区域有限个点集将区域划分为相连的三角面网络，区域中的任 意点落在三角面的顶点、边上或三角形内。如果点不在顶点上，该点的高程值通 常通过线性插值的方法得到（在边上用边的两个顶点的高程，在三角形内则用三 个顶点的高程） 。所以 TIN 是一个三维空间的分段线性模型，在整个区域内连续 但不可微。 TIN 的数据存储方式不仅要存储每个点的高程，还要存储其平面坐标、节点 连接的拓扑关系，三角形及邻接三角形等关系。TIN 模型在概念上类似于多边形 网络的矢量拓扑结构，只是 TIN 模型不需要定义“岛”和“洞”的拓扑关系。 有许多种表达 TIN 拓扑结构的存储方式，一个简单的记录方式是对于每一 个三角形、边和节点都对应一个记录，三角形的记录包括三个指向它三个边的记 录的指针；边的记录有四个指针字段，包括两个指向相邻三角形记录的指针和它 的两个顶点的记录的指针；也可以直接对每个三角形记录其顶点和相邻三角形， 见图 4-1 三角网存储方式。每个节点包括三个坐标值的字段，分别存储 X，X，Z 坐标。 这种拓扑网络结构的特点是对于给定一个三角形查询其三个顶点高程和相 邻三角形所用的时间是定长的，在沿直线计算地形剖面线时具有较高的效率。当 然可以在此结构的基础上增加其它变化，以提高某些特殊运算的效率，例如在顶 点的记录里增加指向其关联的边的指针。 图 4-1 三角网存储方式 图 4-1 三角网存储方式 不规则三角网数字高程由连续的三角面组成， 三角面的形状和大小取决于不 规则分布的测点，或节点的位置和密度。对于不规则分布的高程点，可以形式化 的描述为平面的一个无序的点集 P，点集中每个点 p，对应于它的高程值。将该 点集转成 TIN，然后采用 Delaunay 三角剖分方法。 4.1.1.1 Delaunay 三角网和 Voronoi 图 在介绍 Delaunay 三角网之前，首先需要知道 Voronoi 图，Voronoi 图又叫 泰森多边形或 Dirichlet 图， 它是由一组由连接两邻点直线的垂直平分线组成的 连续多边形组成。N 个在平面上有区别的点，按照最邻近原则划分平面；每个点 与它的最近邻区域相关联。Delaunay 三角形是由与相邻 Voronoi 多边形共享一 条边的相关点连接而成的三角形。Delaunay 三角形的外接圆圆心是与三角形相 关的Voronoi多边形的一个顶点。 Voronoi三角形是Delaunay图的偶图， 如图 4-2 所示。 对于给定的初始点集 P，有多种三角网剖分方式，其中 Delaunay 三角网具 有以下特征 （1）Delaunay 三角网是唯一的； （2）三角网的外边界构成了点集 P 的凸多边形“外壳”； （3）没有任何点在三角形的外接圆内部，反之，如果一个三角网满足此条 件，那么它就是 Delaunay 三角网。 图 4-2 实线为 Delaunay 三角形，虚线为 Voronoi 图 图 4-2 实线为 Delaunay 三角形，虚线为 Voronoi 图 （4）如果将三角网中的每个三角形的最小角进行升序排列，则 Delaunay 三 角网的排列得到的数值最大，从这个意义上讲，Delaunay 三角网是“最接近于 规则化的“的三角网。 Delaunay 三角形网的特征又可以表达为以下特性 （1） 在 Delaunay 三角形网中任一三角形的外接圆范围内不会有其它点存在 并与其通视，即空圆特性； （2）在构网时，总是选择最邻近的点形成三角形并且不与约束线段相交； （3）形成的三角形网总是具有最优的形状特征，任意两个相邻三角形形成 的凸四边形的对角线如果可以互换的话， 那么两个三角形 6 个内角中最小的角度 不会变大； （4）不论从区域何处开始构网，最终都将得到一致的结果，即构网具有唯 一性。 Delaunay 三角形产生的基本准则任何一个 Delaunay 三角形的外接圆的内 部不能包含其他任何点[Delaunay 1934]。Lawson[1972]提出了最大化最小角原 则，每两个相邻的三角形构成凸四边形的对角线，在相互交换后，六个内角的最 小 角 不 再 增 大 。 Lawson[1977] 提 出 了 一 个 局 部 优 化 过 程 （ LOP, local Optimization Procedure）方法。先求出包含新插入点 p 的外接圆的三角形，这 种三角形称为影响三角形（Influence Triangulation） 。删除影响三角形的公共 边（图 b 中粗线） ，将 p 与全部影响三角形的顶点连接，完成 p 点在原 Delaunay 三角形中的插入。 4.1.1.2 Delaunay 三角形网的通用算法－逐点插入算法 基于散点建立数字地面模型，常采用在 d 维的欧几里得空间 Ed 中构造 Delaunay 三角形网的通用算法逐点插入算法，具体算法过程如下 （1）遍历所有散点，求出点集的包容盒，得到作为点集凸壳的初始三角形 并放入三角形链表。 图 4-3 狄洛尼三角形中插入点 图 4-3 狄洛尼三角形中插入点 （2）将点集中的散点依次插入，在三角形链表中找出其外接圆包含插入点 的三角形（称为该点的影响三角形） ，删除影响三角形的公共边，将插入点同影 响三角形的全部顶点连接起来，从而完成一个点在 Delaunay 三角形链表中的插 入。 （3）根据优化准则对局部新形成的三角形进行优化（如互换对角线等） 。将 形成的三角形放入 Delaunay 三角形链表。 （4）循环执行上述第 2 步，直到所有散点插入完毕。 上述基于散点的构网算法理论严密、唯一性好，网格满足空圆特性，较为理 想。由其逐点插入的构网过程可知，在完成构网后，增加新点时，无需对所有的 点进行重新构网，只需对新点的影响三角形范围进行局部联网，且局部联网的方 法简单易行。同样，点的删除、移动也可快速动态地进行。但在实际应用当中， 这种构网算法不易引入地面的地性线和特征线，当点集较大时构网速度也较慢， 如果点集范围是非凸区域或者存在内环，则会产生非法三角形。 为了克服基于散点构网算法的上述缺点，特别是为了提高算法效率，可以对 网格中三角形的空圆特性稍加放松，亦即采用基于边的构网方法，其算法简述如 下 （1）根据已有的地性线和特征线，形成控制边链表。 （2）以控制边链表中一线段为基边，从点集中找出同该基边两端点距离和 最小的点，以该点为顶点，以该基边为边，向外扩展一个三角形（仅满足空椭圆 特性）并放入三角形链表。 （3）按照上述第 2 步，对控制边链表所有的线段进行循环，分别向外扩展。 （4）依次将新形成的三角形的边作为基边，形成新的控制边链表，按照上 述第 2 步，对控制边链表所有的线段进行循环，再次向外扩展，直到所有三角形 不能再向外扩展为止。 4.1.1.3 其它 Delaunay 优化算法 （1）Watson 的算法 Watson 的算法一开始就要形成一个包含所有给定点区域的超级三角形。起 初， 超级三角形是不完全的 如下图所示. 然后，通过算法继续递增地在已经 存在的三角化格网内插入新的点。 查找所有其外接圆包含新的点并且被删除掉给 出已知的插入多边形。这样就产生了新的三角网。这个过一直将持续到用完所有 的插入点以及所有拥有超级三角形顶点的三角形被删除掉。 a．初始的三角化 b．点的插入 a．初始的三角化 b．点的插入 c．元素的删除 d．最终的格网 图 4-4 Watson 的算法 c．元素的删除 d．最终的格网 图 4-4 Watson 的算法 这是 Delaunay 三角化最简单和最广泛应用的算法。 （2）Lawson 的算法或对角线交换算法 如果将一个点加入到一个已经存在的三角网中， 那么就形成了所有新的三角 形的外接圆。如果任何 相邻的区域位于任何三角形的外接圆内，那么用三角形 和它的相邻的区域便可形成一个四边形。 四边形的对角线被转换为一个新的三角 网。这个过程一直持续到不再有的错误三角形和不在需要变换。 （3）限制性的 Delaunay 三角化（Constrained Delaunay Triangulation） CDT 如果指定了平面的点集以及一系列不交错的边缘，那么 CDT 将产生如下的 网格 a、在三角网中包含了所有给定的边缘 b、尽可能使其与 Delaunay 三角网相近 这种算法已经在 N 个给定点的区域在 ONlogN 的时间复杂性下经过了测 试。如果指定点集和不交叉边缘，那么给定点集的 CDT 具有所有新的边缘所具有 的特性，即存在一个这样的圆 a、边缘的终点都落在圆内 b、如果有任一个节点在圆内，那么至少边缘的一个节点是不可见的，例如， 如果线是从点到边缘的两个节点来画的， 那么至少有一条直线会与网格的边缘相 交。 因此，如果没有指定边缘，那么 CDT 与 Delaunay 是一样的。CDT 是一种分 而取之的方法，因为给定的数据是依据任意的尺寸来筛选的，并且将区域划分成 若干带。这个过程将一直重复，直到得到整个区域的 CDT 为止。 （4）扫描线算法（Sweepline algorithm） 在 Voronoi 多边形中，扫描线被定义为一条移动的水平线与 Voronoi 边界 的交点。外接圆是通过找出连接节点的直线的三条平分线的交点来定位的。但是 这样得到 Voronoi 多边形不是太有力。一对一的几何变换在双曲线中绘制出一 条平分线或者一条垂直的半条。 任一两节点的平分线， 命名为 p 和 q, 描述如下 转换将 x 描述为这样，y 为初始值，点与节点的距离为 p。描述的 Voronoi 区域 Rp*和 Rq*是被画出的二分线限制的，为 Cpq,，它是一个双曲线或是一个直 的半条线。只有在扫描线到达所有的两个形成的节点，才能形生二分线。预期的 顶点是两个邻近的二分线的交点。 4.1.2 块段模型 4.1.2 块段模型 4.1.2.1 规则块段模型 规则块段模型包括常见的矩形块段组成的二维或三维模型， 在此主要指三维 块段模型。它是六十年代初为描述浸染状矿床而发展起来的一种建模方法，多用 于大型露天矿床， 这种方法的产生是为了适应当时采用地质条件下方法进行储量 计算工作。其实质是采用一系列的相同尺寸的立方块去表达矿体。 这种技术有几个明显的优点1）在一个资源有限的计算环境中，编制程序 时可以采用隐含的定位技术，即采用一系列的三维数组来存放各个块段的信息， 如品味，岩性等，数组的下标分别对应于块段的行、列、层号，从而可以节省 存储空间与计算时间；2）由于采用相同尺寸的块段，在进行地质统计学品位估 值时，克立格方程组的系数矩阵[K]只需计算一次，可以大大减少计算量，同时 也方便进行储量的统计和计算。 但是这种技术使用起来很不灵活，为了解决运算效率问题，通常要旋转模型 的坐标系，使之与矿体的走向一致。在遇到弧形或多个矿体的时候，有时可能不 止构筑一个模型，需要不同的旋转坐标，并相应的带来切断面的数量问题，它们 加大了处理时的开销。另一方面，这种方法不能精确描述矿体的边界，用一个 2.5m矩形单元模拟一断面为25m10m梯形采场， 可能发生大致50％的边界误差， 即采场中矿石可能比构模过程中得到的要多或少 10％。为了比较精确的描述矿 体边界， 必须相对减小块段尺寸， 但是， 很小块段的克立格估算值时过分平滑的， 品位估值可能不准确，且大大增加了计算工作量与存储空间。 矿床建模包括两方面的内容矿体边界的确定、边界内品位信息的估算。针 对规则块段模型在矿床建模问题上存在的矛盾，结合计算机软硬件水平，提出一 种改革方案，即变块模型。 4.1.2.2 变块模型 变块模型又叫做不规则块段模型。这种建模方法与规则块段模型一样，也是 在整个建模范围内用立方块来模拟矿体，不同的是，这些立方体尺寸不在是固定 的，当对最初的由尺寸相同的方块组成的模型施加了地质解释之后，就沿着边界 或地质界面建立“子块段”，以保证更加逼近详细的地质要求。 在七十年代，“SEAMSYS”控制数据软件包和“OBMRTZ”咨询系统就采用了 可变尺寸块段模型。八十年代末，MICROMINE 系统采用更高级的技术，在矿体的 任何地方可随意变换块段大小，而且也不在存储无用的信息。这样既达到精确描 述矿体边界或地质构造的目的，又减少了存储空间。虽然在这个系统中，管理模 型的程序变得比较复杂，但那是程序开发人员的事，对应用单位来说是大为有利 的。 事实上， 采用变块模型虽然能够比较精确的描述矿体边界， 减少了存储空间。 但是它并未从根本上改变边界确定与品位估值的矛盾。 原因还在于小块段品位克 立格估值问题上。 因此，尽管至此为止的建模方法不仅能够估计品位，而且可以处理详细的地 质结构和模拟表明，但是这并不十分尽人意。建立矿床模型的目的不仅是为了进 行储量计算，更重要的是要在其上进行采矿设计、编制计划等，这在一定程度上 似乎更关心矿体的几何形态，这就需要线框模型的构建。 4.1.3 线框块段嵌套模型 4.1.3 线框块段嵌套模型 线框块段嵌套模型实现了“块段模块”与“线框模型”相结合的混合构模 方法。即在处理矿岩质量信息时采用“块段模型”，在计算矿岩体积、重量、模 拟开采、输出工程图时采用“线框模型”，将“块段模型”的优化计算优势和 “线框模型”可以精确处理矿体边界的特点完美地结合在一起， 解决了目前广泛 采用的用单一“块段模块”或“线框模型”方法建模无法避免的缺点。 混合构模 技术在国际上是最先进的，目前仅美国、英国、澳大利亚等少数国家实现。 4.2 地质变量估值方法 4.2 地质变量估值方法 为了利用有限的地质及生产勘探样品数据， 确定矿床内部相关变量 （如品位、 岩体质量指标、岩性）的分布，必须选用合理的地质变量估值方法。 4.2.1 区域化变量 4.2.1 区域化变量 区域化变量概念的提出是随着统计学理论和方法在地学领域和矿床储量计 算中的应用和发展而提出的，并由此产生了地质统计学这一数学地质分支。区域 化变量一方面具有一般随机变量（如设备生产能力）取值的随机性，另一方面， 这种变量在空间上具有一定的相关性，即一定范围内的变量之间互相存在着影 响。在地学领域，许多变量具有这种特征，如矿床的品位、岩体质量指标和参数 等等。 目前对于岩性以外的地质变量 （区域化变量） 的估值方法主要有趋势面法、 距离幂次反比法、样条函数法、地质统计学方法。本节主要对趋势面法和距离幂 次反比法进行简单的介绍， 由于地质统计学方法是目前国外储量计算与管理的主 要方法，国内自 1977 年引入该理论后，逐渐获得国家储委的认可，成为品位估 计及储量估算的主要方法，将单列一节详细介绍。 4.2.1.1 趋势面法 当某变量体现出在空间上连续变化的趋势时， 可以用一个平滑的数学平面加 以描述，即用已知采样点数据拟合出一个平滑的数学平面方程，再根据该方程计 算无测量值处的数据。这种只根据采样点的属性数据与地理坐标的关系，进行多 元回归分析得到平滑数学平面方程的方法，称为趋势面法，又称全局多项式内插 法Global PolynomialInterpolation，GPI。 趋势面法用数学的方法把观测值划分为 2 部分趋势部分和偏差部分。趋势 部分反映区域性的总的变化，受大范围的系统因素控制；偏差反映局部范围的变 化特点，受局部因素和随机因素控制。设以zyxg.,表示所研究的变量在空间上 的分布特征，其中zyx,,为空间上的坐标。根据前述原理，趋势分析的一个任 务是将 g观测值变换分解为 2 个部分 ezyxTzyxg,,,, （4-1） 式中zyxT,,称为趋势值；e 为剩余值。 在研究区域内的样品点，对于已知的样品数据 i g 、 i x 、 i y 、 i z 通常可用多项 式回归作为回归方程求得它们的趋势值和剩余值， 因为任何函数在一个适当的范 围内总可以用多项式来逼近， 而且通过调整多项式的值可使所求的回归方程适合 问题的需要。理论上来说，多项式次数愈高，则趋势值愈接近观测值，在具体工 作中必须根据实际需要恰当地选择趋势面的次数，一般来说很少超过 5 或 6 次。 4.2.1.2 距离幂次反比法 距离幂次反比法（Inverse Distance Weighted, IDW）是最常用的空间内插 方法之一，是一种与空间距离有关的插值方法，在计算插值点取值时按距离越近 权重值越大的原则， 用若干临近点的线性加权来拟合估计点的值。 具体表达如下。 ∑ ∑ n i p i n i i p i d g d g 1 1 1 1 （4-2） 式中，g是估计值， i g是第nii,, 2 , 1L个样本， i d是距离，p是距离的幂， 它的大小显著影响着估值的结果。距离幂次反比法的幂次不同，则有不同的适用 范围和估值效果，当幂次等于 2 时，该方法就称“距离平方反比法” 。 4.2.2 地质统计学方法 4.2.2 地质统计学方法 地质统计学是根据英文单词 Geostatistics 字面意思翻译过来的， 它是关于 取自地球的大量时间的收集、 分析、 解释和表达的一个数学分支。 按其基本原理， 可定义其为地质统计学是以区域化变量理论为基础，以变异函数为主要工具， 研究在空间分布上既有随机性、又有结构性的自然现象的科学。 早在本世纪初，传统的统计学方法就被用于进行数据分析。在地质矿产方面 最初也是利用传统的统计学作为分析数据的工具，但到了四十年代后期，H.S. 西舍尔（Sichel）判明南非的金矿品位呈对数分布，于是奠定了地质统计学的开 端。 五十年代初，南非工程师 D.G.克立格（Krige）根据多年对南非金铀砾岩型 金矿储量计算的经验，提出一个论点“可以预计，一个矿山总体中的金品位的 相对变化要大于该矿山某一部分中的金品位的相对变化。”换句话说，以较近距 离采集到的样品很可能比以较远距离采集到的样品具有更近似的品位。 这一论点 是描述在多维空间内定义的数值特征的空间统计学据以建立的基础。 到了六十年代， 开始认识到需要把样品值之间的相似性作为样品间距离的函 数来加以模拟，并且得到了半变异函数。随后，法国学者 G.Matheron 将 D.G.krige 等人的成果理论化、系统化，提出“区域化变量”的概念，并于 1962 年发表了应用地质统计学论 ，该著作的出版标志着地质统计学作为一门新兴 的边缘学科而诞生，地质统计学开始进入学术界。在法国的枫丹白露成立了地质 统计学中心（Entre De Geostatistiues） ，吸收培养了一大批学员，不仅为地质 统计学的理论研究而且为地质统计学的传播起到了巨大的作用。 七十年代，地质统计学的研究突飞猛进。这期间，无论是论文发表的篇数、 理论突破的频度，还是在世界各地表现出来的对地质统计学的极大关心程度，都 反映了地质统计学的发展达到了前所未有的繁荣景象。至八十年代，地质统计学 的发展出现了一个全面慢化的过程，这主要是地质统计学的理论已逐步趋于成 熟。 目前，地质统计学主要分为两大学派以法国的 G.Matheron 为首的“枫丹 白露学派”，主要研究领域为参数地质统计学；以美国的 A.G.Gournel 为代表的 “斯坦福地质统计学派”，主要研究领域为非参数地质统计学。前者形成了简单 地质统计学（Simple Kriging） 、普通地质统计学（Ordinary Kriging）等线性 地质统计学理论及泛克立格（Universal Kriging） 、析取克立格（Disjunctive Krigig）等非线性地质统计学理论。他们都要求样品品位的分布参数满足若干条 件；后者则提出了指示克立格法概率克立格法（Probabilitric Kriging）以及 指示条件模拟（Indicator Condition Simulation）等一套理论及方法。 鉴于地质统计学非凡的估值效果，西方国家矿产法规定，矿山申请贷款时， 必须向银行当局提高地质统计学所计算出的矿储量。 地质统计学正式引入我国始于 1977 年。当时，由随美中贸易全国委员会矿 业代表团来华访问的美国福录尔采矿公司的 H.M.Parker 将地质统计学的基本概 念系统的介绍给我国的有关数学地质专家，立即引起冶金部的极大关注。1980 年，中国金属学会冶金地质学会在桂林召开的数学地质遥感地质学术会议上，成 立了“地质统计学协作小组”。该小组担负着冶金系统地质统计学工作的协调、 理论方法的研究、程序的编写以及方法的推广应用和普及工作。为了更有效的推 广应用地质统计学，候景儒、黄竞先等翻译出版了矿产地质统计学 ，同时， 还编著了地质统计学及其在矿产储量计算中的应用 。由于地质统计学工作者 的不懈努力，地质统计学在我国很快的发展起来，1989 年 11 月 1 日，在江西九 江召开的第一届全国地质统计学学术讨论会展示出地质统计学在我国的蓬勃发 展。相关期刊杂志上，发表的有关地质统计学理论研究及地质统计学由于的论文 逐年增加，地质统计学的发展在其领域中取得了显著的成绩。 地质统计学方法，通常也称为克里格方法，这种方法主要具有 2 个方面的意 义定量地评价区域化变量的空间结构性、对未取样位置处的区域化变量进行估 值。 4.2.2.1 变异（或协方差）函数计算方法 4.2.2.1 变异（或协方差）函数计算方法 区域化变量的空间结构性主要通过试验变异函数的计算和理论变异函数的 拟合而获得。常用的试验变异函数计算方法有 一、传统变异函数计算方法 2 1 1 2 N h ii i hZxZxh N h γ −⎡⎤ ⎣⎦∑ （4-3） 二、传统交叉变异函数计算方法 1 1 2 N h ZYiiii i hZxZxhY xY xh N h γ −−⎡⎤⎡⎤ ⎣⎦⎣⎦∑ （4-4） Z x和 Y x为同一区域中两个相关的区域化变量。在协同克立格中，当 Z x和 Y x是相同的区域化变量时，上式变成传统变异函数。 三、非遍历协方差函数 1 1 N h iixy i C hx ym m N h − ∑ （4-5） 11 11 , N hN h xiyi ii mx my N hN h ∑∑ （4-6） i x和 i y为不同变量时，上式变为协同协方差函数。 四、非遍历的相关函数 xy C h e h σ σ − （4-7） x σ和 y σ 分别为点对头和尾的标准，当x和y为不同变量时，即为协同非遍 历协方差函数。 五、一般相关变异函数 2 2 GR y x h h m m γ γ ⎛⎞ ⎜⎟ ⎝⎠ （4-8） 符号与上式相同。 六、两两相关变异函数 2 1 2 2 N h ii RR i ii xy h N hxy γ ⎡⎤− ⎢⎥ ⎣⎦ ∑ （4-9） 经验表明，当数据分布正向偏移时， （五）与（六）能够缓冲数据的稀疏性 和分离性，他们有时比其它方法更能揭示空间结构和各向异性。 七、对数变异函数 2 1 1 lnln 2 N h Lii i hZxZxh N h γ −⎡⎤ ⎣⎦∑ （4-10） 当区域化变量服从对数正态分布或变量值较小时，可以使用这种方法。 八、罗德函数（Rodogram） 1 1 2 N h Rii i hxy N h γ − ∑ （4-11） 九、麦当娜函数（Maogram） 1 1 2 N h Mii i hxy N h γ − ∑ （4-12） （八）和（九）一般用在处理大尺度结构，而不用在有块金效应的方差模型 中。 十、指数变异函数 1 , 0 xZ I x Z xZ ≤⎧ ⎨ ⎩ （4-13） 其中，Z为门坎值（cutoff） ，此为指示克立格的基础。 大多数试验变异函数具有如图 4-5 所示特征。 图 4-5 试验变异函数及其拟合曲线与参数 图 4-5 试验变异函数及其拟合曲线与参数 试验变异函数计算完毕后，需要对其理论变异函数参数进行拟合，理论变异 函数模型有球状模型、指数模型、高斯模型、线性模型等，如图 4-6。其中以球 状模型为多。 图 4-6 理论变异函数模型（依次为球状、指数、高斯、线性） 图 4-6 理论变异函数模型（依次为球状、指数、高斯、线性） 球状模型的数学公式为 ⎪ ⎩ ⎪ ⎨ ⎧ ≥ 0.8，矿体具有随机性 变化。 4.2.2.2 克立格估值方法 4.2.2.2 克立格估值方法 一、简单克立格法（sk） 该方法要求区域变量 Z 满足二阶平稳假设，即 A． []常数） mxZE （4-16） B． {} 2 ,Cov Z xZ xhE Z xZ xhm⋅−⎡⎤ ⎣⎦ （4-17） 估计值 1 n ii i ZmZmλ ∗ − ∑ （4-18） SK 方程组 1 ,,,1,2,..., n iijj i Cov Z ZCov ZZjnλ ∑ （4-19） SK 估计方差 2 111 ,2,,, nnn skijiiijiji iij Cov Z ZCov Z ZCov Z ZCov Z Zδλλλ ⎡⎤ −− ⎣⎦ ∑∑ ∑ （4-20） 二、普通克立格法（OK） 实际应用中，二阶平稳假设很难满足（有时协方差函数不存在） ，放宽要求 在本征假设下 A． 0,,E Z xZ xhxh−∀ ∀⎡⎤ ⎣⎦ （4-21） B．增量 Z xZ xh−⎡⎤ ⎣⎦的方差函数存在且平稳（不依赖于 x,即 2 ,,Var Z xZ xhE Z xZ xhxh−−∀ ∀⎡⎤⎡⎤ ⎣⎦⎣⎦ （4-22） 估计值 1 , n ii i ZZλ ∗ ∑n为信息点数目； OK 方程组 j 1 1 ,, 1,2,..., 1 n jiji n i i Z ZZ Z in λ γμγ λ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∑ ∑ （4-23） 克立格方差 2 1 ,, n ii i Z ZZ Zδλ γγμ − ∑ ， μ为拉格朗日常数； （4-24） 简单克立格法不用计算变异函数，直接用协方差函数求解 SK 方程，而普通 克立格必须先算出变异函数，然后用公式 0Cov hCovhγ− 将其转变为 协方差函数求 OK 方程。另外，普通克立格发有一个无偏约束条件 1 1 n i i λ ∑ ，对数 据分布要求较松，应用范围更广。 三、泛克立格法（UK） 在生产实践中，不少区域化变量存在漂移，连本征假设也满足不了，只能继 续放宽条件 A．假设 Z x有非平稳的一阶矩和二阶矩 2 2, E Z xZ ym xm y DZ xZ yx yγ ⎧−−⎡⎤ ⎪⎣⎦ ⎨ −⎡⎤ ⎪ ⎣⎦⎩ （4-25） B．假设 Z x可分解为漂移 m x与剩余 R x之和，即 Z xm xR x （4-26） C．假设漂移可以表示为1k 个的单项式函数 0,1,2,..., i fxik 的线性组 合 0 k ii i m xa fx ∑ （4-27） D．假设在一定小范围内，平稳实验变异函数可以代替非平稳的实验变异函 数，并在变程内接受 ,, ab Cov ZZCov Z Zabγ−− （4-28） 估计值 1 n ii i ZZλ ∗ ∑ （4-29） UK 方程 10 1 ,,1,2,..., 0,1,..., j j nk l iijivj ii n ll jvv j C v vbC v vjn bblk λμ λ ⎧ − ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∑∑ ∑ （4-30） 其中 11 , j j ll vlvl vv j bfx dx bfx dx vv ∫∫ （4-31） 估计方差 2 10 ,, nk l iilv il C v vC v vbδλμ − ∑∑ （4-32） 泛克立格法不但考虑了数据的空间自相关性，而且考虑了空间趋势性。以上 三种克立格法均为线性地质统计学。 四、析取克立格（IK） 未知的估计值是单个信息可预测函数的和 1 n DKii i ZfZ ∗ ∑ （4-33） 同样要求在无偏最优条件下确定估计值。当点对 ba ZZ ,的二元分布形式已 知时，估计方差才能被计算。Matheron 假设区域化变量服从标准高斯一元和二 元分布，由此得变异函数模型（Hermitian 模型） 。 上式 ii fZ 可用 Hermitian 展开表达 , 0 iii kki k fZZλ η ∞ ∑ （4-34） 其中， ki Zη 为 Hermitian 多项式， , i k λ 是待定权。 泛克立格要求线性非平稳假设，而析取克立格假设区域化变量服从高斯分 布，并要求对分布进行截尾处理。当 iiii ZZfλ时，DK 变 UK，也就是说，泛克 立格可看做是析取克立格的子集。 五、指示克立格（IK） 设区域为 D，针对每个样品点x，在xD∈上定义一个 Z x的阶梯函数 1 , 0 I x Z ⎧ ⎨ ⎩ ZxZ ZxZ ≤ （4-35） Z为门坎值（cutoff） 。任意子区域 A AD∈ 中低于门坎值Z的样品所占区 域A的比例表示为 1 ,, A A ZI x Z dx A φ ∫ （4-36） ,A Zφ称为废石函数，而其补集称为矿石回收率。在边界品位I的条件下， 随机函数 ,I x Z 服从二项分布，其期望值是 {},PrE I x Zob Z xZ≤⎡⎤ ⎣⎦ （4-37） 变异函数为 21 ,,,0,, 2 III h zE I xh ZI x ZCZCh Zγ−−⎡⎤ ⎣⎦ （4-38） ,A Zφ 的估计值表示为 1 ,, N a a A ZI a Zφλ ∗ ∑ （4-39） a代表信息点，IK 方程组和 OK 方程组相同，实际计算时门坎值可选一系列 值，变异函数用中间门坎值计算即可。 指示克立格方程组 1 1 ,;; ; 1,2,..., 1 n iIiji j n i i x xzx A z in λγμγ λ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∑ ∑ （4-40） 指示克立格方差 2 1 , ;, ; n IIKIii i x A zA A zσλ γγμ − ∑ （4-41） 矿产储量计算中，由于特异值及不同矿化母体的存在会影响克立格估值效 果，而指示克立格可以很好的解决上述问题。在非人为地除去以上因素的基础上 用指示克立格法进行有效估值。另外，指示克立格可以对低于或高于某个门坎值 的块段所占的比例进行估计。 六、概率克立格法（PK） 下述公式中，符号与指示克立格法相同。 11 ,, NN aaa aa A ZI a ZZφλη ∗ ∑∑ （4-42） 将信息点 a Z变成有序数据，是原始数据从小到大排序，若 a Z的秩为 1kkn≤≤，则将 a Z转化成k n，相应数据变为 1 2 ,,...,, a n k n n