基于ANUSPLIN的黄土丘陵区气象要素插值方法对比与评估

肖 旭， 郑 诚， 丁成琴， 范陈哲， 白悦江， 林隆超，闫 婷， 高 宇， 史海静,5

（1.陕西省延安市气象局，陕西 延安 716000；2.西北农林科技大学草业与草原学院，陕西 杨凌 712100；3.西北农林科技大学水土保持研究所，陕西 杨凌 712100；4.西北农林科技大学资源环境学院，陕西 杨凌 712100；5.中国科学院水利部水土保持研究所，陕西 杨凌 712100）

气象因子是影响农业生产、人类活动以及生态环境的重要要素，尤其是在干旱半干旱地区。气象要素在农业和环境等相关领域研究较多，主要应用于模型预测和作物生产等方面［1］。关于气象的研究，最初人们依赖于气象站的雨量计监测，但由于气象站点数目较少，空间不完全覆盖，且时效性较差，在完整捕捉气候事件和描述空间变化等方面能力不足。随着空间处理和分析技术的发展，空间插值的方法被用来描述气象因子的空间分布，一定程度上实现了气象空间覆盖，提高了人们对气象空间变化的认识。

目前，插值方法较多，且在系统模拟和反演领域都有不同程度的应用，常见的方法有局部多项式法(Local polynomial interpolation)、样 条 函 数 法(Spline)、克里金插值法(Kriging)和反向距离权重法(Inverse Distance Weight, 简称IDW)等［2-4］。ANUSPLIN 是由Hutchinson［5］开发用于不同空间和时间尺度水文气象插值的方法，其主要特点是可包含线性协变量例如地形变量等，广泛应用于高山丘陵等地形变化起伏较大的地区［6-7］。此外，ANUSPLIN 具有简单的操作流程，不需要对其参数进行单独的预先校准，生成的温度降水时间序列曲面具有较高的精度和可靠性［8］。近年来，ANUSPLIN 插值方法已被国内外学者广泛应用于气候数据的插值［9-10］。黄土丘陵区位于黄土高原中部，是我国生态恢复的重点区域之一，该区沟壑纵横，地形起伏多变，温度降水等气象因子受地形影响较大［11］，传统的插值方法多通过气象站空间距离进行估算，而黄土丘陵区多以山地为主，地形变化足以影响空间温度和降水的差异性分布［12］，单从距离的远近获得气象数据并不能真实的反映该地区的气象变化，从而无法满足植被恢复和重建工作的需求。因此，针对地形复杂的丘陵区，探索包含不同分辨率地形变化的气象插值方法对于进一步完善气象空间数据和相关研究具有重要意义。

本文采用ANUSPLIN 空间插值方法，对延河流域2010—2020年的温度降水进行空间插值，同时选取25 m、90 m、1 km 三个不同精度的数字高程模型(Digital Elevation Model，DEM)作为协变量，对比分析各个尺度温度降水插值结果的精度，评估ANUSPLIN 插值方法在黄土丘陵沟壑区的适用性，并选择区域温度降水最优数据集，补充丘陵区气象插值精度不够的短板，为准确获取地形起伏较大的丘陵区气象空间插值数据提供可参考的理论依据和方法支持。

1 数据与方法

1.1 研究区概况

本研究以典型黄土丘陵区的延河流域为研究区域(36°23′～37°17′N，108°45′～110°28′E)，该区覆盖面积为7687 km2，海拔454～1765 m，该流域从东南向西北依次经过延长、延安、安塞和志丹，是典型的丘陵区（图1）。气候全年较干，夏季多雨，冬季降水较少，处于暖温带大陆性季风半湿润气候向温带半干旱气候的过渡区，全年降雨量496 mm，年平均气温8.8 ℃。从东南向西北，降雨、温度呈现出明显的梯度变化特征。土壤以黄绵土为主, 植被有森林(刺槐和辽东栎)、灌木丛和草地等。

1.2 数据来源与处理

1.2.1 气象数据 基础气象数据来自延河流域及周边站点。数据获取时间为2010—2021年，数据类型为日值数据，包含日平均温度、最高温度、最低温度、日降水量等气候要素，本研究通过多年（11 a）的年均温度和年均降水进行插值。由于是长时间序列数据，部分数据存在异常或缺失等，因此，数据经过筛选、剔除异常数据，最终保留105个气象站（图1）。

1.2.2 DEM 数据 延河流域25 m 分辨率DEM 是通过对日本ALOS 卫星搭载的PALSAR 传感器25 m DEM 数据进行处理得到的。ALOS 25 m DEM 可从NASA 的ASF DAAC (https://search.asf.alaska.edu/#/)免费下载。90 m 和1 km 分辨率DEM 数据来自CGIAR-CSI SRTM中国区域数据。

1.3 研究方法

1.3.1 ANUSPLIN 模型基本原理 ANUSPLIN 是一种使用薄板平滑样条对噪声多变量数据进行透明分析和插值的工具。薄板平滑样条函数法实际上是利用一个平面去拟合所有的气象站点，使得该曲面能通过站点构成“样条”，得到逼近所有控制点弯曲最小的光滑曲面，即利用M 次多项式，对曲线进行分段拟合，用局部的薄板来绘制分段连续的曲线。TPS 是对样条函数法的曲面扩展，插值过程中利用最优的光滑参数实现逼真度和光滑度最佳平衡，保证了精度可靠，同时插值曲面光滑连续［13］。此外，ANUSPLIN 还允许转换自变量和因变量，并允许处理具有缺失数据值的数据集。 当对因变量应用变换时，ANUSPLIN 允许拟合曲面的反向变换，计算相应的标准误差，并纠正这些变换引起的小偏差。研究发现，这在将表面拟合到降水数据和其他自然为正或非负的数据时特别方便［14］。

薄板样条估计z（xi）由一个适当光滑的函数计算，该函数最小为:

式中：yi为观测i的观测数据值；ƒ(xi)为被估计的样条函数；λ为一个正数，称为平滑参数；Jm（ƒ）是用m阶偏导数定义的函数ƒ的粗糙度的度量［m阶偏导数称为SPLINA 输入中的样条的阶(它也被称为粗糙度阶）］［10］。

局部薄盘光滑样条（Partial thin plate smoothing splines）的理论统计模型公式：

式中：Zi为位于空间i点的因变量；T为迭代次数；ƒ(xi)为估算关于xi的未知光滑函数，xi为独立变量；yi为p维独立协变量；b为yi的p维系数；ei为随机误差［10］。

上式中，函数ƒ和系数b是通过最小二乘估计来确定的：

式中：wi为权重的已知局部相对变异系数；Jm(ƒ)为函数ƒ(xi)的粗糙度测度函数，定义为ƒ的偏导（称为样条次数，也叫粗糙次数）；λ为正的光滑参数，在数据保真度与曲面的粗糙度之间起平衡作用。在AUNSPLIN中通常用广义交叉验证GCV的最小化以及最大似然法GML的最小化来确定［15-16］。

1.3.2 插值精度评价指标 本研究选取数据的15%作为验证集，选取过程采用随机抽样，将随机站点数据与插值结果进行比较，以评估气象模拟的准确性和偏差［17］。三种分辨率的插值精度比较适合使用交叉验证法。本研究插值的效果采用均方根误差RMSE来评估，表达式如下：

式中：Zx和Zy分别代表气候数据实测值和插值结果；n为计算序列长度。RMSE 的值越小表示气象模拟效果越好。

此外，将实测值与预测值进行相关性分析，用确定性系数R2表示一元多项式回归方程拟合度的高低。R2取值分布在0～1之间，R2值越大，表示模拟效果越好。

2 结果与分析

2.1 插值结果

基于ANUSPLIN 分别获取了研究区2010—2021 年年均温和累计降水的插值曲面见图2，从图中可以看出，相较于1 km 的插值，90 m 与25 m 分辨率的插值结果细节表现得更加突出，能反映局部地形特征，也有较好的平滑度。

由图2 还可以看出，东部延长县逐日年均温度较高，西部靖边、志丹以及中部的安塞区逐日年均温度较低，三种分辨率的逐日年均温度整体呈“东高西低”的趋势。25 m 分辨率的温度插值范围为7.725～13.329 ℃（图2a），90 m 分 辨 率 为7.786～13.354 ℃（图2b），1 km分辨率为7.872～13.313 ℃（图2c）。

从降雨插值可以看出，低降雨的空间分布主要集中在延河流域西部的靖边地区、中部的宝塔区以及东部延长的少部分区域。25 m 分辨率的逐日降雨插值范围为1.062～1.847 mm，90 m 分辨率的插值范围为1.046～1.852 mm，两者最低降雨量相差0.016 mm，上限相差0.05 mm，较为接近，但1 km 分辨率的插值结果与25 m、90 m 的插值结果相差过大，下限为1.117 mm，比90 m 的结果高0.071 mm、比25 m 的插值结果高0.055 mm；上限为1.73 mm，比两者都低，与25 m 的插值上限相差0.117 mm，整体插值结果范围缩小。

从插值的温度降雨空间图中可以看出，细节特征比较突出，且有明显的“马赛克”纹路，与过去的Cokriging 插值相比，没有出现较大的斑块与“牛眼”，说明ANUSPLIN 有比较好的平滑度，且很好的拟合了地形特征，这与刘林等［18］的研究结论一致。此外，插值结果显示延河流域的中部、西北部降雨较低，东部降雨量较高，这在刘春利等［19-20］的研究结果中得以验证。且土壤水分变化也与之一致［21］。温度插值特征也符合以往气象站公布的数据规律。延河流域的西北部为草原地带［22］，整体呈降水低，温度低，东部为森林草原地带与森林带，降水高，温度高，与插值结果一致，说明ANUSPLIN 模型能较好的适应黄土高原丘陵沟壑区。

2.2 插值结果比较

本研究以25 m、90 m、1 km 三种分辨率的海拔数据为协变量对黄土丘陵沟壑区的温度降雨进行插值，相较于前期研究中采用的1 km×1 km 或500 m×500 m 分辨率［23］，本研究在一定程度上提高了精度，减小了插值误差。从图3 可以得出，对于不同DEM分辨率下温度插值结果，根据计算的均方根误差（RMSE），25 m 分辨率的插值数据精度较好，90 m 分辨率插值数据精度其次，1 km插值数据精度较差。而在确定性系数（R2）方面，25 m 分辨率的R2为0.65，拟合效果较好，1 km分辨率拟合效果其次，90 m分辨率拟合效果稍差。综上所述，在三种不同精度的DEM分辨率模拟情景下，温度插值精度排序为：25 m＞90 m＞1 km。贾洋等［24］的研究表明，DEM 分辨率越高，ANUSPLIN 模型插值精度也越高，温度插值结果也遵循这一规律，高精度的地形数据可以有效反映空间差异，也可提高温度空间插值的精确度。

图3 基于不同分辨率DEM温度降水观测值与预测值分析结果Fig.3 Analysis results of observed and predicted values of temperature rainfall based on DEM with different resolutions

降雨插值结果精度检验按照年均温来检验，空间数据取365 d。对于降雨插值结果的RMSE，90 m<25 m<1 km，即90 m 分辨率的插值精度最高，1 km 的插值精度较低。同时，90 m 分辨率的插值结果与实测结果拟合时，90 m 分辨率也是三种分辨率中最大的，R2=0.833，拟合效果良好，其次是25 m 分辨率，1 km 分辨率的拟合效果最差。这与温度插值规律并不相同，笔者推测出现这一结果的原因可能是由于黄土丘陵区山地较多［25-27］，降水容易在沟地汇集，导致站点数据存在误差。尤其是黄土丘陵区普遍坡长在70～200 m之间，90 m的空间精度与地形变化精度恰好一致，确定性检验仍需后续工作进一步验证。但是，该结果也说明在丘陵区设置气象站点要考虑地形变化。

高分辨率数据集(包括遥感数据和长时间序列的气候数据)是目前黄土丘陵区相关研究迫切需要的。且近年来干旱区气候变化较大［28］，实时气象数据对于气候-植被研究和人类活动具有重要价值。这一气象数据和空间分布研究对于干旱区气候变化研究非常宝贵，同时高精度分辨率的数据集可以为该区域生态恢复和气候变化研究提供基础数据支撑［29］。

3 结论与讨论

本研究选取黄土高原典型丘陵区延河流域作为研究区，以研究区内及周边105 个气象站2010—2021 年逐日气象数据为基础资料，以1 km、90 m、25 m 三种分辨率的数字高程数据为协变量，利用ANUSPLIN 插值方法对研究区不同地形分辨率气象数据进行空间插值，通过气象站站点交叉验证方法，对比分析各分辨率数字高程数据的插值精度水平。研究发现数字高程数据的分辨率影响气象空间插值的结果，具体结果如下：

（1）温度插值显示，延河流域东部延长地区温度较高、西部温度较低符合以往气象站数据规律［30］；降水量空间插值显示中部与西北部较低，与该区域蒸散发时空变化趋势基本特征一致［31］，延安市周边本研究降雨量较少，其主要原因可能是人为干扰导致气象站点数据存在误差，以往研究也表明，延河流域的延安市是人类活动较为频繁的区域［32］。在赵美亮等［33］在西北地区水文空间分布研究中表明，土地利用改变会影响水文要素分布。总的来说，ANUSPLIN 模型对黄土丘陵沟壑区温度降雨空间插值有较好的适应性。

（2）在三种不同的DEM 分辨率模拟场景下，温度插值精度排序为：25 m＞90 m＞1 km，符合数字高程数据分辨率越高，插值精度越高［34］；降雨插值精度排序为：90 m>25 m>1 km,不完全符合以上结论，表明降雨数据空间插值存在最佳分辨率的数字高程数据，最佳数字高程数据与该区域地形特征有关。在野外观测中发现，黄土丘陵区山脉坡长多60～300 m之间，与该结果具有一致性，坡度和坡长等地形数据可能更适合描述降雨数据信息，而不是25 m 的海拔数据。

本研究基于ANUSPLIN 插值方法，结合数字高程数据进行多种空间尺度的空间精度验证，该结果与观测结果拟合较好，说明了结果具有可靠性，所得结论能够为丘陵区插值的应用提供参考。然而，本研究仅考虑了海拔数据作为协变量，然而坡向、坡位等地形因子对插值结果影响也不可忽视，将多种协变量考虑在ANUSPLIN 插值方法中，可以为气象插值精度优化方法提供更有意义的参考依据。