地质统计学讲义.doc

地 质 统 计 学 讲 义 地 质 统 计 学 讲 义 张树泉 教授 北 京 科 技 大 学 第1章 地质统计学的发展历史和现状 1.1地质统计学的发展历史 地质统计学是根据英文单词Geostatistics的字面意思翻译过来的，从词源学上讲，按照韦氏（N .Webster）大词典对于“geo”（地球、土地）和“Statistics”（统计学）两词的释义，地质统计学（Geostatistics）的定义便是“关于取自地球的大量数据的收集、分析、解释和表达的一个数学分支”。就矿山地质统计学的内容范围来说，这一定义是十分恰当的。地质统计学包含经典统计学与空间统计学，其重点是地球状况，也就是说着重于地质特征的分析。按其基本原理可定义为地质统计学是以区域化变量理论为基础，以变异函数为主要工具，研究那些在空间分布上既有随机性，又有结构性的自然现象的科学。 早在上世纪10年代里，传统的统计学方法就已用于分析地质数据。在地质矿产方面最初也是利用传统的统计学作为分析数据的工具，直到上世纪40年代后期，当南非统计学家H.S西奇尔（Sichel）判明南非各金矿的样品品位呈对数正态分布以后，才真正确立了地质统计学的开端。 1951年，南非的矿山工程D.G.克立格（Daniel Krige）在H.S西奇尔研究的基础上提出一个论点“可以预计，一个矿山总体中的金品位的相对变化要大于该矿山某一部分中的金品位的相对变化”。换句话说，以较近距离采集的样品很可能比以较远距离采集的样品具有更近似的品位。这一论点是描述在多维空间内定义的数值特征的空间统计学据以建立的基础。 到上世纪60年代，才认识到需要把样品值之间的相似性作为样品间距离的函数来加以模拟，并且得出了半变异函数。法国概率统计学家马特隆（Matheron）创立了一个理论框架，为克立格作出的经验论点提供了精确而简明的数学阐释。马特隆创造了一个新名词“克立格法”（Kriging），藉以表彰克立格在矿床的地质统计学评价工作中所起到的先驱作用。即1962年，马特隆在克立格和西奇尔研究的基础上，将他们的成果理论化、系统化，并首先提出了区域化变量（Regionalized variable）的概念，为了更好地研究具有随机性及结构性的自然现象，提出了地质统计学（Geostatistics）一词，发表了应用地质统计学，该著作的出版标志着地质统计学作为一门新兴边缘学科而诞生。地质统计学开始进入了学术界。在法国枫丹白露成立了地质统计学中心（Centre de Geostatistiques），培养了一大批学员，不仅为地质统计学的研究而且为它的传播起到了巨大的作用。 70年代，地质统计学的发展突飞猛进。在此期间，从理论突破的频度、论文发表的篇数、以及世界各地对地质统计学所表现的极大关心程度，都说明地质统计学达到了前所未有的发展阶段。80年代，地质统计学的发展出现了一个全面慢化的过程，是因为它的理论已趋于成熟，90年代初，地质统计学已经被广泛地承认是矿床评价的必要部分，在我国已经认可用地质统计学对矿床进行评价的地质报告。目前条件模拟技术广泛应用于石油、采矿、水文、和环境保护等领域中。 1.2 地质统计学的发展现状 1.2.1 地质统计学在国外的发展状况 地质统计学经过40多年的发展,已成为能表征和估计各种自然资源的工程学科,并由法国、南非及一些法语国家推广到几乎全世界。作为一门年轻的边缘学科仍正处于蓬勃发展的阶段。地质统计学研究的最新进展如下 （1）世界上，地质统计学研究可分为两派以法国G.Matheron为代表的“枫丹白露地质统计学派”，主要研究参数地质统计学；和以美国A.G.Journel为代表的“斯坦福地质统计学派”，主要研究非参数地质统计学。前者继续开展以正态假设为基础的析取克立格法及条件模拟的研究，同时把主成份分析和协同克立格法结合起来，提出多元地质统计学的基本思想，形成了简单克立格（Simple Kriging）、普通克立格（Ordinary Kriging）、泛克立格（Universal Kriging）以及析取克立格（Disjunctive Kriging）等一套理论和方法。上述方法的所有计算均依赖于实际样品数据，并求得区域化变量理论模型的若干参数，故称为“参数地质统计学”（Parametric Geostatistics）；后者发展了无须对数据分布做任何假设的指示克立格（Indicator Kriging）、概率克立格（Probabilitric Kriging）以及指示条件模拟（Indicator Condition Simulation）等一套理论和方法，同时考虑如何使用软数据的问题。两个学派都在研究地质统计学中的稳健性问题。 （2）不同的数学方法不断地引进克立格法，多学科渗透形成各种新的克立格法。国外已将地质统计学与贝叶斯统计（Bayes Statistics）数学形态学、模糊数学、自回归、时间序列分析、分数维等相结合研究出许多新的方法[7，8]。例如对于含有特异值（如特高品位），接近高斯分布的数据，将稳健统计学的思想应用于求变异函数，提出了稳健克立格方法；应用条件概率估计空间分布用以解决可回采储量的估计和预报区间的确定问题；将多元区域化变量引入克立格法，利用两个以上具有相关性的变量协同对某一变量进行估值，产生了协同克立格法;将多元区域化变量引入指示克立格法便产生协同指示克立格；将多元统计思想引进克立格产生因子克立格等。在20世纪80年代末至90年代初，地质统计学工作者将分形理论与地质统计学相结合，形成了分形地质统计学 。 （3）条件模拟可保持变量的空间自相关性不变，使观测点上的模拟值等于实测值, 域化变量，而且有可能实现三维空间模拟。目前条件模拟的应用已不局限于评价可采矿石与局部评价问题，对工程位置部署提供指导将是其应用的一个方向，从对矿体的条件模拟向模拟采矿过程，到模拟矿山开发后的各个阶段中可采储量的特征变化，用于评价可采矿石量与品位的局部变化。条件模拟还可用于矿山和选矿厂的管理，并在煤田、石油地质勘探与开发领域也得到应用。主要条件模拟方法有G.马特隆的“转向带法”和A.G儒尔耐耳的“非高斯分布快速模拟法”等。 （4）地质统计学不断应用于新的领域，如矿山地质与矿产经济、石油地质与煤田地质勘探与开发、水文工程地质、环境科学的研究、海洋渔业资源动态管理、海洋地球物理、农林科学、图像处理、制图技术、公共卫生、财政分析以及材料科学的研究等。Oliver和Badr（1994）利用地质统计学技术对三个地区（Hereford,Buxton,Nottingham）的氡（222Rn）进行空间结构及其与各种癌病的相关性分析,结果表明，氡含量高的地区癌病的发生率也高，同时也查明氡含量主要与地质作用有关[27]。 （5）研制出一批高水平的地质统计学方法计算程序软件。在地质统计学的理论及方法基础上开发了许多成熟的应用软件。如美国开发的矿床建模软件包（Deposit Modeling System），功能上可覆盖矿山地质设计的全过程；而MICL（英国矿业计算机有限公司）开发的DATMINE软件包，则集地、测、采于一体；法国巴黎高等矿院地质统计学研究中心研制出两种大型软件系统ISATIS系统及HERESIM系统；澳大利亚的MICROMINE软件,SURPAC软件，加拿大的GEOSTAT软件，CAMET软件和GLS软件系统等。 1.2.2 地质统计学在国内的发展状况 地质统计学是在1977年由美国福禄尔采矿金属有限公司（Flour Mining Meta Incorporation）H.M.Parker博士随美中贸易全国委员会矿业代表团来华访问，传入我国，继而得到进一步的发展。 在美国学者H.M.Parker博士将地质统计学的基本概念和内容系统介绍给我国的数学地质及勘探、矿山设计人员后，我国的有关的学术专业团体的学术活动开展得非常活跃，1980年4月，中国金属学会冶金地质学会在广西桂林召开的第一届遥感地质数学地质学术会议上，正式成立了冶金地质系统的“地质统计学协作小组”。随后几年的各种地质矿业学会上，地质统计学的论文不断增多，与此同时，地质统计学的普及工作相继展开，地矿、冶金、石油、核工业和煤炭等行业，为普及这门学科，举办了各种学术讲座。我国先后派出了许多专家学者到国外学习深造地质统计学，原地矿部、冶金部和有色总公司等先后邀请了国外著名的地质统计学家A.G.儒尔耐耳、D.G克立格、J.M伦杜等来等来华讲学进行学术交流。与此相随，国内出版了关于地质统计学的重要书籍和论文地质统计学，地质统计学及其在矿产储量计算中的应用，矿业地质统计学，线形地质统计学，数学地质的方法与应用，地质统计学及其在储量计算中的应用等。 1989年11月召开的全国第一届地质统计学学术讨论会，地质统计学在我国的发展进入了一个新的阶段，理论研究更加深入，涉及的方法原理更加广泛，整体理论水平与国际水平接近。不仅如此，而且为适应生产的需要，在应用方面有了实质性进展，用地质统计学方法提交的地质勘探成果来开发矿山，还有的开发研制出适于国内生产需要的软件系统，如北京科技大学的克立格法程序系统（3DOK，3DCOK，2DIK）；德兴铜矿应用地质统计学软件，用于指导矿山设计、生产和规划；西安石油学院的克立格绘图系统KMS；原地矿部信息院研制的KPX2.0固体矿产勘查评价自动化系统；武警黄金指挥部1993年用地质统计学方法计算了陕西洛南驾鹿金矿区8801号矿体储量[20]，并提交地质勘探报告；有色总公司北京矿产地质研究所应用地质统计学方法计算了煎茶岭金矿、云南曼家寨锡、锌矿床、查干布拉根银铅锌矿等多个矿床的储量；北京有色设计院用地质统计学方法计算河南金堆城铜矿等矿床储量等等。 第2章区域化变量理论 所谓区域化变量是指以空间点X的三个直角坐标（Xu,Xv,Xw）为自变量的随机场Z（Xu,Xv.Xw）Z（X）。当对它进行了一次观测后，就得到了它的一个现实Z（X），它是一个普通的三元实值函数或空间点函数。区域化变量的两重性表现在观测前把它看成是随机场（依赖于坐标（Xu,Xv,Xw）），观测后把它看成一个空间点函数（即在具体的坐标上有一个具体值）。 G.马特隆定义区域化变量是一种在空间上具有数值的实函数，它在空间的每一个点取一个确定的数值，即当由一个点移到下一个点时，函数值是变化的。 从地质及矿业角度来看，区域化变量具有如下性质 （1）空间局限性即它被限制在一个特定的空间（如一个矿体内）；该空间称为区域化的几何域；区域化变量是按几何支撑定义的。 （2）连续性不同的区域化变量具有不同的连续性，这种连续性是通过相邻样品之间的变异函数来描述的。 （3）异向性当区域化变量在各个方向上具有相同的性质时称各向同性，否则称各向异性。 （4）相关性一定范围内、一定程度上的空间相关性，当超出这一范围后相关性减弱以至消失。 （5）对于任一区域化变量而言，特殊的变异性可以叠加在一般规律之上。 2.1变异函数与协方差 变异函数 研究区域化变量使用变异函数，在一维条件下变异函数定义如下当空间点X在一维X轴上变化时，把区域化变量在x与xh 处的值与的差的方差之半定义为区域化变量在x轴方向上的变异函数，并记为 （21） 在二阶平稳假设下有 h 于是式 （21） 改写成下式 （22） 从式（22）知变异函数依赖于两个自变量x和h，当变异函数与位置x无关，而只依赖于分隔两个样品点之间的距离h时，则就可改写为 （23） 应当说明的是有时把定义为变异函数，则就是半变异函数了，而把直接定义为变异函数时，决不会影响它的性质。 协方差 如前所述，当随机函数中只有一个自变量x时称为随机过程，而随机过程在时刻及处的两个随机变量及的二阶中心混合矩定义为随机过程的协方差函数 （24） 当随机函数依赖于多个自变量时，称为随机场。而随机场在空间点 x和xh 处的两个随机变量和的二阶中心混合矩定义为随机场的自协方差函数 25） 协方差函数一般依赖于空间点x和向量h。当式（25）中h0时，则协方差函数变为，即等于先验方差函数，当其不依赖于x时，简称方差，从而有 （26） 2.2平稳假设及内蕴假设 在地质统计学研究中是用变异函数表示矿化范围内区域化变量的空间结构性的，要用式（22）计算变异函数时，必须要有，这一对区域化变量的若干实现，而在实际工作中（尤其是地质、采矿工作中）只有一对这样的现实，即在 x，xh点只能测得一对数据（因为不可能恰在同一样点上取得第二个样品），也就是说，区域化变量的取值是唯一的，不能重复的。为了克服这个困难，提出了如下的平稳假设及内蕴假设。 平稳假设（stationary assumption） 设一随机函数Z，其空间分布律不因平移而改变，即若对任一向量h，关系式成立时，则该随机函数Z为平稳性随机函数，确切地说，无论位移向量h多大，两个k维向量的随机变量{ ，，, }和{ ，，, }有相同的分布律。通俗地说，在一个均匀的矿化带内，与之间的相关性不依赖于它们在矿化带内的特定位置。这种平稳假设至少要求的各阶矩均存在，且平稳，而在实际工作中却很难满足。在线性地质统计学研究中，我们只需假设其1、2阶矩存在且平稳就够了，因而提出二阶平稳或弱平稳假设。 当区域化变量满足下列两个条件时，称该区域化变量满足二阶平稳 （1） 在整个研究区内，区域化变量的期望存在且等于常数 常数 x （2） 在整个研究区内，区域化变量的空间协方差函数存在且平稳 x,h 当h0时，上式变成 x 即它有有限先验方差。 上述各式中及表示协方差，表示方差。 协方差平稳意味着方差及变异函数平稳，从而有关系式 （27） 内蕴假设 在实际工作中，有时协方差函数不存在，因而没有有限先验方差，即不能满足上述的二阶平稳假设，例如一些自然现象和随机函数，它们具有无限离散性，即无协方差及先验方差，但却有变异函数，这时，我们可以放宽条件，如只考虑品位的增量而不考虑品位本身，这就是内蕴假设的基本思想，当区域化变量的增量满足下列两个条件时，称该区域化变量满足内蕴假设 （1） 在整个研究区域内，随机函数的增量的数学期望为0 x, h （2） 对于所有矢量的增量的方差函数存在且平稳，即 x, h 即要求的变异函数存在且平稳。 内蕴假设可以理解为随机函数的增量只依赖于分隔它们的向量h（模和方向）而不依赖于具体位置x，这样，被向量h分隔的每一对数据 可以看成是一对随机变量的一个不同现实，而变异函数的估计量是 （28） 式中，是被向量h相分隔的试验数据对的数目。 如果随机函数只在有限大小的邻域（例如以a为半径的范围）内是平稳的（或内蕴的），则称该随机函数服从准平稳（或准内蕴）假设，准平稳或准内蕴假设是一种折衷方案，它既考虑到某现象相似性的尺度（scale），也顾及到有效数据的多少。实际工作中，可以通过缩小准平稳带的范围b而得到平稳性，而结构函数（协方差或变异函数）只能用于一个限定的距离|h| ≤b，例如界限b为估计邻域的直径，也可以是一个均匀带的范围，当|h|b时，区域化变量和就不能认为同属于一个均匀带，这时，结构函数或只是局部平稳的，所以，我们把只限于|h|≤b范围内的二阶平稳称为准平稳，把只限于|h| ≤b范围内的内蕴称为准内蕴。上述的这种概念在文本后边将要讨论的克立格估计中十分重要，因此我们可以用这种想法确定适当大小的移动邻域，在该邻域内，随机函数的数学期望和协方差（或变异函数）是平稳的，而且在该邻域内的有效数据足以进行统计推断。显然，平稳假设或内蕴假设可以理解为一种相对的概念。 2.3估计方差 任一估计方法，由于估计时所用样品与被估块段的大小并非严格相等，从而使被估块段的实际值与估计值不同，即产生估计误差，一个储量计算方法的可靠程度就是根据该方法所包含的误差大小来衡量的。最好的估计方法是误差最小的方法。 设有一矿床被分成大小相等的以（）为中心的N个块段，令每一块段的实际品位为（），而用某种方法估计出的估计品位为（），这时就有估计误差 即 （） 可以证明，若 是二阶平稳的话，也是二阶平稳的，因而 （常数） 而且有有限方差，且平稳 当然，我们总是希望估计的平均值与实际值的平均值相同，即 m-m 0 换句话说，不希望有系统误差（所谓无偏性）。 此外，我们总是希望上述大多数误差的绝对值要小一些，并且在某一确定值周围波动，即估计误差的分布具有较小离散性 令为一个二阶平稳的随机函数，其期望为m，协方差为或变异函数为，且只依赖于向量h。 经过推导，估计方差的计算公式如下 （29） 上式可用平均变异函数表示如下 （210） 式（29），（210）中及分别代表当矢量的两个端点各自独立地扫过待估域V及信息域v的协方差函数平均值及变异函数平均值；及分别代表当矢量的两个端点各自独立地扫过任两个信息域v的协方差函数平均值及变异函数平均值；及分别代表当矢量的两个端点各自独立地在待估域V扫过时的协方差函数平均值及变异函数平均值。 当估计量是加权平均值时，估计方差的公式可表示如下 （211） 式中的，表示信息域，V表示待估域，，为，的权系数。或 （212） 2.4离差方差 令V是以点x为中心的开采面，并将其分成以为中心的N个大小相等的生产单元（）。 现在让我们再把v离散成若干个y点，其品位为，则每个以点为中心的单元的平均品位是 以x为中心的开采面V的平均品位是 显然，这N个品位值对它的平均值的离散程度可用其方差表示 （213） 当x固定，则与均为随机变量，而也是一个随机变量，从而可以讨论它的数学期望。至此，我们可以定义离差方差如下在区域化变量（点品位）满足二阶平稳假设条件下，把随机变量的数学期望定义为在开采面V内N个生产单元v的离差方差，记为 （214） 经过推导得到离差方差的通式 （215） 或 （216） 2.5变异函数及结构分析 为表征一个矿床金属品位等特征量的变化，经典统计学通常采用均值、方差等一类参数，这些统计量只能概括该矿床中金属品位等特征量的全貌，却无法反映局部范围和特定方向上地质特征的变化。地质统计学引入变异函数这一工具，它能够反映区域化变量的空间变化特征相关性和随机性，特别是透过随机性反映区域化变量的结构性，故变异函数又称结构函数。 变异函数（variogram）及变异曲线 我们可以把一个矿床看成是空间中的一个域V图2.1，V内的许多值则可以看成是V内一个点至另一个点的变量值，如图中x，xh为沿x方向被矢量h分割的两个点，其观测值分别为及，该两者的差值就是一个有明确物理意义的结构信息，因而可以看成是一个变量。 图2.1 域V内的变量值 区域化变量在空间相距h的任意两点x和xh处的值与差的方差之半定义为区域化变量的变异函数，记为 由上式可以看出，是依赖于x和h两个自变量的，与位置x无关，而只依赖于分隔两个样品点之间的距离h时，则可把变异函数写为 （217） 在实践中，样品的数目总是有限的，把有限实测样品值构制的变异函数称为实验变异函数（experimental variogram），记为 （218） 是理论变异函数值的估计值。 变异函数一般用变异曲线来表示。其常用的模型有球状模型、高斯模型及指数模型[10]。 1） 球状模型公式为 （219） 球状模型在地质采矿中应用最为广泛，如图2.2所示。 图2.2 变异函数曲线 该模型在h0处，作球状模型曲线的切线与总基台的交点的横坐标为，其中a值叫变程（图2.3）。 图2.3 球状模型的变程 2） 高斯模型其通式为 （220） 必须注意式中a不是变程，由于当时，，即当时，，所以，该模型的变程为。 3） 指数模型一般公式为 （221） 由于当h3a时，，即当h3a时，，所以该模型的变程约为3a。 图2.2是一个理想化的变异曲线图。图中的称为块金效应（nugget effect），它表示h很小时两点间品位的变化；a称为变程（range），当时，任意两点间的观测值有相关性，这个相关性随h的变大而减小，当时就不再具相关性，a的大小反映了研究对象（如矿体）中某一区域化变量（如品位）的变化程度，从另一个意义看，a反映了影响范围，例如可以用在范围a以内的信息值对待估域进行估计。C称为总基台值，它反映某区域化变量在研究范围内变异的强度，它是最大滞后距的可迁性变异函数的极限值。 2.6结构分析 在实际工作中区域化变量的变化很复杂，它可能在不同的方向上有不同的变化性,或者在同一方向包含着不同尺度上的多层次的变化性，因此无法用一种理论模型来拟合它，为了全面地了解区域化变量的变异性，就必须进行结构分析。所谓结构分析就是构造一个变异函数模型，对全部有效结构信息作定量化的概括，以表征区域化变量的主要特征。结构分析的主要方法是套合结构，就是把分别出现在不同距离h上和不同方向上同时起作用的变异性组合起来。套合结构可以表示为多个变异函数之和,每一个变异函数代表一种特定尺度上的变异性，其表达式为 （222） 在几个方向上研究区域化变量时，当一个矿化现象在各个方向上性质相同时称各向 同性，反之称各向异性，它表现为变异函数在不同方向上的差异。各向异性按性质可划 分为几何异向性和带状异向性两种。 图2.4 几何异向性 几何异向性，当区域化变量在不同方向上表现出变异程度相同而连续性不同时称为几何异向性，由于这种异向性可以通过简单的几何图形变换化为各向同性而得名。几何异向性具有相同的基台值C（设C00）而变程a不同（见图2.4）。 a b 图2.5 带状异向性 带状异向性，当区域化变量在不同方向上变异性之差，不能用简单的几何变换得到时，就称为“带状异向性”，这时，其具有不同的基台值C，而变程a可以相同或不同（见图2.5a，图2.5b）。 2.7 交叉验证 在进行结构分析后，得到结构模型，如何验证该结构模型是否正确反映矿床实际结构，可以采用交叉验证的办法。具体实施的方法是，依次拿掉一个已知值，然后用该结构模型和待估域周围的已知样品去估该已知值。然后把这些真值和估计值进行比较，不断改变结构参数，当平均相对误差（真值减估计值差的平均值）趋近于“0”，并且实际误差方差（真值和估计值的误差的方差）和理论克立格估计方差之比趋近于“1”时，所选的结构模型最好。 第3章 克立格法 3.1 概述 地质统计学主要是在结构分析的基础上采用各种克立格法来估计和解决实际问题，根据研究目的和条件不同，有各种各样的克立格法相继产生，如当区域化变量满足二阶平稳（或内蕴）假设时，可用普通克立格法；在非平稳条件下采用泛克立格法；为了计算局部可回采储量可用析取克立格法；当区域化变量服从对数正态分布时，可用对数正态克立格法；对有多个变量的协同区域化现象可用协同克立格法；对有特异值的数据可用指示克立格法等。最基本，应用最为广泛的是普通克立格法。 对于任何一种估计方法，都不能要求计算的平均品位估计值和它的实际值完全一样，也就是说，偏差 将是不可避免的，然而，我们要求一种估计方法应当满足下面两点 （1）所有估计块段的实际值与其估计值之间的偏差平均为0，即估计误差的期望应该等于0 （31） 我们称这种估计是无偏的。无偏是指平均说来品位的任何过高或过低的估计，以及由此引起矿石储量的过高或过低都是危险的，因此应该尽量避免。 （2）块段的估计品位与实际品位之间的单个偏差应该尽可能小，即误差平方的期望值（估计方差） （32） 应该尽可能的小。 因此，最合理的估计方法应当提供一个无偏估计且估计方差为最小的估计值。最常用的估计方法是用样品的加权平均求估计值的，也就是说，对于任一待估块段V的真实值的估计值是通过该待估块段影响范围内n个有效样品值（）的线性组合得到的 （33） 式中 是加权因子，是各样品在估计时的影响大小，而估计方法的好坏就取决于如何计算或选择权因子。 对于给定的待估块段V和用来进行估计的一组信息{，}，我们要求出一组权系数（），若使估计方差为最小，则块段V的估计值就能在最小的可能置信区间内产生，而给出最佳、线性、无偏估计的权系数的方法就是克立格法（Kringing）。克立格法是一种最佳局部估计方法，它以最小的估计方差给出块段平均值的无偏线性估计量，即所谓克立格估计量（Kringing estimator）。 所谓局部估计，就是在一个有限的估计邻域内求出某待估块段的最佳估计量，这个估计邻域应该小于矿床的准平稳（均匀）带，因此，所谓最佳局部估计就是要找出一个准平稳带内待估块段的平均品位的最佳估计量。克立格方法就是把矿体划分成许多小块段（待估块段），根据待估块段周围有限邻域内的信息逐块估计，因此，克立格是一个加权滑动平均法，而全体矿体的总体估计是通过对其逐个块段的局部估计的组合而得到。 广义地说，克立格法是一种求最优、线性、无偏内插值估计量的方法（Best Linear Unbiased Estimator，简记BLUE），而具体地说，克立格法就是在考虑了信息样品的形状、大小及其与待估块段相互间的空间分布位置等几何特征以及品位的空间结构之后，为了达到线性、无偏和最小估计方差的估计，而对每一样品值分别赋予一定的权系数，最后进行加权平均来估计块段品位的方法。 显然，克立格法最重要的工作是第一，列出并求解克立格方程组，以便求出诸克立格权系数；第二，求出这种估计的最小估计方差克立格方差。 设是被研究的定义在点支撑上的区域化变量，且假定服从二阶平稳假设，即有期望，及中心化协方差函数 ；或变异函数。现在要求对中心位于的域的平均值 进行估计（图3.1），而在待估域V的周围有一组信息值{，}，在二 阶平稳下，它们的期望，则待估域V的实际值的估计值是这n个有效 数据（）的线性组合为式（33）。 我们的目的是求出式中的n个权系数（），以便保证估计量 无偏，且估计方差最小。由这样的权系数计算出的估计量称为的克立格估计量，而最小估计方差称为克立格方差（Kringing variance），记为。 图3.1 用7个信息样估计待估域V 3.2 普通克立格法 前已述及，对式（33），目的是求出式中的n个权系数（），使为的无偏估计量，且其估计方差最小。 （１） 无偏条件 为了避免系统误差，我们要求估计无偏，即。在二阶平稳条件下,而，要使，即就必须使 （34） 式（34）就是无偏条件。 （２） 普通克立格方程组（ordinary Kringing systems ） 估计方差的计算公式 要求出在无偏条件下，使估计方差 达到极小的诸权系数（）是个求条件极值的问题，即把最优估值问题理解为无偏条件约束下，求目标为估计方差为最小的估值。 为了便于求解，可将极小化为无约束的拉格朗日乘数法求极值的问题，即将约束条件也引入目标函数之中，令 F是n个权系数（）和μ的（n1）元函数，-2μ是拉格朗日乘数（目的是为了后面得到的方程较为简单）。求出F对n个（）和μ的偏导数并令其为0 （） 最后得到普通克立格方程组 （） （35） 它是一个n1个未知数（n个和一个），n1个方程的方程组。在内蕴假设下，式（3--5）可用表示如下 （） （36） （３） 普通克立格方差（Kringing variance） 从式（35）得 将上式代入估计方差公式 由于用上式计算的估计方差为最小估计方差，亦称克立格方差，记为 （37） 若用表示，则式（36）改写为 （3--8） 式中的符号表示如下为分隔矢量h的两个端点分别独立地在域V内移动时，求出区域化变量全部协方差的平均值 为矢量h的两个端点分别独立地在V及中移动，求出区域化变量的全部协方差的平均值 为矢量h的两个端点分别独立地在信息域，中移动，求出区域化变量全部协方差的平均值 （４） 矩阵形式 为了便于书写，克立格方程组与克立格方差均可用矩阵形式表示如下 （39） 或 （310） （3--11） 式中 [K]称为普通克立格矩阵，其中由于有 故[K]为一对称矩阵。 类似地，也可以用表示如下 （312） （313） 式中 3.3泛克里格法（Universal Kriging） 普通克里格法要求区域化变量是二阶平稳的，至少是准平稳或准内蕴的，这时，至少在估计邻域内满足成立。然而，在生产实践中，有些区域化