气温数据的时空变异结构分析

郑 兴 文，舒 红，胡 泓 达

（1.贵州电力设计研究地理信息中心，贵州 贵阳 550002；2.武汉大学测绘遥感信息工程国家重点实验室，湖北 武汉 430079）

0 引言

自然环境中大量的现象可被视为时空随机场的实现［1，2］，可以建立基于随机场的时空模型来分析这些现象或观测数据。但是，目前地质统计学中多是纯空间模型，单纯使用空间模型分析时空数据会造成时间维有用信息的丢失。变异函数（亦称为半变异函数或是变差函数）是地质统计学中重要的组成部分，不仅因为它是许多地质统计学计算（如估计方差、离散方差、正则化变量的变异函数计算等）的基础，还由于它能反映和刻画区域化变量的许多性质［3］。实现空间变异函数模型到时空变异函数模型的扩展，是将变异函数中的空间距离（hs，hs∈Rd，d是物理空间维）扩展为时空距离（h，h∈Rd＋1，d＋1表示物理空间维加一个时间维）的过程。De Ilaoc等指出变异函数从空间到时空的扩展并不是一个简单的过程［4］，要解决时空量纲不一致（或时间和空间“距离”单位的统一）问题和确保时空变异函数模型（或时空协方差）的最优拟合问题等。针对这些问题，一些学者对时空变异函数进行了理论探讨。根据时空的结合方式将时空模型分为时空可分离模型和时空不可分离模型，时空可分离模型［5，6］只是简单地将时间协方差和空间协方差相乘或相加作为时空协方差，这种模型实现简单，但损失了精细的时空结构信息或丢失了时空交互信息，比如积模型中任意两空间点对上的两个时间序列交叉协方差函数相同，而实际上时空协方差函数应该随空间滞后距离呈正比：Cst（h1，ht）∝Cst（h2，ht）。通常，时空不可分离模型［7］使用Bochner理论和选择合适的谱密度构造，模型考虑了数据内在的丰富的时空结构信息，但构造方法复杂，很难获得模型。特别地，De Cesare等提出一种时空变异函数的积和模型［8，9］，其比一般不可分离模型更加灵活。

本文依据积和模型理论给出了时空变异函数的两种实现方法，给出数据时空平稳化的预处理方法。将空间和时间变异函数积和成时空变异函数并实现可视化效果，通过这种时空变异函数分析东北气象数据的时空结构特性。

1 时空变异函数

假设Z＝｛Z（s，t），s、t∈Rd＋1｝（d＋1表示欧式空间维加时间维且一般d≤3）表示为时空随机场，设定时空随机场中两位置的时空距离h＝（hs，ht），hs为矢量，代表样本空间距离同时也包含方向信息，ht为样本时间“距离”。当Z（s，t）满足二阶平稳时可定义其协方差函数为：

显然，协方差函数只与距离有关，与空间和时间位置无关。在随机变量Z（s，t）的期望不变（时空平稳性）且协方差矩阵［Cst］n×n正定的条件下，下列函数为有效变异函数且可表示为：

式中：σ2为Z（s，t）的方差，正定条件即任意ai∈ℜ，i＝1，…，n，和任意正整数n，Cst必须符合：

采用一类积和式变异函数拟合时空地理数据的时空变异结构，协方差函数和变异函数如下：

其中，Cst为时空协方差，Cs为空间协方差，Ct为时间协方差，γst、γs、γt分别是对应的时间、空间、时空变异函数，而Cst（0，0）、Cs（0）、Ct（0）分别是对应的基台值。这里假定γst（0，0）＝γs（0）＝γt（0）＝0，模型中的系数k1、k2、k3满足k1＞0、k2≥0、k3≥0，并通过正定条件由下式决定，推导过程见文献［8］：

式（5）中γs（hs）、γt（ht）的估值通过γst（hs，0）、γst（0，ht）的估值获得，两者通过式（5）和上述假定条件存在如下关系：

其中时间变异函数的滞后向量为rt，空间容差为δs，距离集合的分组数为：

2 时空变异函数的计算

依据时空变异函数的理论模型运用R语言设计时空变异函数，由纯空间和纯时间变异函数通过积和模式构造时空变异函数，本文运用两种方案对时空数据进行纯时间和纯空间变异函数构造，并通过模拟k1值确定Cst（0，0）将纯空间和纯时间变异函数积和为时空模型。

时空数据纯空间和纯时间变异函数构造方案一：纯空间变异函数先对相同测站的时间序列观测数据求平均，然后对不同测站的观测平均值进行变异函数计算。纯时间变异函数采用相同时态不同测站观测属性数据求平均，然后对不同时态的空间平均属性值进行变异函数计算。运用R中的变异函数自动拟合获取时空和空间变异函数公式。时空数据纯空间和纯时间变异函数构造方案二：纯空间变异函数采用对空间测站，不同时态的属性值分别计算变异函数值，获得变异函数块金值、偏基台值、变程等参数，然后对各时态计算的变异函数参数求平均，用块金值、偏基台值、变程的平均值生成最终空间变异函数。纯时间变异函数采用对应时态，不同空间测站时间序列属性数据分别计算变异函数，获得块金值、偏基台值、变程，然后对各测站点属性数据的变异函数参数求平均，用块金值、偏基台值、变程平均值生成时间变异函数。

3 实验分析

3.1 区域及数据介绍

东北三省地处119.70°～132.97°E，38.90°～52.97°N，北靠西伯利亚东部，深受寒冷干燥的冬季风影响，属寒带大陆性季风气候。年气温一般在－5～10℃，冬季较长。周围环山，中部以平原为主，山脉减弱了鄂霍次克海气流的直接侵入，朝鲜半岛和日本列岛阻挡了其与太平洋的直接接触，降低了海洋调节气候的作用。在全国综合自然区划中，东北区属于寒温带、温带、暖温带的湿润、半湿润地区［10］。

本实验数据通过中国气象科学数据共享服务网获取，共有观测站点87个，由于黑龙江肇州、吉林罗子河等6个站点数据严重缺失，实验只采用81个观测站1960年1月－2011年12月的月平均气温数据。

3.2 数据预处理

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

式中随机项R（t）满足二阶平稳，以大连站气温观测时间序列分解为例（图1），将观测气温时间序列（Observed）分解为趋势项（Trend）、季节项（Seasonal）、随机项（Random），可以看出自1960年1月到1989年12月30年的月平均气温值存在明显的周期性，周期为一年，即冬季气温低，夏季气温高，月平均气温年较差在30℃左右，具有明显的季节效应。时间序列存在一个总体升高的趋势，这可能与全球温室效应气温升高有关。

为保持数据连续性和整体变化趋势，时间序列分解后的趋势项和随机项保留，趋势项不明显，对数据平稳性影响不大。针对上述气温时间序列去除周期项后进行自相关法平稳性检测（图2），可见自回归系数随时间延迟阶数增大呈类似周期性余弦衰减，自相关系数很快衰减为0，说明时间序列去周期项后平稳［12］。

变量在空间上也要满足二阶平稳，将每一测站的时间序列经过去季节效应后计算移动平均值，据此探索空间气温变化趋势。用空间剖面图法分析空间站点在东西和南北方向的分布趋势（图3），月平均气温在x轴（东西）方向呈中间略低两侧略高，在y轴（南北）方向呈南高北低趋势，随维度变化明显，南北最大较差为15℃左右。用空间散点图展示空间分布趋势（图4），采用3次多项式value＝ax＋by＋c（x＊y）＋d（x2）＋e（y2）＋f（x3）＋g（y3）进行拟合（x，y为坐标信息，value为月均值温度），趋势面拟合适应度R2检测值为0.73，具有较好的拟合效果，去除趋势后空间数据呈平稳分布（图5）。

图5 月平均气温站点数据去空间趋势后拟合面Fig.5 Fitting surface for monthly average temperature data removed spatial trending

先对各站时间序列进行平稳性处理，再对整个空间数据进行平稳处理，最终达到构造时空变异函数的前提条件，即随机变量满足时空上的二阶平稳。

3.3 时空变异函数实现过程

据时空变异函数的实现方案一对时间和空间数据求均值，根据均值计算纯空间和纯时间变异函数（图6），拟合变异函数采用R库函数包automap中的autofitVariogram（formula，input＿data，mode＝c（"Sph"，"Exp"，"Gau"，"Ste"），kappa＝c（0.05，seq（0.2，2，0.1），5，10））来实现，formula本文采用"z～1"形式是表示为构造普通Kriging法的变异函数，z代表空间或时间月温度均值，其他形式表示运用其他种Kriging（如"z～x＋y"表示泛Kriging等），mode是可供选择的拟合模型，常见的有椭球模型（"Sph"）、指数模型（"Exp"）、高斯模型（"Gau"）和Ste Matern模型（"Ste"）等［13］，kappa是拟合平滑度参数。空间变异函数描述空间相关性随空间滞后的变化情况，随空间距离变大变异函数值升高后趋于平稳表示变量之间相关性变小直到不相关。从图6A中可见变程47 km内温度是相关的，变程以外温度不相关。时间变异函数描述时间变量随时间延迟的相关程度，时间延迟是一维的，这与空间变异函数不同，从图6B可见变程2.8个月内认为时间变量相关，块金值0.97反映了时间变量随机性的大小。

空间和时间变异函数的构造方案二是“先拟合，再对参数求均值”，空间和时间变异函数采用参数的均值构造。最终获得的空间变异函数和时间变异函数的偏基台值（Psill）分别为30.37、1.60，块金值（Nugget）分 别 为 0、0.80，变 程 （Range）分别 为2 035.15、3.25。

对比发现方案一和方案二空间变异函数的参数值具有较大差异，时间变异函数基本相同，分析原因发现第二种方案在确定空间变异函数中对应不同时间拟合变异函数（一共拟合360个）时有8个存在异常现象，异常值严重影响了平均值，因此先去除异常值再进行平均，获得相应参数见表1，结果与方案一基本相同。但方案二所用时间是方案一所用时间的几十倍，而且方案一针对一组数据进行拟合可自动选择最优拟合模型，方案二获得均值参数后要根据过程拟合中的模型选择最后变异函数模型。

空间和时间变异函数构造完成后，根据积和模型构造时空变异函数和时空协方差函数。时空变异函数：r1＋r2－k1＊r1＊r2（r1：纯空间变异函数，r2：纯时间变异函数）；时空协方差函数：k1＊（Snugget＋Spsill）＊（Tnugget＋Tpsill）＋k2＊（Snugget＋Spsill）＋k3＊（Tnugget＋Tpsill）－（r1＋r2－k1＊r1＊r2）（Snugget、Spsill为空间变异函数块金值和偏基台值，Tnugget、Tpsill为时间变异函数块金值和偏基台值）。k1、k2、k3的值满足式（6）要求，k1的值决定了时空变异函数基台值的大小。假设：

比较两种方法构建的时空变异函数：方法二中块金值比方法一小，方法二基台值、变程值比方法一大，说明方法二刻画的时空结构变异性要大，表达的时空结构性强，这与下文中变异程度结果相吻合。

4 结果分析

根据两种生成纯空间和纯时间变异函数方案分别构造时空变异函数和协方差函数（图7、图8），并获取基台值、块金值等相关参数（表2）。

表2 两种方案时空变异函数参数Table 2 Parameters of spatial-temporal variation function computed by two methods

空间变异程度是由随机性因素引起的空间异质性占系统总变异的比例，当变异程度在区间25%～75%时认为变量为中等相关性，变异程度越小表示空间结构性越好。空间变异程度推广到时空变异函数表示为随机性因素引起的时空异质性占系统总变异的比例，当变异程度小于25%认为时空变量具有强时空相关性，说明具有很好的时空结构，变异程度在25%～75%之间认为时空变量具有中等时空相关性，大于75%时空相关性弱。从表3可知方案二变异程度较方案一小但都在25%～75%之间，表示东北月平均气温具有中等的时空相关性，时空结构较好。用方案二拟合的变异函数表述的东北月平均气温时空相关性和时空结构性略强于第一种方案。

可从图6和表1获知方案一和方案二中纯空间和纯时间变异函数的变异程度分别是：方案一，纯空间变异函数变异程度为0，纯时间变异函数变异程度为57%；方案二，纯空间变异函数变异程度为0，纯时间变异函数变异程度为33%。可见东北三省月平均气温在空间上具有较强相关性，空间结构性强，而时间上相关性为中等，时间结构性较好。但时空积和模型的变异程度方案一方案二分别为69.7%和63%，说明时空结构性为中等。时空结构性下降主要可能是时间结构性不强造成的。

无论空间变量还是时间变量都认为在变程内是相关的，超出变程变量不相关。由图7、图8的时空变异函数可知空间变程不随时间延迟而变化，类似的时间变程也不随空间滞后而变化，这反映了时空具有相对稳定的结构特性。通过空间变程和时间变程可以获得东北三省月平均气温认为在空间48 km内、时间3个月内具有一定的相关性，超出此范围则认为东北月平均气温不存在相关性。

5 结论

本文基于一种积和模型的时空变异函数理论研究其实现过程，通过两种构造时空数据纯空间和纯时间变异函数方案分析其对时空变异函数构造的影响，发现第二种方案描述的东北月平均气温时空结构性更好，而方案一实现更为简单。提出一种使时空数据平稳的方法，思路是先对测站点上时间序列去周期性，再去除空间趋势面，以保证时空数据的平稳性。通过预处理的时空数据构造空间和时间变异函数，最终依据积和模型拟合获取时空变异函数并可视化。分析发现东北月平均气温的时空结构性较好，空间结构性好于时间结构性，并且在空间48 km内、时间3个月内（不考虑周期因素）东北月平均气温是相关的，超出此范围不相关。

本文讨论了时空随机场基础上的时空变异函数的积和模型，给出了时空随机场变异函数的计算方法，实验分析了所给出时空变异函数计算方法的可行性，为基于时空变异函数的各种时空估计模型的实现提供了计算基础。

［1］ EYNON B P，SWITZER P.The variability of rainfall acidity［J］.Canadian Journal of Statistics，1983，11（1）：11－23.

［2］ NHU LE D，PETKAU A J.The variability of rainfall acidity revisited［J］.Canadian Journal of Statistics，1988，16（1）：15－38.

［3］ 王仁铎，胡光道.线性地质统计学［M］.北京：地质出版社，1988：38－40

［4］ DE IACO S，MYERS D E，POSA D.Nonseparable space-time covariance models：Some parametric families［J］.Mathematical Geology，2002，34（1）：23－42.

［5］ DE CESARE L，MYERS D，POSA D.Spatial-Temproal Modeling of SO2in Milan District［M］.Geostatistic Wollongong′96.Dordrecht，1996.

［6］ ROUHANI S，HALL T J.Space-time kriging of groundwater data［A］.Geostatistics［C］.Springer Netherlands，1989.639－650.

［7］ CRESSIE N，HUANG H C.Classes of nonseparable，spatiotemporal stationary covariance functions［J］.Journal of the A-merican Statistical Association，1999，94（448）：1330－1339.

［8］ DE CESARE L，MYERS D，POSA D.Estimating and modeling space-time correlation structures［J］.Statistics ＆Probability Letters，2001，51（1）：9－14.

［9］ DE CESARE L，MYERS D E，POSA D.Product-sum covariance for space-time modeling：An environmental application［J］.Environmetrics，2001，12（1）：11－23.

［10］ 周琳.东北气候［M］.北京：气象出版社，1991.1－3.

［11］ 李莎，舒红，徐正全.东北三省月降水量的时空克里金插值研究［J］.水文，2011，31（3）：31－35.

［12］ 李莎，舒红，徐正全.利用时空Kriging进行气温插值研究［J］.武汉大学学报（信息科学版），2012，37（2）：237－241.

［13］ 张仁铎.空间变异理论及应用［M］.北京：科学出版社，2005. g g g