第九讲地统计分析方法

第九讲地统计分析方法第一页，共五十四页，编辑于2023年，星期一2基础知识地统计（Geostatistics）又称地质统计，是在法国著名统计学家G.Matheron大量理论研究的基础上逐渐形成的一门新的统计学分支。它是以区域化变量为基础，借助变异函数，研究既具有随机性又具有结构性，或空间相关性和依赖性的自然现象的一门科学。凡是与空间数据的结构性和随机性，或空间相关性和依赖性，或空间格局与变异有关的研究，并对这些数据进行最优无偏内插估计，或模拟这些数据的离散性、波动性时，皆可应用地统计学的理论与方法。地统计分析的核心就是通过对采样数据的分析、对采样区地理特征的认识选择合适的空间内插方法创建表面。第二页，共五十四页，编辑于2023年，星期一3前提假设

随机过程地统计学认为研究区域中的所有样本值都是随机过程的结果，即所有样本值都不是相互独立的，它们是遵循一定的内在规律的。因此地统计学就是要揭示这种内在规律，并进行预测。正态分布若不符合正态分布的假设，应对数据进行变换，转为符合正态分布的形式，并尽量选取可逆的变换形式。

第三页，共五十四页，编辑于2023年，星期一4平稳性包括两种平稳性：一类是均值平稳；另一类是与协方差函数有关的二阶平稳和与半变异函数有关的内蕴平稳。均值平稳，即假设均值是不变的并且与位置无关；二阶平稳是假设具有相同的距离和方向的任意两点的协方差是相同的，协方差只与这两点的值相关而与它们的位置无关；内蕴平稳假设是指具有相同距离和方向的任意两点的方差（即变异函数）是相同的。二阶平稳和内蕴平稳都是为了获得基本重复规律而作的基本假设，通过协方差函数和变异函数可以进行预测和估计预测结果的不确定性。第四页，共五十四页，编辑于2023年，星期一5区域化变量当一个变量呈现一定的空间分布时，称之为区域化变量，它反映了区域内的某种特征或现象。区域化变量与一般的随机变量不同之处在于，一般的随机变量取值符合一定的概率分布，而区域化变量根据区域内位置的不同而取不同的值。而当区域化变量在区域内确定位置取值时，表现为一般的随机变量，也就是说，它是与位置有关的随机变量。区域化变量具有两个显著特征：即随机性和结构性。第五页，共五十四页，编辑于2023年，星期一

地统计学是以区域化变量理论为基础，以变异函数为主要工具，研究那些在空间分布上既有随机性又有结构性，或空间相关和依赖性的自然现象的科学。

协方差函数和变异函数是以区域化变量理论为基础建立起来的地统计学的两个最基本的函数。地统计学的主要方法之一，克立格法就是建立在变异函数理论和结构分析基础之上的。

第六页，共五十四页，编辑于2023年，星期一（一）协方差函数

协方差函数的概念

区域化随机变量之间的差异，可以用空间协方差来表示。在概率论中,随机向量X与Y的协方差被定义为

E(.)为期望值.区域化变量在空间点x和x+h处的两个随机变量和的二阶混合中心矩定义为Z(x)的自协方差函数，即（4.2.2）(4.2.1)第七页，共五十四页，编辑于2023年，星期一协方差函数的计算公式

式中：h为两样本点空间分隔距离或距离滞后；为在空间位置处的实测值；是在处距离偏离h的实测值[i=1，2，…，]，是分隔距离为h时的样本点对（paris）总数，和分别为和的样本平均数,即(4.2.3)(4.2.4)(4.2.5)第八页，共五十四页，编辑于2023年，星期一

若==m（常数），则上式可以改写为式中：m为样本平均数，可由一般算术平均数公式求得，即

(4.2.6)第九页，共五十四页，编辑于2023年，星期一（二）变异函数

变异函数的概念

变异函数variograms），又称变差函数、变异矩，是地统计分析所特有的基本工具。在一维条件下变异函数定义为，当空间点x在一维x轴上变化时，区域化变量Z(x)在点x和x+h处的值Z(x)与Z(x+h)差的方差的一半为区域化变量Z(x)在x轴方向上的变异函数，记为γ(h)，即

(4.2.7)方差等于平方均值减去均值的平方第十页，共五十四页，编辑于2023年，星期一

在二阶平稳假设条件下，对任意的h有因此，公式可以改写为

从上式可知，变异函数依赖于两个自变量x和h，当变异函数仅仅依赖于距离h而与位置x无关时，可改写成，即

(4.2.9)(4.2.8)第十一页，共五十四页，编辑于2023年，星期一变异函数的性质设Z(x)是区域化变量，在满足二阶平稳假设条件下，变异函数式具有如下性质：

(1)=0，即在h=0处，变异函数为0；

(2)=，即关于直线h=0是对称的，它是一个偶函数；

(3)≥0，即只能大于或等于0；第十二页，共五十四页，编辑于2023年，星期一变异函数的计算公式

设是系统某属性Z在空间位置x处的值，为一区域化随机变量，并满足二阶平稳假设，h为两样本点空间分隔距离，和分别是区域化变量在空间位置和处的实测值[i=1,2,…,N(h)]，那么，变异函数的离散计算公式为(4.2.10)第十三页，共五十四页，编辑于2023年，星期一

这样对不同的空间分隔距离h，计算出相应的和值。如果分别以h为横坐标，或为纵坐标，画出协方差函数和变异函数曲线图，就可以直接展示区域化变量Z(x)的空间变异特点。可见，变异函数能同时描述区域化变量的随机性和结构性，从而在数学上对区域化变量进行严格分析，是空间变异规律分析和空间结构分析的有效工具。第十四页，共五十四页，编辑于2023年，星期一例如：假设某地区降水量Z(x)（单位：mm）是二维区域化随机变量，满足二阶平稳假设，其观测值的空间正方形网格数据如图4.2.1所示（点与点之间的距离为h=1km）。试计算其南北方向及西北和东南方向的变异函数。第十五页，共五十四页，编辑于2023年，星期一图4.2.1空间正方形网格数据（点间距h=1km）

从图4.2.1可以看出，空间上有些点，由于某种原因没有采集到。如果没有缺失值，可直接对正方形网格数据结构计算变异函数；在有缺失值的情况下，也可以计算变异函数。只要“跳过”缺失点位置即可（图4.2.2）。第十六页，共五十四页，编辑于2023年，星期一

首先计算南北方向上的变异函数值，由变异函数的计算公式可得

=385/72=5.35图4.2.2缺失值情况下样本数对的组成和计算过程

☉为缺失值

第十七页，共五十四页，编辑于2023年，星期一

同样计算出最后，得到南北方向和西北—东南方向上的变异函数计算结果见下表。同样可以计算东西方向上的变异函数。

方向

南北

方向

西北—东南

h12345h1.412.824.245.657.07N(h)

362721135N(h)

322113825.359.2617.5525.6922.907.0612.9530.8558.1350.00第十八页，共五十四页，编辑于2023年，星期一变异函数的参数

变异函数有４个非常重要的参数，即基台值（sill）、变程（range）或称空间依赖范围（rangeofspatialdependence）、块金值（nugget）或称区域不连续性值（localizeddiscontinuity）和分维数（fractaldimension）。前3个参数可以直接从变异函数图中得到。它们决定变异函数的形状与结构。变异函数的形状反映自然现象空间分布结构或空间相关的类型，同时还能给出这种空间相关的范围。第十九页，共五十四页，编辑于2023年，星期一

当变异函数随着间隔距离h的增大，从非零值达到一个相对稳定的常数时，该常数称为基台值C0+C。当间隔距离h=0时，γ(0)=C0，该值称为块金值或块金方差（nuggetvariance）。基台值是系统或系统属性中最大的变异，变异函数达到基台值时的间隔距离a称为变程。变程表示在h≥a以后，区域化变量Z(x)空间相关性消失。块金值表示区域化变量在小于抽样尺度时非连续变异，由区域化变量的属性或测量误差决定。

第二十页，共五十四页，编辑于2023年，星期一

上述３个参数可从变异函数曲线图直接得到，或通过估计曲线回归参数得到。第4个参数，即分维数用于表示变异函数的特性，由变异函数和间隔距离h之间的关系确定分维数D为双对数直线回归方程中的斜率，它是一个无量纲数。分维数D的大小，表示变异函数曲线的曲率，可以作为随机变异的量度。但该随机分维数D与形状分维数有本质的不同。第二十一页，共五十四页，编辑于2023年，星期一变异函数的理论模型

地统计学将变异函数理论模型分为3大类：第1类是有基台值模型，包括球状模型、指数模型、高斯模型、线性有基台值模型和纯块金效应模型；第2类是无基台值模型，包括幂函数模型、线性无基台值模型、抛物线模型；第3类是孔穴效应模型。下面有代表性地介绍几种常见的变异函数理论模型。第二十二页，共五十四页，编辑于2023年，星期一

①纯块金效应模型:其一般公式为式中：c0>0，为先验方差。该模型相当于区域化变量为随机分布，样本点间的协方差函数对于所有距离h均等于0，变量的空间相关不存在。

(4.2.11)第二十三页，共五十四页，编辑于2023年，星期一②球状模型:其一般公式为

式中：c0为块金（效应）常数;c为拱高;c0+c为基台值;a为变程。当c0=0，c=1时，称为标准球状模型。球状模型是地统计分析中应用最广泛的理论模型，许多区域化变量的理论模型都可以用该模型去拟合。

(4.2.12)第二十四页，共五十四页，编辑于2023年，星期一

③指数模型:其一般公式为式中：c0和c意义与前相同，但a不是变程。当h=3α时，，即，从而指数模型的变程约为。当c0=0，c=1时，称为标准指数模型。(4.2.13)第二十五页，共五十四页，编辑于2023年，星期一④高斯模型:其一般公式为式中：c0和c意义与前相同，a也不是变程。当时，，即，因此高斯模型的变程约为。当时，称为标准高斯函数模型。(4.2.14)第二十六页，共五十四页，编辑于2023年，星期一⑤幂函数模型:其一般公式为式中：θ为幂指数。当θ变化时，这种模型可以反映在原点附近的各种性状。但是θ必须小于2，若，则函数就不再是一个条件非负定函数了，也就是说它已经不能成为变异函数了。

(4.2.15)第二十七页，共五十四页，编辑于2023年，星期一⑥对数模型:其一般公式为显然，当，这与变异函数的性质不符。因此，对数模型不能描述点支撑上的区域化变量的结构。(4.2.16)第二十八页，共五十四页，编辑于2023年，星期一

⑦线性有基台值模型:其一般公式为

式中:该模型的变程为a，基台值为。

⑧线性无基台值模型:其一般公式为

从式中可以看出，该模型没有基台值，也没有变程。

(4.2.18)(4.2.17)第二十九页，共五十四页，编辑于2023年，星期一例如:某地区降水量是一个区域化变量，其变异函数的实测值及距离h的关系见下表，下面我们试用回归分析方法建立其球状变异函数模型。实测值γ(h)距离h实测值γ(h)距离h2.10.69.24.94.31.110.35.15.72.210.56.26.52.510.97.57.83.111.29.58.83.812.49.8第三十页，共五十四页，编辑于2023年，星期一

从上面的介绍和讨论，我们知道，球状变异函数的一般形式为当时，有第三十一页，共五十四页，编辑于2023年，星期一

如果记，则可以得到线性模型

根据表中的数据，对上式进行最小二乘拟合，得到

(4.2.20)

计算可知，上式的显著性检验参数F=114.054，R2=0.962，可见模型的拟合效果是很好的。(4.2.19)第三十二页，共五十四页，编辑于2023年，星期一

比较(4.2.20)式与(4.2.19)式，并做简单计算可知：c0=2.048，c=1.154，a=8.353，所以，球状变异函数模型为(4.2.21)第三十三页，共五十四页，编辑于2023年，星期一(三)克立格插值方法

克立格（Kriging）插值法，又称空间局部估计或空间局部插值法，是地统计学的主要内容之一。克立格法是建立在变异函数理论及结构分析基础之上的，它是在有限区域内对区域化变量的取值进行无偏最优估计的一种方法。克立格法适用的条件是，如果变异函数和相关分析的结果表明区域化变量存在空间相关性。

其实质是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏、最优估计。

第三十四页，共五十四页，编辑于2023年，星期一

克立格插值（Ｋriginginterpolation)是根据变异函数模型而发展起来的一系列地统计的空间插值方法，包括：普通克立格法（ordinaryＫriging）;

泛克立格法（universalＫriging）;

指示克立格法（indicatorＫriging）;

析取克立格法（disjunctiveＫriging）;

协同克立格法（cokriging）等。下面仅对普通克立格法作一些简单介绍。

第三十五页，共五十四页，编辑于2023年，星期一

首先假设区域化变量满足二阶平稳假设和本征假设，其数学期望为m，协方差函数及变异函数存在。即假设在待估计点（x）的临域内共有n个实测点，即x1，x2，…，xn，其样本值为。那么，普通克里格法的插值公式为

(4.2.22)第三十六页，共五十四页，编辑于2023年，星期一

其中为权重系数，表示各空间样本点处的观测值对估计值的贡献程度。可见，克立格插值的关键就是计算权重系数。显然，权重系数的求取必须满足两个条件：一是使的估计是无偏的，即偏差的数学期望为零；二是最优的，即使估计值和实际值之差的平方和最小。为此，需要满足以下两个条件：第三十七页，共五十四页，编辑于2023年，星期一(1)无偏性。要使成为的无偏估计量，即

当时，也就是当时，则有这时，为的无偏估计量。（2）最优性。在满足无偏性条件下，估计方差为(4.2.23)第三十八页，共五十四页，编辑于2023年，星期一

使用协方差函数表达，它可以进一步写为

(4.2.24)

为使估计方差最小，根据拉格朗日乘数原理，令

(4.2.25)

求F对和的偏导数，并令其为0，得克立格方程组

(4.2.26)

第三十九页，共五十四页，编辑于2023年，星期一

(4.2.27)(4.2.28)整理后得

解线性方程组(4.2.27)式,求出权重系数λi和拉格朗日系数μ,代入公式(4.2.24),可得克立格估计方差第四十页，共五十四页，编辑于2023年，星期一

在变异函数存在的条件下，根据协方差与变异函数的关系：

或，也可以用变异函数表示普通克立格方程组和克立格估计方差，即

(4.2.29)

解线性方程组(4.2.27)式，求出权重系数和拉格朗日乘数μ，代入公式(4.2.24)，可得克立格估计方差，即

(4.2.30)第四十一页，共五十四页，编辑于2023年，星期一上述过程也可用矩阵形式表示，令

则普通克立格方程组为

(4.2.31)解方程组(4.2.31)式，可得

(4.2.32)其估计方差为

(4.2.33)

第四十二页，共五十四页，编辑于2023年，星期一

也可以将克立格方程组和估计方差用变异函数写成上述矩阵形式。令

在以上的介绍中，区域化变量的数学期望可以是已知或未知的。如果m是已知常数，称为简单克立格法；如果m是未知常数，称为普通克立格法。不管是哪一种方法，均可根据方法计算权重系数和克立格估计量。

(4.2.34)(4.2.35)(4.2.36)第四十三页，共五十四页，编辑于2023年，星期一

以图4.2.1为例，4个观测点x1，x2，x3，x4的观测值分别为Z(x1)=37、Z(x2)=42、Z(x3)=36、Z(x4)=35，如果假设降水量的变异函数是向同性（即变异函数在各个方向的变化都相同）的二维球状模型，其具体形式为（4.2.21）式。现在，我们用普通克立格法估计观测点x0的降水量值Z(x0)。根据普通克立格法的基本原理，我们知道，Z(x0)估计的基本公式应该是

第四十四页，共五十四页，编辑于2023年，星期一

根据公式(4.2.32)，可知

(4.2.37)

根据协方差与变异函数的关系以及（4.2.21）式，可得协方差函数

第四十五页，共五十四页，编辑于2023年，星期一

当时，根据克立格矩阵的对称性，当时，，由此计算可得

第四十六页，共五十四页，编辑于2023年，星期一第四十七页，共五十四页，编辑于2023年，星期一

将以上计算结果代入克立格方程组（4.2.31），得

第四十八页，共五十四页，编辑于2023年，星期一

即克立格权重系数分别为：λ1=0.287，λ2=0.210,λ3=0.202,λ4=0.301，μ=-0.473，所以观测点的降水量的克立格估计值为：根据普通克立格法的基本原理，我们知道，Z(x0)估计的基本公式应该是

37.25（mm）。

克立格估计方差为

第四十九页，共五十四页，编辑于2023年，星期一二、应用实例

年降水量和蒸发量，既服从地带性规律，同时又受随机性因素的影响，因此它们是典型的区域化变量。我们以甘肃省53个气象台站多年平均降水量和蒸发量数据（见教材表3.1.2）为实测值，拟合了年降水量和蒸发量的半变异函数理论模型，并采用普通克立格法和双变量协同克里格法，做了空间插值计算，结论如下。