基于Kriging的加密自动气象站要素场插值与改进

王兴　张琳焓　王丽娟　张潇潇

摘 要：为了解决加密自动气象站资料的时空分辨率高，对资料中各个要素场插值计算时耗长，难以应用到对实时性要求高的业务系统中的问题，提出一种改进的Kriging算法，针对加密自动站的空间分布和气象学特点，预先对自动站分区，并对要素场插值的算法流程进行优化，提高其计算速度，以满足相关业务对气象数据实时性的要求。实验结果证明，本方法插值结果更接近实际观测值，且计算速度提升了约20倍。

关键词：克里金；插值；加密自动气象站；模式初始场

中图分类号：TP391.4 文献标识码：A

1 引 言（Introduction）

我国是世界上受气象灾害影响较大的地区之一，随着社会、经济的发展，气象灾害的影响越来越广，造成的损失越来越大，所受到的关注也越来越多。据不完全统计，近年来我国每年因短时强降水、雷雨大风等极端天气所造成的经济损失高达数千亿元，受影响人口超过4亿人次[1]。因此，国家对气象服务，特别是对灾害性天气的预报、预警以及防灾减灾服务的要求越来越高。为了提高气象监测、预测的能力，国家投入了大量资金推动和促进大气监测自动化的进程[2]，自2000年开始，至今已投入业务化运行的加密自动气象站（以下简称自动站或站点）数量超过4万部，且仍在扩充建设。

自动站是一种适用于中尺度区域地面气象观测的专业监测设备，可将温度、湿度、气压、风向、风速等多种气象要素的实况监测数据实时发送到气象业务部门。在我国中东部及沿海地区，约每5—10公里就有一套自动站系统，监测数据发布的频次一般为5分钟，密集时可达到每分钟一次。除了提供高时空分辨率的气象实况观测资料外，自动站数据还可应用于中尺度数值天气预报模式和灾害性天气预报预警等业务系统中[3，4]，用以提供更加精细、准确、客观的模式初始场资料，进而提高数值模式预报的水平。

由于自动站在空间分布上是若干个离散的点，而模式预报等气象业务系统所需的初始场资料主要是等间距的格点场，因此，自动站观测数据需要经过插值处理。常用的插值算法有双线性插值[5]、反距离权重插值[6]、样条插值，以及Kriging（克里金）插值[7-10]等，其中反距离权重插值因其算法简单、易于实现，且考虑了距离对权重的影响，应用较为广泛[11-13]。然而，各个气象要素在空间上的分布并非均匀、平滑的，其受到地形、高程以及大气运动的影响非常大。Kriging算法也称为空间局部估计法，是在变异函数理论和结构分析的基础上，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法[14]。该算法的插值结果往往更接近于实际值，从而能够更加真实地反映各个气象要素的空间分布情况。然而，Kriging算法的运算非常复杂，计算耗时长，对于更新频次高、观测要素多的自动站数据，难以满足现代气象业务系统的实时性要求。

本文针对自动站数据在空间分布上的特点，设计了一种方法来降低Kriging插值的运算量和计算的时耗，以满足自动站数据插值相关应用实时业务化的要求。

2 算法流程设计（The process design of algorithms）

2.1 自动站的分块处理

在插值算法中，每一个待插值点都需要搜寻若干个与其位置相邻的实测点，通行的方法就是对所有实测点进行逐点遍历。以自动站数据为例，假设有40 000个有效自动站数据，目标格点数据的空间分辨率为10km×10km，则能够覆盖我国所有地区的矩形所含格点总数约为：

G1=经向跨度×纬向跨度×单位经向公里数×单位纬向公里数/目标格点数据经向分辨率/目标格点数据纬向分辨率≈|35-73|×|53-3|×110×110/10/10=375 100。

如果对这G1个格点逐一搜寻40 000个自动站中最邻近的n个用以插值计算，那么其运算规模达到1.5×1010。考虑到每次搜寻相邻的自动站时，绝大部分站点因相距待插值点甚远，没有遍历的必要，本文则使用如下方法改进。以1×1的经纬网格为单位，将遍布于全国各地的自动站分别映射到相应地理位置的经纬网格中。即创建一个二维数组G[X，Y]，对于任意一个待插值格点g[i，j]，g[i，j]∈G[x，y]，其中X为我国领土东西径向跨度的值取整加1，Y为南北纬向跨度的值取整加1，x∈[0，X]，y∈[0，Y]。每计算一个格点时，只查找以其为中心，上下左右最相邻的9个单位经纬网格中所有自动站的记录。如此将每个待插值格点遍历自动站的数量大幅降低。

2.2 变异函数的计算

变异函数是指区域化变量Z（x）和Z（x+h）的差值的平方的数学期望，一般标记为，其中，h为这两个变量之间的距离。为了便于算法的工程化，假定与自变量x无关，即变异函数的值不随Z（x）的位置而变化，仅与h有关，亦即函数满足二阶平稳假设[15]，则变异函数的计算可简化为如下式：

（1）

式（1）中，E表示数学期望，n为样本数据的对数，xi表示样本点的值，习惯上将称为变异函数，相应地称为半变异函数，它们反映的是Z（x）与Z（x+h）随距离变异的程度，本节以下内容即是对公式1计算过程的描述。

由于自动站在全国的分布并不规则，总体上呈东南部多西北部少，同时，各气象要素在空间分布上具有一定的气象学规律，如发生天气过程或冷暖锋过境时，其所影响的地区的温度、气压等要素在空间上可显现一定的梯度变化特征和变化趋势方向性，针对这些特征，本文在计算变异函数时，分为以下几步进行：

（1）按方位角将自动站分组。将每次参与插值计算的9个单位经纬网格中所有自动站，以当前待插值格点为圆心，正北方向为0，逐45划分为8个扇形区域，则每个自动站必在且仅在一个扇形区域中。

（2）按距离将自动站分组。将上一步按扇形区域分组后的自动站再按各自与圆心间的距离分组。单位距离按实际情况可调整。距离的计算方法为：

（2）

其中，g[i，j]表示待插值格点，Z[x，y]表示自动站，lat和lon分别表示格点或自动站的纬度和经度。如图1所示，每个圆点代表1个自动站，每个自动站都落在某一扇区中，每个小方格覆盖的地理范围为一个经纬度网格。其中，，d为单位距离的整数倍。

图1 自动站按角度和距离分组示意图

Fig.1 A schematic diagram of AWS grouped

by angle and distance

（3）计算的主要流程如下：

a.确定变异函数的方向数，即划分的扇形区域的数量有N个。

b.确定距离步长为L，步数为M。

c.假定参与一次插值的自动站数量不超过K个。

d.定义三维数组Value[M，N，K]，用于存放自动站记录，由于每条记录包含若干个气象要素信息，因此Value的每一个元素项又是一个一维数组。

e.定义数组Sum[M，N]，作为计数器，Value[M，N，K]每新增一条记录，相应的Sum[M，N]加1。

f.选定一个待插值格点g[i，j]，遍历所在地理位置相对应的9个单位经纬网格中的自动站，直到所有自动站信息都记录到Value[M，N，K]和Sum[M，N]中。

g.开始计算变异函数的值，计算方法为Value数组元素中各个气象要素各自的累加值除以相应Sum元素的值，即求取逐个扇区中，各个自动站要素的平均值。

2.3 变异函数模型和拟合

由2.2节可得到半变异函数的分布图，如实验部分图2所示。在此基础上，进行半变异函数的曲线拟合。该步骤是一个建模的过程，类似于回归分析，目的是将分布图拟合成连续的曲线。常用的变异函数模型有球面模型、指数模型、高斯模型以及线性模型等[15，16]，由于变异函数模型的选择很大程度上需要经过大量试验，反复选择调整，且只能保证某一特定数据条件下的相对最优解，本文球面模型为例进行描述。

球面模型的公式为：

（3）

其中，表示变异函数，c0为块金效应常数，c为偏基台值，c0+c为基台值，a为变程。因为当h=0和h>a时，为常数，因此只要分析的情况。

不妨记，b0=c0，b1=3c/2a，b2=-c/2a3，x1=h，x2=h3，，则可将：

变形为线性方程的形式：

（4）

设有一组数据h和，亦即一组x1，x2和y，将其代入式（4），得到线性方程组：

（5）

为方便求解，将该非齐次线性方程组转换为AX=b的矩阵形式，其中Xi1=hi，Xi2=hi3，，i=1，2，…，n。

由此，将式（5）的求解转化为线性规划的求解问题，有多种有效的解法[17]，此过程不在本文中赘述。

通过上述步骤的计算，即可得到变异函数的各项参数，参数的确立为下一步Kriging插值计算提供了关键性支撑。有关上述变异函数拟合效果的检验，在下文的实验部分进行描述。

2.4 Kriging插值的计算

Kriging插值的基本公式为：

（6）

从形式上看，式（6）与反距离权重插值非常相似，每个待插值点的结果都是若干个相邻实测点的加权平均值。但Kriging方法基于包含自相关的统计模型，其权重不仅取决于待插值点与实测点间的距离，还取决于基于实测点的整体空间排列，即实测点在空间分布上的影响。因此，Kriging方法的式（6）中，取决于实测点与待插值点的距离，以及待插值点周围的实测点之间空间关系的拟合模型。

Kriging插值计算的流程如下：

a.创建一个二维数组G[X，Y]，对于任意一个待插值格点g[i，j]，g[i，j]∈G[X，Y]。

b.依次选择一个g[i，j]，根据2.1节所述的分块搜索策略，确定参与本次插值计算的若干自动站，作为实测点。

c.根据2.2和2.3节所述变异函数模型及函数参数计算方法，求解方程组，得到一组最优的权重系数。

d.将得到的代入式（6），求得待插值点g[i，j]的值。需要说明的是，此处的g[i，j]为一个数值，而实际应用中，由于自动站监测数据是由多个气象要素组成，因此，严格意义上说，g[i，j]的求解是对多个气象要素分别进行。

e.重复步骤b—d，直到所有g[i，j]的值全部求出。

f.输出G[X，Y]即为预期的格点场形式的数据。此数据可用于天气数值预报模式等业务系统中作为初始场资料。

3 实验结果与分析（The experimental results and

analysis）

实验数据选用2015年9月12日05时整点加密自动气象站数据，经过初步质量检验，剔除缺测和异常数值后，该数据的有效站点数为28 343个。

3.1 变异函数拟合的检验

随机选择靠近南京地区的一个网格点，其对应地理坐标为东经118.78，北纬32.05，以该点为圆心（即待插值点），西北方向至正北方向的45夹角为选定的方向，以6km为单位距离，计算此区域内各扇区中的自动站实测温度的半变异函数，结果如图2所示。

图2 实验半变异函数图

Fig.2 A schematic diagram of experimental

semivariogram function

图2中横坐标为单位距离数，纵坐标为半变异函数值，从该图函数值的分布特征可以看出，使用对数函数模型可以得到相对较好的拟合效果，该模型的一般公式为：

（7）

根据式（7），由当前实验数据推算得到，c=0.423，c0=-0.189，拟合出的曲线如图3所示。

图3 实验半变异函数拟合结果

Fig.3 A fitting result of experimental

semivariogram function

3.2 计算速度的检验

在对本文算法进行编码时，还考虑了提取循环语句中反复相同计算的部分，尽可能降低算法的时间复杂度。例如，将每个格点找邻近自动站的代码，改进为先计算出各个1×1的单位经纬网格周边所包括的邻近自动站，将此结果记录下来，而当每个格点找站时，直接读取其所在单位经纬网格预先遍历好的邻近自动站，从而大幅减少无意义的重复遍历和计算。本实验所使用服务器的硬件环境为：Intel Xeon E5-2697 V2 CPU，64GB DDR3内存，SAS接口15 000转/分硬盘。通过在程序中加入记时器，得到每秒处理的格点数量如图4所示。

图4 计算速度对比

Fig.4 Computing speed comparison

图4中，横坐标为计算流逝的时间，单位为秒，纵坐标为每一秒时间内插值计算出的格点数。可以看出，优化前每秒大约只能完成150个格点的计算，以此速度计算完所有G1个格点至少需要42分钟，而采用本文优化方法后，每秒可完成约3 200个格点的计算，所有格点的插值计算历时约115秒。图中还可以看出，在最开始的4秒“优化后”算法完成的插值格点数量为0，这是因为此阶段程序在进行自动站分块、单位经度网格遍历自动站等预处理工作。由此得出优化前后计算速度的对比结果，即优化后的Kriging插值速度提升了约20倍。

3.3 插值结果的检验

为了检验插值结果的正确性，本文使用交叉检验法[18]，随机选取若干自动站实测温度数据，分别与这些自动站最邻近的格点插值温度进行比对，结果如图5所示。

图5 站点实测值与格点插值结果比对

Fig.5 The AWS measured values with the

interpolation results

图5中，深色线条是自动站实测温度，浅色线条是格点插值出的温度。横坐标标注的是各个自动站的站名，纵坐标是温度值，单位为摄氏度。为更加直观地反映Kinging与其他插值方法的优劣，使用反距离权重算法对样本数据进行插值，再对上述相同站点作交叉检验，这两种插值算法的输出结果与自动站实测温度值的绝对误差如图6所示。

图6 两种插值结果与实测值的绝对误差比对

Fig.6 Absolute errors of two kinds of interpolation

results and the measured values

从图6可以看出，Kriging插值的效果普遍好于反距离权重法，Kriging插值结果与实测温度值的绝对误差一般在0—0.45℃，该波动主要有两方面原因。一是由于各个格点与比对的自动站之间距离不同，通常情况下，两者距离越近，插值结果越接近实测值，反之亦然。二是由于各个自动站的海拔高度不同，一般认为在大气非逆温的情况下，近地面每升高100米，温度下降0.65℃，如在青藏高原地区，相邻的自动站之间可能有数百米甚至1km的高度差，在不考虑温度—高度修正的情况下，这对插值的正确性造成了较大的影响。

4 结论（Conclusion）

本文利用Kriging插值算法在空间相关性方面的优势，针对自动站气象要素在空间分布上的规律特征，提出了一种改进的Kriging插值方法，阐述其计算的算法流程。通过实验测试表明，改进后的方法在计算速度上得到约20倍的提升，插值结果与实况监测数据对照，其准确性优于反距离权重插值方法。特别是该方法在速度方面的大幅提升，为自动站数据作为精细化实时初始场资料的一部分应用于数值模式天气预报等业务系统成为了可能。

参考文献（References）

[1] 刘彤，闫天池.我国的主要气象灾害及其经济损失[J].自然灾

害学报，2011，20（2）：90-91.

[2] 吴洪.临近预报系统SWAN培训手册[M].中国气象局气象干

部培训学院，北京，2012.

[3] 刘高平，叶金印.基于分布式架构的县级气象业务系统设计与

实现[J].软件工程师，2015（04）：36-38.

[4] 华韵子，林红.基于RIA的WebGIS在天气预警系统中的应用

研究[J].软件工程师，2015（09）：37-39.

[5] 邓林华，等.基于双线性插值的图像实时消旋系统的设计与实

现[J].计算机与现代化，2010（06）：64-66.

[6] 段平，等.自适应的IDW插值方法及其在气温场中的应用[J].

地理研究，2014（08）：1417-1426.

[7] Krige，D.G.Some Basic Considerations in the Application of

Geostatistics to the Valuation of Ore in South African Gold

Mines[J].Journal of the South African Institute of Mining and

Metallurgy，1976.

[8] 李迪，等.基于Kriging与IDW融合算法的氧化带DTM拟合方

法[J].金属矿山，2013（10）：93-96；100.

[9] 范玉洁，等.降雨资料Kriging与IDW插值对比分析—以漓江

流域为例[J].水文，2014（06）：61-66.

[10] 徐爱萍，胡力，舒红.空间克里金插值的时空扩展与实现[J].计

算机应用，2011（01）：273-276.

[11] 王兴，等.基于OpenCL的雷达外推算法改进与优化[J].计算

机与现代化，2014（08）：81-86；90.

[12] Ying-jie XIA，Li KUANG，Xiu-mei.Accelerating geospatial

analysis on GPUs using CUDA[J].Journal of Zhejiang

University-Science（Computers & Electronics），2011（12）：990-999.

[13] 周丰，郝泽嘉，郭怀成.香港东部近海水质时空分布模式[J].环

境科学学报，2007（09）：1517-1524.

[14] 翟进乾.克里金（kriging）插值方法在煤层分布检测中的应用

研究[D].太原理工大学，2008.

[15] 占文俊，等.一种InSAR大气相位建模与估计方法[J].地球物

理学报，2015（07）：2320-2329.

[16] Q.ZHU，H.S.LIN.Comparing Ordinary Kriging and Regression

Kriging for Soil Properties in Contrasting Landscapes[J].

Pedosphere，2010（05）：594-606.

[17] 冯桂莲.有限元数值解法在MATLAB中的实现及可视化[J].

软件工程师，2015，18（01）：61-64.

[18] Yongping Zhao，Kangkang Wang.Fast cross validation for

regularized extreme learning machine[J].Journal of Systems

Engineering and Electronics，2014（05）：895-900.

作者简介：

王 兴（1983-），男，博士生，讲师.研究领域：气象信息技术

与安全.

张琳焓（1993-），女，本科生.研究领域：大气科学.

王丽娟（1983-），女，硕士，讲师.研究领域：短期天气预报.

张潇潇（1987-），女，硕士，讲师.研究领域：气象信息技术与

教育.