基于时空变异函数的Kriging插值及实现_李莎
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Computer Engineering and Applications 计算机工程与应用
2011， 47 （23）
25
基于时空变异函数的 Kriging 插值及Βιβλιοθήκη Baidu现
2 李 莎 1， ， 舒 红 1， 董 林1 2 LI Sha1， ， SHU Hong1， DONG Lin1
1.武汉大学 测绘遥感信息工程国家重点实验室， 武汉 430079 2.湖北第二师范学院 机械与电气工程系， 武汉 430205 1.State Key Lab of Info. Engineering in Surveying， Mapping and Remote Sensing， Wuhan University， Wuhan 430079， China 2.Department of Mechanical and Electrical Engineering， Hubei University of Education， Wuhan 430205， China LI Sha， SHU Hong， DONG Lin.Research and realization of Kriging interpolation based on spatial-temporal variogram. Computer Engineering and Applications， 2011， 47 （23） ： 25-26. Abstract：The Kriging method is generally used in spatial variable interpolation， but not directly in spatial-temporal variable. It needs to be extended to space-time.The spatial-temporal Kriging， with R language， is applied for the spatial-temporal interpolation research and realization of monthly average temperature.The seasonal part has been removed from original temperature data by time series decomposition.A kind of product-sum variogram in space-time， describing the spatial-temporal correlation， is constructed based on pure spatial variogram and pure temporal one.The realization steps with R software are given. Spatial-temporal Kriging extended from ordinary Kriging is used in the temperature data.The experimental results show that this Kriging method based on spatial-temporal variogram has satisfied accuracy， which supplies an effective approach for interpolation and estimation of spatial-temporal variables. Key words：spatial-temporal correlation； variogram； Kriging interpolation； R language； temperature 摘 要： Kriging （克里金） 算法通常用于对空间变量进行插值， 但不能直接应用于时空变量， 它需要进行时空扩展。以月平均气温
random seasonal trend observed 20 0 －20 5.5 4.0 0 420 0 －4 1995 2000 年份 2005
} rt<-function （ht） { if （hs>range_t） return （psill_t+nugget_t）
3 （psill_t+nugget_t） * （1.5*ht/range_t－0.5* （hs/range_t） ）
26
2011， 47 （23）
Computer Engineering and Applications 计算机工程与应用 （如 “Sph” “ ，Gau” “ ，Exp” 等） ， psill， range 和 nugget 分别是 model 所指变异函数的基台值、 变程和块金值。然而 psill， range 和 nugget 这三个参数是用户预先设定的估计值， 并不代表模型实 际的拟合参数， 实际值被返回至数组 tem.fit 中。 （2） 根据第 （1） 步结果分别构建空间变异函数 rs （hs） 和时 k2 、 k3 的值。本 间变异函数 rt （ht） ， 并由式 （5） 计算出系数 k1 、 实验中， rs （hs） 采用高斯模型拟合， rt （ht） 采用球形模型拟合， 分别由以下 R 程序段得到：
图1
哈尔滨站点时序分解图
3 时空变异函数 3.1 时空变异函数建模
假设 Z (s t) 是定义在 Rk ´ T 上的时空随机过程， 其中 Rk 代表 k 维的欧式空间， T 代表时间， （si， t） i=1， 2， …， n 为时空场 i ， 中的任意样本点位置， 并且 h 为样本点间的时空间隔距离。由 于空间域与时间域的量纲不同且难以统一， 前者以米、 公里等 为计量单位， 后者以小时、 天等为计量单位， 因此时空距离不 能简单的将时空域看成三维空间来进行计算。设定时空距离 h= （hs， h） 当Z （s， t） 满足二阶平稳时， 可定义其协方差函数为： t ，
Z *(s0 t0) = å λi Z (si ti)
i=1 n
（6）
Z *(s0 t0) 为时空点 (s0 t0) 处的估计值，λi 为邻近观测值 式中： Z (si ti) 的加权系数， 引入拉格朗日系数 μ 进行推导可得： | γ11 γ1n 1 || λ1 | | γ01 | | ||| || | | ||| || = | | | （7） |γ γ 1 |||| λ n ||| | γ0n | nn | n1 || | | | | 1 1 0 || μ | | 1 |
1
引言
一个随机函数在空间位置上的测量值可以看成是随机函 数在这些位置上的一次实现值[1]， Kriging 通常用于在没有测量 值的地方对随机函数进行估计或插值， 从而获得连续的空间 数据面。Kriging 插值的重要前提是研究随机变量的空间相关 结构， 即变量性质随距离的相关关系， 构建有效的变异函数。 然而， 众多自然现象中的变量不仅具有空间特性， 同时还 具有时间特性， 因此在对这种时空变量进行插值之前， 除了分 析变量的空间相关特性之外， 还应充分考虑其时间相关特 性。所以原本用于空间 Kriging 的变异函数不能直接用于时空 变量的数据插值中， 必须进行时空扩展， 获得有效的时空变异 函数模型[2-6]。 本文的研究目的是在一类积和时空变异函数的基础上实 现对时空变量的 Kriging 插值。以气温数据为研究对象， 以R
rs<-function （hs） {
2 nugget_s+psill_s* （1－exp （－ （hs/range_s） ） ）
分解后各项的值被返回至数组 x 中， 其中 x$trend 为趋势项， x$seasonal 为季节项， x$random 为随机项。哈尔滨站点的时序 分解图如图 1 所示。将站点原始观测值减去季节项后保存至 文件 residual.csv， 用作后续的变异函数建模和 Kriging 插值。
C st (h s h t) = k1C s (h s)C t (h t) + k 2C s (h s) + k3C t (h t) γ st (h s h t) = (k1C t (0) + k 2)γ s (h s) + (k1C s (0) + k3)γ t (h t) - k1γ s (h s)γ t (h t)
基金项目： 国家高技术研究发展计划 （863） （the National High-Tech Research and Development Plan of China under Grant No.2008AA12Z201） ； 半干旱气候变化教育部重点实验室 （兰州大学） 开放基金。 作者简介： 李莎 （1983—） ， 女， 博士研究生， 讲师， 主要研究领域为地质统计、 时空信息分析； 舒红 （1969—） ， 男， 博士， 教授； 董林 （1984—） ， 男， 博士 研究生。E-mail： ls_li8202@126.com 收稿日期： 2011-03-08； 修回日期： 2011-05-24
[7]
实际采样时， 由于受外界环境等不确定因素的影响， 往往 使得观测数据在某些位置出现缺失； 此外由于受到人力、 物力 的限制， 使得采样值密度不够。这种情况下通常需要对数据 进行必要的插值。Kriging 插值是一种线性的无偏估计方法， 即 BLUE （Best Linear Unbiased Estimator） ， 常用的有简单 Kriging （SK） 、 普通 Kriging （OK） 、 对数正态 Kriging （LK） 、 泛 Kriging （UK） 等。本文选用普通 Kriging 方法实现数据的时空 插值。如下：
C (h s h t) = Cov(Z (s + h s t + h t) - Z (s t))
}
（3） 将 rs （hs） 和 rt （ht） 代入式 （4） 得到时空变异函数模型 rst （hs， ht） 。
4
时空 Kriging 插值
（1）
同时， 变异函数为： γ(h s h t) = 1 E(Z (s + h s t + h t) - Z (s t))2 = σ 2 - C (h s h t) （2） 2 式中 σ 2 为 Z (s t) 的方差。在满足相应的正定条件 [1] 下， 变异 函数是有效的。 实验中往往计算采样点的样本变异函数并选用时空变异 函数模型进行拟合。采用一类积和式变异函数 来拟合月均 气温的时空变异结构， 如下：
数据为例， 运用时空 Kriging 方法结合 R 统计语言进行时空插值研究及其实现。通过时序分解去除气温数据中季节变化项， 在分 别得到空间变异函数和时间变异函数的基础上构建一类积和式时空变异函数来描述变量的时空相关结构， 并给出基于 R 语言的 具体实现步骤。将普通 Kriging 方法进行时空扩展， 应用于气温数据的时空插值中。验证结果表明， 基于时空变异函数的 Kriging 方法能提供较高精度的插值效果， 这为时空变量的插值预测提供了有效的途径。 关键词： 时空相关； 变异函数； 克里金插值； R 语言； 气温 DOI： 10.3778/j.issn.1002-8331.2011.23.007 文章编号： 1002-8331 （2011） 23-0025-02 文献标识码： A 中图分类号： TP39
语言编程实现数据的时序分解、 变异函数建模和时空插值。 实验数据取自中国气象科学数据共享服务网上黑龙江省 31 个 气象站点 1991 年 1 月至 2008 年 12 月共 216 个月的月均气温数据。
2
时序分解
月均气温是一组时间序列值， 具有显著的季节变化， 这势 必影响插值的结果。时间序列一般包括三部分： 趋势项 T 、 季 节项 S 和随机项 R 。其中随机项 R 满足二阶平稳。为了提高 插值的精度， 原始气温数据先就各站点分别作去季节项处 理。定义向量 observed 存放站点 216 个月的原始气温观测数 据， 时序分解可通过如下 R 程序段得到。 observed<-ts （ observed， start=c （ 1991， 1 ） ， end=c （ 2008， 12 ） ， frequency=12） x<-decompose （observed）
2011， 47 （23）
25
基于时空变异函数的 Kriging 插值及Βιβλιοθήκη Baidu现
2 李 莎 1， ， 舒 红 1， 董 林1 2 LI Sha1， ， SHU Hong1， DONG Lin1
1.武汉大学 测绘遥感信息工程国家重点实验室， 武汉 430079 2.湖北第二师范学院 机械与电气工程系， 武汉 430205 1.State Key Lab of Info. Engineering in Surveying， Mapping and Remote Sensing， Wuhan University， Wuhan 430079， China 2.Department of Mechanical and Electrical Engineering， Hubei University of Education， Wuhan 430205， China LI Sha， SHU Hong， DONG Lin.Research and realization of Kriging interpolation based on spatial-temporal variogram. Computer Engineering and Applications， 2011， 47 （23） ： 25-26. Abstract：The Kriging method is generally used in spatial variable interpolation， but not directly in spatial-temporal variable. It needs to be extended to space-time.The spatial-temporal Kriging， with R language， is applied for the spatial-temporal interpolation research and realization of monthly average temperature.The seasonal part has been removed from original temperature data by time series decomposition.A kind of product-sum variogram in space-time， describing the spatial-temporal correlation， is constructed based on pure spatial variogram and pure temporal one.The realization steps with R software are given. Spatial-temporal Kriging extended from ordinary Kriging is used in the temperature data.The experimental results show that this Kriging method based on spatial-temporal variogram has satisfied accuracy， which supplies an effective approach for interpolation and estimation of spatial-temporal variables. Key words：spatial-temporal correlation； variogram； Kriging interpolation； R language； temperature 摘 要： Kriging （克里金） 算法通常用于对空间变量进行插值， 但不能直接应用于时空变量， 它需要进行时空扩展。以月平均气温
random seasonal trend observed 20 0 －20 5.5 4.0 0 420 0 －4 1995 2000 年份 2005
} rt<-function （ht） { if （hs>range_t） return （psill_t+nugget_t）
3 （psill_t+nugget_t） * （1.5*ht/range_t－0.5* （hs/range_t） ）
26
2011， 47 （23）
Computer Engineering and Applications 计算机工程与应用 （如 “Sph” “ ，Gau” “ ，Exp” 等） ， psill， range 和 nugget 分别是 model 所指变异函数的基台值、 变程和块金值。然而 psill， range 和 nugget 这三个参数是用户预先设定的估计值， 并不代表模型实 际的拟合参数， 实际值被返回至数组 tem.fit 中。 （2） 根据第 （1） 步结果分别构建空间变异函数 rs （hs） 和时 k2 、 k3 的值。本 间变异函数 rt （ht） ， 并由式 （5） 计算出系数 k1 、 实验中， rs （hs） 采用高斯模型拟合， rt （ht） 采用球形模型拟合， 分别由以下 R 程序段得到：
图1
哈尔滨站点时序分解图
3 时空变异函数 3.1 时空变异函数建模
假设 Z (s t) 是定义在 Rk ´ T 上的时空随机过程， 其中 Rk 代表 k 维的欧式空间， T 代表时间， （si， t） i=1， 2， …， n 为时空场 i ， 中的任意样本点位置， 并且 h 为样本点间的时空间隔距离。由 于空间域与时间域的量纲不同且难以统一， 前者以米、 公里等 为计量单位， 后者以小时、 天等为计量单位， 因此时空距离不 能简单的将时空域看成三维空间来进行计算。设定时空距离 h= （hs， h） 当Z （s， t） 满足二阶平稳时， 可定义其协方差函数为： t ，
Z *(s0 t0) = å λi Z (si ti)
i=1 n
（6）
Z *(s0 t0) 为时空点 (s0 t0) 处的估计值，λi 为邻近观测值 式中： Z (si ti) 的加权系数， 引入拉格朗日系数 μ 进行推导可得： | γ11 γ1n 1 || λ1 | | γ01 | | ||| || | | ||| || = | | | （7） |γ γ 1 |||| λ n ||| | γ0n | nn | n1 || | | | | 1 1 0 || μ | | 1 |
1
引言
一个随机函数在空间位置上的测量值可以看成是随机函 数在这些位置上的一次实现值[1]， Kriging 通常用于在没有测量 值的地方对随机函数进行估计或插值， 从而获得连续的空间 数据面。Kriging 插值的重要前提是研究随机变量的空间相关 结构， 即变量性质随距离的相关关系， 构建有效的变异函数。 然而， 众多自然现象中的变量不仅具有空间特性， 同时还 具有时间特性， 因此在对这种时空变量进行插值之前， 除了分 析变量的空间相关特性之外， 还应充分考虑其时间相关特 性。所以原本用于空间 Kriging 的变异函数不能直接用于时空 变量的数据插值中， 必须进行时空扩展， 获得有效的时空变异 函数模型[2-6]。 本文的研究目的是在一类积和时空变异函数的基础上实 现对时空变量的 Kriging 插值。以气温数据为研究对象， 以R
rs<-function （hs） {
2 nugget_s+psill_s* （1－exp （－ （hs/range_s） ） ）
分解后各项的值被返回至数组 x 中， 其中 x$trend 为趋势项， x$seasonal 为季节项， x$random 为随机项。哈尔滨站点的时序 分解图如图 1 所示。将站点原始观测值减去季节项后保存至 文件 residual.csv， 用作后续的变异函数建模和 Kriging 插值。
C st (h s h t) = k1C s (h s)C t (h t) + k 2C s (h s) + k3C t (h t) γ st (h s h t) = (k1C t (0) + k 2)γ s (h s) + (k1C s (0) + k3)γ t (h t) - k1γ s (h s)γ t (h t)
基金项目： 国家高技术研究发展计划 （863） （the National High-Tech Research and Development Plan of China under Grant No.2008AA12Z201） ； 半干旱气候变化教育部重点实验室 （兰州大学） 开放基金。 作者简介： 李莎 （1983—） ， 女， 博士研究生， 讲师， 主要研究领域为地质统计、 时空信息分析； 舒红 （1969—） ， 男， 博士， 教授； 董林 （1984—） ， 男， 博士 研究生。E-mail： ls_li8202@126.com 收稿日期： 2011-03-08； 修回日期： 2011-05-24
[7]
实际采样时， 由于受外界环境等不确定因素的影响， 往往 使得观测数据在某些位置出现缺失； 此外由于受到人力、 物力 的限制， 使得采样值密度不够。这种情况下通常需要对数据 进行必要的插值。Kriging 插值是一种线性的无偏估计方法， 即 BLUE （Best Linear Unbiased Estimator） ， 常用的有简单 Kriging （SK） 、 普通 Kriging （OK） 、 对数正态 Kriging （LK） 、 泛 Kriging （UK） 等。本文选用普通 Kriging 方法实现数据的时空 插值。如下：
C (h s h t) = Cov(Z (s + h s t + h t) - Z (s t))
}
（3） 将 rs （hs） 和 rt （ht） 代入式 （4） 得到时空变异函数模型 rst （hs， ht） 。
4
时空 Kriging 插值
（1）
同时， 变异函数为： γ(h s h t) = 1 E(Z (s + h s t + h t) - Z (s t))2 = σ 2 - C (h s h t) （2） 2 式中 σ 2 为 Z (s t) 的方差。在满足相应的正定条件 [1] 下， 变异 函数是有效的。 实验中往往计算采样点的样本变异函数并选用时空变异 函数模型进行拟合。采用一类积和式变异函数 来拟合月均 气温的时空变异结构， 如下：
数据为例， 运用时空 Kriging 方法结合 R 统计语言进行时空插值研究及其实现。通过时序分解去除气温数据中季节变化项， 在分 别得到空间变异函数和时间变异函数的基础上构建一类积和式时空变异函数来描述变量的时空相关结构， 并给出基于 R 语言的 具体实现步骤。将普通 Kriging 方法进行时空扩展， 应用于气温数据的时空插值中。验证结果表明， 基于时空变异函数的 Kriging 方法能提供较高精度的插值效果， 这为时空变量的插值预测提供了有效的途径。 关键词： 时空相关； 变异函数； 克里金插值； R 语言； 气温 DOI： 10.3778/j.issn.1002-8331.2011.23.007 文章编号： 1002-8331 （2011） 23-0025-02 文献标识码： A 中图分类号： TP39
语言编程实现数据的时序分解、 变异函数建模和时空插值。 实验数据取自中国气象科学数据共享服务网上黑龙江省 31 个 气象站点 1991 年 1 月至 2008 年 12 月共 216 个月的月均气温数据。
2
时序分解
月均气温是一组时间序列值， 具有显著的季节变化， 这势 必影响插值的结果。时间序列一般包括三部分： 趋势项 T 、 季 节项 S 和随机项 R 。其中随机项 R 满足二阶平稳。为了提高 插值的精度， 原始气温数据先就各站点分别作去季节项处 理。定义向量 observed 存放站点 216 个月的原始气温观测数 据， 时序分解可通过如下 R 程序段得到。 observed<-ts （ observed， start=c （ 1991， 1 ） ， end=c （ 2008， 12 ） ， frequency=12） x<-decompose （observed）