四川夏季降水量空间插值方法的比较

李 想，李国平

（1.成都信息工程大学大气科学学院，四川 成都 610225；2.气象灾害预报预警与评估省部共建协同创新中心，江苏 南京 210044）

引言

降水是天气预警预报、气候预测预估的最重要的气象要素之一［1-2］。同时获取精度更高的降水资料，可以更准确认识降水的空间分布［3-6］。

更加精确的空间降水数据，除了完善观测设备和观测站网之外，也依赖于更高质量的空间插值方法。插值技术的不断成熟改进了气象台站有限和空间分布不均的不足，但实际业务中没有一种插值方法是普适性的［7］，一般根据研究目的和区域地理特征采用适宜的插值方法。为寻求最佳插值方法，已有研究对同一区域的降水数据采用多种空间插值方法进行比较，例如克里金插值法、样条函数法、趋势面法等方法［8-10］。目前国内主要根据中国高密度气象站点的信息采用不同插值法［11］。例如中国地面降水0.5°×0.5°格点数据集（V2.0）（1961—2013）利用薄盘样条法，并且引入数字高程模型（digital elevation model，DEM）以尽可能地消除中国区域独特地形条件下高程对降水空间插值精度的影响。同时由于协同克里金方法能够引入其余影响因子，也经常被用于业务中［12-14］。例如安徽省的降水插值中效果最佳的是引入经纬度的协同克里金方法［15］。考虑DEM数据的协同克里金方法提高了湖南省的降水插值精度［16］。综合来看，协同克里金方法在降水空间插值上具有较佳的表现。

山区的降水空间插值更为复杂。传统插值依赖于样本空间密度，缺少气象要素的演变机理，忽略了数据的空间分布特征和地形高度等影响因素［17］。对于地形复杂且站点分布不均的山区，降水空间插值的结果往往不太理想［18］。四川山区受热带季风、副热带季风以及青藏高原环流系统的影响，夏季暴雨频发，极易造成泥石流、山体滑坡等次生灾害。因此，获取四川山区夏季精确的降水数据对防震减灾具有重要意义。本文重点针对传统插值效果较差的四川地区，分区域寻找各区域最优的降水影响因子组合并且对不同方法的插值结果进行比较，找出最优插值方法。

1 研究区地理概况、资料和方法

1.1 研究区地理概况

四川省（97°21′E—108°12′E，26°03′N—34°19′N）位于我国西南地区内陆，地处长江上游，总面积约4.86×105km2，地跨青藏高原、横断山脉、四川盆地等地貌单元，地势西高东低。西部是青藏高原东南缘和横断山脉的一部分，海拔为3000～4500 m；东部由盆地底部和盆地边缘组成，底部，海拔300～700 m，由成都平原、眉山-峨眉平原组成；边缘地区以山地为主，海拔为1500～3000 m（图1）。

图1 四川省自动气象站分布Fig.1 Distribution of automatic weather stations in Sichuan Province

1.2 资料

所用资料：（1）四川省气象局共157个自动气象站2010—2019年6—8月逐小时降水资料，用于计算10 a的多年平均夏季降水量；（2）DEM数据来源于先进星载热发射和反射辐射仪全球数字高程模型（advanced spaceborne thermal emission and reflection radiometer global digital elevation model，ASTER GDEM）30 m空间分辨率数据；（3）NDVI数据来源于LAADS（Level-1 and Atmosphere Archive & Distribution System）DAAC（Distributed Active Archive Center）提供的空间分辨率为250 m的MOD13Q1第6版产品，此次使用的数据是的植被指数（VI）值，有两个主要植被层，本文使用的是第一个归一化差异植被指数（normalized difference vegetation index，NDVI）。

1.3 研究方法

首先对四川地区进行区域划分，并对各因子与研究区和各分区的多年平均夏季降水量进行相关性分析，找出每个区域主要的影响因子，再将这些因子进行组合，与降水量进行多元线性回归分析，找出每个区具有最佳拟合效果的因子组合；其次分别对各区进行插值，并对结果进行交叉检验，找出适合各区域的最佳插值方法，同时根据插值结果，给出各区域降水主导的地理影响因子。

共使用6种插值方法，分别是反距离加权（inverse distance weighted，IDW）、径向基函数（radical basis function，RBF）、普通克里金（ordinary Kriging，OK）、协同克里金金（CoKriging，CoK）、局部多项式（local polynomial interpolation，LPI）和经验贝叶斯克里金（empirical Bayesian Kriging，EBK）插值。其中，IDW方法运算快、效率高，但外推能力差，适用于站点分布尽可能均匀且布满整个插值区域的样本数据集［19-20］。RBF工作量小且精度相对较高，适用于样本数据集大，地形平缓的情况［21-22］。OK法又称空间局部插值法，该方法计算速度较慢，适用于区域化变量存在空间相关性的区域［23-24］。OK是应用最广的克里金插值方法，而CoK是OK的扩展形式，它将主变量的空间自相关性和主辅变量之间的交互相关性结合起来，用于无偏最优估值中［19，25］。LPI是一种局部加权最小二乘拟合法，多用于解释局部变异现象、建立平滑表面和确定变量的小范围变异［26-27］。EBK是一种地统计插值方法，它与OK方法不同，是通过估计基础版变异函数来说明所引入的误差，因此大大降低了预测的标准误差［21］。

运用ArcGIS中的交叉验证法对各种插值法进行误差分析。其基本思想是将原始数据随机地分为训练集和验证集，先对训练集进行训练，再利用验证集来测试训练得到的结果，以此作为评价指标［28-29］。选取平均误差、平均绝对误差、均方根误差和综合相对误差作为评判标准。

2 分区及因子选择

2.1 聚类分析

聚类分析是研究多要素（多个变量）的客观分类方法。它的聚类原则是根据某些相似性的指标进行聚类，把对象的个体（样品）进行联合，用分裂或添加的方法进行聚类或串组，故也称串组分析。根据四川各县区的经纬度以及海拔高度，运用聚类分析的7种方法（最远邻元素、最近邻元素、质心连接、组内连接、组间连接、中位数以及快速聚类）并均采用平方欧氏距离进行筛选比对，筛选条件包括最终分得区域数目合适、分得同一区域所在位置相对聚集、区域站点数目合适并平均、分得区域地形特征统一。最终选定组间连接方法将四川地区分为4个区域，分别是区域1（南部地区）、区域2（东北部地区）、区域3（西北部地区）、区域4（中东部地区）（图2）。

图2 四川省分区结果Fig.2 The results after the division in Sichuan Province

2.2 相关性分析

将四川全区及各区多年平均夏季降水量与各因子进行相关分析，相关系数如表1所示。可以看出，四川各区多年平均夏季降水量与6种因子均相关，其中坡向因子与全区和各区相关性最差且未通过α=0.05的显著性检验。全区多年平均夏季降水量与海拔相关系性最好，为-0.461，且通过α=0.01的显著性检验。四川包含多种复杂地貌，使得夏季降水整体上受海拔影响很深，其次是坡度，如迎风坡、背风坡、陡坡、缓坡等对降水的影响占比较高。区域1与经度相关性最好，且通过α=0.05的显著性检验，其次是坡度、海拔、NDVI，但都没有通过α=0.05的显著性检验。区域2与海拔和坡度的相关性通过α=0.01的显著性检验，这是由于区域2西侧小部分区域为盆地边缘地区，地形起伏较大，地形会影响到当地的降水。区域3与NDVI的相关性最好，通过α=0.01的显著性检验，说明在川西高原地区，植被覆盖率对山区降水具有一定的影响。区域4与经度、纬度相关性最好，均通过α=0.01的显著性检验。

表1 四川全区及各区多年平均夏季降水量与各因子相关系数Tab.1 The correlation coefficients between multi-year average summer precipitation in Sichuan and various factors in the whole region and each district

综合全区和4个分区来看，坡向与全区与各分区的相关性都很低且没有通过显著性检验，所以在因子组合中去除坡向这一因子。

2.3 三因子组合的CoK

CoK只能加入3个辅助变量，从5个因子（经度、纬度、海拔、NDVI、坡度）中任意选3个因子作为一个组合，总共有10种组合。以每3个因子作为自变量，多年平均夏季降水量作为因变量，进行多元线性回归分析，结果如表2所示。对于全区而言，R2最高的三种组合方式为AlNdSl、LaAlSl和LaAlNd组合，但拟合效果均不佳，其中AlNdSl组合与多年夏季平均降水量拟合度最高（R2=0.293），R2最高的三种组合方式中均有海拔因子，说明海拔对四川夏季降水影响很大。对于区域1来说，不同组合差距悬殊，R2最高的三种组合方式为LoNdSl、lLoLaSl和LoAlSl，其中最佳拟合组合为LoNdSl，R2=0.556。区域2，10种组合得出的拟合效果差距不大，R2最高的三种组合方式为LoLaAl、LoNdSl和LoLaSl，其中最佳拟合方式是LoLaAl的组合，R2=0.397。综合来看，区域3拟合效果较好，R2最高的三种组合方式为LoAlNd、AlNdSl和LaAlNd，其中R2最高达0.625，为LoAlNd的组合。区域4总体而言拟合效果不如前三个区域的拟合效果好，10组组合中，R2最高的三种组合方式为LoLaNd、LoLaSl和LoLaAl，其中拟合效果较佳的是LoLaNd组合，R2=0.267，其次是Lo-LaSl和LoLaAl，3组组合都含有经度和纬度，说明四川中东部地区夏季降水受地理位置影响比较大。

表2 全区及各区因子组合与多年平均夏季降水量相关性R2Tab.2 The correlation R2 between the combination factors and the multi-year average summer precipitation in each area and the whole area

3 空间插值结果

图3为全区多年平均夏季降水插值结果比较。从四川全区来看，插值结果整体呈现为西北小、东南大的空间分布，即盆地区域降水多，高原地区降水少。川中东部地区出现条状的大值区，为西北—东南向的椭圆形，该区域西部、西北部和南部的高山形成喇叭状地形，使得东来的太平洋东南暖湿气流与盆地周边山地下沉的冷湿气流交汇于此处，形成著名的“华西雨屏”现象［30］。不同的空间插值方法在全区具有一致性，但在局地存在一定差异，OK、3种组合的CoK插值结果表现为川东北地区的大值区与川南北部地区的部分大值区相连，而其他插值结果明显形成了断裂并且在盆地及其东部地区出现小范围的“牛眼”现象，这种现象是IDW常出现的一种现象，这些“牛眼”能突出站点的特征，并在一定程度上提高了插值精度。OK插值以及3种组合的CoK插值差异肉眼几乎难以分辨，必须后续对插值结果进行更进一步的交叉检验来判断几种方法的优劣。

图3 全区多年平均夏季降水插值结果比较（单位：mm）（a）OK，（b）IDW，（c）RBF，（d）AlNdSl组合CoK，（e）LaAlSl组合CoK，（f）LaAlNd组合CoK，（g）LPI，（h）EBKFig.3 Comparison of interpolated results of multi-year average summer precipitation in the whole region（Unit：mm）（a）OK，（b）IDW，（c）RBF，（d）AlNdSl CoKriging interpolation，（e）LaAlSl CoKriging interpolation，（f）LaAlNd CoKriging interpolation，（g）LPI，（h）EBK

图4为区域1多年平均夏季降水插值结果比较。区域1大值区主要集中在其北部地区，且有深入其中部地区的趋势，东南地区也有小范围大值区，而东、西边缘地区值较小。IDW、RBF以及EBK插值结果仍有部分“牛眼”现象。除OK、LPI和LoLa-Sl组合CoK外，其余插值均在东南边缘地区有一大值区。LPI、LoLaSl组合CoK以及OK、EBK插值结果表现为北部大值区较为分散，其余插值则较为集中。

图4 区域1多年平均夏季降水插值结果比较（单位：mm）（a）OK，（b）IDW，（c）RBF，（d）LoNdSl组合CoK，（e）LoLaSl组合CoK，（f）LoAlSl组合CoK，（g）LPI，（h）EBKFig.4 Comparison of interpolated results of multi-year average summer precipitation in District 1（Unit：mm）（a）OK，（b）IDW，（c）RBF，（d）LoNdSl CoKriging g interpolation，（e）LoLaSl CoKriging interpolation，（f）LoAlSl CoKriging interpolation，（g）LPI，（h）EBK

区域2插值结果表现为整体从西北到东南呈小、大、小分布。高值区均呈带状，主要位于盆地北部边缘地区，低值区不同插值方法各不相同但均在南部地区有一个大范围的较低值区域。除3种组合的CoK以及OK插值外，其他插值在阿坝州的茂县、绵阳的北川地区均有明显的小值区（图5）。

图5 区域2多年平均夏季降水插值结果比较（单位：mm）（a）OK，（b）IDW，（c）RBF，（d）LoLaAl组合CoK，（e）LoNdSl组合CoK（f）LoLaSl组合CoK，（g）LPI，（h）EBKFig.5 Comparison of interpolated results of multi-year average summer precipitation in District 2（Unit：mm）（a）OK，（b）IDW，（c）RBF，（d）LoLaAl CoKriging interpolation，（e）LoNdSl CoKriging interpolation，（f）LoLaSl CoKriging interpolation，（g）LPI，（h）EBK

区域3海拔总体较高，降水主要集中在其东南部的九龙一带。RBF、EBK与AlNdSl组合CoK插值结果表现为在阿坝州马尔康和壤塘地区出现相连带状区域，而OK、IDW以及LPI插值则为较为分散的点状。LoAlNd组合CoK插值在甘孜州石渠县的值明显高于其余几种插值，并且LoAlNd和LaAlNd组合的CoK插值在大部分地区都为大值区域（图6）。

图6 区域3多年平均夏季降水插值结果比较（单位：mm）（a）OK，（b）IDW，（c）RBF，（d）LoAlNd组合CoK，（e）AlNdSl组合CoK，（f）LaAlNd组合CoK，（g）LPI，（h）EBKFig.6 Comparison of interpolated results of multi-year average summer precipitation in District 3（Unit：mm）（a）OK，（b）IDW，（c）RBF，（d）LoAlNd CoKriging interpolation，（e）AlNdSl CoKriging interpolation，（f）LaAlNd CoKriging interpolation，（g）LPI，（h）EBK

区域4各种插值方法的插值结果差异不大，高值多集中在其西北部，低值分散在东部地区。除3种组合的CoK以外，其他插值方法在其西北部插值结果呈现出块状明显、边缘清晰的特征。LPI在宜宾筠连附近出现大值区，EBK、OK、IDW、RBF在筠连的大值区不太明显，而3种组合CoK在此地没有表现出大值区域（图7）。

图7 区域4多年平均夏季降水插值结果比较（单位：mm）（a）OK，（b）IDW，（c）RBF，（d）LoLaNd组合CoK，（e）LoLaSl组合CoK，（f）LoLaAl组合CoK，（g）LPI，（h）EBKFig.7 Comparison of interpolated results of multi-year average summer precipitation in District 3（Unit：mm）（a）OK，（b）IDW，（c）RBF，（d）LoLaNd CoKriging interpolation，（e）LoLaSl CoKriging interpolation，（f）LoLaAl CoKriging interpolation，（g）LPI，（h）EBK（Unit：mm）

总体而言，同区域不同插值方法得出的结果差异不大，插值结果大体走向与增减几乎一致，要得到更精确的对照分析，必须对结果进行交叉检验。

4 交叉检验结果

表3为全区与各区多年平均夏季降水交叉检验结果。对于全区而言最好的插值方法是EBK，综合相对误差（每个区域内点的相对误差求得的平均值）只有9.11%，其次是RBF和OK。其原因是在大量样点数据中，RBF能减少误差，而OK在全区更能发挥优势。全区3种组合的CoK均没有达到预期效果，原因可能是全区的降水是多种因子复合影响的结果，对于仅有的3种因子的协同插值达不到预期的效果。但在3种组合CoK中，AlNdSl组合方式最优，说明对于全区来说，海拔、NDVI和坡度对于降水的影响较大。

区域1的插值结果为几个区中最佳，其中最佳插值方法是LoNdSl组合CoK。RBF和IDW插值平均误差虽小，但从标准均方根误差来看略差于其余插值方法，说明这两种插值方法的误差极值效应与CoK相比较大。3种组合的CoK效果均较好，说明此区域在插值方法中引入最能影响此区域的环境因子能够很大程度上提高插值结果的精确性。区域1的降水主要受经度和坡度影响，结合区域1地势来看，区域1处于山地，降水受迎风坡或背风坡的影响较大。

区域2插值结果最佳的是EBK，其次是RBF，这是由于区域2站点多而密集，有利于RBF的插值。3种组合的CoK在区域2的插值效果不好，其原因可能是多方面的：一是该区域大部分为平原地区，地势开阔，海拔平均较低，能影响该区域的因子较为单一，用3种因子组合的CoK方法反而会降低插值精度；二是此区域站点较多，对于RBF和EBK来说更为适合。

区域3插值效果最佳的是AlNdSl组合CoK，其次是RBF。区域3处于川西高原地区，海拔较高且起伏较多，另外区域3降水与NDVI相关性很高，超过了0.5且通过了α=0.01的显著性检验，所以引入AlNdSl的协同克里金在区域3的插值效果最佳。

区域4插值效果最佳的是RBF和EBK，相对误差均在10%左右，3种组合的CoK插值效果一般，但比OK插值效果佳，说明此区域的降水确实由多种因子相互影响，仅3种因子的协同插值不足以体现出该区域的降水分布。若将此区域更细致地划分，会得到更好的插值结果。

综合比较全区与4个区域所有插值方法，RBF和EBK的综合相对误差相对于其他插值方法更小，其原因是选用四川全区时站点密集，有利于RBF插值，并且EBK自带减少误差的子程序；其次为OK和Cok，在选取因子合适的区域，Cok插值可以达到更佳的插值效果，而在影响因子众多或者地形单一的区域，RBF和EBK更胜一筹。总体而言IDW插值效果最差，分区之后所有插值方法的插值精度要高于分区前的精度。

5 结论

基于四川省157个自动气象站点近10 a（2010—2019）夏季降水数据，采用聚类分析进行分区后通过相关性分析和多元回归分析筛选出各区域降水量的地理影响因子。使用协同克里金插值方法的同时，采用传统插值方法进行对比并对插值结果进行交叉检验，结论如下：

（1）影响全区降水的主要地理因子为海拔、坡度、经度和纬度，主导的影响因子组合是海拔、NDVI和坡度；分区细化之后，对南部地区降水影响的主要因子是经度、坡度和海拔，主导的影响因子组合是经度、NDVI和坡度；对东北部地区而言，降水主要因子为坡度，其次是海拔，其主导的因子组合方式是经度、纬度和海拔；西北地区的降水主要因子是NDVI和经度，其主导的因子组合方式为经度、海拔和NDVI；中东部地区的主要降水影响因子是经度和纬度，其主导的因子组合方式是经度、纬度和NDVI。

（2）全区而言，最佳的插值方式是经验贝叶斯克里金插值；对四川南部地区而言，最佳插值方法是经度、坡度和NDVI因子组合的协同克里金插值；对四川东北部地区而言，最佳插值方法是经验贝叶斯函数；对四川西北部地区而言，最佳的插值方式是海拔、NDVI和坡度为主导因子组合的协同克里金插值函数；对四川中东部地区而言，最佳的插值方式是径向基函数插值方法。

（3）分区后插值精度高于分区前插值精度，在所选区域降水影响因素数目适中时，选用协同克里金插值效果更佳。所选区域降水影响因素数目单一或众多的情况下，选用经验贝叶斯克里金插值或径向基函数的效果更佳。