1、LOGO第六章空间统计学分析 空间统计分析方法由来 空间统计分析方法组成空间统计分析方法由分析空间变异与结构的半变异函数和用以空间局部估计的克立格插值法两个主要部分组成，是GIS空间分析的一个重要技术手段。由于空间现象之间存在不同方向、不同距离成分等相互作用，使得传统的数理统计方法无法很好地解决空间样本点的选取、空间估值和两组以上空间数据的关系等问题，因此，空间统计分析方法应运而生。空间统计分析方法的基本原理6.1空间自相关6.2空间局部估计6.3确定性插值法6.4探索性空间数据分析6.5ContentContents s6.1.1 空间统计分析的概念20世纪60年代，法国统计学家Mather
2、on G通过大量理论研究，形成了一门新的统计学分支，即空间统计学。空间统计学是以区域化变量理论为基础，以变异函数为主要工具，研究具有地理空间信息特性的事物或现象的空间相互作用及变化规律的学科。 自相关 空间统计分析的重要任务空间统计分析方法假设研究区中所有的值都是非独立的，相互之间存在相关性。在空间或时间范畴内，这种相关性被称为自相关。揭示空间数据的相关规律和利用相关规律进行未知点预测。由于空间统计分析包含这两个显著的任务，所以涉及两次使用样点数据，第一次用作估计空间自相关，第二次用作未知点预测。6.1.2 空间统计分析中的理论假设 1. 区域化变量区域化变量的两重性表现在观测前把它看成是随机
3、场，依赖于坐标（Xu，Xv，Xw），观测后是一个普通的空间三元函数值或一个空间点函数。 区域化变量是一种在空间上具有数值的实函数，它具有以下属性：空间局限性连续性各向异性区域化变量被限制于一定空间范围，这称为几何域。在几何域内，区域化变量的属性最为明显；在几何域外，不明显。不同的区域化变量具有不同程度的连续性，用区域化变量的半变异函数来描述。当区域化变量在各个方向上具有相同性质时称各向同性，否则称为各向异性。其它属性： 区域化变量在一定范围内呈一定程度的空间相关，当超出这一范围之后，相关性变弱甚至消失。 对于任一区域化变量，特殊的变异性可以叠加在一般的规律之上。 2. 协方差函数 CovZ(t
4、1)，Z(t2)=EZ(t1)EZ(t1)Z(t2)EZ(t2) (6.1) 在随机函数中，当只有一个自变量x时称为随机过程，随机过程Z(t)在时间t1和t2处的随机变量Z(t1)、Z(t2)的二阶混合中心矩定义为随机过程的协方差函数记为CovZ(t1)，Z(t2)，即 CovZ(x)，Z(x+h)=EZ(x)Z(x+h)EZ(x)EZ(x+h) (6.2) Cov(x，x+0)=EZ(x)2EZ(x)2 (6.3)当随机函数依赖于多个自变量时，Z(x)=Z（Xu，Xv，Xw）称为随机场，而随机场Z(x)在空间点x和x+h处的两个随机变量Z(x)和Z(x+h)的二阶混合中心矩定义为随机场Z(x
5、)的自协方差函数，即随机场Z(x)的自协方差函数亦称为协方差函数，一般地，协方差函数依赖于空间点x和向量h。当h=0时，协方差函数变为 3. 变异函数 (x，h)=1/2*VarZ(x)Z(x+h)2 =1/2*EZ(x)Z(x+h)21/2*EZ(x)EZ(x+h)2 (6.4)变异函数在一维条件下，当空间点x在一维x轴上变化时，区域变量Z(x)在点x和x+h处的值Z(x)与Z(x+h)差的方差一半定义为区域变量Z(x)在x轴上的变异函数，记为(x，h)，即 在二阶平稳假设条件下对任意h有 EZ(x+h)=EZ(x) 因此，式（6.4）可改写为 (x，h)=1/2*EZ(x)Z(x+h)2
6、(6.5) 从式（6.5）可知，变异函数依赖于x和h，当变异函数仅依赖于h，与x无关时，变异函数(x，h)可改写成(h)，即 (h)=1/2*EZ(x)Z(X+h)2 (6.6) 4. 平稳性假设及内蕴假设 （1）平稳性假设 设某一随机函数Z(x)，其空间分布律不因平移而改变，即若对任一向量h，关系式 G(z1,z2,x1,x2,)=G(z1,z2,x1+h,x2+h,)成立时，则该随机函数为平稳性随机函数。确切的说，无论位移向量h多大，两个k维向量的随机变量Z(x1),Z(x2),Z(xk)和Z(x1+h),Z(x2+h),Z(xk+h)有相同的分布律。 当区域化变量满足下列两个条件时，称该
7、区域化变量满足二阶平稳： 在整个研究区内，区域化变量Z(x)的数学期望对任意x存在且等于常数，即EZ(x)=m(常数)，任意x。 在整个研究区内，区域化变量的空间协方差函数对任意x和h存在且平稳，即 CovZ(x)，Z(x+h)=EZ(x)Z(x+h)-m2=C(h)，任意x,h （2）内蕴假设 一些自然现象和随机函数具有无限离散性，这时区域化变量Z(x)的增量Z(x)-Z(x+h)满足下列两个条件时，就称该区域化变量满足内蕴假设：在整个研究区内随机函数Z(x)的增量的数学期望为0，即 EZ(x)-Z(x+h)=0， 任意x,h 对于所有矢量的增量的方差函数存在且平稳 VarZ(x)-Z(x+
8、h)=EZ(x)-Z(x+h)2 =2(x,h)=2(h)，任意X,h 即要求Z(x)的半变异函数存在且平稳。 内蕴假设可以理解为： 随机函数Z(x)的增量Z(x)-Z(x+h)只依赖于分隔它们的向量h，而不依赖于具体位置x，这样，被向量h分割的每一对数据Z(x)，Z(x+h)可以看成是一对随机变量Z(x1)，Z(x2)的一个不同现实，而半变异函数(h)的估计量*(h)为 (h)=1/2N(h)*Z(xi)-Z(xi+h)2 式中，N(h)是被向量h相分隔的试验数据对的数目。 准平稳假设如果随机函数只在有限大小的邻域内是平稳的，则称该随机函数服从准平稳假设。准平稳（或准内蕴）假设是一种折中方案
9、，它既考虑到某现象相似性的尺度，也顾及到有效数据的多少。6.2.1 空间自相关理论 空间自相关性在空间统计分析中，相关分析可以检测两种现象的变化是否存在相关性，若所分析的统计量为不同观察对象的同一属性变量，则称之为自相关。通过检测一个位置上的变异是否依赖于邻近位置的变异来判断该变异是否存在空间自相关性根据变异的性质可以将变异分为三种类型：绝对型变异，等级型变异和连续型变异 空间自相关 各向同/异性空间自相关是针对同一个属性变量而言的，当某一测样点属性值高，而其相邻点同一属性值也高时，为空间正相关；反之，为空间负相关。当空间自相关仅与两点间距离有关时，称为各向同性；否则为各向异性。6.2.2 空
10、间自相关分析方法空间自相关方法按功能大致分为两类： 全域型自相关、区域型自相关最为常用的计算空间自相关方法是：Morans I、Gearys C、Getis、Join count以及空间自相关系数图等 1. Morans I法若在区域内有n个空间单元，每个空间单元皆有一个观察值X，空间单元i与空间单元j的空间关系构成Wij的空间相邻矩阵，以1表示i和j相邻，以0表示i和j不相邻。其简单定义为 Wijnn 其中，Wij为表示区位相邻矩阵，Wij=1表示区位相邻，Wij=0则表示区位不相邻。 Moran Index值是应用较广泛的一种空间自相关性判定指标，其计算式为 式中， ， 。Wij表示区位相
11、邻矩阵；Cij表示属性相似矩阵；Xi和Xj分别为i和j空间单元属性数据值，Wij=1代表空间单元相邻，Wij=0代表不相邻，ij，Wii=0。ninjniiijjninjiijninjijninjijijXXnWXXXXWSWCWI11121111211)(1)( )(（6.16） niiXXnS122)(1)(XXXXCjiij 若母体为随机分配，常采用统计验证的方式进一步判定Moran Index的期望值和变异数。I的期望值为 其变异数为其中， ； ； ； )()()(IVarIEIIZ) 1(1)(nIE（6.17） 2202021220212)() 3)(2)(1(62)(3) 33(
12、)(IEnnnWWnWWnnKWnWWnnnIVar ninjijWW110ninjjiijWWW1121)(niiiWWW122)(21214njjniinXXnXXk；Wi和和Wi为相关权重矩阵为相关权重矩阵i及及j行的总和。行的总和。 I值结果一定介于-1到1之间； I0为正相关，数值越大表示空间分布的相关性越大，即空间上聚集分布的现象越明显； I0（正相关） I 0（负相关） 图6.1 空间自相关正负结果示意图 根据各空间间隔自相关值的计算，Morans I公式可改写为 其中，d代表空间间隔；Wij代表区位相邻矩阵。d=1代表空间单元是相邻的；d=2定义为与间隔一个的空间单元相接邻，而
13、与原来的空间单元不相邻。 niijininjijninjijXXXXXXdWdWndI121111)()()()()(（6.19） 区域空间自相关的定义为 其中，Ii为Local Moran Index，Wij为区位相邻矩阵。（6.20） njjiijiXXXXWI1)( 2. Gearys Contiguity Ratio C法 与Morans I类似，其表达式为 ninjijninjiijniiWyyWyynC11112_12_2)()(1（6.21） C = 1，表示不相关；0 C 1表示负相关。 3. Getis统计法 Anselin曾归纳各种空间聚集的研究方法，该方法经常表达为 其中
14、，Wij代表i与j的空间关系，即类似上述空间相邻权重矩阵Wij；而yij则是i与j的观察式。 jijijyw（6.22） 全域型Getis 其中，wij(d)为距离d内的空间相邻权重矩阵。 若i与j相邻，wij(d)=1；若i与j不相邻，wij(d)=0。 区域型Getis 可量测每一个i在距离d的范围内，与每个j的相关程度。 njjnjjiijxxxdwdG11)()(ij （6.23） njjnjjijixxdwdG11)()(ij （6.24） 4. 空间自相关系数图分析法（以某地区为例）空间自相关空间间隔135791113151719-1.00-0.501.000.500.001.50
15、图6.2 某地区某种空间自相关系数图 （1）图中有两处隆起处，代表微视尺度及宏观尺度上，存在显著的聚集分布现象，但聚集现象不存在于中观尺度上。（2）空间间隔为2时，空间自相关值有波峰，即在空间间隔为2时，其空间分布有最大的自相关性。 空间局部估计 常见的克立格插值模型有： 插值一般分为两步：空间局部估计也称空间局部插值，它是利用在地表不同位置采集的样点生成一个连续表面。普通克立格、简单克立格、泛克立格、概率克立格、指示克立格、析取克立格及协同克立格等（1）样点空间结构量化分析；（2）对未知点进行预测 6.3.1 6.3.1 半变异函数分析半变异函数分析 1. 半变异函数及其性质 半变异函数是一
16、个关于数据点的半变异值与数据点间距离的函数，设区域化变量Z(xi)和Z(xi+h)分别是Z(x)在空间位置xi和xi+h上的观测值(i=1,2,N(h)，则半变异函数可由下式进行估计 其中，N(h)是分隔距离为h的样本量。 )(12)()()(21)(hNiiihxZxZhNh（6.25） 主要的理论半变异函数模型有： 常见的分析半变异模型包括：球状模型、指数模型、高斯模型、幂函数模型和抛物线模型等有基台值模型、无基台值模型、孔穴效应模型 2. 影响半变异函数的主要因素 （1）样点间的距离和支撑的大小为了使建立的半变异函数模型能准确地反映各种尺度上的变化特征，要确定采样的最小尺度。在采样之前，
17、首先需要在满足精度的前提下确定最佳的采样尺度。用块段取样时，要考虑支撑的大小，一般采用正则化变量消除其影响。（2）样本数量的大小（3）异常值的影响样本数量在对地统计学中主要指计算实际半变异函数值时的点对数目实际取样工作中点对数目不能无限，一般要求在变程a以内各距离上的点对数目不应小于20对。在小尺度距离上相对要多一些，大尺度距离相对少一些。如果异常值比较多，块金值C0要增大，随机成分的影响加强，而空间自相关的影响消弱。对于半变异函数的模型来讲，块金效应值C0越小越好。 （4）比例效应的影响如果平均值和标准差之间存在明显的线性关系，则比例效应存在，反之亦然。当样品方差随着平均值的增加而增加时，称
18、正比例效应，反之亦然。比例效应的存在会使实际半变异函数值产生畸变，消除比例效应的方法主要是通过对原始数据取对数，或者通过相对半变异函数的求解。 （5）漂移的影响 当漂移存在时，半变异函数值不再是半变异函数的无偏估计。要消除漂移对半变异函数的影响，主要通过建立合适的漂移形式，即EZ(x)=m(x)中，m(x)的函数式，它使半变异函数曲线真实地符合实际半变异函数值。 3. 半变异模型的合并 假设数据中有两个独特的结构，只用单一模型无法表达，就可以用两个单独的模型来模拟这个半变异图，然后将它们合并为一个模型。 4. 半变异模型的步长分组与步长大小的选择 图6.5 12个样点两两形成的样点对示意图在所
19、有样点中两两之间均能形成样点对，如上图。要在半变异云图上画出所有样点对是无法操作的。应设法将样点对按照它们之间的距离和方向进行分组，这个分组过程称为步长分组。2143图6.6 样点对的步长分组示意图在步长分组过程中将样点对按相同距离和方向进行分组，这样每一个点都具有统一的原点，这个特性使理论半变异图具有对称性。上图中，连线1和2具有非常相似的距离和方向。 步长大小的选择 当用规则格网取样时，格网间距通常可以用来确定步长大小；如果数据是通过不规则取样的，步长大小乘以步长数应等于样点间最大距离的0.5倍。 5. 空间数据变化的方向效应 如果在各个方向上Z(x)的变异性相同或相近，称Z(x)为各向同
20、性。反之，称为各向异性。 在结构分析中，半变异函数的变程a在不同方向上的大小反映各向同性或各向异性，如下图所示。C1C2 0 a1 a2 h1(h)2(h)(h)图6.7 半变异函数的各向异性曲线6.3.2 克立格插值法概述 克立格插值法是建立在半变异函数理论分析基础上的，是对有限区域内的区域化变量取值进行无偏最优估计的一种方法。 对于任意待估计点的估计值Z(x0)均可以通过待估测点范围内的n个观测样本值Z(xi)（=1, 2, , n）的线性组合得到，即 其中，入i为权重系数，其和等于1，Z(xi)为观测样本值，它们位于区域内xi位置。)()( 10iniixZxZ（6.26） 由于克立格法
21、是一种无偏最优估计，入i的确定应满足 利用拉格朗日定理，由式（6.27）和式（6.28）可推导出入i与半方差之间的矩阵方程0)()( (00 xZxZ（6.28） min)()( (200 xZxZ（6.27） BA（6.29） 其中， 由式（6.29）代入式（6.26）计算内插估计值Z(x0)0111111212222111211nnnnnnrrrrrrrrrA102010nrrrBn216.3.3 常见克立格模型 1. 普通克立格模型 当区域化变量Z(x)的数学期望EZ(x)=m为未知常数时，常采用普通克立格法进行局部估计。普通克立格模型为 在运用普通克立格法进行局部估计时，设待估块段为V
22、，中心为x，其平均值为ZV，则)()(ssZ（6.31） VVxxZVZd)(1VVmxxZVZEd)(1（6.32） （6.33） 在待估块段V的领域内，存在一组n个已知样点xi（i=1, 2, , n），其观测值为Z(xi)，其数学期望也为m。令ZV#为ZV的线性估计量，由n个已知的样点观测值Z(xi)构成的线性组合，即niiiVxZZ1#)(（6.34） 在满足下面两个条件时，ZV#为ZV的线性无偏、最优估计量 （1）无偏性条件 当 时， ，ZV#为ZV 的无偏估计量。 （2）最优性条件 在满足无偏性条件下，估计方差 为 在无偏性条件下，使估计方差最小，则ZV#为 ZV的无偏、最优估计量
23、。 nii11mZEZEVV#2E212#2 )(niiiVVVExZZEZZE（6.35） 2. 简单克立格模型 简单克立格插值模型可以表示为 简单克立格法可以使用半变异函数或协方差函数进行分析，可进行变换和剔除趋势，也可进行测量误差分析。)()(ssZ（6.36） 3. 泛克立格模型 泛克立格法是在漂移的形式EZ(x)=m(x)和非平稳随机函数Z(x)的协方差已知的情况下，一种考虑到有漂移的无偏线性估计量的空间统计方法，其模型可以表示为)()()(sssZ（6.37） 4. 指示克立格模型 指示克立格法的模型可表示为 其中，I(s)是一个二进制变量。应用二进制变量后，指示克立格法的预测精度
24、将超过普通克立格法。)()(ssI（6.38） 5. 析取克立格模型 析取克立格法的模型表达为 将指示克立格法的指示函数进行一般化处理便得到析取克立格法的指示函数表达式 )()(1ssZf（6.39） niiisZfsZg10)()(（6.41） 6. 协同克立格模型 普通协同克立格法的模型下式所示 协同克立格法应用过程中引用了协同变量，以求预测的结果更好。)()()()(222111ssZssZ（6.42） 6.3.4 克立格模型应用条件XYYXYXs22)(（6.44） 这是一个二阶多项式趋势面方程，由空间坐标(x, y)经线性回归分析获得。如果趋势方程中的回归系数是未知的，便形成了泛克立
25、格模型；如果在任何时候趋势是已知的，会形成简单克立格模型；基于多个变量的克立格模型便形成了协同克立格模型；如果在协同克立格模型中使用的是未经任何变换的Z(s)，便形成了概率克立格模型。 确定性插值法确定性插值法是使用数学函数进行插值，以研究区域内部的相似性或者以平滑度为基础，由已知样点来创建预测表面的插值方法。确定性插值法分为全局性插值法和局部性插值法，又分为精确性插值方法和非精确性插值方法。6.4.1 反距离加权插值法采样点 图6.15 用反距离加权法内插的表面反距离加权插值法是基于相近相似原理反距离加权法使用区域内已知的样点值来预测除样点外的任何位置的值 反距离加权插值法的一般公式为 确定
26、权重的计算公式为 p为参数，可以通过求均方根预测误差的最小值确定其最佳值。 iNiisZSZ10)(（6.45） Nipipiidd100/Nii11（6.46） 6.4.2 全局多项式内插法（a）（b）（c） 图6.16 全局多项式插值法全局多项式插值就像把一张纸插入到那些取值大小不同的样点之间（如上图）由采样点值拟合的全局多项式表面起伏变化平缓，它能够捕捉到数据集中潜在的粗糙数据。 全局多项式插值法适用的情况有：（1）当一个研究区域的表面变化缓慢，可以采用全局多项式插值法对该研究区进行表面插值；（2）检验长期变化的、全局性趋势的影响时一般采用全局多项式插值法6.4.3 局部多项式插值法局部
27、多项式插值法是将一个复杂的表面进行分解，并用每个小平面的中心值来预测研究区中每一点的值，从而拟合出更为准确、真实表面的一种插值方法。局部多项式插值法适于用特定的多项式方程对指定的相邻区域内的所有点进行插值当数据集中含有短程变异时，局部多项式插值表面则能更好地描述这些短程变异。6.4.4 径向基函数插值法 径向基函数法 径向基函数包括五种不同的基本函数：径向基函数法是人工神经网络方法中的一种。由径向基函数生成的表面不仅能够反映整体变化趋势而且可以反映局部变化。平面样条函数、张力样条函数、规则样条函数、高次曲面函数和反高次曲面样条函数采样点 图6.18 用径向基函数法内插的表面径向基函数法就如同将
28、一个橡胶膜插入并经过各个已知样点，同时又使表面的总曲率最小，如上图。径向基函数适用于对大量点数据进行插值计算从而获得平滑表面 探索性空间数据分析(ESDA)ESDA是指利用统计学原理和图形图表相结合对空间信息的性质进行分析、鉴别，用以引导确定性模型的结构和解法。ESDA与EDA区别在于它考虑了数据的空间特性，在方法上它将数据分解为一般趋势和叠加于其上的局部变化两部分。然后用一定的数学函数去拟合由样本点产生的经验变率函数，进行诸如克立格内插等空间操作。6.5.1 探索性空间数据分析的基本理论 1. 基本思想探索性数据分析的基本思想是： 让数据说话，即先分析数据再建立模型；不局限于方法的理论根据，
29、以一种比较松散的、非正式的方式分析数据。 探索性数据分析探索性数据分析的整个操作步骤大体可以划分为两大阶段：探索阶段和证实阶段探索性数据分析提供了各种详细考察一组数据的方法，证实性数据分析估计观察到的模式或效应的再现性。整个探索性数据分析的过程有四个主题，即耐抗性、残差、数据转换以及图示。 （1）耐抗性 耐抗性即对数据的不良表现不敏感 （2）残差 残差是从原始数据中减去概括性统计量后所剩余的部分，其公式为：残差=原数据-拟合值。 （3）数据转换 （4）图示数据转换涉及到用什么样的尺度能够帮助简化对该数据的分析探索性数据分析中最常用到的数据转换方法来自一个被称为“指数变换”的函数族通过显示图示满
30、足分析者观察数据、拟合值、诊断指标及残差表现等的需要，从而揭示数据意料之外的特性及表现。 2. 基本内容 探索性空间数据分析的内容包括以下几个方面： （1）检查数据是否有错误 （2）获得数据的分布特征 （3）对数据规律的初步考察6.5.2 探索性空间数据分析的数学方法 1. 直方图直方图是一种适用于对大量样点数据进行整理加工，找出其统计规律，以便对其总体分布特征进行推断的方法。直方图方法中两个重要的参数：频率分布和概括性的统计指标 （1）位置指标 （2）分布指标 （3）形状指标 位置指标提供该分布中心及其他部分的位置信息均值周围点的分布是描述频率分布的另一个特征，数据的方差是所有观测值与均值的
31、平方离差的均值。偏斜系数用来描述分布的对称性；峰度是基于某种分布曲线的拖尾规模，用来描述这种分布产生离群值可能性大小的指标。 2. Q-Q概率图Q-Q概率图是根据变量分布的分位数对所指定的理论分布分位数绘制的图形，它是一种用来检验样点数据分布的统计图。Q-Q概率图分为两种类型：正态概率图（图6.21）和反趋势正态概率图（图6.22）100 200 300 400 500 600 700 800 900 100200300400500600700800观测值 图6.21 正态概率图正态分布的期望值100 200 300 400 500 600 700 800 900 20406080100观测值
32、0-20 图6.22 反趋势正态概率图与正态直线的偏差 3. 趋势分析 在进行趋势分析时，将三维视图（如图6.23所示）中Z轴数据值分别投影到X、Z平面和Y、Z平面作散点图。投影后的曲线是平直的，表明没有趋势；呈上升趋势的曲线模式，则表明数据中存在全局趋势。X轴Y轴Z轴 图6.23 趋势分布图 4. 半变异/协方差函数云图 半变异/协方差函数表示的是数据集中所有样点对的理论半变异值和协方差，把它们用两点间距离函数来表示，并以此函数作图，称半变异协方差函数云图（如图6.24）。(h)0 0.23 0.46 0.69 0.92 1.15 1.38 1.61 1.84h图6.24 典型的半变异函数云
33、图1.601.280.64 5. 正交协方差函数云图正交协方差函数表示的是两个数据集中所有样点对的理论正交协方差，把它们用两点间距离的函数来表示，并以此函数作图，称正交协方差函数云图。正交协方差函数云图可用来检验两个数据集的空间关联性的局部特征，并可查找两个数据集间关联性的空间变化情况。6.5.3 探索性空间数据分析的应用 探索性空间数据分析探索性空间数据分析是一种图形化的数据分析方式通过它可以完成诸如观察数据的分布、寻找离群值、进行全局趋势分析以及检测空间的自相关和方向变异等任务 1. 检验数据分布 2. 寻找离群值 3. 全局趋势分析直方图中通过观察均值和中值能确定分布的中心寻找离群值的方法有，用直方图查找离群值以及用半变异协方差函数云图识别离群值。预测表面主要由两部分组成，即确定的全局趋势和随机的短程变异。假设全局趋势是确定的，而短程变异是随机的。最终的预测表面是确定表面和随机表面的总和。LOGO
