第5章克里格法

1、一一克里金法概述克里金法概述二二线性克里金法线性克里金法1.1. 简单克里金简单克里金2.2. 普通克里金普通克里金3.3. 泛克里金法泛克里金法三三非线性克里金法非线性克里金法1.1. 对数正态克里金法对数正态克里金法2.2. 指示克里金法指示克里金法3.3. 析取克里金法析取克里金法四四协同克里金法协同克里金法1 1、克里金法概念及种类、克里金法概念及种类概念：概念：又称为又称为空间局部估计空间局部估计或或空间局部插值法空间局部插值法,克里金法是建立在变异函数理论克里金法是建立在变异函数理论及结构分析基础上，在有限区域内对区域化变量的取值进行线性无偏最优估及结构分析基础上，在有限区域内对区

2、域化变量的取值进行线性无偏最优估计的一种方法。计的一种方法。 主要类型：主要类型： 简单克里金法简单克里金法 普通克里金法普通克里金法 Ordinary Kriging 泛克里金法泛克里金法 Universal Kriging 对数正态克里金法对数正态克里金法 Logistic Normal Kriging 指示克里金法指示克里金法 Indicator Kriging 概率克里金概率克里金 Probability Kriging 析取克里金法析取克里金法 Disjuctive Kriging 协同克里金法协同克里金法 Co-KrigingniiivxZxZ1*)()( 设设x为研究区域内任一点

3、为研究区域内任一点0)()(*xZxZEvv待估点的估计待估点的估计值值克里金估计量克里金估计量权重系数权重系数待估点待估点影响范围内影响范围内的的有效样本值有效样本值（1）无偏估计）无偏估计（2）最优估计）最优估计的求解imin)()()()(2*xZxZExZxZVarvvvv显然，估计的好坏显然，估计的好坏取决于权重系数取决于权重系数i（1）数据检查（2）模型拟合（3）模型诊断（4）模型比较 当区域化变量Z(x)的EZ(x)=m已知，则称为简单克里金法 若Z(x)的EZ(x)未知，则称为普通克里金法设区域化变量设区域化变量Z(x)满足二阶平稳假设，其数学期望为常数满足二阶平稳假设，其数学

4、期望为常数m，协方差函数，协方差函数C(h)和和变异函数变异函数 ( (h h) )存在且平稳。存在且平稳。 现现要估计中心点在要估计中心点在x0 的待估块段的待估块段V 的均值的均值ZV(x)， ZV(x)表达式表达式为为 由于由于 EZ(x)=m已知已知令令 Y(x)=Z(x)-m则则 EY(x)=EZ(x)-m= EZ(x)-m=0待估块段新待估值待估块段新待估值dxxZvxZvv)(1)(设在待估块段V附近有n个样点xi(i=1,2,n),其观测值为Z(xi) (i=1,2,n)，则观测值新变量为：Y(xi)=Z(xi)-mY(V)的估计值Yv*是Y(xi) (i=1,2,n)的线性组

5、合，则目标：找出一组权重系数目标：找出一组权重系数 ，使得，使得Yv*成为成为Y(V) 的线性、无偏、最优估计量的线性、无偏、最优估计量), 2 , 1(nii则估计Z(V)的问题转化为估计Y(V)的问题在满足以下两个条件时，Yv*是Y(V)的线性、无偏、最优估计量。（1）无偏性）无偏性 由于由于 所以所以 则则 Yv* 不需要任何条件即是不需要任何条件即是Y(V)的无偏估计量。（2）最优性）最优性 在满足无偏条件下，可推导估计方差公式为：在满足无偏条件下，可推导估计方差公式为：为使估计方差最小，需对上式求为使估计方差最小，需对上式求i的偏导数并令其为0整理得整理得简单简单克里金克里金方程组：

6、方程组：用矩阵表示为：用矩阵表示为：将将简单简单克里金克里金方程组方程组表达式带入估计方差表达式表达式带入估计方差表达式得得简单简单克里金克里金估计估计方差表达式：方差表达式：从简单从简单克里金克里金方程组的方程组的n个方程中便可求得个方程中便可求得n个权重系数个权重系数i，则，则YV(x)的简单的简单克里金克里金估计量为：估计量为：简单简单克里金克里金法的估计精度在很大程度上依赖于法的估计精度在很大程度上依赖于m值的准确度，但是通常情值的准确度，但是通常情况下很难正确估计况下很难正确估计m值，从而导致简单值，从而导致简单克里金克里金估计精度降低。估计精度降低。 简单克里金法计算示例：简单克里

7、金法计算示例：设某一区域气温数据满足二阶平稳假设，协方差函数和变异函数存在，所有采设某一区域气温数据满足二阶平稳假设，协方差函数和变异函数存在，所有采样数据的均值为样数据的均值为16.08度，并将均值作为此区域化变量的数学期望值，将所有度，并将均值作为此区域化变量的数学期望值，将所有采样数据剔除数学期望值后拟合的变异函数模型为球状模型，如下所示。采样数据剔除数学期望值后拟合的变异函数模型为球状模型，如下所示。现用简单克里金方法根据五个已知点的气温数据来估算0点处的气温值设区域化变量设区域化变量Z(x)满足二阶平稳假设，其数学期望为满足二阶平稳假设，其数学期望为m，为未知常数，协方，为未知常数，

8、协方差函数差函数C(h)和变异函数和变异函数(h)存在且平稳。现要估计中心点在存在且平稳。现要估计中心点在x0的待估块段的待估块段V的的均值，即均值，即设待估块段设待估块段V附近有附近有n个样点个样点xi(i =1,2,n),其观测值为其观测值为Z(xi) (i =1,2,n)，待，待估块段估块段V的真值是估计邻域内的真值是估计邻域内n个信息值的线性组合，即个信息值的线性组合，即现要求出权重系数现要求出权重系数i(i =1,2,n)，使，使Z*V(x)为为ZV(x)的无偏估计量，且估计方的无偏估计量，且估计方差最小。差最小。（1）无偏性条件）无偏性条件 由于由于若要满足无偏性条件，需若要满足无

9、偏性条件，需 ，则无偏性条件为：，则无偏性条件为： 即在权系数之和为即在权系数之和为1的条件下估计量是无偏的。的条件下估计量是无偏的。（2）最优性条件）最优性条件 即估计方差最小条件，在满足无偏性条件下，有如下估计方差公式即估计方差最小条件，在满足无偏性条件下，有如下估计方差公式 要求出在满足无偏性条件要求出在满足无偏性条件 下使得估计方差最小的权系数下使得估计方差最小的权系数i(i =1,2,n)， 这这是个求条件极值问题。是个求条件极值问题。根据拉格朗日乘数法原理，建立拉格朗日函数F。求出函数F对n个权系数i的偏导数，并令其为0，和无偏性条件联立建立方程组。整理得普通克里金方程组将解出的将

10、解出的i(i =1,2,n)带入估计量带入估计量公式得到普通公式得到普通克里金克里金估计量：估计量：从普通从普通克里金克里金方程组可得：方程组可得：将此式带入估计方差公式得将此式带入估计方差公式得普通克里金估计方差，记为 ：普通克里金方程组和普通克里金估计方差也可用变异函数(h)表示。在在Z(x)满足二阶平稳条件时，可采满足二阶平稳条件时，可采用协方差或变异函数表达的普通用协方差或变异函数表达的普通克克里金里金方程组及方程组及克里金克里金估计方差计算估计方差计算式进行求解计算；但在本证假设条式进行求解计算；但在本证假设条件下，则只可采用变异函数的表达件下，则只可采用变异函数的表达式进行求解计算

11、。式进行求解计算。为了书写简便和便于计算，普通为了书写简便和便于计算，普通克克里金里金方程组和普通方程组和普通克里金克里金估计方差估计方差均可用矩阵形式表示。均可用矩阵形式表示。协方差函数表达的普通克里金方程组展开得引入矩阵或普通克里金方程组用矩阵形式表达为： 或权重系数 或普通克里金估计方差用矩阵表达为： 或普通克里金计算示例：普通克里金计算示例：设某一区域气温数据满足二阶平稳假设，协方差函数和变异函数存在，拟合的变异函数模型为球状模型，如下所示。数据如下，点的空间分布如图所示。现用普通克里金方法根据已知五个点的气温数据估算0点处的气温值。普通克里金法要求区域化变量Z(x)是二阶平稳或本征的

12、，至少是准二阶平稳或准本征的。在此条件下，至少在估计邻域内有EZ(x)=m（常数）。然而实际中，许多区域化变量Z(x)在估计邻域内是非平稳的，即EZ(x)=m(x)，m(x)称为漂移，这时就不能用普通克里金方法进行估计了，而是要采用泛克里金法进行估计。所谓泛克里金法，就是在漂移的形式EZ(x)=m(x)，和非平稳随机函数Z(x)的协方差函数C(h)或变异函数(h)为已知的条件下，一种考虑到有漂移的无偏线性估计量的地统计学方法，这种方法属于线性非平稳地统计学范畴。漂移漂移：非非平稳区域化变量平稳区域化变量Z(x)的数学的数学期望，在任一点期望，在任一点x上的漂移就是该点上的漂移就是该点上区域化变

13、量上区域化变量Z(x)的数学期望的数学期望。漂移经常用邻域模型来研究。可表达为：在给定的以点x为中心的邻域内的任一点，其漂移m(x)可用如下函数表示。式中，fl(x)为一已知函数；al为未知系数m(x)通常采用多项式形式，在二维条件下，漂移可看成坐标x,y的函数。涨落：涨落：对于有漂移的区域化变量对于有漂移的区域化变量Z(x)，假设可分解为漂移和涨落两，假设可分解为漂移和涨落两部分，部分，式中，m(x) = EZ(x)为点x处的漂移，R(x)称为涨落。1）基本假设假设Z(x)的增量Z(x)Z(y)具有非平稳的数学期望m(x)m(y)和非平稳的方差函数，即假设下式存在：2）协方差函数和变异函数当

14、Z(x)=m(x)+R(x)时，Z(x)的协方差函数C(x,y)为：Z(x)的变异函数(x,y)为：设Z(x)为一非平稳区域化变量，其数学期望为m(x)，协方差函数为C(x,y)且已知，则设Z(x)的漂移m(x)可表示为如下k+1个单项式fl(x)(l=0,1,2,k)的线性组合。已知n个样品点xi(i =1,2,n)，其观测值为Z(xi) (i =1,2,n)，现要用这些样品点估计邻域内任一点x的值Z(x)，Z(x)的泛克里金估计量为：为使Z*(x)为Z(x)的无偏最优估计量，需在以下两个条件下求解权重系数i(i =1,2,n)。1）无偏性条件）无偏性条件若要满足无偏性条件，需则即对任一组系

15、数a0，a1,ak等式均成立，需成立。这k+1个子式称为无偏性条件。2）最优性条件）最优性条件在满足无偏性条件下，用Z*(x)估计Z(x)的泛克里金估计方差为：将无偏性条件带入得要求出在满足无偏性的条件下使得估计方差最小的权系数i(i =1,2,n)，需根据拉格朗日乘数法原理，建立拉格朗日函数F。求出函数F对n个权系数i的偏导数，并令其为0，和无偏性条件联立建立如下方程组。整理得估计Z (x)的泛克里金方程组：泛克里金方程组可用矩阵表示为：其中从泛克里金方程组可得以下两等式：将等式带入估计方差公式可得泛克里金方差，记为：用变异函数(h)表示如下： 设某一区域气温是非平稳的区域化变量，在南北方向

16、（空间坐标的设某一区域气温是非平稳的区域化变量，在南北方向（空间坐标的y方向）上方向）上存在线性漂移，即存在线性漂移，即 。若已知其涨落满足二阶平稳假设，并且拟合的协方差函。若已知其涨落满足二阶平稳假设，并且拟合的协方差函数模型为球状模型，如下所示。数模型为球状模型，如下所示。现用表现用表5 1所示数据，利用泛所示数据，利用泛克里金克里金法根据已知五个点的气温数据来估算法根据已知五个点的气温数据来估算0点点处的气温值。处的气温值。1、对数正态对数正态克里金克里金法法如果区域化变量经对数变换后是正态分布或近正态分布，如果区域化变量经对数变换后是正态分布或近正态分布，则则对区域化变量进对区域化变量

17、进行精确估计的地统计学方法称为对数正态克立格法。行精确估计的地统计学方法称为对数正态克立格法。设区域化变量设区域化变量Z(x)服从对数正态分布，在待估点周围有服从对数正态分布，在待估点周围有n个样点个样点xi(i =1,2,n),其观测值为其观测值为Z(xi) (i =1,2,n)，区域化变量经对数变换后新变量为：，区域化变量经对数变换后新变量为：Y(x)= lnZ(x)，Y(x)为正态分布。假定为正态分布。假定Y(x)满足二阶平稳假设，数学期望为满足二阶平稳假设，数学期望为m，协方，协方差函数差函数C(h)和变异函数和变异函数(h)存在且平稳。存在且平稳。基于对数变换后的采样点数据基于对数变

18、换后的采样点数据Y(xi) (i =1,2,n)，计算实验变异函数并进行变，计算实验变异函数并进行变异函数模型的拟合和选择，然后利用简单克立格或普通克立格估计待估点异函数模型的拟合和选择，然后利用简单克立格或普通克立格估计待估点x处处的值的值Y*(x)。由于估计值由于估计值Y(x)是对数变换后的数值，因此是对数变换后的数值，因此对估计所得对估计所得Y*(x)需进行反变换。需进行反变换。实际研究中常常会需要获取研究区内研究对象大于某一给定阈值的概率分布，实际研究中常常会需要获取研究区内研究对象大于某一给定阈值的概率分布，即要获知研究区内任一点即要获知研究区内任一点x处随机变量处随机变量Z(x)的

19、概率分布。的概率分布。还会碰到采样数据中存在特异值的问题。还会碰到采样数据中存在特异值的问题。（特异值是指那些比全部数值的均值特异值是指那些比全部数值的均值或中位数高的多的数值，其既非分析误差所致，也非采样方法等人为误差引起或中位数高的多的数值，其既非分析误差所致，也非采样方法等人为误差引起，而是实际存在于所研究的总体之中，而是实际存在于所研究的总体之中）。指示克立格法就是为解决上述问题而发展起来的一种非参数地统计学方法。指示克立格法就是为解决上述问题而发展起来的一种非参数地统计学方法。指示克立格法不必去掉重要而实际存在的高值数据的条件下处理各种不同现象指示克立格法不必去掉重要而实际存在的高值

20、数据的条件下处理各种不同现象，并能够给出某点，并能够给出某点x处随机变量处随机变量Z(x)的概率分布。的概率分布。设一区域化变量Z(x)，对于任意给定的阈值z，引入指示函数I(x, z)，表达式如下：指示克立格法步骤如下：（1）确定一阈值，根据指示函数将原数据转换为0或1；（2）利用转换的数据计算指示变异函数，并进行拟合；（3）建立指示克立格方程组，计算待估点值。若把指示函数看做一普通区域化变量，也可直接由简单或普通克立格方法来计算待估点的值。若选择多个阈值则需重复以上步骤。析取克立格析取克立格法法：假设假设已知任意区域化变量（已知任意区域化变量（Z , Z ）及（）及（Z0, Z ）二维概）

21、二维概率分布条件下，对待估点的值或待估点值超过给定阈值的概率进行估率分布条件下，对待估点的值或待估点值超过给定阈值的概率进行估计的一种非线性地统计方法计的一种非线性地统计方法。估值估值步骤步骤：设设区域化变量区域化变量Z(x)在待估点在待估点x0周围有周围有n个样点个样点xi(i =1,2,n),其其观测观测 值值为为Z(xi) (i =1,2,n)， 将原始数据转换为标准正态数据将原始数据转换为标准正态数据 对每个新变量对每个新变量Y(xi)(i=1,2,n) 计算埃尔米特多项式的值计算埃尔米特多项式的值。 计算埃尔米特多项式系数，用埃尔米特多项式来拟合正态变形函数计算埃尔米特多项式系数，用

22、埃尔米特多项式来拟合正态变形函数。 计算待估点析取克立格值计算待估点析取克立格值1、协同区域化变量理论协同区域化变量理论协同克立格法协同克立格法：是多元地统计学研是多元地统计学研究的基本方法，建立在协同区域化究的基本方法，建立在协同区域化变量理论基础之上，利用多个区域变量理论基础之上，利用多个区域化变量之间的互相关性，通过建立化变量之间的互相关性，通过建立交叉协方差函数和交叉变异函数模交叉协方差函数和交叉变异函数模型型，用易于观测和控制的变量对不用易于观测和控制的变量对不易观测的变量进行局部估计。易观测的变量进行局部估计。协同区域化：协同区域化：在统计意义及空间在统计意义及空间位置上均具有某种

23、程度相关性，并位置上均具有某种程度相关性，并且定义于同一空间域中的区域化变且定义于同一空间域中的区域化变量。量。协同区域化变量可用一组K个相关的区域化变量 表示。观测前它是K维区域化变量的向量，即一个随机场，观测后，协同区域化变量是一个空间点函数，可以把 看成是上述K维向量的一个实现。满足二阶平稳假设的协同区域化变量应满足二阶平稳假设的协同区域化变量应满足：满足：（1）每一个协同区域化变量的数学期望存在且平稳：（2）交叉协方差函数存在，且平稳：满足内蕴假设的协同区域化变量应满足内蕴假设的协同区域化变量应满足：满足：（1）每一个协同区域化变量增量的数学期望为0：（2）对于协同区域化变量，交叉变异

24、函数存在且平稳。即（1）交叉协方差函数和交叉变异函数性质交叉协方差函数和交叉变异函数性质1)当当k=k 时时，交叉协方差函数转化为协方差函数，交叉变异函数转化为变异函，交叉协方差函数转化为协方差函数，交叉变异函数转化为变异函数。数。即即2)交叉交叉变异函数变异函数性质性质 交叉变异函数关于交叉变异函数关于k和和k 对称，对称，即即 交叉变异函数关于交叉变异函数关于h和和-h对称，即对称，即 在普通克立格法中变异函数总是大于等于在普通克立格法中变异函数总是大于等于0，但交叉变异函数可以有负值。，但交叉变异函数可以有负值。3)交叉交叉协方差协方差函数函数性质性质交叉协方差函数交叉协方差函数关于关于

25、h和和-h不对称不对称，即即 ，但，但 当当h0 时时， k和和k 顺序顺序不能随意颠倒，不能随意颠倒，即即当当h=0 时时，交叉，交叉协方差转化为直接协方差协方差转化为直接协方差。4）交叉协方差函数和交叉变异函数具有以下关系：5）同一点两个变量点对点协同区域化变量的相关系数为：（2）交叉协方差函数和交叉变异函数计算公式交叉协方差函数和交叉变异函数计算公式设在点x和x+h处，分别测得两个区域化变量的观测值 Zk(x)、Zk(x) 、 Zk(x+h)、Zk(x+h) ，则交叉协方差函数计算公式为:交叉变异函数计算公式为：（3）交叉交叉协方差函数和交叉变异函数计算协方差函数和交叉变异函数计算示例示

26、例采用表5 1、图5 1所示的气温和海拔高度数据，以h=0,h1,4为例，交叉协方差和交叉变异计算过程如下：（1）协同克立协同克立金金估计量估计量1)无偏性条件2)最优性条件对对F求偏导数并令其为零，得协同克立格线性方程组：求偏导数并令其为零，得协同克立格线性方程组：根据协同克立格方程组，协同克立格方差为：根据协同克立格方程组，协同克立格方差为：若有多个变量，则求解若有多个变量，则求解 的协同克立格方程组为：的协同克立格方程组为：协同克立格方差为：协同克立格方差为：要使协同克立格方程组具有唯一解的条件是：要使协同克立格方程组具有唯一解的条件是：设某一区域有两个协同区域化变量：气温设某一区域有两

27、个协同区域化变量：气温u、海拔高程、海拔高程v，均满足二阶平稳假，均满足二阶平稳假设和内蕴假设，其中气温是所要估计的主变量。现在估计邻域内共有设和内蕴假设，其中气温是所要估计的主变量。现在估计邻域内共有5个信个信息样品，如图息样品，如图5 1所示。假定所示。假定1、2、3号点上有气温值，而号点上有气温值，而5个点上均有次要个点上均有次要变量海拔高程值，数据如表变量海拔高程值，数据如表5 6所示。现拟利用协同克立格法估计所示。现拟利用协同克立格法估计0号点的号点的气温值。气温值。估值过程中，气温协方差估值过程中，气温协方差C(ui,uj)根据简单克立根据简单克立金金法计算示例的球状模型法计算示例的球状模型计算，海拔高程协方差计算，海拔高程协方差C(vi,vj)根据式根据