空间分析空间插值与地统计ppt课件

第八讲空间插值与地统计 一阶方法倒距离权重 IDW 插值趋势面分析二阶方法区域化变量协方差函数 变异函数 变异函数模型交叉验证地统计 克里金 Kriging 方法概述 理解不同的克里金方法几种不同的克里金方法普通克里金 OK 简单克里金 SK 泛克里金 UK 指示克里金 IK 协克里金 Ck 1 一 一阶方法 倒距离权重 IDW 插值趋势面分析 2 IDW插值方法假定每个输入点都有着局部影响 这种影响随着距离的增加而减弱 步骤 计算未知点到所有点的距离 计算每个点的权重 权重是距离倒数的函数 计算结果 1 倒距离权重 IDW 插值 1 4最大 6最小 示例 IDW插值 求图中0点的值 IDW插值的一般模型对所有或选定的i进行计算典型地 的取值是1或2除以距离权重的和 保证了权重加起来等于1 例 倒距离加权 IDW 插值结果 IDW插值的缺点 IDW不能得到大于样本最大值或小于样本最小值的估计 对于高程表面 这 抹平 了峰和谷 除非它们的高点和低点是样本的一部分 因为估计值为均值 得到的表面将不通过样本点 8 IDW插值在ArcGIS中的实现 9 IDW插值在ArcGIS中的实现 10 Idaho州降雨量等直线图 IDW插值在ArcGIS中的实现 11 2 趋势面分析 趋势面分析 trendsurfaceanalysis 用数学模型来模拟 或拟合 地理数据的空间分布及其区域性变化趋势的方法 12 确定性插值 函数拟合最普通的形式 多项式e gy ax2 bx c局部 local 分段地全局 global 13 趋势面分析是一种整体内插法 该方法假设一般趋势与局部变化无关 并利用曲面方程来模拟待定点附近地形表面的一般趋势 通常使用的是1次 2次 3次趋势面 过高次的趋势面不利于反映空间趋势 并可能存在趋势面的 畸变 其中 2次趋势面可用待定点附近的6个数据点来计算方程式系数 14 DeterministicSolutions FirstOrderPolynomialInterpolation PredictedModelMeasured SecondOrder third fourth etc PolynomialInterpolation LocalPolynomialInterpolation RadialBasisFunction Spline Interpolation 15 趋势面的性质与特点是一种光滑的数学曲面 能集中地代表地理数据在大范围内的空间分布变化趋势 与实际地理曲面不同 它只是实际曲面的一种近似 实际曲面包括趋势面和剩余曲面两部分 即 实际曲面 趋势面 剩余曲面 16 设Zi xi yi 表示某一地理特征值在空间上的分布 其中 xi yi 为平面上点的坐标 任一观测点Zi可分解为两个部分 即 趋势面 剩余面 17 一阶趋势面FirstOrderTrendSurface 一阶趋势面残差ResidualsfromFirstOrderTrendSurface 18 趋势面参数的确定 最小二乘法 使每一个观测值与趋势值的残差平方和最小 即按建立多元线性方程的方法 使Q对系数b0 b1 bn求偏导 并令这些偏导数等于零 得趋势面的正规方程组 解正规方程组 即可求出系数 从而得到趋势面方程 19 因为任何函数在一定范围内总可以用多项式来逼近 并可调整多项式的次数来满足趋势面分析的需要 一般来说 多项式的次数越高则趋势值越接近于观测值 而剩余值越小 20 多项式趋势面的数学模型 21 22 二阶趋势面SecondOrderTrendSurface 1 x y x2 y2 xy T 三阶趋势面ThirdOrderTrendSurface 1 x y x2 y2 xy x3 y3 x2y xy2 T 23 趋势面的具体计算方法与步骤 原始数据列表 等间隔选取纵横坐标网 将原始数据点入坐标 按多元线性回归分析方法求出趋势面的正规方程组 解出参数 从趋势值等值线图中 获得地理要素的区域性变化规律 用F分布对趋势面进行拟合程度检验 24 即用双三次多项式拟合趋势面 双三次多项式 样条函数 插值 25 该曲面模型有16个待定系数 Cij i j 0 1 2 3 通常用4个数据点 规则格网的4个顶点 的4个函数值组成的4 4方程组求解 如图 这4个函数值是高程Z x方向斜率R y方向斜率S 以及扭矩T 26 其中Z保证曲面通过格网的4个数据点 R S T保证曲面在这4个数据点处光滑连续 双三次多项式 样条函数 内插法是规则格网插值的常用方法之一 这种方法通过一系列曲面片段来拼接地形表面 最终得到一个1阶 2阶连续的表面 该方法属于局部插值 计算负担中等 对于平滑表面拟合效果最好 对于起伏的表面拟合效果最差 27 趋势面分析在ArcGIS中的实现 28 趋势面分析在ArcGIS中的实现 29 趋势面分析在ArcGIS中的实现 Idaho州降雨量等直线图 30 二 二阶方法 IDW插值和趋势面方法的缺陷 IDW 距离权重函数的选择和 邻居 的定义是其 致命缺陷 Achilles heel 趋势面分析 Inasense trendsurfaceanalysisletsthedataspeakforthemselves whereasIDWinterpolationforcesasetstructureontothem 31 二 二阶方法 区域化变量协方差函数 半变异函数 32 地统计学是以区域化变量理论为基础 以变异函数为主要工具 研究那些在空间分布上既有随机性又有结构性 或空间相关和依赖性现象的学科 协方差函数和变异函数是以区域化变量理论为基础建立起来的地统计学的两个最基本函数 地统计学的主要方法之一 克里金方法 Kriging 就是建立在变异函数理论和结构分析基础之上的 33 当一个变量的取值与其空间位置有关时 就称为区域化变量 regionalizedvariable 区域化变量常常反映某种空间现象的特征 用它来描述的现象称之为区域化现象 区域化变量 亦称区域化随机变量 Matheron 1963 将它定义为以空间点x的三个直角坐标为自变量的随机场区域化变量具有两个最显著 也是最重要的特征 随机性和结构性 34 随机变量 随机函数 随机过程 随机场 区域化变量 与时间有关的随机函数 带有多个 2个以上 自变量的随机函数 以空间点的三个直角坐标为自变量 35 http cg ensmp fr Presentation Matheron Matheron en shtml ProfessorGeorgesMatheron 1930 2000 8 7 法国数学家和地质学家 36 区域化变量的功能 由于区域化变量是一种随机函数 因而能同时反映空间变量的结构性和随机性 一方面 当空间点x固定后 Z x 就是一个随机变量 这体现了其随机性 另一方面 在空间两个不同点x与x h处的区域化变量值具有某种程度的相关性 这体现了其结构性 37 区域化变量的组成部分 数据点 结构性可以用均值和常数趋势表示空间相关数据通常呈现正空间相关性随机性测量误差 其他误差 38 distance elevation 结构性 随机性 实际值 39 协方差函数与变异函数 数学期望 方差和协方差数学期望 一阶原点矩方差 二阶中心矩协方差 二阶混合中心矩 40 协方差函数类似地 当Z x 是区域化变量时 对于任意两点si和sj 空间随机过程的协方差函数为 相关系数和方差分别定义为 41 42 若 则过程是二阶平稳的 即均值与方差独立于空间位置并在研究区域上是常数 于是有 称为协方差图或过程的协方差函数 称为相关图或相关函数 显然 协方差函数仅依赖于向量差h 当h 0时 43 若独立性仅是距离的函数 与方向无关 则空间过程是各向同性 isotropy 的 协方差函数就只依赖于距离向量h 若在给定距离和方向上 不同位置数值差异的均值和方差为常数 则下式成立 半变异函数 也称半方差函数 semi variogram 44 半变异函数 semi variogram 区域化变量的基本研究工具 半变异函数就是区域化变量增量平方的数学期望之半 区域化变量在i i h点的值 步长为h的样品对数 步长 h 在一定方向上 距离为h的矢量 45 方差 Variance 变异函数 Variogram 对比 方差 变异函数 46 h a Vetordistanceh 建立经验半变异函数 Semivariogram 测度空间变异 对于间隔距离为h的每一对Z x 和Z x h 测度它们之间差的平方 47 半变异云图 半变异云图 数据集中所有点对的 半变异函数 距离图 如果空间相关性存在 互相靠近的已知点趋向于具有较小的半方差互相远离的已知点趋向于具有较大的半方差 48 半变异云图 semivariogramcloud 49 分组 Binning 如果数据很大 样点对的数目将迅速增加并且变得难以操作 因此 可以将样点对分组 也就是步长 h 分组 如 可将距离在0 1之间的作为第一组 h在1 2之间的作为第二组 依此类推 Binningisaprocessthataveragessemivariancedatabydistanceanddirection hey itsweighted 50 Binnedsemivariogram 51 模型拟合 半变异函数必须拟合一个数学函数或模型 以用于估计任何给定距离的半方差 52 理论变异函数模型 实践中 常用的是变异函数图 notrelatedanymore 变程范围内才有结构性变化 有规律的变化 反映随机性大小 主要来源于区域化变量Z x 在小于抽样尺度h时所具有的内部变异 另外还有抽样分析误差 变异函数是一个单调不减函数 当h超过某一个范围 例如变程 变异函数不再增大 而是趋于一个极限值 即为基台值 实际上等于区域化变量的先验方差 即 即基台值与块金值之差 表示数据中存在空间相关性引起的方差变化范围 53 Thesillisthevalueatwhichthesemivariogramlevelsoff itsasymptoticvalue Therangeisthedistanceatwhichthesemivariogramlevelsoff thespatialextentofstructureinthedata Thenuggetisthesemivarianceatadistance0 0 they intercept Asemivariogramisaplotofthestructurefunctionthat likeautocorrelation describestherelationshipbetweenmeasurementstakensomedistanceapart Semivariogramsdefinetherangeordistanceoverwhichspatialdependenceexists 54 55 ArcGIS GeostatisticalAnalystprovidesthefollowingfunctionstochoosefromtomodeltheempiricalsemivariogram CircularSphericalTetrasphericalPentasphericalExponentialGaussianRationalQuadraticHoleEffectK BesselJ BesselStable 理论半方差模型的类型 56 例 假设某地区降水量Z x 单位 mm 是二维区域化随机变量 满足二阶平稳假设 其观测值的空间正方形网格数据如图1所示 点与点之间的距离为h 1km 试计算其南北方向及西北和东南方向的变异函数 57 图1空间正方形网格数据 点间距h 100m 从图1可以看出 空间上有些点 由于某种原因没有采集到 如果没有缺失值 可直接对正方形网格数据结构计算变异函数 在有缺失值的情况下 也可以计算变异函数 只要 跳过 缺失点位置即可 图2 58 首先计算南北方向上的变异函数值 由变异函数的计算公式可得 385 72 5 35 图2缺失值情况下样本数对的组成和计算过程 为缺失值 59 60 同样计算出最后 得到南北方向和西北 东南方向上的变异函数 下表 同样可以计算东西方向上的变异函数 61 semivariogram南北向 semivariogram西北 东南向 62 3 交叉验证 CrossValidation 对每种插值方法重复下面的步骤 实现对不同插值方法的比较 从数据集中除去一个已知点的测量值 用剩余的点估计除去点的值 比较原始值和估计值 计算出估计值的预测误差 针对每个已知点 进行上述步骤 然后评价不同插值方法的精确度 常用的评价指标是均方根 RMS 63 交叉验证 计算每个点的RMS 计算某种插值方法的平均RMS 1 2 3 选择某种插值方法 4 5 64 Kriging内插和外推背后的理论由法国数学家GeorgesMatheron基于DanielG Krige的硕士论文而发展 Krige是距离加权平均的开创者 用于当时南非的金矿勘探 TheEnglishverbistokrigeandthemostcommonnouniskriging 三 克里金插值方法 克里金 kriging 插值 根据随机场中待测值邻居位置的观测值对待测值进行插值的一组地统计技术 e g 把高程z作为地理位置的一个函数 65 克里金 Kriging 插值法 又称空间局部估计或空间局部插值法 是地统计学的主要内容之一 克里金法建立在变异函数理论及结构分析基础之上 它是在有限区域内对区域化变量进行无偏最优估计的一种方法 其实质是利用区域化变量的原始数据和变异函数的结构特点 对未采样点的区域化变量值进行线性无偏最优估计 克里金法的适用条件 如果变异函数和相关分析的结果表明区域化变量存在空间相关性 66 假设 所有的随机误差都具有二阶平稳性 即 随机误差的均值为零 并且任意两个随机误差之间的协方差只与两者之间的距离和方向有关 而与它们的具体位置无关 克里金 Kriging 方法 67 克里金方法的基本形式 确定趋势项 空间相关随机误差项 理解不同的克里金模型 对误差项的假设 期望值为0 并且和之间的自相关不取决于s点的位置 而取决于位移量h 为确保自相关方差有解 必须允许某两点间的自相关可以相等 如 下面有箭头相连的两对位置点假设具有相同的自相关性 68 趋势值可以被简单地赋予一个常量 即 在任何位置处 理解不同的克里金模型 如果在任何时候趋势已知 无论趋势是否是常量 都形成简单克里金模型 趋势也可以表示为 若趋势中的系数未知 就是泛克里金模型 如果未知 就是普通克里金模型 69 理解不同的克里金模型 可以对Z s 进行一下变换 例如 可以把它换为一个指示变量 若Z s 低于一定的阈值 如空气中的溴氧浓度值0 12ppm 则变为0 若Z s 高于一定的阈值 则变为1 然后对高于阈值的情况进行预测 便形成了指示克里金模型 现在看一下方程的左边 70 普通克里金 OrdinaryKriging OK 简单克里金 SimpleKriging SK 泛克里金 UniversalKriging UK 指示克里金 IndicatorKriging IK 71 协克里金 Cokriging CK 假设和是未知常数 感兴趣的变量是Z1 同时利用Z1的自相关以及与Z2之间的交叉相关 72 TypesofKriging Simplekrigingisoptimalestimationofarandomfield e g T x withaknownmean m x andaknowncovariancePTT x x Ordinarykrigingisoptimalestimationofarandomfield e g T x withanunknownconstantorlinearlytrendingmean butaknownsemivariogramgTT x x Universalkrigingisoptimalestimationofarandomfield e g T x withanunknownpolynomialtrendingmean butaknownsemivariogramgTT x x 73 运用克里金 Kriging 方法进行插值的过程可以分为两步 第一步 变异函数分析 即 样点的空间结构量化分析 是指对样点数据拟合一个空间独立模型 第二步 对未知点值进行预测 利用第一步拟合的变异函数 样点数据的空间分布以及样点数据值对未知点进行预测 74 普通克里格 OK 模型为 其中 s x y 为空间位置 如 观测点 1 5 Z 1 5 100 未知均值 75 首先假设区域化变量满足二阶平稳假设 其数学期望为m 协方差函数及变异函数存在 即 假设在待估计点 x 的邻域内共有n个实测点 即x1 x2 xn 其样本值为Z xi 那么 普通克里格 OK 的插值公式为 第i个观测点的权重 预测点 观测点 76 这个简单的模型类似于倒距离权重法 IDW 不同的是IDW的权重只与距预测点的距离有关 而普通克里格 OK 方法的权重则取决于半变异图 距预测点的距离和预测点周围观测值的空间关系 77 最优 Best 要求尽可能小 从而预测值将尽可能的接近未知值线性 Linear weightedlinearcombinationsofthedata无偏性 Unbiased 一些点的预测值比真值高 一些比真值低 平均起来预测值和真值的差为0 要求权重系数值和为1 Estimation 普通克里金 OK 的目标 BLUE 78 在无偏估计的限制条件下 通过最小化可得到克里金方程 包含所有样点对的半变异值 基于点i和j之间距离的半变异值 包含每一个观测点与预测点之间的半变异值 基于第i个观测点与预测点之间距离的半变异值 待定 79 普通克里金 OK 的计算过程 计算样点对之间的距离 方差 理论半变异值 平均理论半变异值拟合模型计算权系数与预测克里金方差 80 计算样点对之间的距离和方差 半变异值 0 5 方差 示例 普通克里金方法 81 如果数据很大 样点对的数目将迅速增加并且变得难以操作 因此 可以将样点对分组 也就是将步长 h 分组 在该例中 首先将距离在1 2之间的作为第一组 h在2 3之间的作为第二组 依此类推 下表 82 拟合模型 为了预测 需要用理论半方差模型拟合经验半方差 以描述其变程 基台值和块金 range sill nugget 四个常用模型 83 Gaussian Linear Spherical Exponential Where c0 nuggetb regressionslopea rangec0 c sill Assumesnosillorrange 84 拟合模型用平均半变异值为纵坐标 步长 抽样间距 为横坐标 做理论半变异图 半变异值 斜率 距离 13 5 h 生成伽玛矩阵 如 样点 1 5 与样点 3 4 的半变异值为 13 5 2 236 30 19 85 伽玛矩阵 最后一行的1和0根据无偏估计的限制条件求得 伽玛矩阵的逆矩阵 86 计算权系数与预测普通克里金的权系数矩阵为 向量g由未知点生成 如用点 1 4 来计算 计算该点与每个观测点的距离 如 1 4 与 1 5 3 4 1 3 4 5 5 1 的距离 预测点 1 4 的g向量计算结果如下表 半变异值 13 5 h 87 计算权系数 正如预期 权系数随距离的增加而减小 但由于将各点的空间分布也考虑进去 其结果比直接的距离权重要好 88 权系数为负值的处理 DeutschClaytonV Correctingfornegativeweightsinordinarykriging J Computers Geosciences 1996 22 7 765 773 89 克里格方差 为了估计预测结果的不确定性 将权系数向量的每一行和g向量的每一行相乘 然后将结果相加 得到预测克里格方差 其平方根就是克里格标准差 如果假设误差是正态分布的 则95 的置信区间为 克里格预测值 1 96 克里格标准差 95 49 109 75 90 这仅是一个小例子 但从中可以看出几个重要特点 Kriging是计算密集型的方法 需要一个合适的软件 虽然一些GIS软件提供半方差估计 建模和Kriging 该邻域大多严肃的工作者使用特殊软件如GSLIB Variowin 或GS 所有的结果依赖于为从样点数据估计的半方差所拟合的模型 以及相关假设 其中包含一些主观决定 多少个距离组合 拟合的模型是什么 基台值和块金值取多少 除了所讨论的普通克里金 还有其他不同类型的克里金方法 比如 当存在均值漂移 drift 时采用UniversalKriging CoKriging扩展到同时考虑两个或更多变量 www ai geostats orgwww gstat org 91 VariogramModelParameters Wenowlookathowparametersofavariogram covariance modelaffectthe OrdinaryKriging OK weightsSill shape nugget range andanisotropy 92 TheEffectofSill Withanyrescalingofthevariogram neithertheKrigingweightsnortheestimatearechangedwhilethevarianceincreasesbythesamefactorusedtoscalethevariogram 93 Sill Sillof10vs 20 Sill 10 Sill 20 94 TheEffectofShape 高斯变异函数模型分配给近样本点更大的权重 屏幕效应 screeneffect 一个样本位于另一个离未知点更近的样本后面 具有小 或负的 权重 如样本5vs 6 高斯模型比指数模型具有更强的屏幕效应 Weightsthatarelessthan0orgreaterthan1canproduceestimateslargerthanthelargestsamplevalueorsmallerthanthesmallest Weightswithin 0 1 produceestimatesonlywithintheminandmaxofsamplevaluesNegativeweightsmayproducenegativeestimates althoughinmostscienceapplicationsvaluesshouldbepositive 95 Shape Exponentialvs Gaussianmodel Exponential Gaussian 96 TheNuggetEffect 块金效应 nuggeteffect 使权重变得更相似 并导致更大的kriging方差 纯块金效应模型表明完全缺乏空间相关性 或在小于最小抽样间距的尺度上才具有空间依赖性 97 Nugget Nugget 0vs 1 2sill 98 TheEffectofRange 变程减小 kriging方差增大 如果变程太小 那么所有样本点离待估点的距离呈现相同远近 于是估计类似于简单平均 权重为1 n Rangeofhvs 0 5h Range 10 Range 20 99 EffectofAnisotropy MoreweightsaregiventosampleslieinthedirectionofmaximumcontinuityWeightsgiventothesamplesinthemaximumspatialcontinuitywouldincreaseastheanisotropyratiobecomeslarger 100 Anisotropy Directionalvariogramsandcovariancefunctions 101 Kriging在ArcGIS中的实现 探索性数据分析 102 Kriging在ArcGIS中的实现 探索性数据分析 103 Kriging在ArcGIS中的实现 探索性数据分析 104 Kriging在ArcGIS中的实现 探索性数据分析 105 Kriging在ArcGIS中的实现 106 Kriging在ArcGIS中的实现 107 Kriging在ArcGIS中的实现 108 Kriging在ArcGIS中的实现 109 Idaho州降雨量等直线图 Kriging在ArcGIS中的实现 110 误差图 Kriging在ArcGIS中的实现 111 Kriging IDW插值 趋势面 5次多项式 112 IDW插值 普通克里格插值 113 地形插值比较 Voronoi图TINIDWKriging样条规则样条 114 SimpleKriging OdinaryKriging UniversalKriging 115 OdinaryKriging 误差图 UniversalKriging SimpleKriging SimpleKriging OdinaryKriging UniversalKriging 116 ArcGIS 普通克里格的四张图 117 118 参考资料 O SullivanD andUnwinD GeographicInformationAnalysis M JohnWiley Sons Inc 2003 p 45 49 p 265 281 刘湘南等 编著 GIS空间分析原理与方法 第二版 M 科学出版社 2008 pp 195 209 王远飞等编著 空间数据分析方法 M 科学出版社 2007 pp 152 175 ESRI UsingArcGIS GeostatisticalAnalysist 翻译本 生态学空间分析原理与技术 M 科学出版社 2008 张仁铎 空间变异理论与应用 M 科学出版社 2005 1 deSmitM J etal 著 杜培军等译 地理空间分析 原理 技术与软件工具 M 电子工业出版社 2008 119 思考题 IDW插值的基本思想是什么 以普通Kriging插值为例 阐述Kriging插值的基本思想 120 Demo packagegstat http www r project org 121 Meusedataset 155samplestakenonasupportof10 x10mfromthetop0 20cmofalluvialsoilsina5x2kmpartthefloodplainoftheMaas Meuse nearStein NL idpointnumberx ycoordinatesEandN inmcadmiumconcentrationinthesoil inmgkg 1copperconcentrationinthesoil inmgkg 1leadconcentrationinthesoil inmgkg 1zincconcentrationinthesoil inmgkg 1elevelevationabovelocalreferencelevel inmomorganicmatterlossonignition inpercentffreqfloodfrequencyclass 1 annual 2 2 5years 3 every5yearssoilsoilclass codedlimehasthelandherebeenlimed 0or1 ForTlanduselanduse codeddist mdistancefrommainRiverMaaschannel inm 122 library gstat data meuse summary meuse data meuse grid g gstat formula log zinc 1 locat