时空CoKriging的变异函数建模

胡丹桂, 舒 红, 胡泓达

(1.武汉大学 测绘遥感信息工程国家重点实验室， 武汉 430079；2.武汉职业技术学院 计算机学院， 武汉 430074)



时空CoKriging的变异函数建模

胡丹桂1，2, 舒 红1*, 胡泓达1

(1.武汉大学 测绘遥感信息工程国家重点实验室， 武汉 430079；2.武汉职业技术学院 计算机学院， 武汉 430074)

在对地观测中，所研究的地学变量不仅具有时间、空间特征，还受其他变量的影响，采用多元时空相关数据，可以提高时空估值的精度.时空CoKriging是多元时空插值中一种常用的方法，建立时空变异函数和协变异函数是时空CoKriging插值的关键一步.以东北三省为试验区，利用该地区气象站观测数据的月平均空气相对湿度为主变量，引入同时间同位置的月平均空气温度作为协变量，对空气相对湿度和空气温度进行时空变异函数和时空协变异函数建模.实验结果表明，采用和度量时空模型的时空变异函数的随机性空间结构建模的实际拟合效果较好.

多元时空变量； 和度量模型； 变异函数； 空气相对湿度

地统计学是以区域化变量理论为基础，以变异函数为工具，研究那些在空间分布上既有随机性又有结构性的自然现象的科学.在自然科学及社会科学领域，有些变量不仅具有空间特征，而且具有时间特征，这时要把所研究的变量看成是时空随机函数.时空插值方法已成为解决时空离散点到连续体上的一种必不可少的映射方法，时空克里金是时空插值法中常用的一种方法[1].Marc G.Genton[2]，Matthew W. Mitchell[3]等对可分离时空协方差函数进行了理论研究和测试，对其优缺点进行了深入的探讨；Cressie and Huang[4]，Chunsheng Ma[5]，E. Porcu[6]，Gneiting[7]等提出了不可分离的，平稳的协方差函数，允许时空交互；徐爱萍[8]，L.De Cesare[9]，D.E.Myers[10]，S. De Iaco[11]等将积和模型应用于时空变异函数模型的建立中.他们虽然考虑到了时间空间特征，但仅局限于单变量.相对照，传统的多元统计分析方法虽然考虑了多变量，却很少考虑到多元信息的空间特征.但是，在诸如水文、石油、土壤、农林、大气和环保等科学研究中，所研究的变量不仅具有时间、空间特征，而且还在时间空间域中受其他相关变量的影响.将克里金法应用于时空多元变量的插值研究中，一方面需要将克里金法进行时空扩展，另一方面将单变量插值的克里金扩展到多变量协同克里金.建立有效的多元时空协方差变异函数模型是研究时空协同克里金插值的关键一步.因此，本文将传统地统计学扩展到多元时空地统计学，构建多元时空变异函数及协变异函数.

以东北三省72个气象站点1996年～2005年的月均空气相对湿度数据为例，以同时间同位置的月均空气温度为协变量，为时空协同克里金插值建立时空变异函数和时空协变异函数.采用了和度量变异函数模型来拟合变量的时空变异结构，将气象要素模拟从空间维扩展至时空维.同时，考虑了气温数据对相对湿度的影响，在拟合相对湿度变异函数的过程中，加入了气温作为协变量，将时空克里金插值扩展到了时空协同克里金插值.

1 研究区域介绍及数据处理

1.1 区域和数据介绍

实验研究的是我国东北三省(黑龙江、吉林、辽宁)，站点地处38.90°～52.97°N，119.70°～132.97°E之间.地面观测台站所观测的资料来自1996年1月～2005年12月东北三省月空气相对湿度.本实验数据通过中国气象科学数据共享服务网获取，共有观测站点87个，由于部分站点数据严重缺失，实验只采用72个观测站数据，站点分布如图1.实验数据为1996年1月～2005年12月的月平均相对湿度和月平均气温值.

由于空气相对湿度与空气温度相关性较好，我国东北地区，冬季干燥、夏季湿润，空气温度对空气湿度有着一定程度的相关影响，故选取空气温度为协变量.在实验中，对空气湿度和空气温度分别进行时空变异函数及两者的协变异函数拟合.

图1 东北三省观测站点分布图Fig.1 Monitoring station distribution in the three provinces of Northeast China

1.2 数据预处理

时空变异函数构造的重要前提是假设时空变量满足二阶平稳.时空空气相对湿度、空气温度数据可以看作为空间站点上的时间序列，时间序列一般包括：周期项、趋势项、随机项.因此空气相对湿度变量的时间序列分解为：

(1)

为保持数据连续性和整体趋势变化，时间序列分解后的趋势项和随机项保留，因为趋势项不明显，对数据平稳性影响不大.采用时序分解中的加法模型，将变量的周期项提取出来，余下的部分用于时空插值实验.用自相关图检法[12]判断去周期数据的平稳性，如图 4图 5所示，空气相对湿度和空气温度他们的时间序列的自相关系数(ACF)均随延迟时期数很快衰减到±0.1以内，图示表明去周期项后数据近似平稳.

图2 塔河空气湿度时间序列综合分析图Fig.2 Air humidity time series decomposition in Tahe

图3 塔河空气温度时间序列综合分析图Fig.3 Air temperature time series decomposition in Tahe

图4 塔河空气相对湿度时间序列随机项自相关法检测平稳性Fig.4 The test of stationarity using autocorrelation method for Tahe Random

图5 塔河空气温度时间序列随机项自相关法检测平稳性Fig.5 The test of stationarity using autocorrelation method for Tahe random

空气湿度和温度变量在空间上也要满足二阶平稳，将每一测站上的时间序列经过去周期项后计算移动平均值，根据移动平均值探索空间空气湿度和温度的变化趋势.用空间散点图展示空间分布趋势，图 6是对湿度空间站点数据趋势面的拟合.通过实验对比，发现不管是湿度还是温度均采用3次多项式(2)进行拟合，具有较好的拟合效果.

(2)

式中，value为月均空气湿度，x，y为空间坐标信息，a，b，c，d，e，f，g为3次多项式通过拟合得到的系数.去除趋势后空间数据分布呈现平稳性，图 7是月平均空气湿度站点数据去空间趋势后拟合面.分别对空气湿度和空气温度各站的时间序列进行平稳性处理再对整个空间数据进行平稳处理，最终满足构造时空变异函数的前提条件，即随机变量近似时空上的二阶平稳.

图6 月平均湿度空间站点数据趋势面拟合Fig.6 Fitting trend surface for monthly average humidity data

图7 月平均湿度站点数据去空间趋势后拟合面Fig.7 Fitting surface for monthly average humidity data removed spatial trending

2 时空CoKriging变异函数建模方法

(3)

时空协同克里金(CoKriging)法的插值公式：

(4)

2.1 时空直接变异函数建模

假设Z={Z(s，t)，s、t∈Rd+1}(d+1表示欧式空间维加时间维)表示为时空随机场，设定时空随机场中两位置的时空距离h=(hs，ht)，hs为矢量，代表样本空间距离同时也包含方向信息，ht为样本时间距离.当Z(s，t)满足二阶平稳时可定义其协方差函数为[9]：

C(hs，ht)=Cov(Z(s+hs，t+ht)，Z(s，t)).

(5)

显然，协方差函数只与距离有关，与空间和时间位置无关.在地统计学中，随机变量Z(s，t)的期望不变(时空二阶平稳性)，协方差矩阵具有对称正定性.理论上，满足下列要求的连续函数可以定义为有效变异函数[9]：

σ2-C(hs，ht),

(6)式中,σ2为Z(s，t)的方差，正定条件即任意ai∈R，i=1，…，n，和任意正整数n，C必须符合：

(7)

1)可分离时空协方差模型

地统计学中协方差函数模型(如球状模型、指数模型和高斯模型等)不能直接用于时空随机变量的统计分析，必需进行时空扩展.如果一个随机场Z的时空协方差函数可完全表示成纯空间和纯时间协方差函数相乘，则该时空协方差函数被认为是可分离的[14].Mitchell Genton， Gumpertz(2005)[3]，Montserrat Fuentes(2004)[15]等文献中提出了可分离型模型.这类模型相对容易构建，且计算机实现插值效率较高，但往往要求一些比较理想的假设[16-17]，损失了精细的时空结构信息或丢失了时空交互信息.

可分离型时空协方差函数可表示为：

C(hs，ht)=Cs(hs)·Ct(ht),

(8)

式中,C(hs，ht)是变量的时空协方差函数，Cs(hs)、Ct(ht)分别是纯空间协方差和纯时间协方差函数.Cs(hs)、Ct(ht)的具体形式有多种选择，或者同为一种模型，比如高斯模型，或者为不同模型，比如Cs(hs)为球形模型，Ct(ht)为指数模型，或是其他满足正定条件的函数.

2)时空积和模型

另一类称为不可分离型时空协方差函数，它是将已知的纯空间协方差与纯时间协方差函数通过加乘、混合、积分等变换得到.积和式变异函数来拟合时空地理数据的时空变异结构，协方差函数和变异函数如下[9]：

C(hs，ht)=k1Cs(hs)Ct(ht)+

k2Cs(hs)+k3Ct(ht),

(9)

γ(hs，ht)=(k1Ct(0)+k2)γs(hs)+

(k1Cs(0)+k3)γt(ht)-k1γs(hs)γt(ht),

(10)

式中，C(hs，ht)为时空协方差，Cs(hs)为空间协方差，Ct(ht)为时间协方差，γ(hs，ht)、γs(hs)、γt(ht)分别是对应的时空、空间、时间变异函数，而C(0，0)、Cs(0)、Ct(0)分别是对应的基台值.这里，假定γ(0，0)=γs(0)=γt(0)=0，满足k1>0、k2≥0、k3≥0，并通过正定条件由下式决定，推导过程见文献[18]：

(11)

3)时空和度量模型

时空和度量模型也是不可分离时空变异函数模型.采用和度量模型变异函数来拟合时空地理数据的时空变异结构，协方差函数可以表示为[19]：

(12)

式中，C(hs，ht)是空间间隔距离为hs，时间间隔距离为ht的协方差函数；

时空和度量变异函数为[19]：

(13)式中，γ(hs，ht)、γs(hs)、γt(ht) 分别为时空变异函数、纯空间变异函数、纯时间变异函数.这个模型的缺点是求出这10个参数的值存在一定的复杂性.

2.2 时空协变异函数建模

在使用时空CoKriging来研究变量的时空变异性时，关键的一步就是决定变量之间的时空协变异函数.它用来描述两变量之间交叉的时空连续性.协变异函数定义为：

(14)

(15)

即在同一个位置上，将两个变量的样本数值相加，所得之和即是新变量在该点的样本值，然后计算新变量的变异函数：

(16)

而这一新变量的变异函数与原变量的变异函数和协变异函数有如下关系：

(17)

因此可得：

(18)

3 时空CoKriging变异函数的拟合结果

分别采用时空可分离模型、时空积和模型、时空和度量模型3种不同的时空结合方法，对两个时空变量空气湿度和空气温度，采用时间去周期空间去趋势后的残差，来拟合时空变异函数.如图 8是空气相对湿度的残差的时空经验变异函数和理论变异函数.图 9是空气温度的残差的时空经验变异函数和理论变异函数.

图8 空气相对湿度残差的时空变异函数Fig.8 Empirical and fitted space-time variograms of the humidity residuals

图9 空气温度残差的时空变异函数Fig.9 Empirical and fitted space-time variograms of the temperature residuals

并采用均方根差RMSD(root-mean-squared-difference)指标来判断一个经验变异函数拟合的好坏.表1分别是3种不同方法得到的模型的RMSD值，从表 1可以看出，针对这两个时空变量，和度量时空模型的拟合效果最好.本实验中，采用和度量时空模型来建立空气湿度的残差和空气温度的残差的时空变异函数.

图10是将空气湿度的残差和温度的残差之和作为一个新变量，这个新变量的经验变异函数和经和度量时空模型拟合的理论变异函数.

表1 经验变异函数与理论变异函数之间的RMSD值Tab.1 RMSD between the empirical and fitted variograms

图10 空气湿度残差与温度残差之和的时空经验变异函数和理论变异函数Fig.10 Empirical and fitted space-time variograms of the sum of humidity residuals and temperature residuals

使用和度量模型中的时空变异函数模型拟合数据，其中球状模型被选为空间和时间变异函数，等式右边第3项度量时空变异模型为指数模型[19].这3个组件都具有块金(C0)，偏基(C1)和变程(r)的参数，度量时空变异模型还有各向异性比这个参数，一共10个参数.表 2是空气相对湿度的残差，使用和度量模型拟合时空变异函数的10个不同参数的值.表 3是空气温度的残差，使用和度量模型拟合时空变异函数的10个不同参数的值.表 4是空气湿度残差与温度残差之和，使用和度量模型拟合时空变异函数的10个不同参数的值.

表2 相对湿度残差和度量变异函数模型拟合参数Tab.2 Fitting parameters of the relative humidity residuals variogram using sum metric model

经过和度量时空模型拟合后，空气湿度残差的理论变异函数γ11为：

γ11(hs，ht)=1.01+7.41×

(19) 表3 空气温度残差和度量变异函数模型拟合参数Tab.3 Fitting parameters of the air temperature residuals variogram using sum metric model

经过和度量时空模型拟合后，空气温度残差的理论变异函数γ22为：

(20) 表4 空气湿度残差与温度残差之和的 和度量变异函数模型拟合参数Tab.4 Fitting parameters of the sum of the relative humidity residuals and air temperature residuals variogram using sum metric model

(21)

所以，根据式(18)，空气湿度的残差和空气温度的残差时空协变异函数γ12为：

(22)

4 结论

有效的变异函数模型是克里金插值的基础，本文的重点对时空协同克里金的变异函数进行建模分析.在拟合变异函数之前，首先对变量进行时间上去周期，空间上去趋势处理，以保证时空变量的平稳性.得到插值结果后，必须先将其加上之前去除的周期项和趋势项，才是最终估计结果.结果表明，空气相对湿度和空气温度在纯空间和纯时间上均符合球状模型，并且分别用3种不同的方法对这两个变量进行时空变异函数建模，发现在可分离模型、积和模型、和度量模型这3种模型中，和度量模型的拟合效果最好.完成时空直接变异函数拟合后，最后根据D.E.Myers[21]提出的方法，进行时空协变异函数建模.

本文研究了多元时空数据进行时空直接变异函数和时空协变异函数的建模，不仅考虑了时间空间信息，而且还考虑了其他协变量.时空变异函数是时空克里金插值模型的权重参数构建的基础，对后续的时空协同克里金插值精度影响非常大.后续研究中，将研究大样本下的时空协同克里金插值结果(估计量)的统计性质(无偏、最优、相合性、线性插值结果的渐近正态性)，发展适应性多元时空统计模型应用于多个地学领域时空数据分析并评价插值精度.

[1] 李 莎， 舒 红， 徐正全. 利用时空Kriging进行气温插值研究[J].武汉大学学报:信息科学版， 2012， 37(2)：237-241.

[2]GentonMG.Separableapproximationsofspace-timecovariancematrices[J].Environmetrics， 2007， 18: 681-695.

[3]MitchellMW，GentonMG，GumpertzML.TestingforSeparabilityofspace-timecovariances[J].Environmetrics， 2005， 16: 819-831.

[4]CressieN，HuangHC.Classesofnonseparablespatio-temporalstationarycovariancefunctions[J].JournaloftheAmericanStatisticalAssociation， 1999， 94:1330-1340．

[5]MaCS.Spatio-temporalstationarycovariancemodels[J].JournalofMultivariateAnalysis， 2003， 86: 97-107.

[6]PorcuE，GregoriP，MateuJ.Nonseparablestationaryanisotropicspace-timecovariancefunctions[J].EnvironmentalResearchandRiskAssessment， 2006， 21: 113-122.

[7]GneitingT.Nonseparable，stationarycovariancefunctionsforspace-timedata[J].JournaloftheAmericanStatisticalAssociation， 2002， 97: 590-600．

[8] 徐爱萍， 圣文顺， 舒 红. 时空积和模型的数据插值与交叉验证[J]. 武汉大学学报：信息科学版， 2012， 37(7)：766-770.

[9]DeCesareL，MyersDE，PosaD.Product-Sumcovarianceforspace-timemodeling:anenvironmentalapplication[J].Environmetrics， 2001， 12: 11-23．

[10]Myers.DE.Space-timecorrelationmodelsandcontaminantplumes[J].Environmetrics， 2002， 13:535-553.

[11]DeIacoS，MyersDE，PosaD.Space-timeanalysisusingageneralproduct-summodel[J].Statistics&ProbabilityLetters， 2001， 52 (1):21-28.

[12] 王 燕.应用时间序列分析 [M].第2版.北京：中国人民大学出版社， 2008: 111-130.

[13]DeIacoS，PalmaM，PosaD.Modelingandpredictionofmultivariatespace-timerandomfields[J].ComputationalStatistics&DataAnalysis， 2005， 48: 525-547.

[14]MateuJ，PorcuE，GregoriP.Recentadvancestomodelanisotropicspace-timedata[J].StatisticalMethodsandApplications， 2008， 17: 209-23.

[15]FuentesM.Testingforseparabilityofspatial-temporalcovariancefunctions[J].JournalofStatisticalPlanningandInference， 2006， 136: 447-466.

[16]SteinM.Space-timecovariancefunctions[J].JournaloftheAmericanStatisticalAssociation， 2005， 100:310-321．

[17]FuentesM，ChenL，DavisJM.Aclassofnonseparableandnonstationaryspatialtemporalcovariancefunctions[J].Environmetrics， 2008， 19: 487-507.

[18]DeCesareL，MyersDE，PosaD.Estimatingandmodelingspace-timecorrelationstructures[J].Statistics&ProbabilityLetters， 2001， 51(1):9-14．

[19]DanielA.GriffithGB,HeuvelinkM.Derivingspace-timevariogramsfromSpace-TimeAutoregressive(STAR)modelspecifications[J].TheInternationalArchivesofthePhotogrammetry，RemoteSensingandSpatialInformationSciences， 2009， 38(2):15-19.

[20]DimitrakopoulosR，LuoX.SpatiotemporalModeling:CovariancesandOrdinaryKrigingSystems[C]//DimitrakopoulosR.GeostatisticsfortheNextCentury.SpringerNetherlands:KluwerAcademicPublishers， 1994(6):88-93．

[21]MyersDE.MatrixformulationofCokriging[J].MathematicalGeology，1982， 14 (3):249-257.

Variogram modeling in space-time CoKriging

HU Dangui1，2, SHU Hong1, HU Hongda1

( 1.State Key Laboratory of Information Engineering in Surveying，Mapping and Remote Sensing， Wuhan University， Wuhan 430079;2. Institute of Computer， Wuhan Polytechnic， Wuhan 430074 )

The variables used in various environmental studies not only have the spatial temporal dependence， but also are affected by other variables. However， the multivariate space-time correlated data can be used to improve the accuracy of space-time prediction and interpolation. The space-time CoKriging is a commonly used method of multivariate space-time interpolation， with modeling a space-time direct variogram and a cross variogram as a key step. Therefore， taking the three provinces in Northeastern China as study area. this paper modeled the space-time direct variograms and cross variograms for air relative humidity and air temperature by using the monthly mean relative air humidity observed in meteorological stations as a main variable and the air temperature in the same position and at the same time as a covariate variable. Inverse distance weighted interpolation was used to get the current meteorological observation data. The experimental results show that the space-time variogram using sum-metric model is well fitted， compared to the deterministic structure modeling spatial distance inversely．

multivariate space-time variables; Sum-metric model; variograms;air relative humidity

2014-11-05.

国家自然科学基金项目(41171313)；苏州市科技计划项目( SYG201319)；地理空间信息工程国家测绘地理信息局重点实验室开放研究基金项目 (201329 )；湖北省自然科学基金项目 (2014CFB725).

1000-1190(2015)04-0596-07

P208< class="emphasis_bold">文献标识码： A

A

*通讯联系人. E-mail: shu_hong@whu.edu.cn.