4地统计学.ppt
地统计学（第五讲）,常见半方差模型 应用实例 克立格插值方法,1、协方差函数,协方差函数的概念 区域化随机变量之间的差异，可以用空间协方差来表示。区域化变量 在空间点x和x+h处的两个随机变量和的二阶混合中心矩定义为Z(x)的自协方差函数，即,协方差函数的计算公式为：式中：h为两样本点空间分隔距离或距离滞后，为 在空间位置 处的实测值，是 在 处距离偏离h的实测值i=1，2，是分隔距离为h时的样本点对（Paris）总数，和 分别为 和 的样本平均数。,协方差函数的计算公式,若=m（常数），则上式可以改写为：式中：m为样本平均数，可由一般算术平均数公式求得，即：,变异函数,变异函数的概念 变异函数（Variograms），又称变差函数、变异矩，是地统计分析所特有的基本工具。在一维条件下变异函数定义为，当空间点x在一维x轴上变化时，区域化变量Z(x)在点x和x+h处的值Z(x)与Z(x+h)差的方差的一半为区域化变量Z(x)在x轴方向上的变异函数，记为(h)，即,在二阶平稳假设条件下，对任意的h有因此，公式可以改写为 从上式可知，变异函数依赖于两个自变量x和h，当变异函数 仅仅依赖于距离h而与位置x无关时，可改写成，即：,变异函数的性质设Z(x)是区域化变量，在满足二阶平稳假设条件下，变异函数式具有如下性质：(1)=0，即在h=0处，变异函数为0；(2)=，即 关于直线h=0是对称的，它是一个偶函数；(3)0，即 只能大于或等于0；,(4)|h|时，c(0)，或=c(0)，即当空间距离增大时，变异函数接近先验方差,变异函数的计算公式,设 是系统某属性Z在空间位置x处的值，为一区域化随机变量，并满足二阶平稳假设，h为两样本点空间分隔距离，和 分别是区域化变量 在空间位置 和 处的实测值i=1,2,N(h)，那么，变异函数 的离散计算公式为,这样对不同的空间分隔距离h，计算出相应的 和 值。如果分别以h为横坐标，或 为纵坐标，画出协方差函数和变异函数曲线图，就可以直接展示区域化变量Z(x)的空间变异特点。可见，变异函数能同时描述区域化变量的随机性和结构性，从而在数学上对区域化变量进行严格分析，是空间变异规律分析和空间结构分析的有效工具。,实例1，假设某地区降水量Z(x)（单位：mm）是二维区域化随机变量，满足二阶平稳假设，其观测值的空间正方形网格数据如图4.2.1所示（点与点之间的距离为h=1km）。试计算其南北方向及西北和东南方向的变异函数。,图4.2.1 空间正方形网格数据（点间距h=1km）,从图4.2.1可以看出，空间上有些点，由于某种原因没有采集到。如果没有缺失值，可直接对正方形网格数据结构计算变异函数；在有缺失值的情况下，也可以计算变异函数。只要“跳过”缺失点位置即可（见图4.2.2）。,首先计算南北方向上的变异函数值，由变异函数的计算公式可得：=385/72=5.35,图4.2.2 缺失值情况下样本数对的组成和计算过程，为缺失值,同样计算出最后，得到南北方向和西北东南上的变异函数计算结果见下表。同样可以计算东西方向上的变异函数。,变异函数的参数,变异函数有四个非常重要的参数，即基台值（Sill）、变程（Range）或称空间依赖范围（Range of Spatial Dependence）、块金值（Nugget）或称区域不连续性值（Localized Discontinuity）和分维数（Fractal Dimension）。前3个参数可以直接从变异函数图中得到。它们决定变异函数的形状与结构。变异函数的形状反映自然现象空间分布结构或空间相关的类型，同时还能给出这种空间相关的范围。,当变异函数随着间隔距离h的增大，从非零值达到一个相对稳定的常数时，该常数称为基台值C0+C，当间隔距离h=0时，(0)=C0，该值称为块金值或块金方差（Nugget Variance）。基台值是系统或系统属性中最大的变异，变异函数达到基台值时的间隔距离a称为变程。变程表示在ha以后，区域化变量Z(x)空间相关性消失。块金值表示区域化变量在小于抽样尺度时非连续变异，由区域化变量的属性或测量误差决定。,C0,C0+C2,(h),上述三个参数可从变异函数曲线图直接得到，或通过估计曲线回归参数得到。第4个参数，即分维数用于表示变异函数的特性，由变异函数 和间隔距离h之间的关系确定：分维数D为双对数直线回归方程中的斜率，它是一个无量纲数。分维数D的大小，表示变异函数曲线的曲率，可以作为随机变异的量度。但该随机分维数D与形状分维数有本质的不同。,变异函数的理论模型,实际上理论变异函数模型 往往是未知的，需要从有效的空间取样数据中去估计。对各种不同的h值可计算出一系列的 值，因此要有一个理论模型去拟合这一系列的 值,变异函数的理论模型,地统计学将变异函数理论模型分为三大类：第一类是有基台值模型，包括球状模型、指数模型、高斯模型、线性有基台值模型和纯块金效应模型；第二类是无基台值模型，包括幂函数模型、线性无基台值模型、抛物线模型；第三类是孔穴效应模型。下面有代表性地介绍几种常见的变异函数理论模型。,（1）纯块金效应模型。其一般公式为：式中：c00，为先验方差。该模型相当于区域化变量为随机分布，样本点间的协方差函数对于所有距离h均等于0，变量的空间相关不存在。,（2）球状模型。其一般公式为：式中：c0为块金（效应）常数，c为拱高，c0+c为基台值，a为变程。当c0=0，c=1时，称为标准球状模型。球状模型是地统计分析中应用最广泛的理论模型，许多区域化变量的理论模型都可以用该模型去拟合。,（3）指数模型。其一般公式为：式中：c0和c意义与前相同，但a不是变程。当h=3a时，即，从而指数模型的变程 约为。当c0=0，c=1时，称为标准指数模型。,（4）高斯模型。其一般公式为：式中：c0和c意义与前相同，a也不是变程。当 时，即，因此高斯模型的变程 约为。当 时，称为标准高斯函数模型。,（5）幂函数模型。其一般公式为：式中：为幂指数。当变化时，这种模型可以反映在原点附近的各种性状。但是必须小于2，若，则函数 就不再是一个条件非负定函数了，也就是说它已经不能成为变异函数了。,（6）对数模型。其一般公式为：显然，当，这与变异函数的性质不符。因此，对数模型不能描述点支撑上的区域化变量的结构。,（7）线性有基台值模型。其一般公式为：式中该模型的变程为a，基台值为。（8）线性无基台值模型。其一般公式为 从式中可以看出，该模型没有基台值，也没有变程。,例如，某地区降水量是一个区域化变量，其变异函数 的实测值及距离h的关系见下表，下面我们试用回归分析方法建立其球状变异函数模型。,实测值(h),距离h,2.1,0.6,9.2,4.9,4.3,1.1,10.3,5.1,5.7,2.2,10.5,6.2,6.5,2.5,10.9,7.5,7.8,3.1,11.2,9.5,一、变异函数理论模型的最优拟合,地统计学中变异函数的理论模型建立与普通统计学中的回归模型的建立相似。为了使得理论模型能够最充分地描述所研究的某一区域化变量的变化规律，需要对理论模型进行检验一般根据变异函数的计算值，选择合适的理论模型来拟合一条最优的理论变异函数曲线问题（一般称为最优拟合）,一、变异函数理论模型的最优拟合,最优拟合的过程实质上是配合最优模型的过程拟合过程一般包括三个步骤第一步：确定曲线类型第二步：参数最优估计第三步：最优曲线的确定,第一步：确定曲线类型,所谓曲线指常见的变异函数理论模型根据h与变化函数之间的散点图确定曲线类型,第一步：确定曲线类型,例如，某地区降水量是一个区域化变量，其变异函数 的实测值及距离h的关系见下表，下面我们试用回归分析方法建立其球状变异函数模型。,实测值(h),距离h,2.1,0.6,9.2,4.9,4.3,1.1,10.3,5.1,5.7,2.2,10.5,6.2,6.5,2.5,10.9,7.5,7.8,3.1,11.2,9.5,从上面的介绍和讨论，我们知道，球状变异函数的一般形式为：当 时，有：,第二步：参数最优估计,如果记，则可以得到线性模型：根据表中的数据，对上式进行最小二乘拟合，得到：计算可知，上式的显著性检验参数F=114.054，R2=0.962，可见模型的拟合效果是很好的。,(4.2.19),(4.2.20),比较(4.2.20)式与(4.2.19)式，并做简单计算可知：c0=2.048，c=1.154，a=8.353，所以，球状变异函数模型为：,(4.2.21),第三步：最优曲线的确定,首先考虑决定数 R2残差平方和RSS变程块金值基台值,二、克立格插值方法,克立格（Kriging）插值法，又称空间局部估计或空间局部插值法，是地统计学的主要内容之一。克立格法是建立在变异函数理论及结构分析基础之上的，它是在有限区域内对区域化变量的取值进行无偏最优估计的一种方法。克立格法适用的条件是，如果变异函数和相关分析的结果表明区域化变量存在空间相关性。其实质是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏、最优估计。,克里格插值（Kriging Interpolation），是根据变异函数模型而发展起来的一系列地统计的空间插值方法，包括：普通克里格法（Ordinary Kriging）、泛克里格法（Universal Kriging）、指示克里格法（Indicator Kriging）析取克里格法（Disjunctive Kriging）、协同克里格法（Cokriging）等。下面仅对普通克立格法作一些简单介绍。,首先假设区域化变量 满足二阶平稳假设和本征假设，其数学期望为m，协方差函数 及变异函数 存在。即：假设在待估计点（x）的临域内共有n个实测点，即 x1，x2，xn，其样本值为。那么，普通克里格法的插值公式为：,其中 为权重系数，表示各空间样本点 处的观测值 对估计值 的贡献程度。可见，克立格插值的关键就就是计算权重系数。显然，权重系数的求取必须满足两个条件：一是使 的估计是无偏的，即偏差的数学期望为零；二是最优的，即使估计值 和实际值 之差的平方和最小。为此，需要满足以下两个条件：,（1）无偏性。要使 成为 的无偏估计量，即。当 时，也就是当 时，则有 这时，为 的无偏估计量。（2）最优性。在满足无偏性条件下，估计方差为：,使用协方差函数表达，它可以进一步写为：为使估计方差最小，根据拉格朗日乘数原理，令：(4.2.25)求F对 和 的偏导数，并令其为0，得克立格方程组：(4.2.26),(4.2.24),整理后得：(4.2.27)解线性方程组(4.2.27)式，求出权重系数 和拉格朗日乘数，代入公式(4.2.24)，经过后可得克立格估计方差，即(4.2.28),在变异函数存在的条件下，根据协方差与变异函数的关系：或，也可以用变异函数表示普通克立格方程组和克立格估计方差，即：(4.2.29)(4.2.30),上述过程也可用矩阵形式表示，令：则普通克立格方程组为：(4.2.31)解方程组(4.2.31)式，可得：(4.2.32)其估计方差为：(4.2.33),也可以将克立格方程组和估计方差，用变异函数写成上述矩阵形式。令：在以上的介绍中，区域化变量 的数学期望 可以是已知或未知的。如果m是已知常数，称为简单克立格法；如果m是未知常数，称为普通克立格法。不管是那一种方法，均可根据方法计算权重系数和克立格估计量。,以图4.2.1为例，四个观测点x1、x2、x3、x4的观测值分别为Z(x1)=37、Z(x2)=42、Z(x3)=36、Z(x4)=35。如果假设降水量的变异函数是各向同性（即变异函数在各个方向是的变化都相同）的二维球状模型，其具体形式为（4.2.21）式。现在，我们用普通克立格法估计观测点x0的降水量值Z(x0)。根据普通克立格法的基本原理，我们知道，Z(x0)估计的基本公式应该是：,根据公式(4.2.32)，可知：根据协方差与变异函数的关系以及（4.2.21）式，可得协方差函数：,(4.2.37),比较(4.2.20)式与(4.2.19)式，并做简单计算可知：c0=2.048，c=1.154，a=8.353，所以，球状变异函数模型为：,(4.2.21),当 时，根据克立格矩阵的对称性，当 时，由此计算可得：,将以上计算结果代入克立格方程组（4.2.31），得：,即克立格权重系数分别为：1=0.287，2=0.210，3=0.202，4=0.301，=-0.473，所以点的降水量的克立格估计值为：根据普通克立格法的基本原理，我们知道，Z(x0)估计的基本公式应该是：37.25（mm）。克立格估计方差为：,应用实例,年降水量和蒸发量，既服从地带性规律，同时又受随机性因素的影响，因此它们是典型的区域化变量。我们以甘肃省53个气象台站多年平均降水量和蒸发量数据为实测值，拟合了年降水量和蒸发量的半变异函数理论模型，并采用普通克里格法和双变量协同克里格法，做了空间插值计算，结论如下。,（一）半变异函数半变异函数模型，是克立格空间插值的前提条件，同时它也决定着空间插值的精度。一般情况下，半变异函数模型是根据半变异函数云图的分布，选择合适的理论模型，按照估计方差最小的原则，运用最小二乘法求得。图4.2.4和图4.2.5分别给出了年降水量和年蒸发量的半变异函数云图。图4.2.4 年降水量的半变异函数云图,图4.2.5 年蒸发量的半变异函数云图 从图4.2.4和图4.2.5可以看出，年降水量和年蒸发量的块金效应都不明显，这是因为样本点是各个气象站点的实测值，空间分辨率可以忽略不计，另外实验误差和人为性误差基本上都很小。我们选择各种不同的半变异函数理论模型，经过多次拟合计算和对比分析，发现指数模型比较好地描述了年降水量的空间变异规律。其变异函数的具体形式如下：,(4.2.38)(4.2.38)式拟合的适度系数为。我们选择各种不同的半变异函数理论模型，经过多次拟合计算和对比分析，发现球状模型比较好地描述了年蒸发量的空间变异规律。其变异函数的具体形式如下：(4.2.39)(4.2.39)式拟合的适度系数为。,空间插值结果基于半变异函数的理论模型(4.2.38)和(4.2.39)，对甘肃省范围内的年降水量和蒸发量，用普通克里格法进行空间插值计算，得到的结果分别如图4.2.4和图4.2.5。结果讨论 从图4.2.6可以看出，在甘肃省范围内，年降水量的空间分布格局总体上是东南多西北少，并且呈现从东南方向到西北方向逐渐过渡，梯度变化明显；山地多，平地少，南北方向从南部队祁连山脉向北部的沙漠戈壁逐渐减少。,空间插值结果,图4.2.6 甘肃省年降水量的普通克里格空间插值结果,结果讨论从图4.2.6可以看出，在甘肃省范围内，年降水量的空间分布格局总体上是东南多西北少，并且呈现从东南方向到西北方向逐渐过渡，梯度变化明显；山地多，平地少，南北方向从南部队祁连山脉向北部的沙漠戈壁逐渐减少。,年降水量的空间变程很大，最多的东南部是最少的西北部的近10倍，其中，甘南东南部玛曲和禄曲、陇南东南部以及平凉灵台东南地区，年降水量达到691.59 786.75mm之间。400mm等降水线靠近兰州附近，而到了西北端，几乎整个酒泉市、嘉峪关市和张掖市的西北部，年降水量只有59.17-102.08mm。,结果讨论,图4.2.7 甘肃省年蒸发量的普通克里格空间插值结果从图4.2.7可以看出，年蒸发量的空间格局，恰好与年降水量的空间格局相反：西北多、东南少，呈现出由西北向东南逐渐减少的变化趋势，梯度变化明显。,年蒸发量的空间变程虽然小于年降水量，但仍然较大，在西北端的酒泉大部分地区以及民勤北部的腾格里沙漠地区，年蒸发量可以达到2931.30-3522.76mm之间，而在甘南玛曲的部分地区，只有1024.54-1179.88mm。兰州正好处于我国干旱和半干旱区的过渡地带，年蒸发量大致介于1389.77-1508.66mm之间。,
