-
摘要: 在对气象数据进行插值的过程中, 如果只考虑数据的空间信息而忽视数据在时间上的关联, 必然影响插值的精度.针对具有时空特性的气象数据, 提出一种将时空Kriging方法与弹性网方法相结合的新方法.该方法主要利用弹性网算法解决时空Kriging算法中的时空变异函数矩阵为病态矩阵而无法求逆的问题, 通过弹性网算法获得变异函数矩阵方程的稀疏解, 从而提高时空插值的精度.在实际观测的气温数据和AQI数据上的仿真实验验证了该方法对气象时空数据插值的准确性.Abstract: In the process of interpolating meteorological data, if we only consider the spatial information of data and ignore the correlation of data in time, it will inevitably affect the accuracy of interpolation. For the meteorological data with time and space characteristics, this paper combines the spatial-temporal Kriging method with elastic net algorithm. This method uses the elastic net algorithm to solve the problems that the spatial-temporal variational function matrix in the spatial-temporal Kriging algorithm is ill-posed and we can not find its pseudoinverse. The elastic net algorithm is used to obtain the sparse solution of the variational function equation to improve the accuracy of spatial temporal interpolation. The simulation experiments on the observed temperature data and AQI data verify the high accuracy of the proposed method for spatial and temporal data interpolation.
-
Key words:
- Spatial-temporal data /
- spatial-temporal Kriging /
- interpolation /
- elastic net algorithm
-
基于数据的分析方法[1-3]广泛应用于复杂系统建模[4]中, 完备的、准确的数据集是高精度建模的基础.然而在实际观测的气象、交通、环境等领域的数据集中, 存在系统误差、随机误差或者数据缺失、数据异常的情况.
对于数据缺失或异常的问题, 最普遍的解决方法是在数据集中挑选具有连续性的数据子集, 但该方法会对已知数据资源造成极大的浪费.同时, 也会减小数据所记录的某些极端事件的周期, 导致在后续的研究中增加对该极端事件的统计概率[5].还有一种方法就是对存在的数据缺失进行合理的推理得到一个完整的数据集, 即数据插值.常采用的插值方法有多项式插值算法[6]、三次卷积插值[7]、反距离加权算法[8]、普通Kriging插值算法[9-10]等.
仅仅从时间、空间上考虑插值会忽略数据在时间上的趋势性和在空间上的关联性, 从而丢失时空数据集所包含的重要信息, 影响插值精度.因此, 对具有时空特性的数据进行插值的过程中, 需要同时考虑时间趋势性和空间关联性.关于时空插值的方法主要是从时间域或者空间域的方法扩展而来.例如, Li等[11]提出两种时空插值模型, 一种是把时间作为和空间等同的一个维度, 扩展为高维的空间插值.另一种是对空间插值函数和时间插值函数作乘法, 称之为时空乘积插值方法. Antonić等[12]利用具有前向反馈的多层神经网络对克罗地亚的七个气候变量进行时空插值, 通过神经网络训练数据找出变量之间的时间趋势性和空间关联性.
除上述方法外, 还有一大类典型的时空方法是从空间Kriging插值算法衍生而来, 应用较为广泛.例如, 徐爱萍等[13]对空间Kriging算法进行时空扩展, 在时空扩展的过程中把时间看作第三维坐标进行考虑.该方法把时空插值问题转化为高维的空间插值问题, 但忽略了时间和空间各向异性, 影响插值精度. Myers[14]分析考虑各向异性的三个原因:在空间上无法映射时间的独特性质; 在进行逼近和插值的过程中有不同的方向; 各向异性基函数具有部分可微性. De Cesare等[15]针对时空相关性结构提出对于变异函数的不可分离模型—一类积和模型, Iaco等[16]提出一种广义积和模型应用到时空变异函数建模当中, 之后又进一步在文献[17]中针对时空变异函数以及协方差函数给出严格正定的积和模型.
时空Kriging算法广泛应用在气象、水文领域中, 文献[18]对米兰地区的二氧化氮的含量进行时空Kriging插值. Bogaert[19]在其论文中对Ordinary Kriging算法, Cokriging算法以及时空Kriging算法采用交叉验证的方法说明时空Kriging的预测效果更佳. Zeng等[20]对大气二氧化碳含量进行时空Kriging插值. McDaniel等[21]把时空Kriging算法应用到对土壤的甲烷含量和二氧化氮的含量的检测中.李莎等[22]对东北三省月降雨量进行时空Kriging插值研究, 同时与空间Kriging插值算法比较, 插值精度比空间Kriging更高.文献[23]针对时空统计问题, 提出Kriging算法和Kalman滤波相结合, 并得到广泛应用.
上述研究极大地推动了时空数据插值研究的发展, 但是在众多的时空Kriging文献中, 关注的重点集中在利用优化或改进时空变异函数模型来提高时空插值精度, 针对时空Kriging算法求解权重系数过程中存在的变异函数矩阵具有病态特性的问题研究较少.变异函数矩阵的病态特性是由于该矩阵的列向量的线性相关性过大, 表示的特征太过于相似以至于容易混淆所产生的.病态矩阵问题在很多领域都存在, 解决这一问题的方法有很多, 例如Lasso回归[24]和岭回归[25]等. Lasso回归通过施加一范数惩罚稀疏化解空间, 岭回归算法通过损失求解的无偏性来换取解的稳定性.弹性网算法[26-28]结合岭回归和Lasso回归的优点, 既能提高解的稳定性, 又能得到稀疏的解, 故本文引入弹性网算法解决由于时空变异函数矩阵病态导致无法求逆的问题.通过求解时空变异函数矩阵方程, 算出时空域样本点对应的权重系数, 进而计算得到较为精确的时空插值结果.
1. 时空Kriging算法
时空变量是空间向时空的扩展, 对于时空变量的分析与空间变量的分析大体类似.定义$ {D^R} $表示空间域度, $ R $表示空间维度(本文以二维空间为例), $ T $表示时间域度, 若存在$ s = (x, y) \in {D^R}, t \in T $, 那么$ z(s, t) $表示时空随机变量.
1.1 时空数据平稳性
时空Kriging方法是建立在时空数据平稳性这一基础之上的[22], 因此在进行时空数据分析前要进行平稳性分析, 平稳性分析的方法包括平稳假设、二阶平稳假设、本征假设.
平稳假设:在时空域内的时空点集内, 从时空域的某一时空位置移动到另一时空位置, 其时空随机变量的数学期望为定值.
二阶平稳假设:在时空域内, 时空随机变量的协方差对于任意的时空位置都存在且平稳.
本征假设:时空随机变量增量的数学期望对任意时空随机变量都存在并且为0.同时, 时空变量增量的方差对任意时空随机变量都存在且平稳.
1.2 时空变异函数
时空变异函数如式(1)所示
$$ \begin{align} &{\gamma ^*}({h_s}, {h_t}) = \frac{1}{{2n({h_s}, {h_t})}} \times\\ &\qquad\sum\limits_{i = 1}^{n({h_s}, {h_t})} {{{(z({s_i}, {t_i})-z({s_i}+{h_s}, {t_i}+{h_t}))}^2}} \end{align} $$ (1) 式(1)中$ {h_s} $、$ {h_t} $为对应的空间间隔距离和时间间隔距离; $ {z({s_i}, {t_i})} $为样本点$ ({s_i}, {t_i}) $的观测值, $ {z({s_i} + {h_s}, {t_i} + {h_t})} $是与$ z({s_i}, {t_i}) $空间距离相隔$ {h_s} $、时间距离相隔$ {h_t} $的样本观测值; $ {n({h_s}, {h_t})} $是时空分隔距离分别为$ {h_s} $、$ {h_t} $的样本点对总数.
与式(1)对应的时空协方差函数如式(2)所示
$$ \begin{equation} Cov[z(s, t), z(s + {h_s}, t + {h_t})] = C({h_s}, {h_t}) \end{equation} $$ (2) 其中$ C $为协方差的简记形式.求取时空域中未知点观测值计算公式如式(3)所示
$$ \begin{equation} {z^*}({s_0}, {t_0}) = \sum\limits_{i = 1}^n {{\lambda _i}z({s_i}, {t_i})} \end{equation} $$ (3) 式(3)中$ {\lambda _i} $为第$ i $个已知样本点对应的权重系数.其中权重系数是能够满足点$ ({s_0}, {t_0}) $处的估计值与真实值的差最小的系数.
由无偏估计的条件$ E(z_0^ * - {z_0}) = 0 $可以推导得
$$ \begin{equation} \sum\limits_{i = 1}^N {{\lambda _i} = } 1 \end{equation} $$ (4) 式(4)中$ N $表示已知样本点的个数.
求解权重系数即求解式(5)所表示的条件函数
$$ \begin{align} {\min \left\{ {Var(z_0^ * ({s_0}, {t_0}) - {z_0}({s_0}, {t_0}))} \right\}} {\sum\limits_{i = 1}^N {{\lambda _i} = 1} } \end{align} $$ (5) 通过引入拉格朗日系数$ \varphi $进行推导, 得到如下矩阵形式.
$$ \begin{align} &\begin{bmatrix} {{\gamma _{11}}({h_s}, {h_t})}&{{\gamma _{12}}({h_s}, {h_t})}& \cdots &{{\gamma _{1n}}({h_s}, {h_t})}&1\\ {{\gamma _{12}}({h_s}, {h_t})}&{{\gamma _{22}}({h_s}, {h_t})}& \cdots &{{\gamma _{2n}}({h_s}, {h_t})}&1\\ \vdots & \vdots & \ddots& \vdots &1\\ {{\gamma _{1n}}({h_s}, {h_t})}&{{\gamma _{2n}}({h_s}, {h_t})}& \cdots &{{\gamma _{nn}}({h_s}, {h_t})}&1\\ 1&1& \cdots &1&0 \end{bmatrix}\cdot\\&\qquad \begin{bmatrix} {{\lambda _1}}\\ {{\lambda _2}}\\ \vdots \\ {{\lambda _n}}\\ \varphi \end{bmatrix} = \begin{bmatrix} {{\gamma _{01}}({h_s}, {h_t})}\\ {{\gamma _{02}}({h_s}, {h_t})}\\ \vdots \\ {{\gamma _{0n}}({h_s}, {h_t})}\\ 1 \end{bmatrix} \end{align} $$ (6) 1.3 时空变异函数模型
目前, 关于时空变异函数模型的研究主要有两大类:可分离模型和不可分离模型.可分离模型较为简单, 主要是对时间变异函数和空间变异函数以加法或者乘法的形式进行组合得到, 不能够很好表达出时间和空间的相关信息; 而不可分离模型构建过程相对复杂, 但能够更加有效地表现出变量在时间和空间的相关关系.其中一类积和式模型是常用的不可分离模型[22], 表示如下.
$$ \begin{align} {\gamma _{st}}({h_s}, {h_t}) = & [({k_1}{C_t}(0) + {k_2}]{\gamma _s}({h_s})+\\& [{k_1}{C_s}(0) + {k_3}]{\gamma _t}({h_t})-\\& {k_1}{\gamma _s}({h_s}){\gamma _t}({h_t}) \end{align} $$ (7) 与式(7)对应的时空协方差变异函数模型如下式所示.
$$ \begin{align} {C_{st}}({h_s}, {h_t}) = & {k_1}{C_s}({h_s}){C_t}({h_t})+\\& {k_2}{C_s}({h_s}) + {k_3}{C_t}({h_t}) \end{align} $$ (8) $$ \begin{equation} \begin{cases} {k_1} = \dfrac{ [{C_s}(0) + {C_t}(0) - {C_{st}}(0, 0)]}{[{C_s}(0) \cdot {C_t}(0)]}\\ {k_2} = \dfrac{ [{C_{st}}(0, 0) - {C_t}(0)]}{{C_s}(0)}\\ {k_3} = \dfrac{ [{C_{st}}(0, 0) - {C_s}(0)]}{{C_t}(0)} \end{cases} \end{equation} $$ (9) 其中, $ {\gamma _{st}}({h_s}, {h_t}) $、$ {\gamma _s}({h_s}) $、$ {\gamma _t}({h_t}) $分别为对应的时空变异函数、空间变异函数和时间变异函数; $ {C_{st}}(0, 0) $、$ {C_s}(0) $、$ {C_t}(0) $分别为时空变异函数、空间变异函数和时间变异函数对应的基台值, 这3个量可用时空变异函数式(1)进行估计.
1.3.1 积和式时空变异函数模型推导
根据时空Kriging算法的性质可以得到$ {\gamma _{st}}(0, 0) = {\gamma _t}(0) = {\gamma _s}(0) = 0 $, 则由式(7)得
$$ \begin{equation} {\gamma _{st}}({h_s}, 0) = [({k_1}{C_t}(0) + {k_2}]{\gamma _s}({h_s}) \end{equation} $$ (10) $$ \begin{equation} {\gamma _{st}}(0, {h_t}) = [{k_1}{C_s}(0) + {k_3}]{\gamma _t}({h_t}) \end{equation} $$ (11) 设$ {\gamma _{st}}({h_s}, 0) = {\gamma _s}({h_s}) $, $ {\gamma _{st}}(0, {h_t}) = {\gamma _t}({h_t}) $, 则由式(10)、式(11)得
$$ \begin{equation} {k_1}{C_t}(0) + {k_2} = 1 \end{equation} $$ (12) $$ \begin{equation} {k_1}{C_s}(0) + {k_3} = 1 \end{equation} $$ (13) 由式(8)可以得到
$$ \begin{equation} {C_{st}}(0, 0) = {k_1}{C_s}(0){C_t}(0) + {k_2}{C_s}(0) + {k_3}{C_t}(0) \end{equation} $$ (14) 由式(12)$ \sim $(14)联立即可解得式(9).
这类积和式模型的优点是$ {\gamma _{st}}({h_s}, {h_t}) $只由空间变异函数$ {\gamma _s}({h_s}) $、时间变异函数$ {\gamma _t}({h_t}) $、以及时空基台值$ {C_{st}}(0, 0) $决定.上述三个量都可以由式(1)进行估计得到.
1.3.2 时空Kriging算法步骤
空间Kriging算法只能对某一时间点上的空间区域观测值进行插值, 而对于任意时间和任意空间点观测值插值, 则需要用到时空Kriging算法, 其算法步骤总结如下.
1) 由式(1)计算时空变异函数值.
2) 计算空间变异函数值和时间变异函数值, 分别得到最佳变异函数拟合曲线.
3) 根据式(9)分别计算得到$ {k_1} $、$ {k_2} $、$ {k_3} $, 根据式(7)得到时空变异函数积和模型.
4) 根据式(6)求解样本点对应的权重系数, 根据式(3)得到时空插值结果.
1.3.3 时空Kriging算法复杂度分析
设$ N $为样本个数, $ p $为待插值点的个数.计算变异函数矩阵需要进行两次循环次数为$ N $的循环, 则其计算复杂度为$ {\rm O}({N^2}) $. $ m $为空间变异函数散点个数, $ n $为时间变异函数散点个数.利用最小二乘算法拟合空间变异函数曲线和时间变异函数曲线, 其计算复杂度为$ {\rm O}({m^3} + {n^3}) $.最后代入式(9), 计算时空模型的系数, 并代入时空模型中得到时空模型表达式, 其计算复杂度为$ {\rm O}(1) $.计算待插值点与已知点之间的变异函数值, 其复杂度为$ {\rm O}(pN) $.对时空变异函数矩阵$ {{\bf{\gamma }}_{(N + 1) \times (N + 1)}} $求逆, 其对应的计算复杂度为$ {\rm O}({(N + 1)^3}) $.求解对应的权重系数, 计算复杂度为$ {\rm O}({(N + 1)^3}) $.对待插值点进行插值, 需要迭代$ p $次, 则其计算复杂度为$ {\rm O}(p{N^3}) $.综上, 时空Kriging的总计算复杂度为$ {\rm O}(2{(N + 1)^3} + p{N^3} + {m^3} + {n^3}) $.
2. 弹性网算法
本文中由式(3)求解时空插值时, 需要得到已知样本点对应的权值系数$ {\lambda _i} $, 即通过式(6)求解.在求变异函数矩阵的条件数时, 发现式(6)中的变异函数值矩阵大多情况下为病态矩阵, 无法直接求逆.
针对该问题本文提出使用弹性网算法(Elastic net algorithm)[26]求解式(6).
弹性网算法是一种回归分析的方法, 该算法以最小二乘回归算法作为基础, 目标函数中施加$ L_2 $范数惩罚项和$ L_1 $范数惩罚项的线性组合.对于式(6)记作为$ {A\pmb{\lambda}} = {\pmb{Y} } $, 则该式的求解过程可以转化为以下形式:
$$ \begin{equation} L({\beta _1}, {\beta _2}, {\pmb{\lambda}}) = {\left| {{\pmb{Y}} - {A\pmb{\lambda}}} \right|^2} + {\beta _2}{\left| {\pmb{\lambda}} \right|^2} + {\beta _1}{\left| {\pmb{\lambda}} \right|_1} \end{equation} $$ (15) 其中
$$ {\left| {\pmb{\lambda}} \right|^2} = \sum\limits_{j = 1}^{n + 1} {\lambda _j^2} $$ $$ \begin{equation} {\left| {\pmb{\lambda}} \right|_1} = \sum\limits_{j = 1}^{n + 1} {\left| {\lambda _j^{}} \right|} \end{equation} $$ (16) $$ {\hat {\pmb \lambda} } = \min \left\{ {L({\beta _1}, {\beta _2}, {\pmb{\lambda}})} \right\} $$ 上述的过程中设$ \alpha = {\beta _2}/\left({{\beta _1} + {\beta _2}} \right) $, 则式(16)可以等价为下式所示
$$ \begin{equation} \begin{cases} {{\hat {\pmb \lambda}} = \min {\rm{\{ }}{{\left| {{\pmb{Y}} - {A\pmb{\lambda}}} \right|}^2}{\rm{\} }}}\\ {(1 - \alpha ){{\left| {\pmb{\lambda}} \right|}^2} + \alpha {{\left| {\pmb{\lambda}} \right|}_1} \le t, \forall t} \end{cases} \end{equation} $$ (17) 式(17)中$ (1 - \alpha){\left| \pmb{\lambda} \right|^2} + \alpha {\left| \pmb{\lambda} \right|_1} $称为弹性网惩罚项, 是岭回归惩罚项和Lasso回归惩罚项的凸组合, 当$ \alpha $为0时, 式(17)等价于岭回归.对于所有的$ \alpha \in \left({0, 1} \right) $, 弹性网惩罚项为严格凸, 所以弹性网有岭回归和Lasso两者的特征.当$ \alpha $为1时, 式(17)等价于Lasso回归, 此时Lasso惩罚项凸但不是严格凸.
在时空Kriging算法中, 式$ {A\pmb{\lambda} = \pmb Y} $由于矩阵$ {A} $一般为病态的, 得到的$ {\pmb{\lambda}} $值可能会因为$ \pmb Y $的一些微小变化而变化巨大.弹性网算法首先通过Lasso惩罚项进行稀疏规则化, 然后通过岭回归惩罚项保证解的稳键性, 进而提高时空Kriging插值的精度.
3. 仿真实例
为了验证本文所描述的时空Kriging算法(Spatial-temporal Kriging, STKriging)的有效性, 采用两组时空气象数据包括一组气温时空数据和一组空气污染指标AQI时空数据进行仿真实验.为验证本文所提的积和变异函数模型的有效性, 将其与考虑空间信息进行插值的空间Kriging算法[9]、考虑时空信息的时空自相关移动平均算法[29] (Spatial-temporal autoregressive and moving average, STARMA)、利用空间权重矩阵衡量空间关系的时空回声状态网络算法(Spatial-temporal echo state network, STESN)、基于时间信息进行预测的回声状态高斯过程算法(Echo state Gaussian process, ESGP)[30]等进行比较.插值效果的衡量指标为均方根误差(Root mean squared error, RMSE), 标准化均方根误差(Normalized root mean square error, NRMSE)以及平均对称绝对百分率误差(Symmetric mean absolute percentage error, SMAPE), 其对应的计算公式如下所示:
$$ \begin{equation} {\rm{RMSE }} = \sqrt {\sum\limits_{i = 1}^n {\frac{{[{{\hat y}_i} - {y_i}]^2}}{{n - 1}}} } \end{equation} $$ (18) $$ \begin{equation} {\rm{SMAPE}} = \frac{1}{n}\sum\limits_{i = 1}^n {\frac{{2\left| {{y_i} - {{\hat y}_i}} \right|}}{{\left| {{{\hat y}_i} + {y_i}} \right|}}} \end{equation} $$ (19) $$ \begin{equation} {\rm{NRMSE = }}\sqrt {\sum\limits_{i = 1}^n {\frac{{{{[{{\hat y}_i} - {y_i}]}^2}}}{{\sum\limits_{i = 1}^n {[{y_i} - \bar y} {]^2}}}} } \end{equation} $$ (20) 其中$ y_i $表示实际观测值, $ \bar y $表示实际观测值的平均值, $ \hat y_i $表示插值的数值, $ n $表示样本数.
3.1 气温数据仿真实例
采取2016-07-01到2017-06-31杭州、东台、安庆、定海、嵊泗、景德镇、衢州、瑞安、徐州、淮安、射阳、蚌埠、南京等13个城市每天15时的气温所构成4 758个样本进行仿真实验.对13个观测点的经纬度进行记录, 将其转化成直角坐标系坐标, 4 758个样本值是13个观测地点的气温按照时间顺序排列而成, 即为
$$ {T_{\text{杭州}}}(1), \cdots, \;{T_{\text{南京}}}(1), \;{T_{\text{杭州}}}(2), \;{T_{\text{东台}}}(2), \cdots, \; $$ $$ {T_{\text{南京}}}(2), \cdots, \;{T_{\text{杭州}}}(366), \cdots, \;{T_{\text{南京}}}(366) $$ 在数据集中取连续的200个样本值, 根据这200个样本值计算时空变异函数并对该时空变异函数值建立积和模型, 利用该模型并结合弹性网算法对后续第201到213个样本值进行时空插值, 根据数据集的排列顺序选择200个连续的样本值使得第201到213共13个时空点为同一天13个观测点的观测值, 如$ \; {T_{\text{杭州}}}(i), \; {T_{\text{东台}}}(i), \cdots, \; {T_{\text{南京}}}(i) $形式.同理, 对5月31日到6月29日共计30天的13个观测点每天15时气温进行时空插值.为了方便和其余三种算法进行对比, 计算误差指标时选取13个观测点中射阳30天的插值结果及其对应的观测值计算.
3.1.1 气温时空变异函数模型
本文变异函数模型为一类积和模型, 如式(7)所示, 积和模型中用到空间变异函数模型和时间变异函数模型. 图 1表示的是空间变异函数值和空间变异函数拟合模型. 图 2表示的是时间变异函数值和时间变异函数拟合模型. 图 1空间变异函数空间距离的确定根据观测点位置相互之间的距离, 取最小的距离作为空间距离, 图 2时间变异函数时间距离一般根据数据之间的最小的跨度确定, 在本例中为1天.
空间变异函数模型数学表达式如下式所示:
$$ \gamma ({h_s}) = - 0.05763 \times \exp ( - 0.003449 \times {h_s}) + 0.05779{\rm{ }} $$ 时间变异函数模型数学表达式如下式所示:
$$ {\rm{ }}\gamma ({h_t}) = - 0.005204 \times \exp ( - 0.4046 \times {h_t}) + 0.005385 $$ 图 3表示根据式(1)计算所得气温时空变异函数值.
通过图 1 $ \sim $ 图 3估计对应的空间变异函数基台值、时间变异函数基台值以及时空变异函数基台值.根据积和模型公式可以得到$ {k_1} = 60 $, $ {k_2} = 0.1 $, $ {k_3} = 0.4 $, 然后根据上文中确定的空间变异函数模型以及时间变异函数模型, 代入式(7)即可得到气温时空变异函数模型.
3.1.2 气温时空插值结果
对时空域中的射阳观测点的连续30天的数据进行时空插值, 根据式(7)所示的积和模型得到对气温数据的时空插值, 得到的插值效果图如图 4所示.从图中可以看出通过时空Kriging算法得到的插值结果与观测值误差在±5℃以内.
图 4 STKriging对射阳市气温数据插值效果Fig. 4 Interpolation curves of STKriging on temperature data of Sheyang City为了验证本文提出的与弹性网算法相结合的时空Kriging算法的优越性, 通过空间Kriging算法、STARMA算法、STESN算法、ESGP算法对射阳共30天的15时气温进行预测, 其中空间Kriging算法分别使用每一天的12个观测点的观测值对射阳15时气温进行插值, 进行30次即可得到30天的预测值. ESGP算法只考虑时间信息, 通过射阳30天之前的观测值进行单步预测30天的气温. STARMA算法以及STESN算法虽然利用空间权重矩阵考虑不同地点彼此之间的空间关系, 但是由于只是简单的利用距离因素衡量空间关联性, 不能较为准确地反映出彼此之间的空间关系.各算法预测误差如表 1所示.通过表 1所示四种算法计算指标即RMSE、NRMSE、SMAPE三种指标可以得出STKriging算法的误差最小, 插值精度最高.需要特别说明的是, STESN模型的SMAPE指标较大, 这是由于在某些点预测值与真实值之和接近于零, 导致分母较大造成的.
表 1 气温数据仿真结果比较Table 1 Experimental results for temperature data3.2 AQI数据仿真实例
采取山东省2015年1月到2016年12月共两年13个观测站点的AQI日值所构成9 516个样本进行仿真实验.对13个观测点的经纬度进行记录, 将其转化成直角坐标系坐标, 9 516个样本是同一时间的13个观测点按照时间顺序排列.在数据集中取连续的200个样本值, 根据这200个样本值计算时空变异函数并对该时空变异函数值建立积和模型, 利用该模型并结合弹性网算法对后续第201到213个样本值进行时空插值, 根据数据集的排列顺序选择200个连续的样本值使得第201到213共13个时空点为同一天13个观测点的观测值.同理, 和气温数据处理方法相同, 为与其他方法进行对比, 本文选取13个观测点中第二个观测点(位于青岛市市南区香港中路)共计30天的AQI插值结果以及对应的实际观测值计算插值误差.
3.2.1 AQI时空变异函数模型
如式(7)所示, 一类积和模型分别包括时间变异函数模型和空间变异函数模型, 下面是对该数据集求得如图 5所示的时间变异函数和对应的时间变异函数拟合曲线.时间变异函数的时间距离为1天.空间变异函数的空间距离根据观测点之间的最小距离确定.
对应的时间变异函数模型的数学表达式如下式所示:
$$ {\rm{ }}\gamma ({h_t}) = - 0.003133 \times \exp ( - 0.322 \times {h_t}) + 0.003105{\rm{ }} $$ 图 6为空间变异函数和对应的空间变异函数拟合曲线.
对应的空间变异函数模型的数学表达式如下式所示:
$$ \gamma ({h_s}) = 0.00146 \times \exp (0.5178 \times {h_s}) $$ 图 7表示根据式(1)计算所得实验AQI时空变异函数值.
通过图 5$ \sim\! $图 7估计对应的时间变异函数基台值、空间变异函数基台值以及时空变异函数基台值, 根据本文所使用的积和模型公式得$ {k_1} = -6.72\times10^{-4} $, $ {k_2} = 1.16 $, $ {k_3} = 1.1905 $.然后根据上文中确定的空间变异函数模型以及时间变异函数模型, 代入式(7)即可得到AQI时空变异函数模型.
3.2.2 AQI时空插值结果
对时空域中第2个观测点(位于青岛市市南区香港中路)连续30天的数据进行时空插值, 根据式(7)所示的积和模型得到对AQI数据进行时空插值, 得到插值效果图如图 8所示.
图 8 STKriging对青岛市AQI数据插值效果Fig. 8 Interpolation curves of STKriging on AQI data of Qingdao City为了验证本文提出的与弹性网算法相结合的时空Kriging算法的优越性, 通过空间Kriging算法、STARMA算法、STESN算法、ESGP算法同样对13个观测站中第二个观测站(位于青岛市市南区香港中路)的共30天的AQI值进行预测, 得到各个算法误差指标如表 2所示.
表 2 AQI数据仿真结果比较Table 2 Experimental results for AQI data通过表 2可以看出, 时空Kriging对时空数据的时空插值精度高于其他算法.空间Kriging算法只考虑了空间信息, ESGP算法只考虑时间信息, STARMA算法以及STESN算法只简单利用距离作为衡量空间关联性的标准, 并不能较为准确地衡量空间关联性.然而, 本文提出的时空Kriging算法同时考虑时间和空间的相关性, 故插值精度更高.
4. 结论
为解决气象数据观测值缺失或者异常的问题, 本文提出一种新型的时空Kriging算法.该算法既考虑时间信息又考虑空间要素, 对数据进行时空插值.将弹性网算法与一般时空Kriging算法相结合, 解决变异函数矩阵为病态矩阵而影响时空插值精度的问题.基于两组实际观测的气象时空数据的仿真实验验证了本文所提方法的有效性, 其比单独考虑空间信息的空间Kriging算法、单纯考虑时间信息的ESGP算法、考虑时空信息的STARMA算法以及STESN算法的精度都要高, 具有良好的应用前景.
-
图 4 STKriging对射阳市气温数据插值效果
Fig. 4 Interpolation curves of STKriging on temperature data of Sheyang City
图 8 STKriging对青岛市AQI数据插值效果
Fig. 8 Interpolation curves of STKriging on AQI data of Qingdao City
-
[1] 胡寿松, 张正道.基于神经网络的非线性时间序列故障预报.自动化学报, 2007, 33(7): 744-748 doi: 10.1360/aas-007-0744Hu Shou-Song, Zhang Zheng-Dao. Fault prediction for nonlinear time series based on neural network. Acta Automatica Sinica, 2007, 33(7): 744-748 doi: 10.1360/aas-007-0744 [2] 伦淑娴, 胡海峰.基于罚函数内点法的泄露积分型回声状态网的参数优化.自动化学报, 2017, 43(7): 1160-1168 doi: 10.16383/j.aas.2017.c160541Lun Shu-Xian, Hu Hai-Feng. Parameter optimization of leaky integrator echo state network with internal-point penalty function method. Acta Automatica Sinica, 2017, 43(7): 1160-1168 doi: 10.16383/j.aas.2017.c160541 [3] 苏本跃, 蒋京, 汤庆丰, 盛敏.基于函数型数据分析方法的人体动态行为识别.自动化学报, 2017, 43(5): 866-876 doi: 10.16383/j.aas.2017.c160120Su Ben-Yue, Jiang Jing, Tang Qing-Feng, Sheng Min. Human dynamic action recognition based on functional data analysis. Acta Automatica Sinica, 2017, 43(5): 866-876 doi: 10.16383/j.aas.2017.c160120 [4] 宋贺达, 周平, 王宏, 柴天佑.高炉炼铁过程多元铁水质量非线性子空间建模及应用.自动化学报, 2016, 42(11): 1664-1679 doi: 10.16383/j.aas.2016.c150819Song He-Da, Zhou Ping, Wang Hong, Chai Tian-You. Nonlinear subspace modeling of multivariate molten iron quality in blast furnace ironmaking and its application. Acta Automatica Sinica, 2016, 42(11): 1664-1679 doi: 10.16383/j.aas.2016.c150819 [5] Di Piazza A, Lo Conti F, Noto L V, Viola F, La Loggia G. Comparative analysis of different techniques for spatial interpolation of rainfall data to create a serially complete monthly time series of precipitation for Sicily, Italy. International Journal of Applied Earth Observation and Geoinformation, 2011, 13(3): 396-408 doi: 10.1016/j.jag.2011.01.005 [6] Ly S, Charles C, Degré A. Geostatistical interpolation of daily rainfall at catchment scale: The use of several variogram models in the Ourthe and Ambleve catchments, Belgium. Hydrology and Earth System Sciences, 2011, 15(7): 2259-2274 doi: 10.5194/hess-15-2259-2011 [7] 周登文, 申晓留.边导向的双三次彩色图像插值.自动化学报, 2012, 38(4): 525-530 doi: 10.3724/SP.J.1004.2012.00525Zhou Deng-Wen, Shen Xiao-Liu. Edge-directed bicubic color image interpolation. Acta Automatica Sinica, 2012, 38(4): 525-530 doi: 10.3724/SP.J.1004.2012.00525 [8] Wang Y, Akeju O V, Zhao T Y. Interpolation of spatially varying but sparsely measured geo-data: A comparative study. Engineering Geology, 2017, 231: 200-217 doi: 10.1016/j.enggeo.2017.10.019 [9] Jeffrey S J, Carter J O, Moodie K B, Beswick A R. Using spatial interpolation to construct a comprehensive archive of Australian climate data. Environmental Modelling & Software, 2001, 16(4): 309-330 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=ccce6481f2b1c81e50171ca87c1defb5 [10] Xiao M Y, Zhang G H, Breitkopf P, Villon P, Zhang W H. Extended co-Kriging interpolation method based on multi-fidelity data. Applied Mathematics and Computation, 2018, 323: 120-131 doi: 10.1016/j.amc.2017.10.055 [11] Li L X, Revesz P. Interpolation methods for spatio-temporal geographic data. Computers, Environment and Urban Systems, 2004, 28(3): 201-227 doi: 10.1016/S0198-9715(03)00018-8 [12] Antonić O, Križan J, Marki A, Bukovec D. Spatio-temporal interpolation of climatic variables over large region of complex terrain using neural networks. Ecological Modelling, 2001, 138(1-3): 255-263 doi: 10.1016/S0304-3800(00)00406-3 [13] 徐爱萍, 胡力, 舒红.空间克里金插值的时空扩展与实现.计算机应用, 2011, 31(1): 273-276 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jsjyy201101074Xu Ai-Ping, Hu Li, Shu Hong. Extension and implementation from spatial-only to spatiotemporal Kriging interpolation. Journal of Computer Applications, 2011, 31(1): 273-276 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=jsjyy201101074 [14] Myers D E. Anisotropic radial basis functions. International Journal of Pure and Applied Mathematics, 2008, 42(2): 197-203 [15] De Cesare L, Myers D E, Posa D. Estimating and modeling space-time correlation structures. Statistics & Probability Letters, 2001, 51(1): 9-14 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=6f7474ad57568e5108834d108c0752a5 [16] De Iaco S, Myers D E, Posa D. Space-time analysis using a general product — Sum model. Statistics & Probability Letters, 2001, 52(1): 21-28 https://www.sciencedirect.com/science/article/pii/S0167715200002005 [17] De Iaco S, Myers D E, Posa D. On strict positive definiteness of product and product — Sum covariance models. Journal of Statistical Planning and Inference, 2011, 141(3): 1132-1140 doi: 10.1016/j.jspi.2010.09.014 [18] De Cesare L, Myers D E, Posa D. Product-sum covariance for space-time modeling: An environmental application. Environmetrics, 2001, 12(1): 11-23 doi: 10.1002/1099-095X(200102)12:1<11::AID-ENV426>3.0.CO;2-P [19] Bogaert P. Comparison of Kriging techniques in a space-time context. Mathematical Geology, 1996, 28(1): 73-86 doi: 10.1007/BF02273524 [20] Zeng Z C, Lei L P, Hou S S, Ru F, Guan X H, Zhang B. A regional gap-filling method based on spatiotemporal variogram model of CO$_2$ columns. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(6): 3594-3603 doi: 10.1109/TGRS.2013.2273807 [21] McDaniel M D, Simpson R G, Malone B P, McBratney A B, Minasny B, Adams M A. Quantifying and predicting spatio-temporal variability of soil CH$_4$ and N$_2$O fluxes from a seemingly homogeneous Australian agricultural field. Agriculture, Ecosystems & Environment, 2017, 240: 182-193 [22] 李莎, 舒红, 徐正全.东北三省月降水量的时空克里金插值研究.水文, 2011, 31(3): 31-35 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=sw201103007Li Sha, Shu Hong, Xu Zheng-Quan. Study on spatial-temporal Kriging interpolation of monthly precipitation in three provinces of northeast China. Journal of China Hydrology, 2011, 31(3): 31-35 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=sw201103007 [23] 杜宇健, 萧德云. Kriging和Kalman技术在时空统计上的应用.见: 2004中国控制与决策学术年会论文集, 2004. 197-200, 205Du Yu-Jian, Xiao De-Yun. Application of Kriging and Kalman filter in spatial temporal statistics. In: Proceeding of the 2004 Chinese Control and Decision Conference, 2004. 197-200, 205 [24] 张桂梅, 孙晓旭, 刘建新, 储珺.基于分数阶微分的TV-$L^1$光流模型的图像配准方法研究.自动化学报, 2017, 43(12): 2213-2224Zhang Gui-Mei, Sun Xiao-Xu, Liu Jian-Xin, Chu Jun. Research on TV-$L^1$ optical flow model for image registration based on fractional-order differentiation. Acta Automatica Sinica, 2017, 43(12): 2213-2224 [25] 张文林, 张连海, 牛铜, 屈丹, 李弼程.基于正则化的本征音说话人自适应方法.自动化学报, 2012, 38(12): 1950-1957 doi: 10.3724/SP.J.1004.2012.01950Zhang Wen-Lin, Zhang Lian-Hai, Niu Tong, Qu Dan, Li Bi-Cheng. Regularization based eigenvoice speaker adaptation method. Acta Automatica Sinica, 2012, 38(12): 1950-1957 doi: 10.3724/SP.J.1004.2012.01950 [26] Liu Q S, Sun Y B, Wang C T, Liu T L, Tao D C. Elastic net hypergraph learning for image clustering and semi-supervised classification. IEEE Transactions on Image Processing, 2017, 26(1): 452-463 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=f4419e755c4bb0d167e2cc8b4990b033 [27] You C, Li C G, Robinson D P, Vidal R. Oracle based active set algorithm for scalable elastic net subspace clustering. In: Proceedings of the 2016 IEEE Conference on Computer Vision and Pattern Recognition. Las Vegas, NV, USA: IEEE, 2016. 3928-3937 [28] Zhang Z, Lai Z H, Xu Y, Shao L, Wu J, Xie G S. Discriminative elastic-net regularized linear regression. IEEE Transactions on Image Processing, 2017, 26(3): 1466-1481 doi: 10.1109/TIP.2017.2651396 [29] Pavlyuk D. Short-term traffic forecasting using multivariate autoregressive models. Procedia Engineering, 2017, 178: 57-66 doi: 10.1016/j.proeng.2017.01.062 [30] Chatzis S P, Demiris Y. Echo state Gaussian process. IEEE Transactions on Neural Networks, 2011, 22(9): 1435-1445 doi: 10.1109/TNN.2011.2162109 期刊类型引用(5)
1. 詹兆康，胡旭光，赵浩然，张思琪，张峻凯，马大中. 基于多变量时空融合网络的风机数据缺失值插补研究. 自动化学报. 2024(06): 1171-1184 . 
本站查看2. 宋阿伟，田红，王胜，刘前，谢五三，唐为安，戴娟，丁小俊，吴蓉. 基于风雨综合指数的安徽省台风灾害房屋风险评估方法. 暴雨灾害. 2024(03): 363-370 . 百度学术3. 许文洁，霍亮，沈涛，陈壮，赵博林. 地面沉降时空形变场的构建方法研究. 测绘科学. 2022(01): 204-211 . 百度学术4. 谭理庆，彭琦，曹阳，杨鑫，唐帅，刘俊. 不同轨道类型LEO卫星轨道拟合及预报精度研究. 全球定位系统. 2022(02): 44-51 . 百度学术5. 周特军，陆佳政，吴传平，李波，谭艳军，刘毓. 输电线路山火直升机喷洒灭火试验与应用. 高电压技术. 2021(08): 2811-2819 . 百度学术其他类型引用(6)
-
下载:
下载:
计量
- 文章访问数: 1424
- HTML全文浏览量: 286
- PDF下载量: 238
- 被引次数: 11
