普通克里金法.pdf

天津城市建设学院学报第 1 4 卷第 1 期2 0 0 8 年3 月 J o u rnal o f T i a n j i n I n s t i t u t e o f Ur b a n C o n s t r u c ti o n V o 1 1 4 No 1 Ma r 2 0 0 8 克里金插值算法在等高线绘制中的应用 李莉 ，胡建平 b (天津城市建设学院 a 土木工程系；b 电子与信息工程系，天津 3 0 0 3 8 4)摘要：介绍了普通克里金插值方法的原理及步骤，利用普通克里金插值法对离散的数字高程 模型(DE M)数据进行了插值加密；在 Wi n d o ws环境下，以 V i s u a l C+6 0为开发平台，以 O p e n G L为工具，建立了基于 DE M 的等高线模型；对克里金插值法与距离反比插值法做 了比 较，结果表明，克里金插值后生成的等高线模型比距离反比插值后生成的等高线模型精度高、误 差 小 关键词：数字高程模型；等高线；克里金插值法；距离反比 中图分类号：T P 3 0 1 6 文献标识码：A 文章编号：1 0 0 6 6 8 5 3(2 0 0 8)0 1 0 0 6 8 0 4 等高线是表示地形最常见的形式，可以从总体上 把握研究对象的空间变化特征，在水文、环境、气象、规划等各个领域有着广泛的应用 利用数字高程模型(D E M)数据可以绘制等高线 要取得 D E M 数据地 图，首先必须在地平面上选取有限的离散点，然后通 过航天或者航空遥感测量技术，取得这些点的高程 值，从而利用这些高程值近似地表示地面形状 对于 那些在 D E M 地形图中没有记录高程值的点，则需要 利用其附近点的高程值和适当的空间内插值函数来 计算其高程值【l】D E M 内插就是根据若干相邻参考 点的高程，求出待定点上的高程值，在数学上属于插 值问题 任何内插方法都是基于原始地层起伏变化的 连续光滑性，或者说临近数据点间有很大的相关性，才可能由临近的数据点内插出待定点的高程】考虑 到被描述对象的随机性和一定范围内呈现出的某种 程度的相关性，克里金插值(K r i g i n g)是合适的内 插方法，而且它还可以求出插值误差，有利于了解插 值结果的可靠程度 1 克里金算法原理 1 1 克里金算法的基本原理 克里金插值法，又称空间局部估计或空间局部插 值法，是地统计学的主要内容之一 克里金插值法在 空间相关范围分析的基础上，用相关范围内的采样点 收稿 日期：2 0 0 7 1 1-0 l：修订日期：2 0 0 7 1 2 2 8 作者简介：李莉(1 9 8 3 一)，女，天津人，天津城市建设学院硕士生 来估计待插点属性值 克里金法是建立在变异函数 理论及结构分析基础之上的，它在有限区域 内对区域 化变量的取值进行无偏、最优估计 变异函数和相关 分析的结果表明，区域化变量存在空间相关性，其实 质是利用区域化变量的原始数据和变异函数的结构 特点，对未采样点的区域化变量的取值进行线性无 偏、最优估计 与其它的插值方法相比，克里金法的 显著特点是使误差的方差最小 克里金插值法主要包 括：简单克里金、普通克里金、泛克里金、指示克里 金、析取克里金和协同克里金等 普通克里金法是一 种较为常用的方法 假设区域化变量Z(x 1 满足二阶 平稳假设和本征 假设，其数学期望为m，协方差函 数c(h 1 及变异函 数 r(h 1 存在 即 E E Z(x)=m c(h)=E z()z(+|I1)一 m y(h)=z()一 z(+)二 设Z(x)是一个二阶平稳的随机函数，它在，1 个 位置取样：z()，z()，z()，点X o 处的 估计 量为 z ()=z()(1)i=1 其中：为权重系数，表示各空间样本点 处的观 维普资讯 湟 壹堡 堂 堂拯 李 莉等：克里 金插值算法在等高 线绘制中的 应用 测值z()对估计值Z (X)的贡献程度 克里金算法 的关键就是计算权重系数，权重系数的计算必须满 足两个条件：使Z (X)是 z()的无偏估计量，即E Z ()=E z()j 当E z ()=，z 时，即 厂 n n n E I z()I=E z()=，z 时，则 有 =1；L i=1 J i=1 i=1 使估计方差最小，也就是使估计值Z f X)和实际值 z()之 差的 平方和最 小 即 =E l Z(x)-Z (X)l=E I z()一 z()I 用 协 方 差 函 数 可 以 表 达 为 L i=1 =c(，)+c(，x,)-2 Z c(，)(2)要使估计方差最小，根据拉格朗 日乘数原理，令：F=一 2(喜 一 1 求 F 对 和 的 偏 导 数，并令偏导数为 0，得克里金方程组 整理后 2 c(x ，1 2 t-0 (3)(i=1，2，n)(4)解线性方程组(4)，求出权重系数 和拉格朗 日 乘数，代入公式(1)和(2)，分别求出估计值和 估计方差 在变异函数存在 的条件下，根据协方差 函数 c(h)及变异函数r(h)的关系式：c(h)=c(O)一 r(h)，用变异函数表示普通克里金方程组和克里金估计方 差，即 =1 i=1 ，)+=，)=1 (i=1 2一，n)(5)=(，)一 (，)+(6)i=l 方程组也可用矩阵表示，令 e l 1 C1 2 C 2 1 C2 2 C n l C n 2 1 1 。C 1 1 C 2 1 ：C n n 1 1 0 =：一，D=6 9 c(，X)c(，X)c(，)1 则普通克里金方程组为 K2=D(7)解方程组(7)，可得=K D 其估计方差为=c(x，1 一 D (8)1 2 克里金插值的计算步骤 克里金插值的前提是：根据空间场的结构，选择 适当的变异函数模型，并求出变异函数 其计算步骤 如下(1)网格化 选择区域的范围及网格的大小(2)计算被估点(即网格节点)坐标(3)根据搜索策略【3(近点距离搜索，方位搜 索)选择合适的参估点 (4)根据已经求出的变异函数，求出方程组的 系数 F C 1 1 C 12 e l 1 K=c 2 1 c 2 2C 2 n 1 1 2C n n 1 1 1 1 0 ：一，D=c(，X)c(，)c(，X)1 (5)解方程组(采用L u 方法)，求权系数 (6)用Z ()=z()求 被 估 点 的 值，其 i=1 中n 是插值点的个数(7)重复(2)(6)步，直到网格节点的值全 部求出(8)输出结果 1 3 计算实例 图1 为空间正方形网格数据，四个观测点，的高程值分别为Z()=3 7、Z(x 2)=4 2、z(x 3)=3 6、z(x 4)=3 5 假设高程值 的变异函数 一 I i 2 lJ=维普资讯 7 0 3 7一 r 3 7一 l(3 7 u 一3 5 3 8 3 7 3 7 3 3 3 4 3 5 3 8 X 0(j 3 5 3 7 (3 6)3 6 3 5 3 6 x 4(3 5)3 6 3 5 3 4 3 3 3 2 2 9 2 8 J 3 8 3 7 3 5 3 0 2 9 3 0 3 2 图 1 空间正方形网格数据(点间距 r-1 k m)是 同 同性(即变 异 函数 在 各个 方 向 的变化 郡 相 I J)的球状模型，则用回归分析法求出其球状变异函数模 型(有关变异函数的求解过程见文献 3 )为 f 0 =0 ()20 4 8+1 1 5 4(三 l_ _ )0 8 5 3 5 (9)用普通克里金法估计观测点x o 的高程值 Z(x o)，Z(x o)估计的基本公式为 z ()=z()(l，2，3，4)当 i=-，时，C 1 l=c 2 2=C 3 3=C 4 4=c(0)=C o+c=2 0 48+1 1 5 4=3 2 0 2 根据克里金矩阵的对称性，当i 时，C j=c (1 一 1)=3 2 0 2 一y l X i X I 由此计算可得=c =C 0 4=3 2 0 2 一 )：3 2 0 2 一 )=3 2 0 2 一【2 0 4 8+1 1 5 4(三2 8 5 3 5 2 5 3 5)=0 8 7 0 、8 q 3=c 31=3 2 0 2 一 f 1 _ 0 5 4 2 c l4=c 4 1=c 0 2=3 2 0 2 一 f 厮1=o。7 1 1 c 23=c 32=3 2 0 2 一 f 1=0 6 0 1 c 4=c 43=3 2 0 2 一 f 1 _ 0 3 8 3 =c 42=3 2 0 2 一 f 1 _ 0 4 6 6 c o1=3 2 0 2 一 f)_ 0 9 5 2 c o3=3 2 0 2 一 f 1 _ 0 5 7 1 天津城市建设学院学报2 0 0 8 年筮 鲞 筮 塑 将以上计算结果代人式(7)得 3 2 0 2 0 8 7 0 0 5 4 2 0 7l l l 0 2 8 7 0 2l 0 0 2 0 2 0 3 0l -0 47 3 0 8 7 0 3 2 0 2 0 6 0l 0 46 6 l 0 5 42 0 6 0l 3 20 2 0 3 8 3 l 0 7l l l 0 4 6 6 l 0 3 8 3 l 3 2 0 2 l l 0 0 9 5 2 0 7l l 0 5 7l 0 8 7 0 l 即克里金权重系数分别为：=0 2 8 7，=0 2 1 0，=0 2 0 2，=0 3 0 1，=0 4 7 3 根据普通 克里 金的 基本原理，点的高程值z()的 估计值为 Z 0 =0 2 8 7$z()+0 2 1 0*Z(x 2)+o 2 0 2$z()+o 3 0 1$z(x 4)=0 2 8 73 7+0 21 042+0 2 023 6+0 3 01 3 5=3 7 2 4 6mm 克里金估计方差为 4 o K =c(，)一 C(X i，)+=i=1 3 2 0 2一f 0 2 8 70 9 5 2十0 21 0 x0 71 1+0 2 0 2 x 0 5 7 1+0 3 0 1 x 0 8 7 O)+0 4 7 3=2 66 9mm 2 应 用实例 基于克里金插值方法，在 Wi n d o w s环境下，以 V i s u a l c+6 0为开发平台，以 O p e n G L为工具建立 了基于 D E M 的等高线模型 建模时分别使用了距离 反比插值算法与克里金插值算法，图 l 一 4 分别是距离 反比插值和克里金插值后生成的二维和三维等高线 模型 图 1 克里金插值后生成的二维等高线 五 加 加 2 娌 “加 鸲$加 维普资讯 天津城市建设学院学报 李莉等：克里金插值算法在等高线绘制中的应用 图2 距离反比插值后生成的二维等高线 图3 克里金插值后生成的三维等高线 图 4 距离反比插值后生成的三维等高线 从图中可以看出，克里金插值后生成的等高线模 型比较光滑，精度更高 克里金方法与距离反比法的 比较如表 1 所示 表 1 克里金方法与距离反比法的比较 性质 插值方法 逼近程度 外推 能力 运算速度 使用范 围 距离 数据分布均 很 差 快 数据分布 反比法 匀时效果好 均匀 克里金 高 很强 慢 均可，使用 方法 广泛 3 结论 7l 在绘制等高线的过程中，选择的插值方法会直接 影响到插值误差 的大小，从而影响所建模型的精 度 因此需要对被描述的数据进行仔细地研究和分 析，最终根据数据的特点，找到一种误差较小的插值 方法 相对来说，克里金插值能较好地反映各种地形 变化，但其计算量很大，尤其在对大面积区域、大数 据量内插时，这是一个不能不考虑 的因素 在实际的 建模过程中，不但要考虑到模型的精度，还要考虑到 系统 的运算速度，最终确定一种性价比最好的方 案 随着克里金插值法广泛应用于降雨分布、石油工 程、地质学等领域，它的局限性逐渐被认识，并在应 用过程中得以完善 参考文献：1 彭黎，金益民 利用 D E M 数据绘制等高线的自 适应 算法 J 微机发展，2 0 0 5，1 5(4)：3 3 3 5 2 李志林 数字高程模型 M 武汉：武汉大学出版社，2 o 0 3 3 颜慧敏 空间插值技术的开发与实现 D 成都：西南 石油学院，2 0 0 5 4 贾志刚 精通 O p e n G L M 北京：电子工业出版社，1 9 9 8 5 张仁铎 空间变异理论及应用 M 北京：科学出版 社，2 0 0 5 6 苏姝，林爱文，刘庆华 普通 K 啦i I l g法在空间内插中 的运用 J 江南大学学报：自然科学版，2 0 0 4，3 (1)：1 8 2 1 (编辑：胡玉敏)Ap p l i c a t i o n o f Kr i g i n g I n t e r p o l a t i o n i n Co n t o u r Cr e a t i o n LI Li HU J i a n p i n g b (a D e p a r t me n t o f C i v i l E n g i n e e ri n g；b De p a r t me n t o f E l e c t r o n i c s a n d I n f o r ma t i o n E n g i n e e r i n g，T I UC，T i a n j in 3 0 0 3 8 4，C h i n a)Ab s t r a c t：Us i n g o r d i n a r y k r i g i n g i n t e r p o l a t i o n me t h o d，the d i s c r e t e DEM d a t a i s i n t e r p ol a t e d I n the W i n d o ws e n v i r o n me n t，wi th Vi s u a l C+6 0 a s d e v e l o p me n t p l a t f o rm an d Op e n GL a s the t o o l，c o n t o u r mod e l b a s ed o n the DEM i s c r e a t ed Kr i g ing i n t e rpo l a t i o n me thod i s c o mp a r ed wi th i n v e r s e d i s t a n c e i n t e rpo l a t i o n me thod Re s u l t s h o ws tha t c o n t o u r mo d e l，wh i c h i s c r e a t e d b y the u s ing o f k r i g i n g int e rpo l a t i o n，i s mo r e p r eci s e r an d h a s l e s s e r r o r Ke y wo r d s：DEM；c o n t o u r；k r i g i n g i n t e rpo l a ti o n；i n v e r s e d i s t a n c e 维普资讯