降雨量空间插值方法比较研究

李巍　范文义　毛学刚　王兰霞

摘要对反距离加权法、样条函数法、普通克里格法及协同克里格法等几种常用的降雨空间插值方法的优缺点进行了分析和比较，考虑到高程对降雨量影响较大，在协同克里格法的基础上将高程作为第2类影响因素引入降雨量的空间插值方法中，并提出了引入高程信息的协同克里格法。将4种方法（反距离加权法、样条函数法、普通克里格法及协同克里格法）用于大兴安岭地区降雨量的插值计算，结果表明，考虑高程信息的协同克里格法的插值效果明显优于其他3种方法。

关键词降雨量；空间插值方法；协同克里格法；高程信息

中图分类号S161.6文献标识码A文章编号0517-6611（2014）12-03667-03

基金项目“十二五”国家科技支撑项目（2011BAD08B01）；黑龙江省青年基金（QC2012C102）；黑龙江省普通高校重点实验室项目（KJKF1203）；黑龙江省教育厅科学技术研究项目（12511490，12531599）。

作者简介李巍（1982- ），女，黑龙江哈尔滨人，讲师，在读博士，从事3S技术研究及3S在土壤侵蚀、土地规划方面研究。*通讯作者，教授，博士后，博士生导师，从事3S技术研究及应用。

气象资料是进行地学研究的一个基本而重要的参数，气象资料可以反映出区域的森林植被分布、土壤侵蚀分布、水文生态的规律[1]。由于经济条件和技术手段的制约，大部分区域内设置的气象站点观测密度均不高，这就使得研究者们无法在研究区内得到连续有序的气象空间数据分布。因此利用现有的气象观测站的数据通过空间插值对观测数据进行补充就尤为重要。笔者利用大兴安岭地区2000～2010年降雨量资料，采用4种常用的降雨空间插值方法对降雨量插值进行比较分析。

1资料与方法

选取大兴安岭地区2000～2010年的多年平均降雨量资料，采用反距离加权法、样条函数法、普通克里格法及协同克里格法4种方法对大兴安岭地区降雨量进行插值分析，并对插值结果进行检验，找出最优插值方法。

1.1插值方法

1.1.1反距离加权法。现阶段，对地学研究最常用的插值方法就是反距離加权插值方法（IDW）[2]，它是从地理学第一定律的角度推出，也叫距离倒数乘方法。1972年，美国国家气象局提出了反距离加权法。反距离加权法就是根据数据点间的空间距离远近加权插字段进行的插值方法。距离中心越近的点，其估算值越受影响，中心点的影响随着离它的距离越远而减小。反距离加权插值的计算公式是：z（x0）=ni=1z（xi）（di0）/ni=11（di0）p，式中，z（x0）是插值点的预估值，z（xi）（i=1，2，3，…，n）是实测样本值，n表示参与插值的实测样本数；di0表示插值点和第i0个站点之间的距离；p是距离的幂，插值的结果被p影响，p越大，内插后的效果越平滑，反之，p越小，内插后的效果越尖锐，一般根据最小平均绝对误差的大小来确定p。在日常使用中，经常出现反距离平方加权法，就是令p=2的反距离加权法。如果p=3，就进行了反向距离立方插值。在ArcGIS中，p的默认值为2。此次降雨量插值时p的值设为2，站点搜索的范围设为周边6个，即n=6。

1.1.2样条函数法。样条函数（SPLINE）通过一个使表面整体曲率最小的数学函数来估计单元值，所得表面较为光滑。样条函数法主要是针对一些特征点，对估计方差进行控制，适用多项式拟合的方法产生平滑的插值曲线。在降雨量空间插值中，样条函数就是用一个最小曲率面来充分逼近各降雨量观测点，进而推算出整个研究区的降雨量分布。样条函数的计算公式为：Z=ni=1λiR（γi）+T（x，y），式中，Z为降雨量的预测值；n为参与此次插值的观测站点数；λi是线性方程组求解确定的系数；γi是预估测点到第i点的距离，R（γi）和T（x，y）的表达式分别为：R（γi）=γ24[ln[γ2π]+c-1]+τ2[k0[γτ]+c]+ln[r2π]2π、T（x，y）=a1+a2x+a3y，式中，τ2为权重系数，γ为一直点与采样点之间的距离，k0是修改后的贝塞尔函数，c为常数，a为线性方程的系数。此次插值所用的样条函数为规则样条函数，指数选择0.1，点点搜索范围为临近的6个观测站点。

1.1.3普通克里格插值法。普通克里金插值法（OK）是区域化变量的线性估计，是由法国的地理学家Matheron和南非矿山工程师Krige提出[3]，普通克里格法最初用于矿山勘探。它假设数据变化呈正态分布，通过对数据的空间分析获取权重值，插值的整个过程相当于在未知区域化变量Z的期望值时对样点进行加权滑动求取平均值的过程。普通克里格方法的公式为：Z（x0）=ni=1λiZ（xi）、ni=1λi=1，式中，Z（x0）为待插值点的估计值，Z（xi）为第i个样本点的实测值，n为参与计算的观测站样本个数，λi为第i个样本点的权重系数。权重λi的选择必须保证Z（x0）能进行无偏估计，且估计的方差小于其他线性组合生成的方差。克里格插值方法的关键问题是选择合适的变异函数，目前经常使用的变异函数包括高斯、考克斯（幂）、指数、球面及线性模型等。此次研究选择考克斯变换模型作为基本变异函数模型，采用离预估点最临近的6个观测站点的数据进行插值。

1.1.4协同克里格插值法。Dirks等发现当站网密度较高时，普通克里金方法的插值效果与其他常用方法相比并无多大优势[4]。Borgam等也曾得到过同样的结论[5]。但随着高程的增加，降雨量有增加的趋势。Hevesi等研究了年平均降雨量与高程的相关性，并得到了其相关系数达0.175的结果[6]。因此，在此采用协同克里格方法（OCK），并将高程作为第2影响因素引入到对降雨量的空间插值中来[7-8]。协同克里格法的公式为：Z（x0）=ni=1λiz（xi）+λ[y（x0）-my+mz]，式中，Z（x0）为x0点的预测值；z（xi）是第i站点的观测值；y（x0）是x0点的高程；n为降雨量观测站点的个数；my和mz为海拔高程和降雨量的全局平局值；λi和λ为协同克里格插值的权重系数。尝试使用集中变异函数，结果显示指数模型的插值效果最好，站点搜索范围仍然为临近6个站点。

1.2检验方法以上4种插值方法利用交叉检验的方法对插值的结果进行对比分析。交叉检验方法是首先假设一个站点的观测值是未知的，用这个站点周围的其他站点的观测值来进行插值获取估计值，再假设另一个站点的观测值未知，用周围其他站点的观测值进行插值获取估计值，以此轮换将17个站点分别作为未知站点，获取17个站点的估计值，然后计算这17个站点的估计值和观测值之间的误差值，根据误差值的大小来判断4种插值方法效果的好坏。

利用平均相對误差指标可以判断出估计值和实测值之间的误差大小，利用平均绝对误差指标可以估算出插值获取的估计值的误差范围，因此，在此利用MRE（平均相对误差）、MAE（平均绝对误差）这2个指标作为误差的检验指标。MRE、MAE的表达式分别是：MRE=1nni=1Zai-ZeiZai、MAE=1nni=1|Zai-Zei|，式中，Zai是第i个站点的观测值；Zei是第i个站点的插值估计值；n为作为交叉检验观测站的数量。

2结果与分析

2.1多年平均降雨量插值为了更加直观地看到各种插值方法的效果，图1给出了大兴安岭地区不同方法得到的2000～2010年的多年平均降雨量插值结果。从图1可知，各种方法得到的插值结果明显不同，总体上大兴安岭地区的年降雨量在337～475 mm，强降水中心在呼玛县南部。

2.2不同插值模型的插值精度分析

2.2.1年尺度插值结果比较。从4种方法对多年年均降雨数据插值结果（图2）来看，在所有的插值方法中，优劣排序为OCK

2.2.2月尺度插值结果比较。利用2007年全年12个月的降雨量数据进行月尺度插值结果对比，从4种方法对2007年的月降雨数据插值结果（图3）来看，6～9月的插值结果的不确定性远远小于11月～次年3月。

2.3降雨量空间插值精度验证分析利用4种方法对多年

3结论

（1）从4种方法对多年年均降雨数据插值结果来看，在所有的插值方法中，优劣排序为OCK

（2）从4种方法对2007年的月降雨数据插值结果来看，6～9月的插值结果的不确定性远远小于11月～次年3月。

（3）基于普通协克里格方法对大兴安岭地区4个及其周边14个站点进行多年平均降雨量插值，得到研究区多年平均降雨空间分布特征，大兴安岭地区的年降雨量在358～462 mm，强降水中心在呼玛县南部，其最大值为462 mm，降水的低值区漠河县西部，年降水量为358 mm。

参考文献

[1] MOHAMED A S.Reliabilty estimation of rainfallrunoff models[D].New York：State University of New York，1999.

[2] LAMN.Spatial interpolation methods：a review[J].The Amercian Cartographer，1983，10（2）：129-149.

[3] 侯景儒，黄竞先.地质统计学的理论与方法[M].北京：地质出版社，1990：69-78.

[4] DIRKS K N，HAY J E，STOW C D，et al.Highresolution studies of rainfall on Norfolk Island.Part：interpolation of rainfall data[J].J Hydrol，1998，208（3/4）：187-193.

[5] BORGAM，VIZZACCARO A.On the interpolation of hydrologic variables：formal equivalence of multiquadratic surface fitting and Kriging[J].J Hydrol，1997，195：160-171.

[6] HEVESI J A，FLINT A L，ISTO J D.Precipitation estimation in mountainous terrain using multivariate geostatistics.part I：structuralanalysis[J].J Appl Meteor，1992，31：661-676.

[7] 王家华，高海余，周叶.克里金地质绘图技术[M].北京：石油工业出版社，1999：157-179.

[8] 朱会义，贾绍凤.降雨信息空间插值的不确定性分析[J].地理科学进展，2004，23（2）：34-42.