李任君，高懋芳，李 强，李百寿

（1. 桂林理工大学测绘地理信息学院，广西桂林541006；2. 中国农业科学院农业资源与农业区划研究所/农业农村部农业遥感重点实验室，北京100081；3. 江苏师范大学地理测绘与城乡规划学院，徐州221000）

0 引言

降水是气象观测的重要因素之一，高精度的降水数据是气象、农业干旱分析的重要基础［1］。气象数据表现为空间内的不均匀性、不连续性，通过站点测得的气象数据，以一定的空间插值算法求解未知点的降水量是获得降水空间数据的重要途径［2-3］。插值过程中，由于降水具有间断性和空间不连续性的特点［4］，影响因素众多，且存在很强的季节变化，与气温空间插值相比，降水数据空间插值不确定性更大，更加复杂困难［5］。在所有插值方法中，基于统计插值技术的克里金法和薄盘光滑样条函数法最为适用，这些技术建模只将空间分布作为观测数据的函数而不需要其先验知识和物理过程，能提高插值的准确度［6-8］。一些研究表明，利用平滑样条法对降水的插值效果明显优于其他方法［9-10］。钱永兰将局部薄盘样条插值结果分别与反向距离权重法和普通克里金法的插值结果进行对比，显示ANUSPLIN 的插值误差最小，适当增加站点数量和提高DEM 精度可进一步提高ANUSPLIN 软件的插值精度［11］。刘正佳利用ANUSPLIN 软件，以高程、距海岸线距离及高程—海岸线距离分别为协变量对多年平均降水量进行空间插值，发现以高程为协变量比以距海岸线距离和高程—海岸线距离为协变量所获得的插值结果更加精确［12］。谭剑波采用克里金和ANUSPLIN 空间插值软件对青藏高原东南缘年均温度和年累计降水进行插值，发现降水插值中将高程作为协变量并不能体现ANUSPLIN 插值的优势，但仍然优于克里金空间插值方法［13］。徐翔利用ANUSPLIN 软件对中国康滇区多年月平均温度和降水进行空间插值，并与克里金法和反距离加权法比较得出，ANUSPLIN 更适用于复杂山地环境下气象要素的空间插值［14］。

虽然目前已有大量的研究在不同时空尺度上分析了我国不同地区和流域降水要素的插值，但是还没有一种插值方法适用于所有的气象要素，应根据研究目的和研究区域的自然地理地形特征选取适宜的插值方法及参数［15］。由于对平原地区的降水插值方法研究甚少，基于此，文章对黄淮海平原458 个气象站点的多年平均降水数据进行处理，将连续365 d 的30 年平均降水数据每16 d 的降水累积量为一个试验样本，共得到23 个降水量数据样本。以分辨率分别为90 m 和1 km 的数字高程模型（DEM）作为协变量，采用三变量局部薄盘样条函数，通过ANUSPLIN 软件进行降水插值。使用预留验证集对该种插值方法在黄淮海地区降水插值的精度进行评价，并将插值结果与克里金插值方法进行对比，以期为相关科学研究提供理论依据和数据支撑。

1 材料与方法

1.1 研究区与数据来源

黄淮海地区是我国重要的粮、棉、油生产基地，对国民经济发展起着举足轻重的作用［16］。该地区属温带大陆性季风气候，冬季寒冷干燥，夏季高温多雨，年平均降水量为500～600 mm，全年降水主要集中在夏季，秋、冬、春3 季均为水分亏缺期，经常出现先涝再旱，甚至连年干旱，是全国受旱面积最大的区域之一，严重制约着该地区农业生产力。

研究从中国气象数据共享网（http：//cdc.cma.gov.cn/）收集了研究区内457 个气象站点 1981 年 1 月至 2010 年 12 月每日（08：00- 次日 08：00）的降水量数据，选取 428 个气象站点作为试验站点输入模型进行插值，预留29 个均匀分布在研究区的气象站点作为验证站点，研究区地理位置以及样点分布见图1。试验中采用的90 m 和1 km 分辨率的数字高程模型数据（Digital Elevation Model，DEM）来源于CGIAR-CSI SRTM 中国区域数据（http：//srtm.csi.cgiar.org/），其具有精度高、应用广泛、获取免费等优点。

图1 研究区DEM 及气象站点分布Fig.1 DEM of the study area and distribution of meteorological stations

1.2 数据处理

应用ANUSPLIN 软件，首先需要将降水数据通过SPSS 软件处理成程序要求的标准格式，包括站点代码、经度、纬度、高程信息和与气象有关的属性数据集，对数据进行检查修正处理，剔除重复站点、缺值、错值，数据格式详见表1，以9 月30 日—10 月16 日样本数据其中5 个站点的数据为例。ANUSPLIN 软件所需要的DEM 数据需要转换成ASCII 形式书写，利用ArcGIS 软件对DEM 数据进行拼接裁剪、投影变换为ALBERS投影，然后利用转换工具转为ASCII 格式文件，同时检查和保证转换过程中数据的正确性。

表1 1981—2010 年平均降水数据格式（rain.dat）Table 1 Data format of average precipitation from 1981 to 2010

1.3 模型检验方法

为了体现ANUSPLIN 对气象数据插值的优势，通过对气象站点进行验证，计算均方根误差（RMSE）来评价ANUSPLIN 和常用的克里金插值方法的精度，公式为：

式（1）中，Zx和Zy分别代表降水量实测值和预测值，n为测站数。

ANUSPLIN 软件提供了用于判别误差来源和插值质量的统计参数，该研究选取广义交叉验证均方根（GCV）、残差均方根（MSR）和方差估计均方根（VAR）等指标进行最优模型选取，具体公式详见ANUSPLIN 的用户指南［17］。

2 结果与分析

2.1 最优模型选取

由数据预处理部分得到ANUSPLIN 软件特定格式的降水数据和DEM 数据。根据数据和ANUSPLIN 的用户指南编写运行ANUSPLIN 的SPLINA 和LAPGRD 这两个程序模块的批处理命令文件。APLINA 程序模块负责生成曲面拟合结果的系数以及误差统计文件，LAPGRD 程序模块通过输入SPLINA 程序模块生成的表面系数文件得到插值表面。通过对前人工作经验的总结［18-19］，该试验选取了适用于降水插值的6 个模型（变量、协变量和样条次数多种组合），以经纬度、高程作为变量对多年平均降水数据进行试验，6 种模型参数设置见表2。ANUSPLIN 在日志文件（Log file and List file）中提供了一系列用于判别误差来源和诊断插值质量的统计参数，如均值、方差、标准差、拟合曲面参数的有效数量估计Signal（信号自由度）、剩余自由度Error、光滑参数RHO、广义交叉验证GCV、期望真实均方误差MSE、均方残差MSR、方差估计VAR 及其平方根等，这些统计参数和模型诊断参数用来判断站点是否充分，数据是否存在错误，是否存在合适的模型可以成功插值，帮助用户轻松选择出最佳空间分布模型［20-21］。利用黄淮海平原428 个站点的降水数据运行6 种模型，对产生的判别误差来源和插值质量的统计参数进行整理选取最优模型，以9 月30 日至10 月16 日的数据为例，误差统计参数见表3。

表2 六种待选薄盘光滑样条函数模型Table 2 Six optional smooth spline function models for thin disks

表3 平均降水数据ANUSPLIN 插值输出统计校验Table 3 Comparisons of the ANUSPLIN output statistics for rainfall interpolation

依据最佳模型判断标准：基于广义交叉验证或最大似然法，统计误差最小，信噪比（SNR，信号自由度与剩余自由度之比）最小，信号自由度小于站点的1/2，SPLINA 生成的误差文件中无*符号。Signal 表示拟合曲面的复杂程度，大于观测站点的一半表示数据存在显著错误或者是存在不适合表面模型的数据，拟合过程找不到最优光滑参数，这些情况在误差统计文件中以*符号标出。可以根据插值产生的误差统计文件，检验并消除原始数据在位置和数值上的错误，发现程序执行过程中的问题，并做出及时的调整［22-23］。对运行SPLINA 与LAPGRD 的6 种模型产生的误差统计数据进行分析发现：高程作为协变量比高程作为独立变量的插值精度要高；其他几种模型的统计误差数据都相差不大，DEM 为90 m 的插值精度普遍高于分辨率为1 km 插值。当DEM 分辨率为90 m 时，以经纬度为变量、高程为协变量，样条次数设为4 的模型精度最高，所有数据都可以找到适合的表面模型，精度明显高于样条次数为2 或者3 的结果，其中模型6 为最优待用模型。

2.2 插值结果

应用ANUSPLIN 软件，以90 m 分辨率的高程数据作为第三变量，选用模型6 为最优待用模型，对1981—2010 连续30 年23 个时段的平均降水数据样本进行插值，在ArcGIS中将输出表面文件转换成tiff 图像格式文件，进行显示操作，结果见图2。

图2 1981 年1 月1 日至2010 年1 月4 日研究区全年降水分布图（单位：mm）Fig.2 Distribution map of precipitation in the study area from January 1,1981 to January 4，2010

从图2 空间分布上可看出西北部地区全年降水量较低，降水主要分布在南部地区，由南到北降水量呈递减趋势分布，5 月份开始，降水区域慢慢向北部扩大，8 月份时全区降水主要分布在北部和中部，南部降水量减少；之后降水分布主要在东部，西部降水减少，容易发生北旱南涝和东涝西旱，和实际的降水空间分布一致；时间分布上可看出黄淮海平原1 月份降水最少，7 月份降水最多，全年降水过分集中在7—8 月份，9 月份开始，全区降水量逐渐减少，发生先涝后旱的概率比较高，经常造成连年干旱。

2.3 插值精度分析

为了更好地体现ANUSPLIN 对气象数据插值的优势，利用常用的克里金插值方法对研究区6 月10—25 日期间的降水数据进行空间插值，如图3 所示，通过ANUSPLIN 获得的插值表面细节更加突出，能反映局部地形特征，使降水随高程的梯度变化更加直观，并且有比较好的平滑度。对冬季（1 月1—16 日）和夏季（6 月26—7 月12 日）的降水插值结果进行验证，以均方根误差（RMSE）为验证几种插值模型的标准，RMSE 通过插值获得的验证站点的降水数据与其观测值计算所得，图4 所示为冬季降水插值，采用ANUSPLIN 插值，DEM 分辨率分别为90 m、1 km 作为高程数据和克里金方法获得的插值结果的RMSE 分别为：0.379 mm，0.376 mm，0.451 mm；图5 为夏季降水插值，DEM 分辨率分别为90 m、1 km 作为高程数据和克里金方法3 种插值结果的RMSE 分别为：4.186 mm，4.237 mm，4.310 mm。

图3 基于克里金和ANUSPLIN 的降水插值栅格（a：克里金插值；b：ANUSPLIN 插值）Fig.3 Precipitation interpolation grid based on Kriging and ANUSPLIN

根据降水插值的误差比较结果得出：（1）无论是夏季降水还是冬季降水，ANUSPLIN的插值精度都高于克里金插值方法，夏季降水插值精度明显低于冬季降水插值；（2）90 m分辨率和1 km 分辨率的DEM 分别作为ANUSPLIN 的协变量输入时，在冬季降水量较少的数据插值结果中差别不是很大，在夏季降水插值中，以90 m 分辨率的DEM 作为高程输入模型获得的插值精度更高，DEM 栅格特征更明显。

图4 冬季降水插值结果和观测值的散点分布（a）ANUSPLIN 插值-90 mDEM；（b）ANUSPLIN 插值-1 kmDEM；（c）克里金插值Fig.4 The scatterplot of true value and predictive value of the precipitation interpolation in winter

图5 夏季降水插值结果和观测值的散点分布（a）ANUSPLIN 插值-90 mDEM（b）ANUSPLIN 插值-1 kmDEM（c）克里金插值Fig.5 The scatterplot of true value and predictive value of the precipitation interpolation in summer

3 结论

由于不同区域地形复杂多样，气候差异性大，气象站数量有限，影响降水数据插值精度的因素众多，插值使用的模型和方法都需要根据研究区域实际情况进行选择。该研究利用国内外广泛应用的专用气象插值软件ANUSPLIN，以高程作为插值变量因子对黄淮海平原连续30 年地降水数据进行插值，得到与遥感数据的时空分辨率相匹配的栅格数据，对插值数据进行验证，得出结论：（1）与克里金方法相比，ANUSPLIN 插值精度更高，可以更细致地描绘出黄淮海平原地形因素对降雨空间分布的影响，分辨率为90 m 的DEM 作为高程输入数据所获得的插值结果比分辨率为1 km 的插值结果精度相对高，但是前者试验时间明显长于后者；（2）夏季降水量集中的6—8 月份的插值精度比降水量较少的冬季降水插值精度低，虽然提高DEM 的分辨率，可以使夏季降水插值精度有所提升，但是并没有改善降雨量的增加与精度误差呈正相关的关系。通过使用ANUSPLIN 对黄淮海平原的降水插值方法的研究，探讨插值精度与DEM 分辨率的关系，获得较高精度的平均降水栅格数据，为黄淮海平原干旱分析并指导当地农作物灌溉生产提供重要的科学依据。