山西省逐日日平均气温插值方法适用性评估＊

卫倩倩，陈霄健，闫佳逸，张 宁

（山西省气象信息中心，山西 太原 030006）

在全球变化研究中，地学模型和气候模型等相关研究的重要参数之一，是以高分辨率、格点化的气候数据作为环境因子[1]。而气温是气候变化研究中最基本指标之一，它不但在气候变化检测研究中有着基础地位，同时又与气候预测、气候模式以及其他一系列的气候变化问题研究有十分紧密的关系[2]。由于经纬度、海陆分布以及地貌特征与下垫面特性的差异等不同，体现在气温的空间分布上，会表现为具有明显的区域特征[3]。通过地面气象站观测到的气温数据在空间上是离散的、孤立的点，只能代表观测站点有限区域内的气温状况，难以反映气候空间变化连续过渡的基本特征，站网密度的增加也不能有效地将站点气温数据连接成面。因此，为获得格点化的气候数据，通常会选择进行空间插值。

空间插值的实质是通过已知的离散样点的数据，按照某种数学关系估算出未知点的数据。近年来，研究人员针对空间插值方法展开了大量的研究[4-6]，常用的空间插值方法主要有克里金法、反距离加权法、泰森多边形法、多项式回归法、样条函数法等。张海平等[7]针对ArcGIS软件中7种常用的空间插值方法，从不同视角进行适用性分类，比较各插值方法的差异，为研究人员如何选择合适的插值模型提供参考。使用薄板局部光滑样条插值方法的ANUSPLIN空间插值软件，相较于传统的插值方法，只将空间分布作为观测数据的函数，同时考虑了相关因素的影响，插值结果也更符合气象要素特点，插值结果的准确度也有所提高。

目前使用ANUSPLIN软件进行气候数据插值的应用较多，对较大时间尺度（月、年）或较小时间尺度（日），或不同区域面积上的不同气象要素的插值研究均有所涉及[1，8-11]。本研究利用ANUSPLIN软件，对山西省2019年国家级地面气象站点的逐日日平均气温进行插值，并与采用普通克里金法、反距离权重法插值的结果进行对比，选取绝对误差、相对误差2个评估指标，分析比较3种插值方法的结果差异，找出更加适合山西省的逐日平均气温空间插值的方法，同时对省级气象观测站日平均气温的缺测数据的插补进行思考，提高其数据完整性。

1 研究区域和数据来源

1.1 研究区域

本文的研究区域为山西省，地貌如图1所示。山西省地处黄河流域中部，因地属太行山以西，故名山西。东与河北省为邻，太行山作为天然屏障；西、南部以黄河为堑，与陕西省、河南省相望；北与内蒙古自治区毗连。地理位置范围在北纬34°34′～40°44′，东经110°14′～114°33′之间，面积15.67万km2。山西境内是典型的由广泛黄土覆盖的山地高原，地势东北高西南低，地貌类型复杂多样，有山地、丘陵、台地、平原，山多川少，山地、丘陵面积占全省总面积的80.1%，大部分地区海拔在1 500 m以上，最高点为五台山主峰叶斗峰，海拔3 061.1 m，为华北最高峰（信息来源于山西省人民政府网站省情概貌）。

图1 山西省地貌图

1.2 数据来源

本文研究所使用的数据为山西省2019年全年共874个站点的日平均气温资料，其中包括109个国家级地面气象站（基准站、基本站、原一般站）和764个省级气象观测站（国家级业务考核站）。同时，为尽可能地减小由于地缘边界站点减少造成的插值误差[12]，选取周边省份的136个国家站数据同时参与插值运算，经纬度范围为北纬34°～41°，东经110°～115°，共计245个国家级站点数据。以上数据均来源于山西省气象信息共享平台，数据内容包括站点的经度、纬度、海拔高度和日平均气温。数据的准确性对插值结果影响较大，因此调取的数据均已通过气象资料业务系统（MDOS）质量控制检查。由于气温的分布与地形、海拔高度等密切相关，在插值过程中引入高程数据作为协变量，所需DEM数字高程数据下载自地理空间数据云网站（www.gscloud.cn），选取SRTMDEMUTM 90 M分辨率数据高程数据产品。

2 研究方法

2.1 空间插值方法

2.1.1 薄盘光滑样条法

薄盘光滑样条（Thin Plate Smoothing Spline，TPS）法是对样条函数法的曲面扩展，常用于不规则分布数据的多变量平滑内插，利用光滑参数来达到数据保真度和拟合曲面光滑度之间的优化平衡[11]。局部薄盘光滑样条是对薄盘光滑样条原型的扩展[12]，它除普通的样条自变量外，允许引入线性协变量子模型，如温度和海拔之间、降水和海岸线的关系等。

局部薄盘光滑样条的理论统计模型为：

式（1）中：zi为位于空间i点的因变量；f为要估算的关于xi的未知光滑函数；xi为d维独立变量；b为yi的p维系数；yi为p维独立协变量；ei为随机误差[13]。

函数f和系数b通过公式（2）的值最小来获得，即由最小二乘法估计来确定：

式（2）中：ρ为光滑参数；Jm（f）为函数f（xi）的粗糙度测算函数，定义为f的m阶偏导，通常由广义交叉验证GCV的最小化、最大似然法GML的最小化、或期望真实平方误差MSE的最小化来确定。

ANUSPLIN软件是基于普通薄盘和局部薄盘样条函数对多变量数据进行内插的工具，能同时进行多个表面的空间插值，包括通过提供全面的统计分析、数据诊断等过程，同时也支持灵活的数据输入和表面查询功能。该软件包含了18个插值模型，研究人员可以根据模型诊断结果，为不同气象要素选择最佳的插值模型。

2.1.2 普通克里金法

普通克里金（Ordinary Kriging，OK）法是从地统计学的克里格法中演化出的一种插值方法，它以区域性变量理论为基础，半变异函数为分析工具，根据待估样点有限邻域内若干已测定的样点数据，对未知样点进行线性无偏、最优估计。它不仅考虑预测点位置与已知数据位置的相互关系，而且还考虑变量的空间相关性，是目前应用最广的克里金插值方法[14-15]。其算法如下：

式（3）（4）中：z（x0）为x0点的预测值；z（xi）为xi点的测量值；λi为测量值对预测值的权重系数。

普通克里金法的优点是可估计测试参数的空间变异分布以及估计参数的方差分布，但计算步骤烦琐，计算量大，且变异函数需要根据经验人为选定[14]。

2.1.3 反距离权重法

反距离权重法（Inverse Distance Weighted，IDW）是一个均分的过程，通过对采样点进行线性加权来决定输出栅格值的大小，基于相似相近原理，加权与距离成反比[14-15]。该方法以插值点和样点的距离为权重进行加权平均，算法如下：

式（5）中：z为估算值；n为参与插值的样点数目；zi为第d i个样点的观测值；di为插值点到第i个样点的距离；p为用于计算距离权重的幂指数，一般p=2。

IDW是根据点数据生成栅格图层的常用方法，但它易受样本点极值的影响，计算中常出现一种孤立点数据明显高于周围数据点的情况，即“牛眼”现象[14]。

2.2 研究步骤

选用研究区内245个国家级站点2019年的日平均气温数据作插值运算，为评估插值结果的准确性，同时对省级气象观测站数据的准确性进行校验，选取765个省级气象观测站作为校验站点。在进行插值前，对数据的完整性再次进行检验，剔除其中日平均气温缺测的站点。

利用ANUSPLIN软件进行TPS插值。首先，将参与插值的国家级站点的日平均气温数据、DEM高程数值进行格式转换，编写软件运行所需的cmd文件，运行splina模块，得到相应的日志文件。对照最佳模型判断标准，通过分析选择，从ANUSPLIN软件的18个待选插值模型中，确定模型9为最佳模型，即选用经度、纬度、高程（m）作为自变量，样条次数为4，进行插值运算，得到输出残差文件、光滑参数文件、表面文件、列表文件、拟合表面系数的误差协方差文件等。其次，根据得到的表面文件、拟合表面系数的误差协方差文件，设定相应的输出格式，编写运行lapgrd模块所需的cmd文件，得到最终的插值表面文件及插值预测误差表面等相关文件。

利用ArcGIS软件进行IDW和OK这2种插值。将参与插值的国家级站点的日平均气温数据，经格式转换，保存为shapefile格式文件，再依次选择IDW和OK模块工具，进行相应的插值运算，得到插值表面文件。

分析对比3种插值方法得到的插值表面文件，同时对使用TPS方法得到的插值预测误差表面文件进行分析。

根据经纬度信息，依次将TPS、IDW、OK这3种方法得到的插值表面数据提取到相应的省级气象观测站点，选用相对误差、绝对相对误差作为评估指标，对3种插值方法的结果进行评估。

日平均气温属于计算值，对于误差较大的站点，再结合气象资料业务系统（MDOS）进行分析，查看是否有因小时值缺测或质控码标识为错误导致的日平均气温计算疑误，对省级气象观测站的数据的准确性进行二次验证，同时也是对插值结果准确性的验证。

3 结果与分析

选用典型月份的每月15日，即01-15、04-15、07-15、10-15这4 d的逐日平均气温数据进行插值结果分析。

3.1 插值表面对比分析

图2中的（a）（b）（c）（d）依次为使用TPS、IDW、OK这3种插值方法得到的01-15、04-15、07-15、10-15这4 d的平均气温插值表面。从图中可以明显看出，使用TPS方法进行插值，得到的插值表面温度分布与地形走势更为一致，同时纹理清晰，分辨率更高、细节更突出，对于一些地势起伏地带的温度变化也有很好地体现，可以更好地反映气温随海拔高度的变化。使用IDW和OK这2种插值方法得到的插值表面，温度变化趋势不明显，不能很好地体现山西境内地势的高低不同对温度产生的影响，其中，使用IDW方法进行插值时，“牛眼”现象尤为明显，但对于小范围的地势起伏还是可以较好地反映出温度随海拔高度的变化。而使用OK插值方法得到的插值表面，对于小范围的地势起伏地带，高低温区别不明显，得到的插值表面趋向于平滑的连续曲面。但三者都对海拔较高、温度偏低的东北部一带和海拔较低、温度也相应偏高的中南部临运盆地一带的平川地区，这种大范围的、地势高低有明显不同的高低温区有很好的指示。

图2 日平均气温插值表面

3.2 TPS法预测误差表面分析

图3中的（a）（b）（c）（d）依次为01-15、04-15、07-15、10-15用TPS插值方法得到的预测误差表面图。从图中可以看出，由于引入了周边省份的国家级站点数据进行插值，边缘地带的插值结果的准确性也得到了保证。插值误差较大的区域主要分布在中部偏西一带，对比站点分布图及DEM数字高程图，该范围内海拔较高，且站点分布间隔相对较远，同时参与插值的周边省份站点也偏少，因此，预测误差相对较高，这点也与前人研究结论相同[1]。

图3 TPS插值的预测误差表面

3.3 评估指标分析

图4和图5分别为用TPS、IDW、OK这3种插值方法得到的01-15和07-15的绝对误差分布图。经统计分析，01-15的插值绝对误差小于等于2℃的比例依次为91.49%、78.40%、79.97%，平均绝对误差依次为0.93℃、1.35℃、1.29℃；07-15的插值绝对误差小于等于2℃的比例依次为92.37%、78.55%、70.66%，平均绝对误差依次为1.06℃、1.49℃、1.75℃。由此可见，用TPS方法插值得到的结果明显优于用IDW和OK方法的插值结果，从图6和图7中，可以得到相同的结论。

图4 01-15绝对误差分布图

图5 07-15绝对误差分布图

图6 01-15相对误差分布图

图7 07-15相对误差分布图

3.4 MDOS入库数据比对分析

选用01-15利用TPS插值方法得到的插值数据与MDOS入库小时数据进行分析判断。

01-15有10个省级气象观测站日平均气温为缺测，依次查询MDOS入库的小时气温数据情况，其中有3个台站仅有一个定时时次数据未正常入库，导致日平均气温缺测。还有2个台站缺测时次也小于6个时次，但根据省级气象观测站的日平均气温计算规则，只要有1个定时时次缺测，则日平均气温缺测，因此，这样会造成一些台站由于个别定时时次的缺测造成日平均气温的缺测，因此，这时可以考虑用插值数据来进行缺测补值，当然也要考虑对插值误差的订正。

4 结论与讨论

利用局部薄盘光滑样条函数法进行日平均气温空间插值，与反距离权重法和普通克里金法相比，由于考虑了温度同海拔高度变化的关系，在插值过程中，引入高程作为协变量，插值结果更符合气温随地势高低的变化规律，精度更高、光滑度更好。但3种插值方法对于大范围的、地势高低有明显不同的高低温区的插值精度区别不大。

选用绝对误差、相对误差2个评价指标，对比分析3种插值方法的插值结果精度，发现用TPS方法插值得到的结果明显优于用IDW和OK方法得到的插值结果。

省级气象观测站的日平均气温计算方法为4次定时值的平均，只要有一次定时缺测，则日平均气温相应缺测。由于传输不稳定或仪器故障等原因，引起的部分时次缺测现象比较常见，可以考虑利用插值数据进行缺测补值，提高气象服务的数据完整性，对于如何进行误差订正，这部分还有待研究。

在下一步的研究当中，将以此研究为基础，选用TPS方法，利用ANUSPLIN软件完成山西省2019年逐日平均气温格点数据集的建立，并逐步扩展到长时间序列的高分辨率的其他气象要素格点数据集的建立，并尽可能选取更多的误差指标及各气象要素相应的特征值进行验证[6，16]，提高插值精度，为气候研究提供参考。