EN中文

地球信息科学学报 >

2017 , Vol. 19 >Issue 8: 1001 - 1010

DOI: https://doi.org/10.3724/SP.J.1047.2017.01001

地球信息科学理论与方法

基于多尺度最小二乘支持向量机优化的克里金插值方法

  • 车磊 , 1, 2 ,
  • 王海起 , 1, * ,
  • 费涛 1 ,
  • 闫滨 1 ,
  • 刘玉 1 ,
  • 桂丽 1 ,
  • 陈冉 1 ,
  • 翟文龙 1
展开
  • 1. 中国石油大学（华东）地球科学与技术学院,青岛 266580
  • 2. 中国电波传播研究所青岛分所,青岛 266107
*通信作者：王海起（1972-）,男,博士,副教授,研究方向为地理信息科学。E-mail: wanghaiqi@upc.edu.cn

作者简介：车 磊（1990-）,男,硕士生,研究方向为地图制图学与地理信息工程。E-mail:

收稿日期: 2016-10-31

  要求修回日期: 2017-06-20

  网络出版日期: 2017-08-20

基金资助

国家自然科学基金项目(41471322)

山东省自然科学基金项目(ZR2012DM010)

收起

Kriging Interpolation Method Optimized by Multi-scale Least Squares Support Vector Machine

  • CHE Lei , 1, 2 ,
  • WANG Haiqi , 1, * ,
  • FEI Tao 1 ,
  • YAN Bin 1 ,
  • LIU Yu 1 ,
  • GUI Li 1 ,
  • CHEN Ran 1 ,
  • ZHAI Wenlong 1
Expand
  • 1. School Of Geosciences, China University of Petroleum (East China), Qingdao 266580, China
  • 2. China Research Institute of Radiowave Propagation Qingdao Branch, Qingdao 266107, China;
*Corresponding author: WANG Haiqi, E-mail:

Received date: 2016-10-31

  Request revised date: 2017-06-20

  Online published: 2017-08-20

Copyright

《地球信息科学学报》编辑部 所有
Fold

摘要

克里金插值方法根据待估位置点、已知样本数据点的位置关系和区域化变量的空间相关性,实现空间加权估计,满足估计的无偏性和最优性。传统方法理论模型形状固定且选择具有人为主观性,无法反映空间数据的变化趋势及其空间多尺度特征。本文为解决上述问题,提出了一种基于多尺度最小二乘支持向量机优化的克里金插值方法,此方法为拟合实验变异函数提供了一种新的思路。从实际样本数据的变化趋势出发,采用最小二乘支持向量机拟合实验变异函数,并利用不同尺度小波核反映不同尺度下的空间变化。最后,实验环节包括模拟和应用,模拟主要验证经多尺度最小二乘支持向量机优化后插值方法的科学有效性以及准确性,应用主要研究青岛市PM2.5浓度时空分布特征,为城市生态科学防护及控制提供理论依据。结果表明,基于多尺度最小二乘支持向量机优化的克里金插值方法能够更好地刻画变异函数,反映不同尺度下的空间变化细节,从而在一定程度上提高插值的精度,是一种可选的克里金插值方法。

关键词： 克里金插值; 最小二乘支持向量机; 变异函数; 多尺度; 小波核函数

本文引用格式

车磊 , 王海起 , 费涛 , 闫滨 , 刘玉 , 桂丽 , 陈冉 , 翟文龙 . 基于多尺度最小二乘支持向量机优化的克里金插值方法[J]. 地球信息科学学报, 2017 , 19(8) : 1001 -1010 . DOI: 10.3724/SP.J.1047.2017.01001

Abstract

Kriging interpolation method realizes spatial weighted estimation that meets the unbiasedness and optimality according to the position relationship between the estimated location sites and the known sample sites and regionalized variable spatial correlation. Traditional theoretical model shape is fixed and chosen with subjectivity, which can't reflect the changing trend and multi-scale spatial characteristics. The choice of scale and the treatment of scale effects also need to be considered. To solve the problems above, we propose a method of kriging interpolation optimized by multi-scale least squares support vector machine (LS-SVM), which provides a new idea for fitting experimental variogram. Starting from the changing trend of the actual sample data, least squares support vector machine fits experimental variogram and the results conform to the spatial changing trend of data itself. Secondly, the wavelet kernel as the LS-SVM kernel function, parameters can be adjusted according to different parts of the nuclear, which is flexible and variable. Finally, the multi-scale wavelet kernel using wavelet multi-resolution characteristics, can reflect the different details of spatial changes, to avoid the single scale LS-SVM ignoring the spatial details of the problem. Followed that, the experiment includes simulation and application. Experimental simulation mainly verifies scientific validity and accuracy by the optimized interpolation of multi-scale least squares support vector machine. Meanwhile, experimental application research of PM2.5 concentrations of temporal and spatial distribution provides the theoretical basis for city ecological protection and controlling. Final results show that kriging interpolation algorithm optimized by multi-scale least squares support vector machine is superior to the traditional method and single scale optimized kriging interpolation algorithm. It would be better to depict the variation function and reflect the different scales of spatial changes in details to further improve the accuracy of the interpolation to some extent, which is an optional kriging interpolation method.

Key words： kriging interpolation; least squares support vector machine; variogram; multi-scale; wavelet kernel function

1 引言

空间插值（Spatial interpolation）是根据已知样本点的属性信息进行插值或推理来生成面状数据或估计待估点的数值[1-2],其本质就是通过建模来拟合生成尽可能逼近要素空间分布特征的函数关系[3]。空间插值方法种类广泛,常用的包括泰森多边形法、反距离加权法、移动拟合法、线性内插法、样条函数法、趋势分析法、克里金插值法等。其中,克里金插值法（Kriging interpolation）又称空间自协方差最佳插值法,这种方法认为空间属性的变化建立于空间位置关系和空间自相关关系的基础上,故定义具有这种变化的属性变量称为区域化变量（Regionalized variable）[4-5]。区域化变量独有的随机性和结构性使探索空间结构和空间变异规律变得有迹可循,变异函数（Variogram）应运而生,其可以对区域化变量的连续性、相关性、尺度性等进行空间描述。克里金插值方法的实质也是利用区域化变量的已知样本数据点变异函数结构特点,对待估位置点区域化变量的取值进行无偏、最优估计[6]。
克里金插值过程中,插值模型的精度取决于模型对空间变异性和空间相关性的反映程度[3]。插值模型需要用理论变异函数拟合实验变异函数,如何选择合适的理论变异模型关系到插值效果的优劣。传统克里金插值方法基于有限空间样本利用现有的理论变异函数模型进行拟合,对于变异函数的选择具有主观性和依赖性,且现有理论模型形状固定,无法反映实际样本的空间变化趋势。其次,空间变化趋势往往具有多尺度的特征,对于尺度的选择和尺度效应的处理也是需要考虑的[7],仅以传统方法拟合的实验变异函数会忽略空间变化的尺度效应。
对于克里金插值方法中对变异函数的优化,许多学者已进行了探究。传统数学统计方面,部分学者提出极大似然法、线性规划法等方法[8-9],此类方法往往会带来这样的问题,即计算量将成倍增加,但优化效果并不明显;人工智能方面,一些专家提出利用粒子群算法、遗传算法等方法[10-12],这类方法在一定程度上效果良好,但可能陷入局部最优困境,且结果优劣很大程度上取决于优化程度的好坏,受优化算法影响较大。同时,上述2类优化方法均未更深层次地考虑空间多尺度特性,从而忽略了空间尺度对空间细节的反映。
为解决上述问题,本文提出一种基于多尺度最小二乘支持向量机（以下简称多尺度LS-SVM）优化的克里金插值方法,此方法为拟合实验变异函数提供了一种新的思路。从已知样本数据点的空间变化趋势出发,通过多尺度LS-SVM拟合实验变异函数。经实验验证,此方法能更好地刻画变异函数,反映不同尺度下的空间变化细节,从而在一定程度上提高插值的精度,是一种可选的克里金插值方法。此方法有以下特点：首先,通过LS-SVM拟合实验变异函数,拟合的结果符合数据本身的空间变化趋势;其次,小波核作为LS-SVM核函数,能根据不同部分调整核参数,灵活多样;最后,多尺度小波核利用小波多分辨率的特性,可以体现空间变化的不同细节,避免了单尺度LS-SVM忽略空间细节的问题。

2 相关理论原理

2.1 普通克里金插值方法

变异函数,又称为变程方差函数,是地质统计学分析所特有的基本工具。它既能描述区域化变量的空间结构性变化,又能描述其随机性变化。
区域化变量z(x)在位置点x处和位置点x+h处值之差的方差一半定义为z(x)在x轴方向上的变 异函数。离散样本的实验变异函数可以通过式（1）计算。
γ h = 1 2 N h ∑ i = 1 N h z x i - z x i + h 2 i = 1 , … , N h （1）
式中：i+1,…, N(h);h为样本点对的距离;N(h)代表样本点对距离为h所有样本点对的个数;z(xi)和z（xi+h)分别是区域化变量z（x）在空间位置xi和xi+h处的真实值。当γ(h)不依赖于它的方向变化,只依赖于距离h时,则称γ(h)各向同性。
为得到任意距离下的变异函数以进行后续的插值,需要拟合实验变异函数,数据量较大时一般需在拟合之前进行分组操作。传统理论拟合模型包括球状模型、指数模型、高斯模型、线性模型等,以下将对普通克里金插值方法的原理进行说明。
克里金插值方法估计需要同时满足无偏性和估计方差最小：
E z x 0 - z * x 0 = 0 （2）
Var z x 0 - z * x 0 = min （3）
求解上述方程,经推导,用变异函数来表示如下,即：
∑ i = 1 k λ i γ x i , x j + μ = γ x 0 , x j j = 1 , … , k ∑ i = 1 k λ i = 1 （4）
式中：k为样本观测站点的总个数;λi为克里金权重系数,表示各空间样本点xi处的区域化变量z(xi)对待估位置点x0的贡献程度;μ为克里金拉格朗日乘子;γ(xi,xj)为空间样本点xi与空间样本点xj距离下的实验变异函数值;γ(x0,xj)为待估位置点x0与空间样本点xj距离下的实验变异函数值。
待估位置点x0处的区域化变量估计值z*（x0）可用式（5）进行估计：
z * x 0 = ∑ i = 1 k λ i z x i （5）
式中：x1,…,xk为已知样本的观测站点;z(x1),…,z（xk）为对应的观测值。
普通克里金插值方法在满足无偏性和估计方差最小这2个原则的条件下,建立了含有约束条件的拉格朗日函数。其中,无偏性体现在约束条件上,而方差最小体现在求极值问题上,拟合后的变异函数计算值直接作用在λi,导致λi的权重发生变化,进而对待估位置点的估值产生影响。

2.2 最小二乘支持向量机

最小二乘支持向量机（Least Squares Support Vector Machine,LS-SVM）用于函数回归时用平方和误差损失函数替代Vapnik的ε不敏感损失函数[13],同时用等式约束替代标准SVM中的不等式约束[14]。
给定N个样本的数据集 x m , y m m = 1 N , x m ∈ R d ,其中第m个输入xm对应的输出为ym,回归函数f(x)的基本形式如式（6）所示。
f x = ω T φ x + b （6）
式中：ω为权系数向量（列向量）;ωT代表其转置向量;φ（x）为输入空间到特征空间的映射函数;b为常数项。
LS-SVM对应的优化问题为：
min ω , b , e J ω , b , e = 1 2 ω T ω + 1 2 γ ∑ m = 1 N e m 2 = 1 2 ω 2 + 1 2 γ ∑ m = 1 N e m 2 （7）
式中：γ称为正则化参数或惩罚参数;em∈R为误差项。
满足等式约束条件：
y m = ω T φ x m + b + e m m = 1 , 2 , ⋯ , N （8）
上述优化问题可由拉格朗日方程求解,并基 于KKT（Karush-Kuhn-Tucker）条件得到最终模型（式（9））。
f x = ω T φ x + b = ∑ m = 1 N α m φ x m , φ x + b = ∑ m = 1 N α m K x m , x + b （9）
式中：αm为拉格朗日乘子;K(xm,x)为核函数。

2.3 小波函数

小波函数是一种时频局部化的工具,是把数据、函数或算子分割成不同频率的成分[15],然后再用分解的方法研究不同尺度下的成分,故非常适合做空间多尺度分析。其中,低频信息作为空间变化的大尺度框架,而高频信息反映空间变化的细节部分。目前应用较多的实数和复数小波函数中,以Morlet小波、Gaussian小波、Mexican hat小波等为主要应用[16],具体形式分别如下：
Morlet小波函数：
φ x = e - x 2 2 cos ω 0 x （10）
Gaussian小波函数：
φ x = x e - x 2 2 （11）
Mexican hat小波函数又称为墨西哥帽函数或者Marr小波：
φ x = 1 - x 2 e - x 2 2 （12）
小波函数作为SVM核函数,需满足2个条件：Mercer和平移不变定理[17]。以Morlet小波函数为例,经推导验证,同时满足这二者,最终得到Morlet小波核函数形式,如式（13）所示。
K x m , x = ∏ m = 1 d cos ω 0 × x m - x a m exp - x m - x 2 2 a m 2 （13）

3 基于多尺度LS-SVM优化的克里金插值方法

3.1 多尺度LS-SVM优化的克里金插值方法原理

多尺度LS-SVM以二尺度LS-SVM为例[18-20],具体实现如下：
二尺度LS-SVM模型在回归问题上,输入变量和输出变量分别为距离和实验变异函数值,故假定拟合的数据集为 { h i , γ h i } i = 1 n ,其中hi∈Rd,在这里d=1,表示第i个样本点对的距离,作为自变量,γ(hi)∈R,表示在距离hi下的实验变异函数值,作为因变量,n代表分组后需要拟合的实验变异函数值总个数。
在二尺度问题中,假定尺度1上小波核函数尺度因子偏大,尺度2上小波核函数尺度因子偏小。在较大尺度因子上,f1(x)用集合 { h i , γ h i } i = 1 n 来回归。回归模型可以用式（14）表示：
f 1 h = ω 1 T φ 1 h + b 1 （14）
在较小尺度因子上,f2(x)用集合 { h i , γ h i - f 1 h } i = 1 n 来回归,回归模型可用式（15）表示。
f 2 h = ω 2 T φ 2 h + b 2 （15）
将二尺度的模型融合为最终的回归模型,可以表示为：
f h = ω 1 T φ 1 h + ω 2 T φ 2 h + b （16）
其中,b=b1+b2
需要说明的是,f1(x)是对原始样本数据回归,f2(x)则是在f1(x)的基础上对残差进行回归,故又可以表示为：
f 1 h = ω 1 T φ 1 h + b （17）
f 2 h = ω 2 T φ 2 h （18）
对应的优化问题如下：
（19）
其相应的约束条件为：
y i = ω 1 T φ 1 h i + b + e 1 i i = 1 , 2 , ⋯ , n （20）
y i - f 1 h i = ω 2 T φ 2 h i + e 2 i i = 1 , 2 , ⋯ , n （21）
对偶问题的拉格朗日函数为：
L ω 1 , ω 2 , b , e 1 , e 2 , α 1 , α 2 = J ω 1 , ω 2 , b , e 1 , e 2 - ∑ i = 1 N α 1 i ω 1 T φ 1 h i + b + e 1 i - y i - ∑ i = 1 N α 2 i ω 2 T φ 2 h i + e 2 i + f 1 h i - y i （22）
根据KKT条件,该拉格朗日函数的最优解条 件为：
∂ L ∂ ω 1 ＝ 0 → ω 1 = ∑ i = 1 N α 1 i φ 1 h i + ∑ i = 1 N α 2 i φ 1 h i ∂ L ∂ ω 2 ＝ 0 → ω 2 = ∑ i = 1 N α 2 i φ 2 h i ∂ L ∂ e 1 i ＝ 0 → α 1 i = γ 1 e 1 i i = 1,2 , ⋯ , n ∂ L ∂ e 2 i ＝ 0 → α 2 i = γ 2 e 2 i i = 1,2 , ⋯ , n ∂ L ∂ b ＝ 0 → ∑ i = 1 N α 1 i + ∑ i = 1 N α 2 i = 0 ∂ L ∂ α 1 i ＝ 0 → ω 1 T φ 1 h j + b + e 1 j = y j j = 1,2 , ⋯ , n ∂ L ∂ α 2 i ＝ 0 → ω 2 T φ 2 h j + e 2 j + ω 1 T φ 1 h j + b = y j j = 1,2 , ⋯ , n （23）
经拉格朗日方程求解,矩阵形式为：
0 1 ̅ T 1 ̅ T 1 ̅ Ω 1 ＋ I γ 1 Ω 1 1 ̅ Ω 1 Ω 1 + Ω 2 ＋ I γ 2 × b α β = 0 y y （24）
式中： α ＝ α 11 , α 12 , ⋯ , α 1 n T , β ＝ α 21 , α 22 , ⋯ , α 2 n T , 1 ̅ ＝ 1,1 , ⋯ , 1 T , y = y 1 , y 2 , ⋯ , y n T ,I为单位矩阵, Ω 1 和 Ω 2 为小波核函数, Ω 1 = K 1 h p , h q = φ 1 h p T φ 1 h q p , q = 1 , ⋯ , n , Ω 2 = K 2 h p , h q = φ 2 h p T φ 2 h q p , q = 1 , ⋯ , n 。
最终,二尺度LS-SVM模型表示为式（25）。
f h = ω 1 T φ 1 h + ω 2 T φ 2 h + b = ∑ i = 1 n α 1 i + α 2 i K 1 h i , h + b 1 + ∑ j = 1 n α 2 j K 2 h j , h + b 2 （25）
式中：K1和K2为不同尺度下的小波核函数。

3.2 基于多尺度LS-SVM优化的克里金插值方法 流程

基于多尺度小波LS-SVM优化的克里金空间插值方法实现流程具体可分为如下4步：
（1）通过离散变异函数公式计算所有样本点对的实验变异函数值;
（2）采用多尺度小波LS-SVM拟合实验变异函数值,得到理论变异函数模型;
（3）建立克里金空间插值方程组,根据理论变异函数模型求解克里金权重系数λi;
（4）根据克里金权重系数λi计算待估位置点的区域化变量估计值;
其中,步骤（2）中,多尺度小波LS-SVM拟合实现变异函数曲线过程,又可细分为以下4步：
（1）计算所有点对的实验变异函数值;
（2）实验变异函数值进行分组操作;
（3）通过多尺度小波LS-SVM拟合分组后的实验变异函数;
（4）得到最终理论变异函数。

4 实验与分析

为检验插值效果,采用交叉验证的方法[21],即假定每个观测站点要素值是未知的,通过其余插值观测点预测来估测,最后计算所有已知位置点的实际观测值和估测值误差,以此评判插值结果,评价指标采用平均绝对误差（MAE）和均方根误差（RMSE）2个评价指标[22]。MAE反映估测的总体误差,RMSE反映样本数据的估测灵敏度和极值效应,二者值越小,表明效果越好。其定义分别如下：
MAE = 1 S ∑ i = 1 S y i - y ^ i （26）
RMSE = 1 S ∑ i = 1 S y i - y ^ i 2 （27）
式中：S代表待估位置点的总个数;yi为待估位置点的区域化变量真实值,i=1,…, S; y ^ i 为待估位置点的区域化变量估计值,i=1,…, S。

4.1 实验模拟

模拟采用Spatial Interpolation Comparison 97数据（简称SIC97数据）,该数据为1986年5月8日瑞士467个降水观测站的降水数据,其中样本观测站点总个数为100,即k=100,待估位置点的总个数为367,即S=367。数据来源：https://www.researchgate.net/profile/Gregoire_Dubois/publication/281292 076_Spatial_Interpolation_Comparison_97_（SIC97）_dataset,通过已知观测点来插值待估点,观测点与待估点分布如图1所示。
Fig. 1 The distribution diagram of SIC97 data

图1 SIC97数据分布示意图

由于样本点数目较多,拟合实验变异函数之前,进行分组操作,分组后需要拟合的实验变异函数值总个数为20,即 n = 20 。图2和图3分别表示分组前的变异函数云图和分组后的变异函数云图,图4表示需要拟合的实验变异函数值。
Fig. 2 Variogram cloud before grouping

图2 分组前的变异函数云图

Fig. 3 Variogram cloud after grouping

图3 分组前的变异函数云图

Tab. 1 Comparisons of Kriging interpolation results based on the theoretical variogram models

表1 各理论变异模型拟合后的克里金插值结果对比

类别 球状模型 指数模型 高斯模型 单尺度
LS-SVM（RBF）		 单尺度
LS-SVM
(wave)		 二尺度
LS-SVM
(wave)
MAE 39.43 40.86 110.46 39.81 39.09 38.54
RMSE 56.21 57.33 154.82 56.78 55.97 55.70
Fig. 4 The fitting experimental variogram values

图4 需要拟合的实验变异函数值

为对比多尺度LS-SVM优化的克里金插值方法的效果,采用传统理论变异函数模型：球状模型、指数模型、高斯模型。此外,也对比了单尺度RBF核函数和Morlet小波核函数LS-SVM理论模型。其中,单尺度和多尺度LS-SVM的惩罚参数和核参数寻优都是通过交叉验证的方法,单尺度RBF核函数LS-SVM的惩罚参数γ=168,核参数σ2=8×103;单尺度Morlet小波核LS-SVM的惩罚参数γ=681,核参数α=328;二尺度Morlet小波核LS-SVM的惩罚参数γ1=3.7×103,γ2=378,核参数γ1=0.03,a2=14。最终拟合的实验变异函数如图5所示。
Fig.5 Variogram fitting curve

图5 变异函数拟合曲线

表1直观地显示出各理论变异模型拟合后的克里金插值结果,从表中对比可以发现,在此实验数据中,传统克里金插值拟合模型中,球状模型最优;单尺度RBF核函数LS-SVM模型优于高斯模型和指数模型,但略低于球状模型;单尺度Morlet核函数LS-SVM模型略优于传统克里金插值拟合最好的球状模型;多尺度Morlet核函数LS-SVM模型优于单尺度LS-SVM模型,在所有拟合模型中插值效果最优。结合图5可更加直观地体现多尺度LS-SVM模型在拟合实验变异函数中的变化,其中,大尺度保证了整体的空间变化趋势,小尺度体现了细节变化,这样最终使空间多尺度特征得以反映,插值精度较传统方法具有一定程度上的提高,尽管提高并不明显,但结合本方法对变异函数拟合的准确性来看,本文为拟合变异函数提供了一种思路,该方法是一种可选的克里金插值方法。

4.2 实验应用

将该插值方法应用于青岛市PM2.5浓度数据的空间插值分析。PM2.5是指环境空气中动力学当量直径小于等于2.5 μm的颗粒物,富含大量有害物质,研究PM2.5数据的时空分布特征对城市环境防护及控制具有实践意义。
4.2.1 数据源
青岛市19个国控空气质量监测站点,包括青岛市区9个站点、黄岛区（包括原胶南市）2个、即墨市2个、平度市2个、胶州市2个、莱西市2个,其分布如图6所示。PM2.5数据选取日均值,时间跨度2015年10月至2016年9月,共计12个月。
Fig. 6 The map of air quality monitoring sites in Qingdao city

图6 青岛市空气质量监测站点分布图

4.2.2 实验分析及对比
通过对青岛市19个国控站点PM2.5数据的统计分析可以发现,无论月均值还是最值,PM2.5浓度12月最高,而8月浓度最低（图7）,统计结果与此文所述也基本相符[23]。故选取特征最明显的12月和8月作为样本数据,以此估测青岛市PM2.5分布特征。
Fig. 7 Monthly PM2.5 concentration values of Qingdao city

图7 青岛市各个月份PM2.5浓度值

实验流程与SIC97数据的流程类似,故不再详细介绍。本实验为验证本文方法得到的插值结果是否符合样本实际的大致趋势,采用普通克里金插值方法（指数模型）得到图8（a）和图8（c）所示的8月和12月青岛市PM2.5分布图,而图8（b）和图8（d）则是通过多尺度LS-SVM优化的克里金插值方法得到的8月和12月青岛市PM2.5分布图。
Fig. 8 Results of PM2.5 interpolation methods in Qingdao

图8 青岛市PM2.5各插值方法结果分布图

对比2种方法得到插值结果整体趋势,多尺度LS-SVM优化的克里金插值方法结果大概符合样本通过普通克里金方法（指数模型）插值出来的结果,空间变化趋势相近,因此认为此方法具有一定可靠性。一个更为保守的说法是,采用多尺度LS-SVM优化的克里金插值方法得到的结果与实验样本经克里金插值方法（指数模型）得到的结果趋势一致,此思路具有一定的借鉴性。而从2种方法得到的插值结果差异上看,本文方法插值结果的空间细节变化特征更明显,且存在特征变化急剧的位置,体现了多尺度LS-SVM模拟拟合变异函数中的多尺度效应以及细节变化效应,而普通克里金方法插值结果空间变化较为平缓,拟合实现变异函数时采用的是传统理论变异函数模型,只需求解固定参数,故一定程度上会忽视空间细节的变化。
Tab.2 The comparison of the interpolation results of each model in December

表2 12月各模型插值结果对比

指标类别 高斯模型 指数模型 球状模型 单尺度
LS-SVM（RBF）		 单尺度
LS-SVM
(wave)		 二尺度
LS-SVM
(wave)
MAE 14.41 12.52 12.54 12.47 12.46 12.42
RMSE 17.66 15.73 15.77 15.57 15.57 15.31
整体的时间变化趋势,PM2.5浓度12月远高于8月,8月PM2.5浓度在13~30 μg/m3的区间范围,而12月PM2.5浓度在70~130 μg/m3的区间范围,二者相差较大。结合图7可以得到结论,青岛市PM2.5季节变化明显,表现为冬季高、夏季低的时间特征。整体的空间变化趋势,从东向西PM2.5浓度由低到高,表现为逐渐上升的趋势。8月青岛市的东南、东北方向PM2.5浓度较低,并呈现逐渐向西加重的变化态势,层次感明显;而12月,呈现2个集中分布的较高浓度地带,同时也表现为东南、东北方向PM2.5浓度偏低,西北方向和中部偏高的态势。因此根据其时空分布特征来重点防护与治理具有重要的研究意义。
最后,为对比多尺度小波LS-SVM优化的克里金插值方法,进行各理论变异模型拟合的结果对比。表2和表3分别代表交叉验证检验得到的12月和8月各理论变异模型拟合后的克里金插值结果误差对比。总结来看,多尺度LS-SVM优化的克里金插值方法产生的平均绝对误差和均方根误差相比传统的克里金方法（高斯模型、指数模型、球状模型）、单尺度LS-SVM优化的克里金插值方法均较小,在一定程度上提高了插值精度,可以作为一种可选的克里金插值方法。
Tab. 3 The comparison of the interpolation results of each model in August

表3 8月各模型插值结果对比

指标类别 高斯模型 指数模型 球状模型 单尺度
LS-SVM（RBF）		 单尺度
LS-SVM
(wave)		 二尺度
LS-SVM
(wave)
MAE 3.62 3.47 3.49 3.56 3.44 3.44
RMSE 4.54 4.20 4.27 4.58 4.23 4.19

5 结语

传统克里金插值方法基于变异函数,根据待估位置点、已知样本数据点的位置关系和区域化变量的空间相关性,实现空间加权估计,但存在以下问题：传统方法理论模型形状固定且选择具有人为主观性;空间数据的变化趋势无法反映;空间数据多尺度特征被忽略。为解决上述问题,本文提出了一种基于多尺度最小二乘支持向量机优化的克里金插值方法,此方法为拟合实验变异函数提供了一种新的思路。从实际样本数据的变化趋势出发,采用最小二乘支持向量机拟合实验变异函数,并利用不同尺度小波核反映不同尺度下的空间变化。最后,实验环节包括模拟和应用,模拟主要验证经多尺度最小二乘支持向量机优化的插值方法的科学有效性以及准确性,应用主要研究青岛市PM2.5时空分布特征,为城市生态科学防护及控制提供理论依据。最终表明,基于多尺度最小二乘支持向量机优化的克里金插值方法能更好地刻画变异函数,反映不同尺度下的空间变化细节,从而在一定程度上提高插值精度,是一种可选的克里金插值方法。本文提出这种方法,一方面是基于对空间数据空间特征（如位置、分布特征等）和属性特征（如空间相关性、空间异质性）等性质的理解;另一方面基于对变异函数理论拟合模型（球状、指数、高斯等）的分析,辅以真实的数据进行实验测试,结果真实地表明该方法能在一定程度上提高插值方法的精度,尽管提高并不明显,但结合本方法对变异函数拟合的准确性来看,综合理论分析与实践,提出针对克里金插值问题的一种创新研究思路。通过实验模拟（经典SIC97插值数据）和应用（青岛市PM2.5数据）,认为具有一定的借鉴性。值得注意的是,利用样本点对拟合实验变异函数时,样本是影响插值准确性的重要因素。当样本的数量不大时,其空间特征可能并不是空间对象性质的真实反映;而样本数量较大时,其拟合程度的好坏会直接作用待估位置点的估计值。因此,不能完全依靠自动化拟合实验变异函数,需要结合专家经验知识来进一步判断,同时如何准确合理地实现分组操作也值得研究。此外,今后还可更细致地对空间进行尺度划分进行研究,从宏观与微观进行适当选取。
致谢：感谢青悦开放环境数据中心（https://data.epmap.org）提供的环境数据处理支持。

The authors have declared that no competing interests exist.

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
Lam N S N. Spatial interpolation methods: A review[J]. The American Cartographer, 1983,10(2):129-150.Two forms of spatial interpolation, the interpolation of point and areal data, are distinguished. Traditionally, point interpolation is applied to isarithmic, that is, contour mapping and areal interpolation to isopleth mapping. Recently, areal interpolation techniques have been used to obtain data for a set of administrative or political districts from another set of districts whose boundaries do not coincide. For point interpolation, the numerous methods may further be classified into exact and approximate. Exact methods include most distance-weighting methods, Kriging, spline interpolation, interpolating polynomials, and finite-difference methods. Approximate methods include power-series trend models, Fourier models, distance-weighted least-squares, and least-squares fitting with splines. Areal interpolation methods, on the other hand, are classified according to whether they preserve volume. Traditional areal interpolation methods which utilize point interpolation procedures are not volume-preserving, whereas the map overlay and pycnophylactic methods are. It is shown that methods possessing the volume-preserving property generally outperform those that do not.

DOI

[2]
王劲峰,武继磊,孙英君,等.空间信息分析技术[J].地理研究,2005,24(3):464-472.在GIS技术日趋成熟和空间数据极大丰富的今天,通过分析空间数据探索空间过程机理正变得日益迫切。空间信息分析技术至少包括以下六个主要方面:(1)空间数据获取和预处理;(2)属性数据空间化和空间尺度转换;(3)空间信息探索分析;(4)地统计;(5)格数据分析;(6)复杂信息反演和预报。本文提出了解决具体应用问题一般的空间数据分析计算、结果解释和反馈程序。认为空间过程的一般共性和作为共同的研究对象,各种不同的方法技术最终可能导致空间数学(spatialmathematics)的产生,同时发展鲁棒的空间分析软件包对于普及空间数学是必要的。

[ Wang J F, Wu J L, Sun Y J, et al.Techniques of spatial data analysis[J]. Geographical Research, 2005,24(3):464-472. ]

[3]
朱会义,刘述林,贾绍凤.自然地理要素空间插值的几个问题[J].地理研究,2004,23(4):425-432.资源管理、灾害管理、生态环境治理以及全球变化研究的需要强化了部分自然地理要素空间插值研究的重要性。这些要素空间插值的核心是建立充分逼近要素空间分布特征的函数方程。对于给定的区域与要素样本值 ,插值函数可以有多种模型形式。各类模型的精度受其理论基础、模型算法、时空尺度效应、样本数据属性等因素的综合影响。通过对国际主要插值研究成果进行分析 ,文章认为各类模型插值精度的差异缘于模型对插值要素空间变异性与空间相关性的反映 ,具体应用中 ,只有对已知样本数据进行变异性与相关性分析才能选出适当的插值方法。

DOI

[ Zhu H Y, Liu S L, Jia S F.Problems of the spatial interpolation of physical geographical elements[J]. Geographical Research, 2004,23(4):425-432. ]

[4]
Matheron G.Principles of geostatistics[J]. Economic Geology, 1963,58(8):1246-1266.

DOI

[5]
孙洪泉. 地质统计学及其应用[M].北京:中国矿业大学出版,1990.

[ Sun H Q.Geostatistics and application: Quantitative geography[M]. Beijing: China University of Mining and Technology Press, 1990. ]

[6]
徐建华. 计量地理学[M].北京:高等教育出版社,2006.

[ Xu J H. Quantitative geography[M]. Beijing: Higher Education Press, 2006. ]

[7]
孟斌,王劲峰.地理数据尺度转换方法研究进展[J].地理学报,2005,60(2):277-288.尺度问题在地理学、生态学和水文科学等众多领域都具有非常重要的地位。近年来,随着对地观测技术和地理信息科学的飞速发展,解决地理数据的尺度转换问题成为目前地理信息科学及相关研究中的热点和难点问题之一。在地理信息科学相关领域中,地图学和遥感科学研究人员从尺度概念的理解到尺度转换的理论和方法都做了大量的研究,对解决地理数据空间特征的尺度转换做出了重要的贡献。在地理数据属性特征的尺度转换研究领域,地理信息科学研究者提出的面域插值方法是解决此问题的主要方法之一。同时,在社会经济科学领域,“小区域统计学”也发展了一套相关的理论和方法,试图解决统计单元的变更问题。文章在全面回顾和比较不同研究领域解决“尺度转换

DOI

[ Meng B, Wang J F.A review on the methodolody of scaling whith geo-data[J]. Geographica Sinica, 2005,60(2):277-288. ]

[8]
王仁铎,胡光道.线性地质统计学[M].北京:地质出版社,1988.

[ Wang R D, Hu G D.Linear geological statistics[M]. Beijing: Geological Publishing House, 1988. ]

[9]
矫希国,刘超.变差函数的参数模拟[J].物探化探计算技术,1996,18(2):157-161.

[ Jiao X G, Liu C.Estimation of variation parameter[J]. Computing Techniques for Geophysical and Geochemical Exploration, 1996,18(2):157-161. ]

[10]
Li M, Li G, Azarm S.A kriging meta-model assisted multi-objective genetic algorithm for design optimization[J]. Journal of Mechanical Design, 2008,130(3): 031401.ABSTRACT The high computational cost of population based optimization methods, such as multi-objective genetic algorithms (MOGAs), has been preventing applications of these meth-ods to realistic engineering design problems. The main challenge is to devise methods that can significantly reduce the number of simulation (objective/constraint functions) calls. We present a new multi-objective design optimization approach in which the Kriging-based metamodeling is embedded within a MOGA. The proposed approach is called Kriging assisted MOGA, or K-MOGA. The key difference between K-MOGA and a conventional MOGA is that in K-MOGA some of the design points are evaluated on-line using Kriging metamodeling instead of the actual simulation model. The decision as to whether the simulation or its Kriging metamodel should be used for evaluating a design point is based on a simple and objective criterion. It is determined whether by using the objective/constraint functions' Kriging metamodels for a design point, its "domination status" in the current generation can be changed. Seven numerical and engineering examples with different degrees of difficulty are used to illustrate applicability of the proposed K-MOGA. The results show that on the average K-MOGA converges to the Pareto frontier with an approximately 50% fewer number of simulation calls compared to a conventional MOGA.

DOI

[11]
张强,许少华,于文涛,等.粒子群算法在克里金三维地质建模中的应用[J].大庆石油学院学报,2011,35(1):85-89.

[ Zhang Q, Xu S H, Yu W T, et al.Application of particle swarm optimization to Kriging three-dimensional geological modeling[J]. Journal of Daqing Petroleum Institute, 2011,35(1):85-89. ]

[12]
严华雯,吴健平.加权最小二乘法改进遗传克里金插值方法研究[J].计算机技术与发展,2012,22(3):92-95.数据内插被广泛应用于地统计分析领域,克里金插值作为其中最为有效的方法之一,其原理是通过建立变异函数理论模型,得到可靠的权重值和拉格朗日系数,构成求解待测点的线性组合.为了有效地提高插值精度,文中利用加权最小二乘法优化遗传算法中的适应度函数,进而改进普通基于遗传算法优化的克里金插值方法.并且在MATLAB中利用外部工具箱确定模型参数,最后通过实例验证,将该方法与普通克里金插值以及遗传克里金插值结果进行对比,发现采用该方法,插值效果较好且误差也较小,证明了通过加权最小二乘法可以有效改进普通遗传克里金插值方法.

DOI

[ Yan H W, Wu J P.Research on genetic algorithm kriging optimized by weight least square[J]. Computer Technology and Development, 2012,22(3):92-95. ]

[13]
Vapnik V.The nature of statistical learning theory[M]. New York: Springer, 1995.

[14]
Suykens J A K, Vandewalle J. Least squares support vector machine classifiers[J]. Neural Processing Letters, 1999,9(3):293-300.In this letter we discuss a least squares version for support vector machine (SVM) classifiers. Due to equality type constraints in the formulation, the solution follows from solving a set of linear equations, instead of quadratic programming for classical SVM's. The approach is illustrated on a two-spiral benchmark classification problem.

DOI

[15]
Daubechies I.Ten lectures on wavelets[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1992.

[16]
张建锋,崔树军,李国敏.常用小波及其时-频特性[J].地学前缘,2012,19(6):248-253.

[ Zhang J F, Cui S J, Li G M.Frequently-used wavelets and their time-frequency properties[J]. Earth Science Frontiers, 2012,19(6):248-253. ]

[17]
周建萍,郑应平,王志萍.基于Morlet小波核多类支持向量机的故障诊断[J].华东电力,2008,36(8):76-80.故障诊断问题实质上是一个模式识别问题,即多分类问题。采用Morlet小波来构造支持向量机（Sup-port Vector Machine,SVM）的核函数,Morlet小波核SVM比普通SVM具有更好的鲁棒性和更强的泛化能力。在一对一算法的基础上实现Morlet小波核多类支持向量机的故障诊断,并将此方法成功应用于电厂汽轮发电机组的故障诊断。实验仿真结果表明Morlet小波核多类SVM故障分类器比BP神经网络训练和测试速度快,且其分类精度在高斯噪声干扰下还保持100%,比BP神经网络高出11.8%。因此该方法能够快速而准确地对电厂汽轮发电机组的故障进行诊断,满足电力系统实时操作的要求。

DOI

[ Zhou J P, Zheng Y P, Wang Z P.Fault diagnosis based on Morlet wavelet kernel multi-class support vector machine[J]. East China Electric Power, 2008,36(8):76-80. ]

[18]
王琴,沈远彤.二尺度最小二乘小波支持向量回归[J].工程地球物理学报,2009,6(4):516-520.

[ Wang Q, Shen Y D.Multi-scale least squares support vector regression[J]. Chinese Journal of Engineering Geophysics, 2009,6(4): 516-520. ]

[19]
张相胜,王蕾,潘丰.多尺度最小二乘小波支持向量机的回归建模[J].计算机工程,2012,38(10):175-177.普通最小二乘支持向量机算法用于多尺度回归建模时精度较低。针对该问题，选取墨西哥草帽小波函数作为最小二乘支持向量机的核函数，设计一种基于小波核的多尺度最小二乘小波支持向量机。在此基础上，通过解二次优化问题求出多尺度回归建模问题的全局最优解，最终得出的多尺度回归模型能够有效地逼近多尺度信号。仿真结果表明，该算法具有较高的精度。

DOI

[ Zhang X S, Wang L, Pan F.Regression modeling of multi-scale least square wavelet support vector machine[J]. Computer Engineering, 2012,38(10):175-177. ]

[20]
王琴,沈远彤.基于压缩感知的多尺度最小二乘支持向量机[J].自动化学报,2016,42(4):631-640.提出一种基于压缩感知(Compressive sensing, CS)和多分辨分析(Multi-resolution analysis, MRA)的多尺度最小二乘支持向量机(Least squares support vector machine, LS-SVM). 首先将多尺度小波函数作为支持向量核, 推导出多尺度最小二乘支持向量机模型, 然后基于压缩感知理论, 利用最小二乘匹配追踪(Least squares orthogonal matching pursuit, LS-OMP)算法对多尺度最小二乘支持向量机的支持向量进行稀疏化, 最后用稀疏的支持向量实现函数回归. 实验结果表明, 本文方法利用不同尺度小波核逼近信号的不同细节, 而且以比较少的支持向量能达到很好的泛化性能, 大大降低了运算成本, 相比普通最小二乘支持向量机, 具有更优越的表现力.

DOI

[ Wang Q, Shen Y D.Multi-scale least squares support vector machine using compressive sensing[J]. Acta Automatica Sinica, 2016,42(4):631-640. ]

[21]
Johnston K, Ver Hoef J M, Krivoruchko K, et al. Using ArcGIS geostatistical analyst[M]. Redlands: Esri, 2001.

[22]
Han S, Schneider S M, Evans R G.Evaluating co-kriging for improving soil nutrient sampling efficiency[J]. Transactions of the ASAE, 2003,46(3):845.ABSTRACT The spatial variability of soil texture and soil nitrate-N, P, and K was studied in two center-pivot irrigated fields (89 ha total). Two soil texture components (clay and silt) were found to be correlated with soil nitrate-N, P, and K, and were used as auxiliary variables in the cokriging procedure to estimate soil nitrate-N, P, and K at unsampled locations. With a sampling density of 2.7 sites/ha (61 61 m grid) as the baseline, removal of 50% of the sampling sites resulted in a normalized mean absolute error (NMAE) of 34.5%, 22.9%, and 15.3% for soil nitrate-N, P, and K, respectively. These numbers reduced to 32.8%, 20.7%, and 12.0% when only 25% of the sampling sites were removed. The study showed that the cokriging technique provided slightly better estimates than the ordinary kriging method for soil P and K at a higher sampling density (>2.1 sites/ha). However, when a variable has a large random variation, such as the soil nitrate-N, cokriging did not provide better estimates than ordinary kriging. The results of this study provide guidelines on the selection of kriging or cokriging in improving the soil nutrient sampling efficiency. nderstanding the spatial distributions of soil properties, particularly soil nitrate-N, P, and K, within a field is important for the development of site-specific fertilizer management strategies. Because of limited resources, intensive soil sampling is impractical. Typically, soil samples are collected from only a small number of sites in the field. The data at the sampling locations are then used to predict values at unsampled locations and thus to generate the spatial distributions. Spatial interpolation techniques are needed to study the spatial distributions of soil properties. If a soil property is autocorrelated in space, kriging methods, based on Mather- on's regionalized variable theory (Matheron, 1963), can be utilized to provide the best linear unbiased estimates (BLUE). Many researchers have applied the kriging tech- niques to evaluate variability of soil physical and chemical

DOI

[23]
杜艳伟,程建光,王静,等.2013年青岛市PM2.5时空分布及来源分析[J].城市环境与城市生态, 2015,28(2):17-19.

[ Du Y W, Cheng J G, Wang J, et al. Temporal and spatial distribution of PM2.5 in Qingdao city and analysis of its sources in 2013[J]. Urban Environment & Urban Ecology, 2015,28(2):17-19. ]

Options
文章导航

/

〈 〉