1、2018/10/14,1,地统计学方法,资源与环境学院 杨勇,2018/10/14,2,华中农业大学 资源与环境学院,设想一下这样的问题,？,这块地的土壤养分情况如何？ 不仅需要知道一个总体情况 而是要知道每个地方的不同含量 方便为那些含量低的地方施肥 该怎么办呢？,2018/10/14,3,华中农业大学 资源与环境学院,方案一,Step1: 密集采样 Step2: 把土样运回实验室 Step3: 晒干，磨碎，化学分析耗时，耗力，耗财 得到的是点状数据 面状连续分布呢？ 未采样地的状况如何呢？,2018/10/14,4,华中农业大学 资源与环境学院,方案二,算法分析,2018/10/14,5,
2、华中农业大学 资源与环境学院,实例：,（a）有机质,（b）全氮,（c）有效磷,2018/10/14,6,华中农业大学 资源与环境学院,1.1 地统计学的发展和概念,一、地统计学发展简史 地统计学（Geostatistics）是20世纪50年代初在南非采矿业中为了计算矿石储量而发展应用起来的，首先被采矿工程师Krige和统计学家Sichel应用于南非的采矿工作中。50年代后期，法国Matheron在此基础上提出了区域化变量理论，形成了地统计学的基本框架。,2018/10/14,7,华中农业大学 资源与环境学院,地统计学发展简史,70年代，计算机的出现，这项技术被引入到地学领域。1975年在罗马举
3、行了关于该学科的第一个国际性会议后，陆续有多个相关国际会议举行。我国的地统计学研究和应用是1977年由侯景儒、黄竞先等首先进行的。现已广泛运用于地质、土壤、农业、气象、海洋、生态、森林和环境治理等方面,2018/10/14,8,华中农业大学 资源与环境学院,二、地统计学的概念,定义：地统计学是以区域化变量理论为基础，以变异函数为主要工具，研究那些在空间分布上既有随机性又有结构性，或空间相关性和依赖性的自然现象的科学。（王政权，1999）,2018/10/14,9,华中农业大学 资源与环境学院,1.2 地统计学的应用（土壤）,土壤属性的空间分布特征是土壤污染治理、土地管理和现代农业的重要依据之一
4、。 土壤是一个形态和过程都相当复杂的自然综合体，成土过程中不同的物理、化学、生物等因素的影响，使得土壤性质具有高度的空间异质性。人类活动进一步加剧了土壤属性的变异性和不确定性。 同时，土壤本身处于一个时刻变化的动态过程，因此，对土壤空间性质进行描述和定律研究相当困难。,2018/10/14,10,华中农业大学 资源与环境学院,1.2 地统计学的应用（土壤）,自上世纪七八十年代地统计学引入土壤学研究中以来，随着学科发展和应用方向的扩展，地统计学方法已经成为土壤学特别是大尺度土壤学研究的一个重要工具。 地统计学在土壤物理性质空间变异中的应用 地统计学在土壤化学性质空间变异中的应用 地统计学在土壤重
5、金属污染空间变异中的应用 地统计学在采样策略中的应用 地统计学在其他特性中的应用,2018/10/14,11,华中农业大学 资源与环境学院,地统计学在土壤物理性质空间变异中的应用,湖北咸宁 据：罗勇，陈家宙，2008,土壤容重空间变异,土壤饱和导水率空间变异,2018/10/14,12,华中农业大学 资源与环境学院,地统计学在土壤化学性质空间变异中的应用,（a）有机质,（b）全氮,（c）有效磷,（d）速效钾,湖北沙洋 据：杨勇，贺立源，2010,2018/10/14,13,华中农业大学 资源与环境学院,地统计学在土壤重金属污染空间变异中的应用,武汉市东湖高新技术开发区 据：张贝，杨勇，2010
6、,2018/10/14,14,华中农业大学 资源与环境学院,1.3 地统计学在土壤科学中的应用展望,地统计学和土壤多源数据的处理 利用多源的相关数据预测目标属性的分布 地统计学和土壤过程的空间建模 利用多源数据模拟土壤发生发展的过程 地统计学和土壤特性的不确定性模拟 土壤属性超过某一阈值的概率 地统计学和土壤过程的时空变异 地统计学与精确农业 土壤综合特性的空间变异性研究 ,2018/10/14,15,华中农业大学 资源与环境学院,样本数据的统计分析和预处理,描述性统计 频数分布：直方图 集中趋势的度量：平均数、中位数、众数 离散型度量：极差、方差 偏度和峰度 数据检验和分布分析 异常值的识别
7、和处理：平均值加标准差法、四倍法 正态分布的检验方法：直方图法、PP、QQ、 数据转换处理：对数转换、平方根转换、反正弦转换 相关分析和回归分析 回归分析 相关分析,2018/10/14,16,华中农业大学 资源与环境学院,区域化变量,当一个变量呈空间分布时，称之为“区域化”。这种变量常常反映某种空间现象的特征，用区域化变量描述的现象称之为区域化现象。如生态学、土壤学和地质学中许多研究的变量都具有空间分布的特点，实质上都是区域化变量。 在研究区域内所有点处的样品数据的实测值就是一个区域化值，其相应的函数z(x)就是一个区域化变量，也是该区域随机模型（函数）Z(x)的一个实现。,2018/10/
8、14,17,华中农业大学 资源与环境学院,平稳假设,1、平稳性：表示当将既定的n个点的点集从研究区域某一处移向另一处时，随机函数的性质保持不变，也称为平移不变性。即随机函数分布的规律性不因位移而改变，是严格平稳的，具有平稳性。,2018/10/14,18,华中农业大学 资源与环境学院,二阶平稳性假设,2、二阶平稳性假设（弱平稳性假设）：随机函数的均值为一常数，且任何两个随机变量之间的协方差依赖于它们之间的距离和方向，而不是它们的确切位置： 条件1：,数学期望：反映随机变量取值的集中特征，是随机变量取得数字的代表数。该条件表示：在整个研究区内，区域化变量的数学期望对任意x存在，且等于常数,201
9、8/10/14,19,华中农业大学 资源与环境学院,二阶平稳性假设,条件2：在整个研究区内，区域化变量的协方差函数对任意x和h存在，且平稳，即：,协方差：两个不同参数之间的方差就是协方差，用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况，即当两个变量是相同的情况。期望值分别为E(X) = 与 E(Y) = 的两个实数随机变量X与Y之间的协方差定义为：COV(X，Y)=E(X-E(X)(Y-E(Y) ，若两个随机变量X和Y相互独立 ，则他们的协方差为0。,2018/10/14,20,华中农业大学 资源与环境学院,本征假设,条件1：条件2：,r(h)称为半方差函数，也叫变异函数 本征假设是
10、地统计学中对随机函数的基本假设 事实上，当作用于大区域时，本征假设的第一个条件很难满足，空间变异的漂移或趋势面可能存在，由于这种漂移，第二个条件也不能满足，但地统计学理论的基础是本征假设，因此，有必要去认识一个随机过程是否是平稳性的,在研究区域内，区域化变量Z(x)的增量的数学期望对任意x和h存在且等于0,在研究区域内，区域化变量的增量Z(x)-Z(x+h)的方差对任意x和h存在且平稳,2018/10/14,21,华中农业大学 资源与环境学院,平稳假设,就严格性而言： 平稳性假设二阶平稳性假设本征假设本征假设是地统计学中对随机函数的基本假设,2018/10/14,22,华中农业大学 资源与环境
11、学院,变异函数和协方差函数,变异函数和协方差函数存在以下关系：,2018/10/14,23,华中农业大学 资源与环境学院,协方差具体计算方法,设Z(x)为区域化随机变量，并满足二阶平稳条件，h为两样本点空间分割距离，Z(xi)和Z(xi+h)分别是Z(x)在空间位置xi和xi+h上的观测值，则协方差函数的计算公式为：,N(h)是分隔距离为h时的样本对数总数,2018/10/14,24,华中农业大学 资源与环境学院,变异函数具体计算方法,公式：,值分别是：4，3，4，5，7，9，7，8，7，7，则：,2018/10/14,25,华中农业大学 资源与环境学院,2018/10/14,26,华中农业大
12、学 资源与环境学院,2018/10/14,27,华中农业大学 资源与环境学院,变异函数散点图,2018/10/14,28,华中农业大学 资源与环境学院,变异函数的理论拟合模型,理论变异函数用来拟合一些列经验变异函数值，供后续进行插值估计时使用。选用理论变异函数模型是，要根据经验半方差图的性状来选取合适的模型,2018/10/14,29,华中农业大学 资源与环境学院,变异函数的理论拟合模型,变异函数的理论模型：有基台值模型无基台值模型,2018/10/14,30,华中农业大学 资源与环境学院,有基台值模型球状模型,C0：块金常数 C0+C ：基台值 C：拱高 a：变程 应用最广的模型,2018/
13、10/14,31,华中农业大学 资源与环境学院,有基台值模型指数模型,C0：块金常数 C0+C ：基台值 C：拱高 3a：变程 当C0=0,C=1时，称为标准指数函数模型,2018/10/14,32,华中农业大学 资源与环境学院,有基台值模型高斯模型,C0：块金常数 C0+C ：基台值 C：拱高：变程 当C0=0,C=1时，称为标准高斯函数模型,2018/10/14,33,华中农业大学 资源与环境学院,三种常用模型比较,0.95,2018/10/14,34,华中农业大学 资源与环境学院,有基台值模型线性有基台值模型,C0：块金常数 C0+C ：基台值 C：拱高 A ：常数，表示直线斜率 当C0
14、=0,C=1时，称为标准指数函数模型,2018/10/14,35,华中农业大学 资源与环境学院,有基台值模型纯块金效应模型,2018/10/14,36,华中农业大学 资源与环境学院,无基台值模型线性无基台值模型,2018/10/14,37,华中农业大学 资源与环境学院,无基台值模型幂函数值模型,2018/10/14,38,华中农业大学 资源与环境学院,无基台值模型对数值模型,2018/10/14,39,华中农业大学 资源与环境学院,套合模型,在实际中，有时区域化随机变量Z(x)的变化相当复杂，往往包含各种尺度及各种层次的变化，反映在变异函数r(h)上，就是单一的模型结构不能将其合理表达，而是多
15、层次的结构相互叠加在一起，地统计学上称为套合。所谓套合结构，就是把分别出现在不同距离h上或不同方向上同时起作用的变异性组合起来，对全部有效的结构信息，作定量化的概括，以表示区域化变量的主要特征。,2018/10/14,40,华中农业大学 资源与环境学院,套合模型,土壤是一个不均与、具有高度空间异质性的复合体，它与土壤母质、气候、水文、地形和生物等因素有关，分析土壤空间变异的因素，可将其变异分为系统变异（土壤形成因素相互作用造成）和随机变异（可以观测到的，但与土壤形成印务无关且不能直接分析的）两大类。如由h分开的两个点x和x+h的土壤某一性质Z(x)和Z(x+h)。当h趋近于0时，可以认为两点间
16、的差异完全是由取样和测定误差造成，当h逐步增大，如h1m，差异可能还要加上诸如水分等因素，当h100m时，在新的变异要考虑地形的作用。,2018/10/14,41,华中农业大学 资源与环境学院,套合模型,当h一定时，变异函数r(h)应包含小于h的所有影响因素，因此，绝大多数变异函数都由下面两个变异函数组成：r(h)=r0(h)+r1(h)，即一个代表纯块金方差，一个代表空间相关的方差。一般情况下，套合模型可以用放映各种不同尺度变化的多个变异函数之和表示，即：,ri(h)可以是相同的或不同的理论模型,2018/10/14,42,华中农业大学 资源与环境学院,套合模型,如，区域化变量Z(x)的变异
17、性由r0(h),r1(h)和r2(h)组成，其中,2018/10/14,43,华中农业大学 资源与环境学院,套合模型,三者组成的套合模型为：,2018/10/14,44,华中农业大学 资源与环境学院,套合模型,2018/10/14,45,华中农业大学 资源与环境学院,最优拟合参数最优估计,变异函数的理论模型主要是曲线模型，将曲线模型经过适当的变换，化为线性模型，然后用最小二乘法原理求未知参数的估计。,2018/10/14,46,华中农业大学 资源与环境学院,基于优化搜索算法的参数拟合,对于结构复杂的变异函数理论模型，特别是套合结构模型，参数复杂，难以用一般的通用方法求解出模型中的参数。但一些智
18、能优化算法，如遗传算法、模拟退火算法、蚁群算法能够使用统一的流程求解出接近最优的参数。,2018/10/14,47,华中农业大学 资源与环境学院,基于遗传算法的变异函数理论模型参数估计,1、多尺度套合模型的规范表达,2018/10/14,48,华中农业大学 资源与环境学院,基于遗传算法的变异函数理论模型参数估计,从上式可以看出，需求解的参数为2n+1个（因为第一个模型总是纯块金模型）。而在实际计算时，可以令 ，这样方便从经验半方差图中识别 ci取值区间。并有以下约束：,2018/10/14,49,华中农业大学 资源与环境学院,基于遗传算法的变异函数理论模型参数估计,编码策略及初始群体产生 假设
19、需要顾及m(m=2n+1)个参数，每个参数的取值范围和估值精度分别是Umin,Umax和Qi，则将m个参数分别以L1,L2,Lm为长度进行二进制编码，其中则每条染色体长度为 ，染色体中每个参数编码对应的解码 公式为：,以这种编码方式随机产生T组染色体,2018/10/14,50,华中农业大学 资源与环境学院,基于遗传算法的变异函数理论模型参数估计,确定个体适应度评价函数,2018/10/14,51,华中农业大学 资源与环境学院,基于遗传算法的变异函数理论模型参数估计,遗传操作 遗传算法主要包括3个基本算子，即选择、交叉和变异，为此，需确定交叉概率Pc和变异概率Pm，3个过程执行以后，将产生新一
20、代种群，并记录适应度最高的染色体,2018/10/14,52,华中农业大学 资源与环境学院,主要空间插值法简介,分类： 确定性方法：基于实测数据的相似性程度或平滑程度，利用数学函数进行插值（如逆距离加权法） 地统计方法：利用实测数据的统计特性来量化其空间自相关程度，生产插值面并评价预测的不确定性,2018/10/14,53,华中农业大学 资源与环境学院,主要空间插值法简介,分类： 整体插值法：利用整个实测数据集来预测局部插值法：在大面积的研究区域上选取较小的空间单元，利用预测点周围的临近样点来进行预测,2018/10/14,54,华中农业大学 资源与环境学院,空间整体插值法,1、全局多项式插值
21、法（趋势面分析法）：即用数学公式表达感兴趣区域上的一种渐变的趋势。 平面： 曲面： 多项式中的参数系数往往用最小二乘法求解。但该方法是不精确的插值方法，很少有实测点刚好在生产的插值面上，而是或高或低于插值面，高低数值相加，之和近似为0。,2018/10/14,55,华中农业大学 资源与环境学院,空间整体插值法,全局多项式插值法的插值结果往往呈条带状（左图），适合于描述那些呈明显趋势分布的属性，不适合描述那些空间分布波动较大（较破碎，右图）的自然属性,2018/10/14,56,华中农业大学 资源与环境学院,空间整体插值法,2、变换函数插值法：根据一个或多个空间参量的经验方程进行整体空间插值，这
22、种经验方程称为变换函数。即用与被预测属性相关的其他属性建立回归方程，进行空间预测：,b0,b1,b2为回归系数，p1，p2为独立空间变量，z(x)为被预测属性,2018/10/14,57,华中农业大学 资源与环境学院,空间局部插值法,1、泰森多边形插值：由一组连续多边形组成，多边形的边界是由相邻两点直线的垂直平分线组成。,特性： （1）每个多边形内仅包含一个离散数据点。 （2）在多边形内的任一点k(x，y)同Pi(xi，yi)之间距离总小于它同其它离散点Pj(xj，yj)之间距离。 （3）泰森多边形的任意一个顶点必有三条边与它连接，这些边是相邻三个泰森多边形两两拼接的公共边。 （4）泰森多边形
23、的任意一个顶点周围存在三个离散点，将其连成三角形后其外接圆的圆心即为该顶点，该三角形称泰森三角形,2018/10/14,58,华中农业大学 资源与环境学院,空间局部插值法,各泰森多边形内的每一点属性均由各多边形内的已知点确定，若求数据域内任意一点数据属性Z（xi，yi），则需首先判断待求点所落入的多边形，然后再由控制该多边形的已知点Z（x，y）推算得到。,2018/10/14,59,华中农业大学 资源与环境学院,空间局部插值法,2、三角测量插值法：将采样点用直线与其相邻点连接成三角形，三角形内部包括任何样点，形成一个包括多个倾斜三角板的多面体（TIN）,未测点只可能在三角形内或三角形边线上，利
24、用线性插值即可求得 缺点是每个预测值只是根据三个实测值得到，且有时会产生突变现象,2018/10/14,60,华中农业大学 资源与环境学院,空间局部插值法,3、逆距离加权法（IDW）：利用被预测区域点周围的实测值来预测未采样点的值，实测点离预测点越近，则对插值的结果影响越大。 其中p为实测值对预测值的影响级，若p=0，则每一个权重是一样的，预测值是所有实测值的平均值，当p增加时，相距较远的点的权重迅速减小，2最为常用。 由于IDW方法只考虑距离进行权重分配，所以临近实测点的贡献往往很大，而造成空间分布的多点中心现象。,2018/10/14,61,华中农业大学 资源与环境学院,空间局部插值法,4
25、、局部多项式插值法（移动内插法）：多项式插值法将整个区域考虑成一个平面或曲面，而局部多项式插值法是在划定的领域内（窗口内）用其中的实测数据来拟合不同次数的多项式。,2018/10/14,62,华中农业大学 资源与环境学院,空间局部插值法,5、简单移动平均法：,6、样条插值法：,2018/10/14,63,华中农业大学 资源与环境学院,空间局部插值法,7、克里格方法：和IDW一样，也是一种局部估计的加权平均，但是它对各实测点权重的确定是通过半方差分析获取的，可分为线性克里格法和非线性克里格法。 （1）普通克里格 （6）概率克里格 （2）简单克里格 （7）贝叶斯克里格 （3）泛克里格 （8）普通协
26、同克里格 （4）指示克里格 （5）析取克里格,2018/10/14,64,华中农业大学 资源与环境学院,克里格法实质上是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏最优估计的一种方法，从数学角度讲就是一种对空间分布的数据求线性最优无偏内插估计量的一种方法。是根据待估样点有限领域内若干已测定的样点数据，在考虑样点形状、大小和空间相互位置关系，它们与待估样点相互空间位置关系，以及变异函数提供的结构信息之后，对该待估样点进行的一种线性无偏最优估计,2018/10/14,65,华中农业大学 资源与环境学院,普通克里格法,假定Z(x)是满足本证假设的一个随机过程
27、，该随机过程有n个观测值z(xi)，要预测未采样点x0处的值，则线性预测值Z*(x0)可以表示如下：Kriging是在使预测无偏并有最小方差的基础上，去确定最优的权重值，满足以下两个条件： （1）无偏性条件 （2）最优条件：,2018/10/14,66,华中农业大学 资源与环境学院,普通克里格法,在本证假设条件下，上左边的式子可以表示为：根据方差最小原则，借助拉格朗日乘子，普通克里格的预测方程组为：预测方差为：,2018/10/14,67,华中农业大学 资源与环境学院,普通克里格法,克里格公式也可以用矩阵的形式表示，对点状克里格，有：,2018/10/14,68,华中农业大学 资源与环境学院,
28、普通克里格法 实例,2018/10/14,69,华中农业大学 资源与环境学院,普通克里格法 实例,2018/10/14,70,华中农业大学 资源与环境学院,2018/10/14,71,华中农业大学 资源与环境学院,2018/10/14,72,华中农业大学 资源与环境学院,2018/10/14,73,华中农业大学 资源与环境学院,简单克里格法,如果我们知道区域随机变量的平均值，那么我们可以利用这种先验知识通过简单克里格法来提高预测的精度，这种克里格预测方法仍然是线性加和，但将随机过程的平均值包括了进去，这种随机过程必须是二阶平稳的，预测公式为：,2018/10/14,74,华中农业大学 资源与环
29、境学院,简单克里格法,权重利用以下公式计算：用矩阵形式表示为： 其中：则： 预测方差为：,2018/10/14,75,华中农业大学 资源与环境学院,协同克里格,协同克里格是利用两个变量之间的互相关性，用其中易于观测的变量对另一变量进行局部估计的方法。协同克里格法比普通克里格法能明显改进估计精度及采样效率。但在实际应用中，协同克里格法要求有一个已知的相关函数，这就需要在很多地点同时采样，测定二个函数间的相互关系。与相关函数一样，这种相互关系也受样本数目多少的影响。,2018/10/14,76,华中农业大学 资源与环境学院,协同克里格,协同克里格法是建立在协同区域化变量理论基础上的，通过建立交叉协
30、方差函数和交叉变异函数模型，然后用协同克里格法对未抽样点的变量进行估值,2018/10/14,77,华中农业大学 资源与环境学院,协同区域化的概念,在实际中，每一种区域化现象都与许多变量有关，同一个区域化现象可以用几个相关变量表示。在地统计中，把某一点上某一性质的观测值，与在统计分析上依赖于相邻一点上的另一性质的观测值，这两种性质之间的相关性称为协同区域化或横相关。描述这种协同区域化现象的变量称为协同区域化变量。,2018/10/14,78,华中农业大学 资源与环境学院,协同区域化的概念,研究协同区域化变量具有许多优点，例如，在土壤空间变异性研究中，有些土壤性质的测定难度和费用较高，而另一些土
31、壤性质的测定相对简单易行。因而，可以用较容易测定的土壤某一性质之值去估计另一种测定难度大、费用高的土壤性质之值的变化。 此外，在空间变异分析中，如果能用一种变量的信息去弥补所遗漏或提供另一变量的信息，这无疑是非常有意义的。协同区域化变量的莅临和方法将提供解决两个变量空间相关和估计的问题。,2018/10/14,79,华中农业大学 资源与环境学院,协同区域化变量理论,设K个协同区域化变量Z1(x),Z2(x),Zk(x)，组成一组K维区域化变量的向量Z1(x),Z2(x),Zk(x)，在观测前，它是一个K维区域化变量，观测后，它可以看成K维向量的一个实现。在二阶平稳假设条件下，协同区域化变量有：
32、 （1）每一个Zk(x)（K=1,2,K）的数学期望存在且平稳，即,2018/10/14,80,华中农业大学 资源与环境学院,协同区域化变量理论,（2）对每对区域化随机变量的交叉协方差函数为：（3）在满足本征假设条件时，区域化变量的增量数学期望为0，则每对区域化协同变量的交叉变异函数存在，为：,2018/10/14,81,华中农业大学 资源与环境学院,交叉协方差函数和交叉变异函数的计算公式,设在点x和点x+h处，分别测得两个变量的观测值Zk(x)，Zk(x)，Zk(x+h)，Zk(x+h)，则交叉协方差函数计算公式为：其中N(h)为样本对数,2018/10/14,82,华中农业大学 资源与环境
33、学院,交叉协方差函数和交叉变异函数的计算公式,交叉变异函数的计算公式为：,2018/10/14,83,华中农业大学 资源与环境学院,协同克里格法,Zvk0的估计量为Zvk0#,是K个协同区域化随机变量全部有效数据（观测值）的线性组合：以2个协同变量来说明克里格估计方程组和协同克里格估计方差,2018/10/14,84,华中农业大学 资源与环境学院,协同克里格线性方程组,设在点x0处的某变量的平均值为u0，在x0附近有两个已观测的协同区域化随机变量ui（i=1,2,n）和vj（j=1,2,m）。则平均值u0的估计值u0#构成的协同克里格线性估计量为：其中ai，bj为协同克里格权重系数，为使u0#
34、为u0的最优无偏线性估计量，必须满足：,2018/10/14,85,华中农业大学 资源与环境学院,协同克里格线性方程组,1 无偏性条件：只有当 时，才能成立，因此，它就是无偏性条件,2018/10/14,86,华中农业大学 资源与环境学院,协同克里格线性方程组,2 最优性条件 在满足无偏性条件下，协同克里格估计方差为：,2018/10/14,87,华中农业大学 资源与环境学院,协同克里格线性方程组,设,2018/10/14,88,华中农业大学 资源与环境学院,则为使协同克里格估计方差最小，令求上式的偏导数并令其为0，得到协同克里格线性方程组,2018/10/14,89,华中农业大学 资源与环境
35、学院,2018/10/14,90,华中农业大学 资源与环境学院,经整理得2个变量的协同克里格线性方程组的一般表达式：,2018/10/14,91,华中农业大学 资源与环境学院,这是一个n+m+2阶线性方程组，解该方程组得到协同克里格权重系数ai和bj，然后代入 中，得到协同克里格线性无偏最优估计量。此时，协同克里格估计方差为：,2018/10/14,92,华中农业大学 资源与环境学院,实例,2个协同区域化随机变量u和V，其中u0为待估样点，在其周围有u1,u2,V1,V2和V3已知样本,2018/10/14,93,华中农业大学 资源与环境学院,实例,根据已知的理论模型计算两个变量的协方差函数C
36、u(h)和Cv(h)及交叉协方差函数Cuv(h),2018/10/14,94,华中农业大学 资源与环境学院,则协同克里格线性方程组为,解上述方程组得a1=0.512，a2=-0.216，b1=0.488，b2=-0.397，b3=0.666，u1=205963 u2=13823，将这些协同克里格权重系数代入协同克里格线性估计量方程得u0的估计值为356，其协同克里格估计方差为681549。,2018/10/14,95,华中农业大学 资源与环境学院,泛克里格法,漂移的概念 在普通的克里格法中，要求区域化变量Z(x)是二阶平稳的或本证的，至少是准平稳或准本证假设条件，在有限的估计领域内Z(x)的数
37、学期望是一个常数，即EZ(x)=m存在。然而在许多情况下，区域化变量在研究区域内是非平稳的，其数学期望不是一个常数，即EZ(x)=m(x)，m(x)在地统计学上被定义为非平稳区域化变量的漂移，其表达式为 m(x)=EZ(x),2018/10/14,96,华中农业大学 资源与环境学院,泛克里格法,在漂移存在的条件下就不能用普通克里格方法进行空间局部估计，而要采用泛克里格法进行估计。漂移一般采用多项式表示：其中fi(x)为一已知多项式函数，ai为未知系数。当漂移为线性时，一维和二维条件下漂移m(x)的形式为：,2018/10/14,97,华中农业大学 资源与环境学院,泛克里格法,当漂移为二次（非线
38、性）时，一维和二维条件下漂移m(x)的形式为由于漂移的存在，泛克里格法在估计某一点Z(x)的估计值Z#(x)时，首先要估计该点上漂移m(x)的估计值m#(x)，这就要求在某种假设条件下确定非平稳区域化变量的协方差函数和变异函数,2018/10/14,98,华中农业大学 资源与环境学院,泛克里格的基本假设,设Z(x)是一个非平稳的区域化变量，可表示为：一般情况下，上述假设很难满足，因此，考虑Z(x)的增量的情况。 设Z(x)的增量具有非平稳数学期望和非平稳方差函数，并可表示为：,2018/10/14,99,华中农业大学 资源与环境学院,泛克里格的基本假设,如果非平稳区域化变量Z(x)可以分解成两
39、部分，一部分是在较大尺度下可以观察到现象变化m(x)，另一部分是在较小尺度下的变化R(x)，即在给定的尺度下，m(x)可以表示为一个多项式，即单项式函数fi(x)的线性组合：其中x为领域内任一点，ai为未知系数，如果包含x0点，则ai应有n+1个,2018/10/14,100,华中农业大学 资源与环境学院,非平稳条件下的协方差和变异函数,当Z(x)=m(x)+R(x)时，Z(x)的协方差函数CZ(x,y)为而Z(x)的变异函数为,2018/10/14,101,华中农业大学 资源与环境学院,即Cz(x,y)=CR(x,y)，Z(x)的协方差等于R(x)的协方差 rz(x,y)=rR(x,y)，Z
40、(x)的变异函数等于R(x)的变异函数如y=x+h，则rz(h)=rR(h)。只要能求出R(x)的变异函数rR(h)，就可以求得Z(x)的变异函数rz(h)。但m(x)一般为未知多项式函数，无法用R(x)=Z(x)-m(x)来计算rR(h)。为解决这个问题，地统计学泛克里格法中先对R(x)的变异函数进行估计，也就是先估计R#(x)=Z(x)-m#(x)来，即对m#(x)进行估计，根据R#(x)的变异函数rR#(h)与理论的rR(h)进行比较，当rR#(h)=rR(h)时，就可以求Z(x)的变异函数rR(h)来，因此，泛克里格方法估计有两个部分，一个是m(x)的估计，另一个是Z(x)的估计,20
41、18/10/14,102,华中农业大学 资源与环境学院,漂移m(x)的泛克里格法估计,在研究区域内，设Z(x)是一个非平稳区域化随机变量，并满足以下条件：,2018/10/14,103,华中农业大学 资源与环境学院,漂移m(x)的泛克里格法估计,设在研究区内有n个已知样点xa(a=1,2,n），其观测值为Z(xa)=Za，目的是用这些已知样点估计区域内任一固定点x处的漂移值m(x)。设m(x)的估计量为m#(x)，它可以表示这几个已知样点数据的线性组合，即,2018/10/14,104,华中农业大学 资源与环境学院,漂移m(x)的泛克里格法估计,为使m#(x)成为m(x)的无偏最优估计，则必须
42、满足： 1 无偏性，即只有当 时，无偏性存在，此时L=0,1,2,n,2018/10/14,105,华中农业大学 资源与环境学院,漂移m(x)的泛克里格法估计,2 最优性 m#(x)估计m(x)的方差可表示为即，上式达到最小,2018/10/14,106,华中农业大学 资源与环境学院,漂移m(x)泛克里格线性方程组,在满足无偏性条件下，用拉格朗日乘数法求方差的最小值得到漂移m(x)的克里格线性方程组可矩阵表达为,2018/10/14,107,华中农业大学 资源与环境学院,其中,2018/10/14,108,华中农业大学 资源与环境学院,这是一个具有n+K+1个未知数的n+K+1个方程组成的方程
43、组，称为漂移m(x)泛克里格线性方程组，泛克里格估计方差为：,2018/10/14,109,华中农业大学 资源与环境学院,Z(x)的泛克里格法估计,在研究区域内设Z(x)是一个非平稳区域化变量，并满足以下条件：,2018/10/14,110,华中农业大学 资源与环境学院,Z(x)的泛克里格法估计,设在研究区内有n个已知样点xa(a=1,2,n），其观测值为Z(xa)=Za，目的是用这些已知样点估计区域内任一固定点x处的漂移值Z(x)。设Z(x)的估计量为Znk#(x)，它可以表示这几个已知样点数据的线性组合，即,2018/10/14,111,华中农业大学 资源与环境学院,Z(x)的泛克里格法估
44、计,为使Z#(x)成为Z(x)的无偏最优估计，则必须满足： 1 无偏性,只有当,2018/10/14,112,华中农业大学 资源与环境学院,Z(x)的泛克里格法估计,2 最优性,上式（方差）达到最小,2018/10/14,113,华中农业大学 资源与环境学院,Z(x)的泛克里格线性方程组,在满足无偏性条件下，根据拉格朗日乘数法求方差的最小值，得到Z(x)的克里格线性方程组：,矩阵形式：,2018/10/14,114,华中农业大学 资源与环境学院,2018/10/14,115,华中农业大学 资源与环境学院,指示克里格法,在有的情况下，我们需要知道某区域化变量Z(x)在某地超过阈值z的概率：指示克
45、里格是实现方法之一,2018/10/14,116,华中农业大学 资源与环境学院,指示克里格法,指示码：仅有0和1两个值，即可表示某种物质存在和不存在，也可表示连续变量是否大于某个阈值如果z(x)是一个随机过程Z(x)的实现，那么 可以被认为是一个指示随机函数 的实现,2018/10/14,117,华中农业大学 资源与环境学院,指示克里格法,通过这种转换，某种程度上使预测结果更接近实际应用。如评价土壤重金属污染问题，通过设定合理的污染浓度阈值z，那么就可以将一个连续的随机变量Z(x)转化为一个指示函数，对这个指示函数而言，1表示没有受到污染，可以被接受，0表示受到污染，不能被接受。,2018/1
46、0/14,118,华中农业大学 资源与环境学院,指示克里格法,也可以定义多个阈值，并且为每个阈值生成一个新的指示变量。如定义S个阈值为zc(s),s=1,2,S,从而得到多个指示变量：这些可以看做是相应的随机函数的实现，,2018/10/14,119,华中农业大学 资源与环境学院,指示半方差函数,其半方差值可以通过转化后的指示数据计算得到：,后续的插值方法和普通克里格法一样，预测位置得到的插值结果为一个0和1之间的数，表示,2018/10/14,120,华中农业大学 资源与环境学院,指示克里格应用实例,以鲁西北禹城地区土壤重金属Pb含量为例，针对其存在特异值与偏态分布的特点，利用指示克里格法研
47、究了土壤Pb含量小于特定阈值的条件概率分布,县级采样尺度396个点（2KM）镇级采样尺度285个点（0.5KM）,2018/10/14,121,华中农业大学 资源与环境学院,指示克里格应用实例,2018/10/14,122,华中农业大学 资源与环境学院,取土壤Pb含量0.10.9分位数共9个值作为阈值T，分别为18.64，24.23，28.40，33.24，37.96，41.32，44.48，47.51，51.19， 计算小于各阈值条件的指示变异函数（步长为500m）,2018/10/14,123,华中农业大学 资源与环境学院,2018/10/14,124,华中农业大学 资源与环境学院,（球状模型）,2018/10/14,125,华中农业大学 资源与环境学院,2018/10/14,126,华中农业大学 资源与环境学院,思考题,1、简述空间属性插值的一般流程 2、克里格方法有哪些，各适用于什么情况？,
