基于3D TPS空间插值的中国降水量空间分布响

杜国明，曾志芳

(1.中山大学地理科学与规划学院 广东省城市化与地理环境空间模拟重点实验室，广东 广州 510275;2.中山大学图书馆，广东 广州 510275)

草地植被与降水关系紧密，降水的减少是造成草地退化的重要因素[1]，对二者关系的研究一直是一个活跃的研究领域[2-3]。尤其我国是一个干旱、缺水严重的国家，人均水资源只有2 200 m3，仅为世界平均水平的四分之一[4]，是全球13个人均水资源最贫乏的国家之一，且分布极不均衡：东南多，西北少。研究降水量的空间分布对草地资源分布、区域水文和水资源分析以及区域水资源管理、旱涝灾害管理、生态环境治理都具有重要意义[5]。受观测台站数量及分布的限制，空间插值是十分必要的。常用插值方法有普通克里金、样条法、反距离加权插值法等[6-7]，但在常用插值方法中通常忽略了高程值的影响，而草地植被、降水量与高程等地形因子有着密切的相关性，这样造成较大的插值误差，尤其对我国地形复杂多样、山区面积广大的区域，在研究降水量等气候要素的空间分布时更应该考虑高程等地形因子的影响。为了弥补这一不足，本研究利用三维薄板样条函数(Three-Dimensional Thin Plate Splines，简称3D TPS)进行中国降水量空间数据的插值研究。

薄板样条函数(TPS)[8-10]作为一种空间插值方法，寻找一个通过所有控制点的弯曲最小(由一个能量函数定义)的光滑曲面，就像一个薄铁板，通过给定的几个“样条”，使铁板表面尽量光滑。TPS具有光滑、连续、弹性好的特点，且不要求测量点呈规格网排列，从而广泛应用于各种曲面的拟合以及复杂曲面的数学表示[11]。3D TPS除了考虑位置坐标(x,y)，还充分考虑高程值z，z不是协变量而是同位置坐标(x,y)一样，以自变量形式存在，这样，3D坐标(x,y,z)共同参与插值过程，这对于地形复杂区域的空间插值研究具有重要意义。

1 三维薄板样条函数插值原理

1.1数据 数据包括中国国界图、中国数字高程模型(DEM)图和中国气象站分布图。其中，中国国界图为矢量图，可分为面状图层和线状图层两种。中国DEM图是在美国地质调查局(USGS)的地球资源观测卫星数据中心生产的全球DEM基础上，结合中国国界图截取生成的DEM，它是GRID格式数据，本研究采用的分辨率为2.5分(即0.041 67度)，大约相当于5 km。中国气象站分布图以矢量点状图shape格式存储，在图形文件中包含每个站点的位置(x,y)信息，在属性数据库中存放年均降水量等属性信息，内含1951―2007年年均降水量值，共有454个气象站。在ArcMap中将中国DEM图和中国气象站点分布图叠加，结果如图1所示。

图1 中国数字高程及气象站空间分布图Fig.1 Digital elevation model and spatial distribution of meteorological stations in China

1.23D TPS原理介绍 假定有n个数据点的观测值r(xi,yi,zi)在三维坐标点(xi,yi,zi)处，则观测值可用如下公式表示[8-13]：

r(xi,yi,zi)=f(xi,yi,zi)+ε(xi,yi,zi),i=1,2,…,n.

(1)

式中，f(xi,yi,zi)是观测值r(xi,yi,zi)对应的估计函数，ε(xi,yi,zi)是估计误差，它是零均值、非相关随机误差；(xi,yi,zi)是三维欧氏空间的坐标。

若要估计函数f最大限度地接近观测值，可以通过最小化公式ξi得出，公式ξi如下：

λJm(f)

(2)

式中，平滑参数λ>0，Jm(f)是三维空间关于估计函数f的平滑项，被定义为：

(3)

f(xi,yi,zi)可通过以下公式得出：

(4)

式中，m是维度，在3D TPS中的m=4。φj是多项式，其中φ1=1,φ2=x,φ3=Y,φ4=z。参数aj和bi可由下式求出：

(5)

式中，a=(a1,a2,…,am)，b=(b1,b2,…,bn)T，v=(v1,v2,…,vn)。E是n×n矩阵，Eij=U(dij)。T是m×n矩阵，Tij=φ(tj)。U是径向基函数，dj是3D坐标点(x,y,z)和(xi,yi,zi)间的欧氏距离，公式表示为：

(6)

dij是(xi,yi,zi)和(xj,yj,zj)间的欧氏距离，可通过下式得出：

(7)

1.3方法与步骤 利用ArcGIS中的ArcMap完成数据预处理及成图，在Visual Basic.Net平台上利用ArcGIS Engine进行二次开发完成数据提取、处理和运算等。

1.3.1数据预处理 根据全球DEM，利用ArcGIS，通过中国国界图掩膜，得到中国DEM。将中国气象站点以shape格式点状图层形式存储，(x,y)坐标存放在图形文件中，年均降水量等属性存放在属性数据库中。

1.3.2提取气象站点的3D坐标(x,y,z)及年均降水量 在Visual Basic.Net平台上利用ArcGIS Engine的二次开发技术，获取2D坐标(x,y)和年均降水量P。高程值z可通过提取属性数据库的属性信息得到，也可通过DEM获取。两种方法之间有一定误差，由于在插值时，用DEM提取的高程值代替气象站点的实测海拔可以提高插值精度[14]，所以本研究采用DEM获取高程值的方法。

1.3.3格网化数据 先均匀生成等间距的格网，再对每个格网点进行空间插值。另外，还可采用中国DEM作为模板，即采用2.5′ × 2.5′ 作为格网单元，且与中国DEM的行列数、起始点坐标相同，每个格网的大小一样。这样，插值时共有1 134×1 479=1 677 186个格网点。

1.3.43D TPS插值 根据已知气象站点的3D坐标(x,y,z)和年均降水量值P，推求未知点——无观测值的每个格网点处的年均降水量值，由于格网点间是均匀分布的，利用ArcGIS Engine的二次开发技术不难求出每个格网点的3D坐标(x′,y′,z′)，再根据公式(2)可以求出每个格网点处的年均降水量P′。

1.3.5生成中国年降水量空间分布图 可用中国DEM作为模板，用年均降水量值P′，替换中国DEM中每个格网点处的高程值，就可以得到中国降水量空间分布图。

1.3.6与其他方法的对比分析 与二维薄板样条法(简称2D TPS)、反距离加权法(Inverse distance weighted，简称IDW)、样条法(Spline)、普通克里金、泛克里金进行比较分析[15-16]。

本研究采用标准误差(SD)作为评价指标，即插值平均误差平方和的算术平方根。设n个估计值与真实值的误差为ε1、ε2、……、εn，则标准误差公式为：

(8)

标准误差的值越小，插值精度越高；反之越低。

2 结果与分析

2.1试验结果 图2是利用3D TPS插值的结果，由于插值过程中充分考虑高程值对降水量的影响，插值效果较好。然而，由于本研究区域较广，降水量变化很大，从东南到西北，降水量总体呈递减趋势。东南沿海区域降水量多且地势较低，西北降水少。对于某些局部，如西藏南部降水量明显增多区域，图2中清晰可鉴；而在图3和图4中不很明显。另外，图2右下角的南海诸岛，由于缺乏气象站点，且岛礁相对大陆面积较小，结果不是很清楚，图中环形区域仅为示意图。

图3是利用2D TPS插值的结果，由于2D TPS仅考虑了地理坐标(x,y)，没有高程值参与插值，在地理位置相近区域，降水量相差不多。因此，许多插值细节不能体现出来。

图4是利用普通克里金插值的结果，该方法仅考虑地理位置、没有考虑高程值对降水量的影响，因此，除了与2D TPS方法有类似的问题外，普通克里金插值还出现了一些斑块状区域以及亮点，在局部区域形成个别“牛眼”现象，说明与周围区域颜色有较大差异。这可能由于采用的气象站数量有限且分布不均匀，同时没有考虑高程值影响造成的。

图2 基于3D TPS插值的中国年均降水量空间分布图Fig.2 Spatial distribution map of annual precipitation in China based on 3D TPS interpolation

图3 基于2D TPS插值的中国年均降水量空间分布图Fig.3 Spatial distribution map of annual precipitation in China based on 2D TPS interpolation

图4 基于普通克里金插值的中国年均降水量空间分布图Fig.4 Spatial distribution map of annual precipitation in China based on ordinary Kriging interpolation

2.2精度评价 为了评价3D TPS插值方法的精度，进一步进行各种插值方法的对比分析。除了TPS外，还有反距离加权插值法(IDW)、样条法，克里金插值再细分为普通克里金和泛克里金两种插值方法。从454个气象站点中，将20个气象站点观测值作为校验值，为了提高可信度，这20个站点是随机选取的，且尽量覆盖整个研究区域；这20个站点不参与插值，仅用来验证插值结果。

3D TPS的SD值最小(SD=113.95)，插值精度最高，2D TPS次之(SD=156.82)，再次是泛克里金(SD=160.97)，与2D TPS相差不多，其他依次是IDW(SD=181.56)和普通克里金(SD=197.19)，精度最低的是样条法(SD=241.49)(表1)。由于20个气象站点的选取是随机的，且不参与插值，仅用于验证，因此，结果的可信度还是较高的。

进一步分析得知，绝对误差最大的分别是元江气象站和五台山气象站，如，元江气象站的绝对误差最小的是由3D TPS插值得到390.81；而绝对误差最大的是样条法，高达707.94。五台山气象站除3D TPS插值得到的绝对误差55.67外，其他方法插值得到绝对误差都在300左右。这两个气象站共同特点是分别与邻域气象站点的高程值有较大差异，如：元江气象站高程为401 m，而周围气象站点高程都在1 000 m以上，大大低于周围气象站的高程值；五台山气象站的高程为2 208 m，而周围气象站的高程值除少数略高于1 000 m(但在1 500 m以下)外，大多低于1 000 m，远高于周围气象站高程值。其他绝对误差较大的气象站，如，安庆气象站等，情况类似。由于3D TPS插值充分考虑了高程值的影响，在地形起伏较大、高程变化明显的区域，插值精度最高，优势较为明显。

表1 各种插值方法的比较Table 1 The comparison of different interpolation methods

最后，将本研究与其它文献研究进行比较。文献[17]应用趋势面模型对我国34个重要城市模拟，得到SD=155.95，仅大于本研究3D TPS模型的SD，即文献[17]应用的趋势面模型插值精度低于3D TPS，而高于本研究中其它5种插值方法。由于文献[17]也考虑了高程值的影响，所以会比本研究中其它5种方法(没有考虑高程影响)插值精度较高。但是，在同样考虑了高程值影响的情况下，3D TPS插值精度更高。

3 结束语

我国是一个多山的国家，地形变化复杂，草地资源分布广阔。在降水量空间插值研究中，高程值往往被忽略，以致造成较大误差，本研究采用3D TPS对降水量进行空间插值，除了考虑地理位置外，还将高程值作为重要的插值因子，充分考虑了高程值对降水量的影响。通过试验，与其他插值方法进行比较分析，得出3D TPS对降水量进行空间插值的精度最高，尤其在高程值变化较大的区域优势更为明显。这为森林、草地植被、水资源、生物量等空间分布的研究提供了重要的参考依据。

[1] 崔庆虎,蒋志刚,刘季科,等.青藏高原草地退化原因述评[J].草业科学,2007,24(5):20-26.

[2] 郭姜宁.森林及草地植被与降水[J].中国草业科学,1987,4(4):58-60.

[3] 刘桂霞,苗玉华,李记开.放牧干扰和年际间降水量变化对地木耳生长速度的影响[J].草业科学,2011,28(9):1649-1652.

[4] 张晓松.我国人均水资源为世界平均水平1/4[J].国土经济,2001(1):46.

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

[6] 储少林,周兆叶,袁雷,等.降水空间插值方法应用研究——以甘肃省为例[J].草业科学,2008,25(6):19-23.

[7] 马轩龙,李春娥,陈全功.基于GIS的气象要素空间插值方法研究[J].草业科学,2008,25(11):13-19.

[8] Duchon J.Splines minimizing rotation-invariant semi-norms in Sobolev spaces[J].Constructive Theory of Functions of Several Variables,1977,571:85-100.

[9] Hutchinson,M.F.Interpolating mean rainfall using thin plate smoothing splines[J].International Journal of Geographic Information Systems,1995,9(4)：385-403.

[10] Wahba G.Spline Models for Observational Data[M].Philadelphia:Society for Industrial and Applied Mathematics,1990:10-99.

[11] 孙海燕,丁咚.薄板样条函数及复杂曲面的数学表示[J].测绘工程,2006,15(2):7-8.

[12] 阎洪.薄板光顺样条插值与中国气候空间模拟[J].地理科学,2004,24(2):163-169.

[13] Sinthanayothin C,Bholsithi W,Tharanon W.Tooth alignment of the dental cast using 3D thin plate spline[A].Proceedings of the 3rd IASTED International Conference on Advances in Computer Science and Technology[C].Anaheim,CA,USA:ACTA Press,2007:255-259.

[14] 游松财,李军.海拔误差影响气温空间插值误差的研究[J].自然资源学报,2005,20(1):140-144.

[15] 徐吉宏,柳小妮,张德罡,等.基于ArcGIS的天祝草地综合顺序分类研究[J].草业科学,2011,28(5):866-870.

[16] 朱猛蒙,孙玉荣,张蓉,等.基于GIS的苜蓿斑蚜区域化预测预报技术初步研究[J].草业学报,2011,20(2):163-169.

[17] 尚宗波,高琼,杨奠安.利用中国气候信息系统研究年降水量空间分布规律[J].生态学报,2001,21(5):689-694.