黄土高原降雨量空间插值精度比较——KRIGING与 TPS法

万龙，马芹，张建军，付艳玲，张晓萍

(西北农林科技大学资源环境学院;中国科学院水利部水土保持研究所;黄土高原土壤侵蚀与旱地农业国家重点实验室:712100，陕西杨凌)

降雨是地表水资源的主要补给来源，认识降雨时空演变规律对科学利用和调控有限地表水资源尤为重要。近年来，随着GIS在多种领域的广泛应用，很多学者对于降雨量空间插值方法的优劣性做了深入研究，认为各种插值方法都有自己的优点，但也有明显的缺陷［1-3］。如泰森多边形方法(Thiessen)适用于站点密集地区，要求地形大致相同，且不考虑高程影响［4］，反距离加权法(IDW)通过权重调整空间插值等值线结构，但该方法也不考虑地形因素对降雨的影响［5］。基于统计插值技术的克里金插值方法(KRIGING)利用半方差函数来反映插值对象的空间相关性，被认为是空间数据的最佳内插法［6］。薄板光滑样条函数方法(Thin Plate smoothing Spine-TPS)是样条函数方法的一种，兼顾了插值曲面的平滑性与精确度［7］。除了插值方法外，站点数量及分布均匀程度、时间尺度变化和栅格单元大小等因素，均会引起空间信息的不确定性［8-9］。

黄土高原地区年降雨量在300～650 mm之间，东南季风、地理位置及地势构架因素的影响，空间上呈现出降雨量由东南向西北递减的趋势，时间上则表现为年际变化大，月际分配极不均匀的特点。受地面强烈切割和起伏的影响，稀疏覆被地表加强了夏季热力对流作用，降雨形式多暴雨和冰雹［10］。在气象站点观测值基础上，如何准确估算站点外区域的气象信息，是科技和管理人员实践中经常面对的问题。KRIGING插值方法是国内外文献中通过站点数据获取空间信息常用的方法［11］。为气象数据曲面拟合而编写的专用软件ANUSPLIN是基于样条插值理论，允许引进多元协变量的线性子模型，近年来在国际上得到了广泛应用［12］。针对黄土高原降雨时空变异性大的特征，哪种方法更能准确反应气候要素的空间信息，少有文献探讨。笔者拟利用黄土高原的降雨数据对这2种方法在多年平均、年和月尺度上进行插值精度比较，以期为土壤侵蚀模型的建立提供参考。

1 研究区概况

内蒙古河口镇—陕西龙门区间(河龙区间)是黄河中游处于山西—陕西交界的区域，面积约11.3万km2，以黄土丘陵沟壑地貌为主，海拔 450～2 079 m，切割深度在100～500 m之间，≥25°陡坡面积比例10% ～40%，地表割裂度30% ～70%。东部以吕梁山为界，西南邻北洛河，西北为毛乌素沙地。属于干旱半干旱气候，年内降雨集中于5—9月，暴雨强度大，因引起剧烈侵蚀，每年输入黄河的泥沙量超过三门峡以上总输沙量的70%［13］。河龙区间植被以草原为主，北部为典型草原和荒漠草原，南部为灌丛草原、草甸草原和森林草原。区间内黄绵土广泛分布，由南向北逐渐粗化。截至1996年底，河龙区间总人口约800万，其中农业人口约697.2万，土地利用不合理，水土流失严重。

2 数据资料及研究方法

2.1 数据资料

黄河中游河龙区间及毗邻区共包括50个气象站点，其中27个站点分布于河龙区间(图1)，站点分布经纬度范围为 E106.5°～113.5°，N35°～40.5°。50个气象站点1981—2000年的月降雨数据从国家气象站获得。各站点月降雨量数据累加，计算多年平均降雨量值。在ArcGIS支持下使用KRIGING方法和采用ANUSPLIN专用软件的TPS方法，对多年平均、年及月降雨量进行插值，利用河龙区间内27个站点数据进行交叉验证和对比，比较2种方法插值精度的优劣。插值设置投影参数为Albers等积投影，第1条纬线为N 25°，第2条纬线为N 47°，中央经线为 E 110°，插值栅格为500 m ×500 m。研究区气象站点平均分布密度为4200 km2/站。

2.2 KRIGING插值法

KRIGING插值法是由法国地理学家Georges Martheron和南非采矿工程师D.G.Krige创立的地质统计学的最佳内插方法［14］，近年来被广泛应用于地学研究中。区域化变量和变异函数是克里金插值的基石。

图1 河龙区间气象站点分布图Fig.1 Location of the meteorological stations in the He-Long Section

当空间点x在一维x轴上变化时，把区域化变量在x与x+h处的Z(x)与Z(x+h)差的方差之半定义为区域化变量Z(x)在x轴方向上的变异函数，即半方差函数，记为γ(h)，在满足二阶平稳假设条件下，有

式中:N(h)是研究区内距离为h的点对数，Z(x)和Z(x+h)分别为x及x+h点的变量值。

图2描述了半方差函数与站点点对距离h之间的关系，随着距离h的增加，观测点的半方差增大，相关性也逐渐降低。当h达到一定程度，即h=a(变程)时，γ(h)逐渐趋近于定值。理论上γ(h)在h=0时γ(h)=0，但是，有时会在原点附近出现不连续的现象。这种现象称为块金效应(Nugget Effect)，这时的γ(h)就叫做块金效应值(C0)，而(C+C0)为基台值。当C0/(C+C0)比较小时，各采样点的加权值变化比较大，数据的空间自相关性较好。随着C0/(C+C0)的增大，空间数据相关性变差，各个采样点的加权值变化减小，数据的空间自相关程度降低。块金效应值等于基台值时，即纯块金效应，说明采样数据没有包含空间相关性的任何量化信息，这种情况的估计结果就完全是采样点的算术平均值，各个采样点的加权值都一样大［15］。

图2 半方差函数变异曲线Fig.2 Semi-variogram variation curve

由于计算得到的实验半方差函数是离散的且变化也往往不规则，在克里金估计中难以直接应用于指导克里金插值，所以要用理论模型去拟合。半方差理论模型常采用的模型有球面模型，模型参数通过实验数据用最小二乘法拟合得到。文中变程设置为各点对的平均距离50 967 m，用离估计点最近的12个站点的数据来进行插值。

拟合半方差重要的用途是确定局部内插需要的权重因子，其值按采样点数据的半方差图的统计分析原理计算，即:

式中:Zi(xi)为已知第i个站点的观测值;λi为第i个站点的权重;^Z(x0)为待估点x0的预测值。λi权重的选择应使 ^Z(x0)是无偏估计。

2.3 TPS方法

薄板样条函数(TPS)插值法，具有连续、光滑的数学特性，是自然样条函数在多维空间的推广。它除普通的样条自变量外允许引入线性协变量子模型，如温度，海拔等。由于黄土高原降雨与海拔无明显的相关关系;故仅采用以经度和纬度为双变量薄板光滑样条函数进行插值，没有引入海拔值作为协变量因子［16］，且文中样条次数采用通常的二次样条函数。

TPS 法插值原理［17］如下:

式中:zi为位于空间i点的因变量;xi为d维样条独立变量矢量;f(xi)为要估算的关于xi的未知光滑函数;yi为p维独立协变量矢量;b为yi的p维系数;ei为具有期望值为0和方差为wiσ2的自变量随机误差，wi为作为权重的已知相对误差方差;σ2为所有数据点上为一常数的误差方差 ，但通常未知;N为观测值个数。

函数f和系数b通过下述最小二乘法估计确定:

式中:Jm(f)为函数f(xi)的粗糙度测度函数，定义为函数f的m阶偏导(在ANUSPLIN中称为样条次数，也叫糙度次数);ρ为正的光滑参数，在数据保真度与曲面的粗糙度之间起平衡作用，本文中由广义交叉验证GCV(generalized cross validation)的最小化来确定［18］。ANUSPLIN软件在日志文件(Log file)中提供了用于判别误差来源和插值质量的统计参数，可用于确定拟合过程中的最优光滑参数。

2.4 误差分析方法

采用交叉验证的方法对数据进行比较，依次移去一个点，用所有其他站点插值，得到被移去站点的预测值，该预测值与实测值的误差为该站点的交叉验证误差。常用的评价插值精度的方法有标准均方根误差(Root Mean Squared Interpolation Error，E)，可以反映估计值与实测值误差的绝对量［19］。为了消除对异常值的敏感程度，同时利用一致性指标(the index of agreement，A)来评价KRIGING方法和TPS方法插值精度的优劣［20］。公式如下:

式中:Oi为第i个站点的实际观测值;Pi为第i个站点的交叉验证预测值;为n个站点的实际观测值的平均值;n为用于参与验证的站点数。A值在0～1之间，1表示预测值与实测值非常吻合，0表示不吻合。

2.5 全域型Moran’s I指数

计算空间自相关的方法有许多种，最为知名也最为常用的是Moran’s I方法。用全域型Moran’s I指数来表示研究区要素空间自相关程度。全域型Moran’s I指数是基于统计学相关系数的变异数和共变量关系推算得来。假设数组 xi全域型的Moran’s I的公式如下:

式中:I为moran’s I指数值，Wij为研究范围内每一个空间单元 i与空间单元 j(i、j=1，2，…，n)的空间相邻权重矩阵，以1表示i与j相邻，以0表示i与j不相邻;n为数组xi的元素个数。

I值标准化Z(I)的公式为

式中:E(I)为I的期望值;Var(I)为I的方差。对I值进行显著性确定时，在5%显著水平下，Z(I)大于1.96时，表示研究范围内某现象的分布有显著的空间自相关性。

3 结果与分析

3.1 多年平均降雨量插值结果比较

受东南季风的影响，河龙区间多年(1981—2000年)平均降雨量由东南向西北递减，其空间变异情况为:平均值424.0 mm，极值比1.66，标准差52.9 mm，离差系数0.125，Moran’s I为0.57，统计值 Z=8.0，空间上具有极显著的自相关性。KRIGING与TPS方法的插值结果均能很好地反映多年平均降雨量的空间变化趋势，具有大致相同的空间变化梯度，TPS方法插值法降雨量的等值线更加平滑和柔和(图3)。

对河龙区间内27个站点进行KRIGING与TPS方法交叉验证，标准均方根误差E分别为21.2和21.7 mm，一致性指标A分别为0.775和0.765。这2种方法对河龙区间多年平均降雨量的插值精度差别很小，E差值为多年平均降雨量的0.12%，A相差0.01。KRIGING方法(图3(a))要略微优于TPS方法(图3(b))。从数据点对的半方差云图(图3(c))可知，数据的空间相关性较好，无块金效应，插值精度较高。

图4显示了KRIGING与TPS插值方法对27个站点多年平均降雨量的交叉验证误差，可见，这2种方法交叉验证误差具有极显著的相关性，近85%的误差集中在-30～30 mm之间。交叉验证最大误差站点均出现在兴县，且均大于40 mm。插值误差较大的区域集中在山西河曲、五寨、兴县、离石、隰县等站点分布较稀疏的地区。

3.2 年尺度插值结果比较

图3 河龙区间1981—2000年多年平均降雨量插值结果Fig.3 Interpolation maps and semi-variogram cloud of the average annual precipitation of 1981—2000

图4 KRIGING与TPS方法对河龙区间多年平均降雨量插值交叉验证误差对比Fig.4 Comparison of the cross-validation error for average annual precipitation interpolation between KRIGING and TPS methods at stations within He-Long Section

分别采用KRIGING与TPS方法对1981—2000年降雨量进行插值，对河龙区间内27个站点进行交叉验证，计算平均E及A值，结果见图5。结果表明，KRIGING和TPS方法对河龙区间1981—2000年降雨量插值结果的标准均方根误差E(图5(a))近95%的年份均介于40～70 mm之间，其均值分别为56.2和57.0 mm，一致性指标A值(图5(b))均主要介于0.4～0.8之间，平均值分别为0.582和0.576。这2种方法的插值精度非常接近，平均E差值为多年均值的0.19%，平均A值相差0.006，KRIGING方法仍略微优于TPS方法。各年份2种方法A值相差基本在±0.05范围内，与多年平均降雨量插值比较，2种方法的年插值精度平均A值均降低约0.19(19%)。

从图5中选取E值较小、A值较大的1987年和E值较大、A值较小的1996年进行具体比较分析。1987年河龙区间暴雨次数较少，降雨强度较低，各站点降雨量差别较小，均在510 mm以下;而1996年一些站点降雨次数及暴雨强度较大，降雨量在各站点差别大，个别站点如兴县、安塞、清涧等地区降雨量甚至超过600 mm。研究区这2年的降雨量空间变异情况统计结果见表1。可知，尽管1987年各站点年降雨极值比、离差系数较1996年大，但1987年各站点年降雨量空间自相关性要较1996年好得多，也就是说，1996年各站点年降雨量空间变异性较1987年大得多。

图5 河龙区间1981—2000年年均降雨量KRIGING和TPS插值方法的交叉验证对比Fig.5 Cross-validation for annual precipitation using KRIGING and TPS methods

表1 河龙区间1987和1996年降雨量空间变异统计Tab.1 Spatial variation of annual precipitation in 1987 and 1996

对于不同空间变异的具体年份，KRIGING和TPS方法插值结果差异较大(图6)。1987年(图6(a))降雨量空间自相关性较好，KRIGING和TPS方法插值精度均较高，交叉验证平均 A值分别为0.711和0.671，而1996年(图6(b))降雨量的空间自相关性较差，2种方法插值的精度均较低，交叉验证平均A值分别为0.407和0.471。与KRIGING方法相比，TPS方法的平滑作用可以改变局部精度而影响到全局效果。如1987年东南部插值的平滑作用，使全局交叉验证A值提高了0.04，E值减少了6.9 mm;但1996年局部平滑使全局A值降低了0.06，E值增加了1.6 mm。

图6 河龙区间1987、1996年降雨量KTIGING和TPS方法插值效果图和半方差云图比较Fig.6 Interpolation maps of annual precipitation in 1987 and 1996 in the He-Long Section using KRIGING and TPS methods

3.3 月尺度插值结果比较

河龙区间降雨量月际分配极不均匀，空间变异非常大。研究时段内，春季4月份约占区间多年平均降雨量的3% ～8%，平均降雨量10～35 mm，从南向北逐渐减少且梯度平缓。夏季7月份约占多年平均降雨量的15% ～30%，平均降雨量60～120 mm，从南向北逐渐减少，有多个强降雨中心。采用KRIGING和TPS方法对降雨较平稳的4月份，以及空间变异大、暴雨强度高的7月份的降雨数据进行27个气象站点插值精度的交叉验证，结果见图7。

可知，对于4月份的插值结果(图7(a)、(b))，KRIGING和TPS方法交叉验证E值范围均在2～18 mm之间，平均为6.2和6.6 mm。A值介于0.40～0.90之间，平均为0.629和0.631。对于7月份的插值结果(图7(c)、(d))，KRIGING和TPS方法 E值范围均在20～60 mm之间，平均为34.3和35.0 mm，A值介于0.10～0.70之间，平均为0.419和0.421。结果表明:无论降雨过程平稳的4月还是空间变异性较大的7月，KRIGING和TPS方法插值精度平均结果差别较小;这2种方法对月雨量插值交叉验证平均一致性指标A值均相差0.002，但在具体年份有较大差别，2种方法A值绝对值差异基本在±0.05范围内波动;在降雨变异较大的7月，2种方法插值的平均精度均低于年值插值结果约0.16(16%)，更低于多年平均插值结果约0.35(35%)，但对于降雨过程较平稳的4月，2种方法的插值精度均高于年值插值精度约0.05(5%)，低于多年平均插值结果约0.14(14%)。

由于7月降雨量插值平均精度较差，各年份间差异较大，故对插值精度较高的1992年7月和精度较低的1989、1996年7月插值表面、半方差云图及站点交叉验证误差进行深入对比，其空间变异统计结果见表2，插值结果比较见图8和图9。

图7 用KRIGING与TPS方法对河龙区间1981—2000年4和7月降雨量插值交叉验证误差对比Fig.7 Comparison of cross-validation in April and July precipitation using KRIGING and TPS methods

表2 河龙区间1989、1992、1996年7月降雨量空间变异统计Tab.2 Spatial variation of July precipitation in 1989，1992 and 1996

1992年(图8(b))7月，要素在空间上呈正自相关，并且自相关性很显著。2种方法插值表面光滑，TPS方法更为平滑，半方差随着站点距离的增大而增大，TPS方法平滑参数大小适中，误差判断GCV较小(2.17)。由于有显著的空间自相关性，所以，2种插值方法的精度较高，2种方法交叉验证预测值与实测值一致性较好(图9)，交叉验证有效系数 ＞0.7。

而1989(图8(a))和1996年(图8(c))7月，要素空间自相关性很差，2种方法插值降雨量空间分布不一致，KRIGING方法插值表面不光滑。1989年，半方差甚至随站点距离增大基本不变，导致KRIGING方法插值产生纯块金效应，交叉验证预测值仅为相邻一些站点的平均计算结果，与实测值没有相关关系(图9)，站点交叉验证A＜0.2。TPS方法对1989和1996年插值平滑参数过大和过小，数据存在短相关，找不到合适的平滑表面，GCV较大，分别为3.01和5.50，交叉验证A＜0.2。

总之，在1989和1996年的7月，由于空间上较低的相关性或呈现负的相关性，2种方法插值均不能正确反映降雨量的实际空间分布状况，尤其在出现暴雨的站点，交叉验证预测值与实际值可相差至100 mm以上，交叉验证的预测值与实际值误差较大，仅为各站点平均结果，基本没有相关关系，2种方法插值结果均不能客观地表现该时期的空间信息。对于这样的空间分布状况，引入协变量可能会改善其精度。

4 结论

图8 河龙区间1989、1992和1996年7月降雨量KTIGING和TPS方法插值图和半方差云图比较Fig.8 Interpolation maps of July precipitation of 1989，1992 and 1996 in the He-Long Section using KRIGING and TPS methods

图9 河龙区间1989、1992和1996年7月降雨量实测值与交叉验证预测值相关关系Fig.9 Correlation between the predicted and actual values of July precipitation in 1989，1992 and 1996 at the stations within the He-Long Section

1)无论多年平均、还是年和月尺度，KRIGING和TPS方法插值结果都能正确反映河龙区间降雨量的空间分布趋势，不同时间尺度上，2种方法的交叉验证区域时段交叉验证平均一致性指标A值相差均在±0.01范围内。不同时间尺度2种方法的面平均插值精度均没有显著性和实质性差异。随着尺度变小，降雨量空间分布变异性增大，2种方法插值精度均显著降低，不同年份间一致性指标A值绝对值差别基本在±0.05范围内波动。仅从数值上可以看出，多年平均和年尺度，KRIGING方法要略优于TPS法。

2)KRIGING和TPS方法在不同时间尺度插值精度均与全域要素自相关程度显著相关。2种方法多年平均降雨量的插值精度均优于4月，其次是年值、7月。与多年平均降雨量插值平均一致性指标A值比较，4月插值精度均降低约14%，年值均降低约19%，而7月均降低约35%。

3)TPS方法插值降雨量的等值线更加平滑和柔和，但河龙区间降雨空间变异性较大，与KRIGING方法相比，TPS方法的局部平滑作用可以改善或削弱全局插值精度。

［1］蔡福，于慧波，矫玲玲，等.降水要素空间插值精度的比较［J］.资源科学，2006，28(6):73-79

［2］李新，程国栋，卢玲.空间内插方法比较［J］.地球科学进展，2000，15(3):260-265

［3］林忠辉，莫兴国，李红轩，等.中国陆地区域气象要素的空间插值［J］.地理学报，2002，57(1):47-56

［4］王淑英，陈守煜.加权平均的权重优选算法及其应用［J］.水利学报，2003(12):1-7

［5］何红艳，郭志华，肖文发.降水空间插值技术的研究进展［J］.生态学杂志，2005，24(10):1187-1191

［6］李军龙，张剑，张丛，等.气象要素空间插值方法的比较分析［J］.草业科学，2006，23(8):6-11

［7］刘志红，Tim R McVicar，Li Lingtao，等.基于 5 变量局部薄盘光滑样条函数的蒸发空间插值［J］.中国水土保持科学，2006，4(6):23-30

［8］许家琦，舒红.降水数据空间插值的时间尺度效应［J］.测绘信息与工程，2009，34(3):29-30

［9］朱会义，贾绍凤.降雨信息空间插值的不确定性分析［J］.地理科学进展，2004，23(2):34-42

［10］张汉雄.黄土高原的暴雨特性及其分布规律［J］.地理学报，1983，38(4):416-425

［11］Bargaoui Z K，Chebbi A.Comparison of two kriging interpolation methods applied to spatiotemporal rainfall［J］.Journal of Hydrology，2009，365(1-2):56-73

［12］Tiat A，Henderson R，Turner R，et al.Thin plate smoothing spline interpolation of daily rainfall for New Zealand using a climatological rainfall surface［J］.International Journal of Climatology，2006，26(14):2097-2115

［13］冉大川，柳林旺，赵力仪，等.黄河中游河口镇至龙门区间水土保持与水沙变化［M］.郑州:黄河水利出版社，2000

［14］王广德，过常龄.“Krige”空间内插技术在地理学中的应用［J］.地理学报，1987，42(4):366-375

［15］刘永社，印兴耀，贺维胜.空间相关分析因素对储层建模中克里金估计结果的影响［J］.石油大学学报，2004，28(2):24-27

［16］刘志红，Tim R McVicar，Li Lingtao，等.基于 ANUSPLIN的时间序列气象要素空间插值［J］.西北农林科技大学学报，2008，36(10):227-234

［17］Hutchinson M F.ANUSPLIN Version 4.3 User Guide［M］.Canberra:The Australia National University，Center for Resource and Environment Studies，2004

［18］Hancock P A，Hutchinson M F.Spatial interpolation of large climate data sets using bivariate thin plate smoothing splines［J］.Environmental Modelling ＆ Software，2006，21(12):1684-1694

［19］Newlands N K，Davidson A，Howard A，et al.Validation and inter-comparison of three methodologies for interpolating daily precipitation and temperature across Canada［J］.Environmetrics，2010，22(2):205-223

［20］张晓萍，张橹，穆兴民，等.黄河中游河口一龙门区间多年平均流域水平衡特征［J］.地理学报，2007，62(7):753-763