1、地统计分析方法地统计分析方法地统计学 地统计学（Geostatistics）又称地质统计学，是在法国著名统计学家G.Matheron大量理论研究的基础上逐渐形成的一门新的统计学分支。地统计学：以区域化变量理论为基础，以变异函数为主要工具，研究在空间分布上既有随机性又有结构性，或空间相关和依赖性的自然现象的科学。地统计学的理论基础是区域化变量理论，两个最基本的函数是协方差函数和变异函数，克立格法是其主要方法之一。地统计学与经典统计学 经典统计学的相同点：经典统计学的相同点：它们都是在大量采样的基础上，通过对样本属性值的频率分布或均值、方差关系及其相应规则的分析，确定其空间分布格局与相关关系。地统
2、计学的优势：地统计学的优势：地统计学既考虑到样本值的大小，又重视样本空间位置及样本间的距离，弥补了经典统计学忽略空间方位的缺陷。1.区域化变量 当一变量呈现为空间分布时，称之为区域化变量。区域化变量，亦称区域化随机变量，G.Matheron（1963）将它定义为以空间点x的三个直角坐标xu，xv，xw为自变量的随机场:Z Z（x x）=Z=Z（x xu u，x xv v，x xw w）区域化变量与一般的随机变量不同之处在于，一般的随机变量取值符合一定的概率分布，而区域化变量根据区域内位置的不同而取不同的值。而当区域化变量在区域内确定位置取值时，表现为一般的随机变量，也就是说，它是与位置有关的随
3、机变量与位置有关的随机变量。区域化变量的特征 区域化变量是一个随机随机变量，它具有局部的、随机的、异常的特征。区域化变量具有一般的或平均的结构结构性质，即变量在点x与偏离空间距离为h的点x+h处的值Z(x)和Z(x+h)具有某种程度的相似性，即自相关性，这种自相关性的程度依赖于两点间的距离h及变量特征。区域化变量还具有空间局限性空间局限性（即这种结构性表现为一定范围内）、不同程度的连续性不同程度的连续性和不同程度的各向不同程度的各向异性异性（即各个方向表现出的自相关性有所区别）等特征。2.协方差函数 空间协方差函数和变异函数能反映区域化变量的随机性和结构性。1 1）协方差函数的概念：）协方差函
4、数的概念：在概率论中，随机变量X与Y的协方差定义为：Cov(X,Y)=E(X-E(X)(Y-E(Y)(3.8.1)Cov(X,Y)=E(X-E(X)(Y-E(Y)(3.8.1)区域化变量在空间点x和x+h处的两个随机变量的二阶混合中心矩定义为Z(x)的自协方差函数，即:Cov(Z(x),Z(x+h)=EZ(x)Z(x+h)-EZ(x)EZ(x+h)(3.8.2)容易看出，该函数依赖于空间点x和向量h。2）协方差函数的计算公式 若若 （常数）（常数）则公式（则公式（3.8.33.8.3）可改写为）可改写为 式中，m为样本平均数，可由一般平均数求得，即：mhxZxZii)()(2)(1)()()(
5、1)(*mhxZxZhNhchNiii(3.8.6)niixZNm1)(13.变异函数1 1）变异函数的概念）变异函数的概念2）变异函数的性质 设Z(x)是区域化变量，在满足二阶平稳假设条件下，变异函数公式（3.8.9）具有如下性质：(0)=0，在h=0时，变异函数为0；(h)=(-h)，即(h)关于直线h=0是对称的，它是一个偶函数；(h)0，(h)表示的方差只能大于或等于0；|h|时，(h)c(0)，或()=c(0)，即当空间距离增大时，变异函数接近先验方差-(h)必须是一个条件非负定函数，由-(xi-xj)构成的变异函数矩阵在条件 时，为非负定的。NiimxZhNc122)()(1)0(
6、nii103）变异函数的计算公式 设Z(x)是系统某属性Z在空间位置x处的值，Z(x)为一区域化随机变量，并满足二阶平稳假设，h为两样本点空间分隔距离，Z(xi)和Z(xi+h)分别是区域化变量Z(x)在空间位置xi和xi+h处的实测值i=1,2,N(h)，那么，变异函数r(h)的离散计算公式为：2)(1)()()(21)(*hNiiihxZxZhNhr(3.8.10)因此，对不同的空间分隔距离h，计算出相应的c(h)和r(h)值。如果分别以h为横坐标，c*(h)或 r*(h)为纵坐标，画出协方差函数和变异函数曲线图，就可以直接展示区域化变量Z(x)的空间变异特点。可见，变异函数能同时描述区域
7、化变量的随机性和结构性，从而在数学上对区域化变量进行严格分析,是空间 变异规律分析和空 间结构分析的有效 工具。)()0()(hCChr例1 假设某地区降水量Z(x)(单位:mm)是二维区域化随机变量，满足二阶平稳假设，其观测值的空间正方形网格数据如图所示(点与点之间的距离为h=1km)。试计算其南北方向及西北和东南方向的变异函数。从上图可以看出，空间上有些点，由于某种原因没有采集到。如果没有缺失值，可直接对正方形网格数据结构计算变异函数；在有缺失值的情况下，也可以计算变异函数。只要“跳过”缺失点位置即可。缺失值情况下样本数对的组成和计算过程为缺失值 首先计算南北方向上的变异函数值，由变异函数
8、的计算公式可得222222)3837()3836()3635()3537()3742()4240(3621)1(35.572/385)3228()3438()3029()2935()3533()3340()4036()2932()3236()3637()3741()4137()3336()3637()3739()3939()3034()3437()3738()3839()3940(3535)3535()3542()4242()3536()3743()4340()3735()3538(222222222222222222222222222222）（同样计算出：同样计算出：同理可以计算东西方向和西
9、北东南方向上的变异函数。最后，得到南北方向和西北东南方向上的变异函数计算结果见下表：9.26)2(17.55)3(25.69)4(22.90)5(4）变异函数的参数 变异函数的四个重要参数：基台值、变程、块金值和分维数。前三个参数可以从变异函数图中得到。基台值基台值当变异函数(h)随着间隔距离h的增大，从非零值达到一个相对稳定的常数时，该常数称为基台值C0+C，它是系统或系统属性中最大的变异；变程变程变异函数(h)达到基台值时的间隔距离a称为变程。表示在ha以后，区域化变量Z(x)空间相关性消失。块金值块金值当间隔距离h=0时，(0)=C(0)，该值称为块金值或块金方差。表示区域化变量在小于抽
10、样尺度时非连续变异，由区域化变量的属性或测量误差决定。分维数分维数表示变异函数的特性，由变异函数(h)和间隔距离h之间的关系确定：hDhr)24()(25）变异函数的理论模型 地统计学将变异函数理论模型分为3大类：有基台值模型，包括球状模型、指数模型、高斯模型、线性有基台值模型和纯块金效应模型；无基台值模型，包括幂函数模型、线性无基台值模型、抛物线模型；孔穴效应模型。纯块金效应模型纯块金效应模型:其一般公式为 式中：c00，为先验方差。该模型相当于区域化变量为随机分布，样本点间的协方差函数对于所有距离h均等于0，变量的空间相关不存在。000)(0hchh 球状模型球状模型:其一般公式为 式中：
11、c0为块金（效应）常数;c为拱高;c0+c为基台值;a为变程。当c0=0，c=1时，称为标准球状模型。球状模型是地统计分析中应用最广泛的理论模型，许多区域化变量的理论模型都可以用该模型去拟合。ahccahahahcchh03300)223(00)(指数模型指数模型:其一般公式为式中：c0和c意义与前相同，但a不是变程。当h=3时，即 ，从而指数模型的变程约为3。当c0=0，c=1时，称为标准指数模型。0)1(00)(0hecchhah195.0113eeahcca0)3(高斯模型高斯模型:其一般公式为式中：c0和c意义与前相同，a也不是变程。当 时 ，即 ，因此高斯模型的变程约为 .当c0=0
12、,c=1时，称为标准高斯函数模型。0)1(00)(220hecchhahah3195.011322eeahcca0)3(a3幂函数模型幂函数模型:其一般公式为式中：为幂指数。当变化时，这种模型可以反映在原点附近的各种性状。但是必须小于2，若2，则函数r(-h)就不再是一个条件非负定函数了，也就是说它已经不能成为变异函数了。20,)(Ahh对数模型对数模型:其一般公式为显然，当h0,logh-，这与变异函数的性质 r(h)0不符。因此，对数模型不能描述点支撑上的区域化变量的结构。hAhlg)(线性有基台值模型线性有基台值模型:其一般公式为 式中:该模型的变程为a，基台值为c0+c。线性无基台值模
13、型线性无基台值模型:其一般公式为 从式中可以看出，该模型没有基台值，也没有变程。ahcahAhhh00c00c)(00)(0hAhhch 例例2:2:某地区降水量是一个区域化变量，其变异函数 的实测值及距离h的关系见下表，下面我们试用回归分析方法建立其球状变异函数模型。)(h从上面的介绍和讨论，我们知道，球状变异函数的一般形式为ahccahahahcchh03300)223(00)(ah 0330)2()23()(hachacch当 时，有 如果记 ，则可以得到线性模型 根据表中的数据，对上式进行最小二乘拟合，得到 计算可知，上式的显著性检验参数F=114.054，R2=0.962，可见模型的
14、拟合效果是很好的。32132100,21,23,),(hxhxacbacbcbhy22110 xbxbby2192007.0731.1048.2xxy比较两式，并做简单计算可知：c0=2.048，c=1.154，a=8.353，所以，球状变异函数模型为535.8202.3535.80)535.821535.823(154.1048.200)(33*hhhhhh(3.8.21)4.克立格（Kriging）法简介 1）克立格法概述 克立格法克立格法，又称空间局部估计或空间局部插值法，建立在变异函数理论及结构分析基础之上，是在有限区域内对区域化变量的取值进行无偏无偏最优估计最优估计的一种方法。克立格
15、法是根据待估样本点（或块段）有限邻域内若干已测定的样本点数据，考虑了样本点的形状、大小和空间相互位置关系，与待估样本点的相互空间位置关系，以及变异函数提供的结构信息，对待估样本点值进行的一种线性无偏最优估计。1）克立格法概述(1)(1)适用条件适用条件 变异函数和相关分析的结果表明区域化变量存在空间相关性。其实质实质是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏、最优估计。(2)(2)克立格法的类型克立格法的类型 普通克立格法；泛克立格法；协同克立格法；对数正态克立格法；指示克立格法；折取克立格法等。克立格法是一簇空间局部插值模型的总称。克立格法是一簇
16、空间局部插值模型的总称。2）克立格估计量 对于研究区域内任一点x的测量值Z(x)，其估计值 的估算公式为 估计量 是实际值Zv(xi)的克立格估计量克立格估计量。其中 为权重系数，是各已知样本Z(xi）在估计 时影响大小的系数，估计 的好坏取决于怎样计算或选择权重系数 。问题的关键在于求各点的权重系数。显然，权重系数的求取必须满足两个条件：两个条件：一是使 的估计是无偏的，即偏差的数学期望为零；二是最优的，即使估计值 和实际值Zv(x)之差的平方和最小。公式为：公式为：niiivxZxZ1*)()()(*xZv)(*xZvi)(*xZv)(*xZvi)(*xZv)(*xZv(3.8.22)mi
17、n)()()()(0)()(2*xxExxVarxxEZZZZZZvvvvvv(3.8.23)(3.8.24)3)普通克立格法 首先假设区域化变量Z(x)满足二阶平稳和本征假设，其数学期望为m，协方差函数c(h)及变异函数r(h)存在。即 假设在待估块段V的邻域 内，有一组n个已知样本 v(xi)(i=1,2,n)，其实测 值为Z(xi)(i=1,2,n)mxZE)(2)()()(mhxZxZEhc2)()(21)(hxZxZEh 克立格法的目标就是求一组权重系数i(i=1,2,n),使得加权平均值 成为待估地段V的平均值ZV（x0）的线性、无偏最优估计量，即克立格估计值。为此，需要满足以下两
18、个条件：niiiVxZZ1*)(无偏性。无偏性。要使 成为ZV的无偏估计量，即 ,当 时，也就是 当 时，则有 这时，为ZV的无偏估计量。最优性。最优性。在满足无偏性条件下，估计方差为mxZExZEniiiniii11)()()(*xZEZEVZV*mZEV*nii11ZV*212*2)()(niiiVVExZxZEZZE 使用协方差函数表达，它可以进一步写为 为使估计方差最小，根据拉格朗日乘数原理，令上式为 求公式F对i和的偏导数，并令其为0，得克克立格方程组立格方程组ninjniiijijiEVvcvvcVVc1112),(2),(),()1(212niiEFniiijinjjiFVvcv
19、vcF110)1(202),(2),(2（3.8.29）niiiEVvcVVc12),(),(整理后得:解上述n+1阶线性方程组,求出权重系数i和拉格朗日系数,可得克立格估计方差niiijinjjVvcvvc111),(),(（3.8.30）（3.8.31）在变异函数存在的条件下，根据协方差与变异函数的关系：c(h)=c(0)-r(h)或r(h)=c(0)-c(h)，也可以用变异函数表示普通克立格方程组和克立格估计方差，即niiijinjjVvvv111),(),(),(),(12VVVvniiiK 上述过程也可用矩阵形式表示，令 则普通克立格方程组为:解方程组式，可得:其估计方差为:1),(
20、),(),(,01111112121212222111211VvcVvcVvcDcccccccccKnnnnnnnnDKDK1DVVcTK),(2用变异函数表示普通克立格方程组和克立格方差例3 在例1中，假设降水量的变异函数为例2中的函数，它是一个各向同性的二维球状模型，已知四个观测点x1，x2，x3，x4的观测值分别为Z(x1)=37(mm)、Z(x2)=42(mm)、Z(x3)=36(mm)、Z(x4)=35(mm),试用普通克立格法内插估计观测点x0的降水量值Z(x0)。根据普通克立格法的基本原理，我们知道，Z(x0)估计的基本公式应该是 )4,3,2,1()(1*0ixZniiiZ 根
21、据公式，可知 根据协方差与变异函数的关系以及式，可得协方差函数1011111111040302011444342413433323124232221141312114321cccccccccccccccccccc535.80535.80)535.821535.823(1 154.1535.8202.3)(33*hhhhhhc 因为i=0,1,2,3,4;j=1,2,3,4,故 当i=j时，c11=c22=c33=c44=c(0)=c0+c=2.048+1.154=3.202 当ij时，由 ，及根据克立格矩阵的对称性，得jijiijxxxxcc202.3)(870.0)535.8)2(21535
22、.8223(154.1048.2202.3)2(202.3)11(202.33322042112ccc571.0)3(202.3952.0)1(202.3466.0)23(202.3383.0)14(202.3601.0)22(202.3711.0)12(202.3542.0)13(202.320320122422422433422322322024114223113ccccccccccccc 将以上计算结果代入克立格方程组，得 所以，克立格权重系数分别为：1=0.287，2=0.210,3=0.202,4=0.301，=-0.473。473.0301.0202.0210.0287.01870
23、.0571.0711.0952.0011111202.3383.0466.0711.01383.0202.3601.0542.01466.0601.0202.3870.01711.0542.0870.0202.314321 观测点x0的降水量的克立格估计值为：克立格估计方差为）（mmxZxZxZxZZ25.3735301.036202.042210.037287.0)(301.0)(202.0)(210.0)(287.04321*01.9304(mm)0.473-0.870)0.3010.5710.2020.7110.2100.952(0.2873.202),(),(412iiiKVvcVVc
