南阳师院第5章-第2节地统计分析方法.ppt

第2节地统计分析方法 一 地统计方法的基本原理 一 区域化变量 二 协方差函数 三 变异函数 四 克立格插值方法二 应用实例 地统计学是以区域化变量理论为基础 以变异函数为主要工具 研究那些在空间分布上既有随机性又有结构性 或空间相关和依赖性的自然现象的科学 协方差函数和变异函数是以区域化变量理论为基础建立起来的地统计学的两个最基本的函数 地统计学的主要方法之一 克立格法就是建立在变异函数理论和结构分析基础之上的 当一个变量呈现为空间分布时 就称之为区域化变量 RegionalizedVariable 这种变量常常反映某种空间现象的特征 用区域化变量来描述的现象称之为区域化现象 区域化变量 亦称区域化随机变量 G Matheron 1963 将它定义为以空间点x的三个直角坐标为自变量的随机场 区域化变量具有两个最显著 而且也是最重要的特征 即随机性和结构性 一 地统计方法的基本原理 一 区域化变量 二 协方差函数 1 协方差函数的概念区域化随机变量之间的差异 可以用空间协方差来表示 区域化变量在空间点x和x h处的两个随机变量和的二阶混合中心矩定义为Z x 的自协方差函数 即 协方差函数的计算公式为 式中 h为两样本点空间分隔距离或距离滞后 为在空间位置处的实测值 是在处距离偏离的实测值 i 1 2 是分隔距离为h时的样本点对 Paris 总数 和分别为和的样本平均数 2 协方差函数的计算公式 若 m 常数 则上式可以改写为 式中 m为样本平均数 可由一般算术平均数公式求得 即 三 变异函数 1 变异函数的概念变异函数 Variograms 又称变差函数 变异矩 是地统计分析所特有的基本工具 在一维条件下变异函数定义为 当空间点x在一维x轴上变化时 区域化变量Z x 在点x和x h处的值Z x 与Z x h 差的方差的一半为区域化变量Z x 在x轴方向上的变异函数 记为 h 即 在二阶平稳假设条件下 对任意的h有因此 公式可以改写为从上式可知 变异函数依赖于两个自变量x和h 当变异函数仅仅依赖于距离h而与位置x无关时 可改写成 即 2 变异函数的性质设Z x 是区域化变量 在满足二阶平稳假设条件下 变异函数式具有如下性质 1 0 即在h 0处 变异函数为0 2 即关于直线h 0是对称的 它是一个偶函数 3 0 即只能大于或等于0 4 h 时 c 0 或 c 0 即当空间距离增大时 变异函数接近先验方差 5 必须是一个条件非负定函数 由 构成的变异函数矩阵在条件时 为非负定的 3 变异函数的计算公式 设是系统某属性Z在空间位置x处的值 为一区域化随机变量 并满足二阶平稳假设 h为两样本点空间分隔距离 和分别是区域化变量在空间位置和处的实测值 i 1 2 N h 那么 变异函数的离散计算公式为 这样对不同的空间分隔距离h 计算出相应的和值 如果分别以h为横坐标 或为纵坐标 画出协方差函数和变异函数曲线图 就可以直接展示区域化变量Z x 的空间变异特点 可见 变异函数能同时描述区域化变量的随机性和结构性 从而在数学上对区域化变量进行严格分析 是空间变异规律分析和空间结构分析的有效工具 例如 假设某地区降水量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缺失值情况下样本数对的组成和计算过程 为缺失值 同样计算出最后 得到南北方向和西北 东南上的变异函数计算结果见下表 同样可以计算东西方向上的变异函数 4 变异函数的参数 变异函数有四个非常重要的参数 即基台值 Sill 变程 Range 或称空间依赖范围 RangeofSpatialDependence 块金值 Nugget 或称区域不连续性值 LocalizedDiscontinuity 和分维数 FractalDimension 前3个参数可以直接从变异函数图中得到 它们决定变异函数的形状与结构 由于数据对 Wz z 经过了标准化 因此界外值可易由2 sigma规则可视化地识别出来 变异函数的形状反映自然现象空间分布结构或空间相关的类型 同时还能给出这种空间相关的范围 当变异函数随着间隔距离h的增大 从非零值达到一个相对稳定的常数时 该常数称为基台值C0 C 当间隔距离h 0时 0 C0 该值称为块金值或块金方差 NuggetVariance 基台值是系统或系统属性中最大的变异 变异函数达到基台值时的间隔距离a称为变程 变程表示在h a以后 区域化变量Z x 空间相关性消失 块金值表示区域化变量在小于抽样尺度时非连续变异 由区域化变量的属性或测量误差决定 C0 C0 C2 h 上述三个参数可从变异函数曲线图直接得到 或通过估计曲线回归参数得到 第4个参数 即分维数用于表示变异函数的特性 由变异函数和间隔距离h之间的关系确定 分维数D为双对数直线回归方程中的斜率 它是一个无量纲数 分维数D的大小 表示变异函数曲线的曲率 可以作为随机变异的量度 但该随机分维数D与形状分维数有本质的不同 5 变异函数的理论模型 地统计学将变异函数理论模型分为三大类 第一类是有基台值模型 包括球状模型 指数模型 高斯模型 线性有基台值模型和纯块金效应模型 第二类是无基台值模型 包括幂函数模型 线性无基台值模型 抛物线模型 第三类是孔穴效应模型 下面有代表性地介绍几种常见的变异函数理论模型 1 纯块金效应模型 其一般公式为 式中 c0 0 为先验方差 该模型相当于区域化变量为随机分布 样本点间的协方差函数对于所有距离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的关系见下表 下面我们试用回归分析方法建立其球状变异函数模型 从上面的介绍和讨论 我们知道 球状变异函数的一般形式为 当时 有 如果记 则可以得到线性模型 4 2 19 根据表中的数据 对上式进行最小二乘拟合 得到 4 2 20 计算可知 上式的显著性检验参数F 114 054 R2 0 962 可见模型的拟合效果是很好的 比较 4 2 20 式与 4 2 19 式 并做简单计算可知 c0 2 048 c 1 154 a 8 353 所以 球状变异函数模型为 4 2 21 四 克立格插值方法 克立格 Kriging 插值法 又称空间局部估计或空间局部插值法 是地统计学的主要内容之一 克立格法是建立在变异函数理论及结构分析基础之上的 它是在有限区域内对区域化变量的取值进行无偏最优估计的一种方法克立格法适用的条件是 如果变异函数和相关分析的结果表明区域化变量存在空间相关性其实质是利用区域化变量的原始数据和变异函数的结构特点 对未采样点的区域化变量的取值进行线性无偏 最优估计 克里格插值 KrigingInterpolation 是根据变异函数模型而发展起来的一系列地统计的空间插值方法 包括 普通克里格法 OrdinaryKriging 泛克里格法 UniversalKriging 指示克里格法 IndicatorKriging 析取克里格法 DisjunctiveKriging 协同克里格法 Cokriging 等 下面仅对普通克立格法作一些简单介绍 首先假设区域化变量满足二阶平稳假设和本征假设 其数学期望为m 协方差函数及变异函数存在 即 假设在待估计点 x 的临域内共有n个实测点 即 其样本值为 那么 普通克里格法的插值公式为 其中为权重系数 表示各空间样本点处的观测值对估计值的贡献程度 可见 克立格插值的关键就就是计算权重系数 显然 权重系数的求取必须满足两个条件 一是使的估计是无偏的 即偏差的数学期望为零 二是最优的 即使估计值和实际值之差的平方和最小 为此 需要满足以下两个条件 1 无偏性 要使成为的无偏估计量 即 当时 也就是当时 则有这时 为的无偏估计量 2 最优性 在满足无偏性条件下 估计方差为 使用协方差函数表达 它可以进一步写为 4 2 24 为使估计方差最小 根据拉格朗日乘数原理 令 4 2 25 求F对和的偏导数 并令其为0 得克立格方程组 4 2 26 整理后得 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 37 根据协方差与变异函数的关系以及 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个气象台站多年平均降水量和蒸发量数据 见教材表3 1 2 为实测值 拟合了年降水量和蒸发量的半变异函数理论模型 并采用普通克里格法和双变量协同克里格法 做了空间插值计算 结论如下 一 半变异函数半变异函数模型 是克立格空间插值的前提条件 同时它也决定着空间插值的精度 一般情况下 半变异函数模型是根据半变异函数云图的分布 选择合适的理论模型 按照估计方差最小的原则 运用最小二乘法求得 图4 2 4和图4 2 5分别给出了年降水量和年蒸发量的半变异函数云图 图4 2 4年降水量的半变异函数云图 同样也可以作出蒸发量的半变异函数运图 可以看出 年降水量和年蒸发量的块金效应都不明显 这是因为样本点是各个气象站点的实测值 空间分辨率可以忽略不计 另外实验误差和人为性误差基本上都很小 我们选择各种不同的半变异函数理论模型 经过多次拟合计算和对比分析 发现指数模型比较好地描述了年降水量的空间变异规律 其变异函数的具体形式如下 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可以看出 在甘肃省范围内 年降水量的空间分布格局总体上是东南多西北少 并且呈现从东南方向到西北方向逐渐过渡 梯度变化明显 山地多 平地少 南北方向从南部队祁连山脉向北部的沙漠戈壁逐渐减少 年降水量的空间变程很大 最多的东南部是最少的西北部的近10倍 其中 甘南东南部玛曲和禄曲 陇南东南部以及平凉灵台东南地区 年降水量达到691 59 786 75mm之间 400mm等降水线靠近兰