第5章克里格法.ppt

第五章 克里金法. .提    纲1 1克里金法概述克里金法概述2 2线性克里金法线性克里金法1.1.简单克里金简单克里金2.2.普通克里金普通克里金3.3.泛克里金法泛克里金法3 3非线性克里金法非线性克里金法1.1.对数正态克里金法对数正态克里金法2.2.指示克里金法指示克里金法3.3.析取克里金法析取克里金法4 4协同克里金法协同克里金法. .一、克里金法概述一、克里金法概述1 1、克里金法概念及种类、克里金法概念及种类概念：概念：又称为又称为空间局部估计空间局部估计或或空间局部插值法空间局部插值法,克里金法是建立在变异函数理论克里金法是建立在变异函数理论及结构分析基础上，在有限区域内对区域化变量的取值进行线性无偏最优估及结构分析基础上，在有限区域内对区域化变量的取值进行线性无偏最优估计的一种方法计的一种方法 主要类型：主要类型： 简单克里金法简单克里金法 普通克里金法普通克里金法 Ordinary Kriging 泛克里金法泛克里金法 Universal Kriging 对数正态克里金法对数正态克里金法 Logistic Normal Kriging 指示克里金法指示克里金法 Indicator Kriging 概率克里金概率克里金 Probability Kriging 析取克里金法析取克里金法 Disjuctive Kriging 协同克里金法协同克里金法 Co-Kriging. .2、克里金估计量•设设x为研究区域内任一点为研究区域内任一点待估点的估计待估点的估计值值克里金估计量克里金估计量权重系数权重系数待估点待估点影响范围内影响范围内的的有效样本值有效样本值（（1）无偏估计）无偏估计（（2）最优估计）最优估计显然，估计的好坏显然，估计的好坏取决于权重系数取决于权重系数λi. .3、克里金法估值过程（1）数据检查（2）模型拟合（3）模型诊断（4）模型比较. .•当区域化变量Z(x)的E[Z(x)]=m已知，则称为简单克里金法•若Z(x)的E[Z(x)]未知，则称为普通克里金法二、线性克里金法. .1、简单克里金法•设区域化变量设区域化变量Z(x)满足二阶平稳假设，其数学期望为常数满足二阶平稳假设，其数学期望为常数m，协方差函数，协方差函数C(h)和变异函数和变异函数γ ( (h h) )存在且平稳。

存在且平稳 现要估计中心点在现要估计中心点在x0 的待估块段的待估块段V 的均值的均值ZV(x)，， ZV(x)表达式为表达式为 •由于由于 E[Z(x)]=m已知已知•令令 Y(x)=Z(x)-m•则则 E[Y(x)]=E[Z(x)-m]= E[Z(x)]-m=0•待估块段新待估值待估块段新待估值. .1、简单克里金法•设在待估块段V附近有n个样点xi(i=1,2,…n),其观测值为Z(xi) (i=1,2,…n)，则观测值新变量为：Y(xi)=Z(xi)-m•Y(V)的估计值Yv*是Y(xi) (i=1,2,…n)的线性组合，则目标：找出一组权重系数目标：找出一组权重系数 ，使得，使得Yv*成为成为Y(V) 的线性、无偏、最优估计量的线性、无偏、最优估计量则估计Z(V)的问题转化为估计Y(V)的问题. .1、简单克里金法•在满足以下两个条件时，Yv*是Y(V)的线性、无偏、最优估计量。

1）无偏性）无偏性 由于由于 所以所以 则则 Yv* 不需要任何条件即是不需要任何条件即是Y(V)的无偏估计量2）最优性）最优性 在满足无偏条件下，可推导估计方差公式为：在满足无偏条件下，可推导估计方差公式为：. .1、简单克里金法•为使估计方差最小，需对上式求为使估计方差最小，需对上式求λi的偏导数并令其为0•整理得简单整理得简单克里金克里金方程组：方程组：•用矩阵表示为：用矩阵表示为：•将简单将简单克里金克里金方程组表达式带入估计方差表达式得方程组表达式带入估计方差表达式得简单简单克里金克里金估计方差表达式：估计方差表达式：. .1、简单克里金法从简单从简单克里金克里金方程组的方程组的n个方程中便可求得个方程中便可求得n个权重系数个权重系数λi，则，则YV(x)的简单的简单克里金克里金估计量为：估计量为：简单简单克里金克里金法的估计精度在很大程度上依赖于法的估计精度在很大程度上依赖于m值的准确度，但是通常情值的准确度，但是通常情况下很难正确估计况下很难正确估计m值，从而导致简单值，从而导致简单克里金克里金估计精度降低。

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

差最小 .2、普通克里金法（（1）无偏性条件）无偏性条件 由于由于若要满足无偏性条件，需若要满足无偏性条件，需 ，则无偏性条件为：，则无偏性条件为： 即在权系数之和为即在权系数之和为1的条件下估计量是无偏的的条件下估计量是无偏的2）最优性条件）最优性条件 即估计方差最小条件，在满足无偏性条件下，有如下估计方差公式即估计方差最小条件，在满足无偏性条件下，有如下估计方差公式 要求出在满足无偏性条件要求出在满足无偏性条件 下使得估计方差最小的权系数下使得估计方差最小的权系数λi(i =1,2,…,n)，， 这这是个求条件极值问题是个求条件极值问题 .2、普通克里金法•根据拉格朗日乘数法原理，建立拉格朗日函数F•求出函数F对n个权系数λi的偏导数，并令其为0，和无偏性条件联立建立方程组•整理得普通克里金方程组. .2、普通克里金法•将解出的将解出的λi(i =1,2,…,n)带入估计量带入估计量公式得到普通公式得到普通克里金克里金估计量：估计量：•从普通从普通克里金克里金方程组可得：方程组可得：•将此式带入估计方差公式得将此式带入估计方差公式得普通克里金估计方差，记为     ：：•普通克里金方程组和普通克里金估计方差也可用变异函数γ(h)表示。

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

在此条件下，至少在估计邻域内有E[Z(x)]=m（常数）然而实际中，许多区域化变量Z(x)在估计邻域内是非平稳的，即E[Z(x)]=m(x)，m(x)称为漂移，这时就不能用普通克里金方法进行估计了，而是要采用泛克里金法进行估计•所谓泛克里金法，就是在漂移的形式E[Z(x)]=m(x)，和非平稳随机函数Z(x)的协方差函数C(h)或变异函数γ(h)为已知的条件下，一种考虑到有漂移的无偏线性估计量的地统计学方法，这种方法属于线性非平稳地统计学范畴 .（1）漂移和涨落•漂移漂移：：非平稳区域化变量非平稳区域化变量Z(x)的数学的数学期望，在任一点期望，在任一点x上的漂移就是该点上的漂移就是该点上区域化变量上区域化变量Z(x)的数学期望的数学期望•漂移经常用邻域模型来研究可表达为：在给定的以点x为中心的邻域内的任一点，其漂移m(x)可用如下函数表示式中，fl(x)为一已知函数；al为未知系数m(x)通常采用多项式形式，在二维条件下，漂移可看成坐标x,y的函数•涨落：涨落：对于有漂移的区域化变量对于有漂移的区域化变量Z(x)，假设可分解为漂移和涨落两，假设可分解为漂移和涨落两部分，部分，•式中，m(x) = E[Z(x)]为点x处的漂移，R(x)称为涨落。

.（（2 2））非平稳区域化变量的协方差函数和变异函数非平稳区域化变量的协方差函数和变异函数•1）基本假设•假设Z(x)的增量[Z(x)－Z(y)]具有非平稳的数学期望[m(x)－m(y)]和非平稳的方差函数，即假设下式存在：•2）协方差函数和变异函数•当Z(x)=m(x)+R(x)时，Z(x)的协方差函数C(x,y)为：•Z(x)的变异函数γ(x,y)为：. .（3）Z(x)的泛克里金法估计•设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) .（3）Z(x)的泛克里金法估计•1）无偏性条件）无偏性条件•若要满足无偏性条件，需•则•即对任一组系数a0，a1,…,ak等式均成立，需•成立。

这k+1个子式称为无偏性条件 .（3）Z(x)的泛克里金法估计•2）最优性条件）最优性条件•在满足无偏性条件下，用Z*(x)估计Z(x)的泛克里金估计方差为：•将无偏性条件带入得•要求出在满足无偏性的条件下使得估计方差最小的权系数λi(i =1,2,…,n)，需根据拉格朗日乘数法原理，建立拉格朗日函数F .（3）Z(x)的泛克里金法估计•求出函数F对n个权系数λi的偏导数，并令其为0，和无偏性条件联立建立如下方程组•整理得估计Z (x)的泛克里金方程组：•泛克里金方程组可用矩阵表示为：•其中. .（3）Z(x)的泛克里金法估计•从泛克里金方程组可得以下两等式：•将等式带入估计方差公式可得泛克里金方差，记为：•用变异函数γ(h)表示如下： . .（4）泛克里金法计算示例•设某一区域气温是非平稳的区域化变量，在南北方向（空间坐标的设某一区域气温是非平稳的区域化变量，在南北方向（空间坐标的y方向）上方向）上存性漂移，即存性漂移，即 若已知其涨落满足二阶平稳假设，并且拟合的协方差函数若已知其涨落满足二阶平稳假设，并且拟合的协方差函数模型为球状模型，如下所示模型为球状模型，如下所示•现用表现用表5‑ ‑1所示数据，利用泛所示数据，利用泛克里金克里金法根据已知五个点的气温数据来估算法根据已知五个点的气温数据来估算0点处点处的气温值。

的气温值 .三、非线性克里金法•1、对数正态对数正态克里金克里金法法•如果区域化变量经对数变换后是正态分布或近正态分布，如果区域化变量经对数变换后是正态分布或近正态分布，则则对区域化变量进对区域化变量进行精确估计的地统计学方法称为对数正态克立格法行精确估计的地统计学方法称为对数正态克立格法•设区域化变量设区域化变量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)存在且平稳存在且平稳•基于对数变换后的采样点数据基于对数变换后的采样点数据Y(xi) (i =1,2,…,n)，计算实验变异函数并进行变，计算实验变异函数并进行变异函数模型的拟合和选择，然后利用简单克立格或普通克立格估计待估点异函数模型的拟合和选择，然后利用简单克立格或普通克立格估计待估点x处处的值的值Y*(x)。

•由于估计值由于估计值Y(x)是对数变换后的数值，因此是对数变换后的数值，因此对估计所得对估计所得Y*(x)需进行反变换需进行反变换 .2、指示克里金法•实际研究中常常会需要获取研究区内研究对象大于某一给定阈值的概率分布，实际研究中常常会需要获取研究区内研究对象大于某一给定阈值的概率分布，即要获知研究区内任一点即要获知研究区内任一点x处随机变量处随机变量Z(x)的概率分布的概率分布•还会碰到采样数据中存在特异值的问题还会碰到采样数据中存在特异值的问题特异值是指那些比全部数值的均值特异值是指那些比全部数值的均值或中位数高的多的数值，其既非分析误差所致，也非采样方法等人为误差引起，或中位数高的多的数值，其既非分析误差所致，也非采样方法等人为误差引起，而是实际存在于所研究的总体之中而是实际存在于所研究的总体之中））•指示克立格法就是为解决上述问题而发展起来的一种非参数地统计学方法指示克立格法就是为解决上述问题而发展起来的一种非参数地统计学方法•指示克立格法不必去掉重要而实际存在的高值数据的条件下处理各种不同现象，指示克立格法不必去掉重要而实际存在的高值数据的条件下处理各种不同现象，并能够给出某点并能够给出某点x处随机变量处随机变量Z(x)的概率分布。

的概率分布 .2、指示克里金法•设一区域化变量Z(x)，对于任意给定的阈值z，引入指示函数I(x, z)，表达式如下：•指示克立格法步骤如下：•（1）确定一阈值，根据指示函数将原数据转换为0或1；•（2）利用转换的数据计算指示变异函数，并进行拟合；•（3）建立指示克立格方程组，计算待估点值若把指示函数看做一普通区域化变量，也可直接由简单或普通克立格方法来计算待估点的值•若选择多个阈值则需重复以上步骤 .3、析取克里金法•析取克立格法析取克立格法：：假设已知任意区域化变量（假设已知任意区域化变量（Z , Z ）及（）及（Z0, Z ）二维概）二维概率分布条件下，对待估点的值或待估点值超过给定阈值的概率进行估计率分布条件下，对待估点的值或待估点值超过给定阈值的概率进行估计的一种非线性地统计方法的一种非线性地统计方法•估值步骤估值步骤：：设区域化变量设区域化变量Z(x)在待估点在待估点x0周围有周围有n个样点个样点xi(i =1,2,…,n),其观测其观测 值为值为Z(xi) (i =1,2,…,n)，，•将原始数据转换为标准正态数据将原始数据转换为标准正态数据•对每个新变量对每个新变量Y(xi)(i=1,2,…,n) 计算埃尔米特多项式的值。

计算埃尔米特多项式的值•计算埃尔米特多项式系数，用埃尔米特多项式来拟合正态变形函数计算埃尔米特多项式系数，用埃尔米特多项式来拟合正态变形函数•计算待估点析取克立格值计算待估点析取克立格值. .四、协同克里金法•1、协同区域化变量理论协同区域化变量理论•协同克立格法协同克立格法：：是多元地统计学研是多元地统计学研究的基本方法，建立在协同区域化究的基本方法，建立在协同区域化变量理论基础之上，利用多个区域变量理论基础之上，利用多个区域化变量之间的互相关性，通过建立化变量之间的互相关性，通过建立交叉协方差函数和交叉变异函数模交叉协方差函数和交叉变异函数模型型，，用易于观测和控制的变量对不用易于观测和控制的变量对不易观测的变量进行局部估计易观测的变量进行局部估计•协同区域化：协同区域化：在统计意义及空间在统计意义及空间位置上均具有某种程度相关性，并位置上均具有某种程度相关性，并且定义于同一空间域中的区域化变且定义于同一空间域中的区域化变量•协同区域化变量可用一组K个相关的区域化变量 表示观测前它是K维区域化变量的向量，即一个随机场，观测后，协同区域化变量是一个空间点函数，可以把 看成是上述K维向量的一个实现。

.1、协同区域化变量理论•满足二阶平稳假设的协同区域化变量应满足二阶平稳假设的协同区域化变量应满足：满足：•（1）每一个协同区域化变量的数学期望存在且平稳：•（2）交叉协方差函数存在，且平稳：•满足内蕴假设的协同区域化变量应满足内蕴假设的协同区域化变量应满足：满足：•（1）每一个协同区域化变量增量的数学期望为0：•（2）对于协同区域化变量，交叉变异函数存在且平稳即. .2、交叉协方差函数和交叉变异函数•（（1））交叉协方差函数和交叉变异函数性质交叉协方差函数和交叉变异函数性质•1)当当k=kˊ 时，交叉协方差函数转化为协方差函数，交叉变异函数转化为变异函时，交叉协方差函数转化为协方差函数，交叉变异函数转化为变异函数即•2)交叉变异函数性质交叉变异函数性质–交叉变异函数关于交叉变异函数关于k和和k ˊ对称，对称，即即–交叉变异函数关于交叉变异函数关于h和和-h对称，即对称，即–在普通克立格法中变异函数总是大于等于在普通克立格法中变异函数总是大于等于0，但交叉变异函数可以有负值但交叉变异函数可以有负值•3)交叉交叉协方差协方差函数函数性质性质•交叉协方差函数交叉协方差函数关于关于h和和-h不对称不对称，，即即 ，但，但 •当当h≠0 时时，， k和和k ˊ顺序顺序不能随意颠倒，不能随意颠倒，即即•当当h=0 时时，交叉，交叉协方差转化为直接协方差协方差转化为直接协方差。

.2、交叉协方差函数和交叉变异函数•4）交叉协方差函数和交叉变异函数具有以下关系：•5）同一点两个变量点对点协同区域化变量的相关系数为：. .2、交叉协方差函数和交叉变异函数•（2）交叉协方差函数和交叉变异函数计算公式交叉协方差函数和交叉变异函数计算公式•设在点x和x+h处，分别测得两个区域化变量的观测值 Zk(x)、Zkˊ(x) 、 Zk(x+h)、Zkˊ(x+h) ，则交叉协方差函数计算公式为:•交叉变异函数计算公式为：. .2、交叉协方差函数和交叉变异函数•（（3））交叉协方差函数和交叉变异函数计算示例交叉协方差函数和交叉变异函数计算示例•采用表5‑1、图5‑1所示的气温和海拔高度数据，以h=0,h1,4为例，交叉协方差和交叉变异计算过程如下：. .3、协同克立金法估值•（1）协同克立协同克立金金估计量估计量. .（（2））协同克立协同克立金法方程组金法方程组. .（（2））协同克立协同克立金法方程组金法方程组•1)无偏性条件•2)最优性条件 . .（（2））协同克立协同克立金法方程组金法方程组•对对F求偏导数并令其为零，得协同克立格线性方程组：求偏导数并令其为零，得协同克立格线性方程组：. .（（2））协同克立协同克立金法方程组金法方程组•根据协同克立格方程组，协同克立格方差为：根据协同克立格方程组，协同克立格方差为：•若有多个变量，则求解若有多个变量，则求解 的协同克立格方程组为：的协同克立格方程组为：•协同克立格方差为：协同克立格方差为：. .（（2））协同克立协同克立金法方程组金法方程组•要使协同克立格方程组具有唯一解的条件是：要使协同克立格方程组具有唯一解的条件是：. .（（3）协同克里金法使用条件）协同克里金法使用条件. .（（4）协同克里金法计算示例）协同克里金法计算示例•设某一区域有两个协同区域化变量：气温设某一区域有两个协同区域化变量：气温u、海拔高程、海拔高程v，均满足二阶平稳假，均满足二阶平稳假设和内蕴假设，其中气温是所要估计的主变量。

现在估计邻域内共有设和内蕴假设，其中气温是所要估计的主变量现在估计邻域内共有5个信息个信息样品，如图样品，如图5‑ ‑1所示假定所示假定1、、2、、3号点上有气温值，而号点上有气温值，而5个点上均有次要变量个点上均有次要变量海拔高程值，数据如表海拔高程值，数据如表5‑ ‑6所示现拟利用协同克立格法估计所示现拟利用协同克立格法估计0号点的气温值号点的气温值 .（（4）协同克里金法计算示例）协同克里金法计算示例•估值过程中，气温协方差估值过程中，气温协方差C(ui,uj)根据简单克立根据简单克立金金法计算示例的球状模型法计算示例的球状模型计算，海拔高程协方差计算，海拔高程协方差C(vi,vj)根据式根据式(5.74)计算，气温海拔交叉协方差计算，气温海拔交叉协方差C(ui,vj)根据式根据式(5.75)计算，海拔气温交叉协方差计算，海拔气温交叉协方差C(vi,uj)根据式根据式(5.76)计算•简单简单克立克立金金法计算示例的球状模型法计算示例的球状模型：：. .复习思考题1.请简述克立格法的概念及其种类请简述克立格法的概念及其种类2.请简述克立格估值的过程请简述克立格估值的过程3.请比较反距离加权法和普通克立格法的异同。

请比较反距离加权法和普通克立格法的异同4.将简单克立格方程组用变异函数进行表达，并写出相应的矩阵表达式，然后利用将简单克立格方程组用变异函数进行表达，并写出相应的矩阵表达式，然后利用变异函数表达的矩阵来计算简单克立格法的计算示例变异函数表达的矩阵来计算简单克立格法的计算示例5.用协方差函数表达的矩阵来计算普通克立格法的计算示例用协方差函数表达的矩阵来计算普通克立格法的计算示例6.若一区域化变量的拟合模型为纯块金效应模型，请推导说明用普通克立格法进行若一区域化变量的拟合模型为纯块金效应模型，请推导说明用普通克立格法进行估值时，参与插值的采样点的权重特点估值时，参与插值的采样点的权重特点7.若普通克立格法的计算示例中，若普通克立格法的计算示例中，1号点的坐标为号点的坐标为(1099904.96, -337837.61)，即与，即与0号点的位置重合，请计算这时号点的位置重合，请计算这时0号点的估计值，并说明克立格插值法的特点号点的估计值，并说明克立格插值法的特点 .复习思考题8、、利用普通克立格法和如下球状模型对下表利用普通克立格法和如下球状模型对下表0号点进行估计号点进行估计9、、根据用变异函数表达的泛克立格方程组写出对应的矩阵表达式，并计算泛克根据用变异函数表达的泛克立格方程组写出对应的矩阵表达式，并计算泛克立格法示例。

立格法示例10、、在在ArcGIS软件中实现协同克立格法插值，并对结果进行验证软件中实现协同克立格法插值，并对结果进行验证。