几种插值方法在微DEM构建中的应用

宋向阳,吴发启

（西北农林科技大学资源环境学院,陕西杨陵 712100）

DEM是地形分析的主要数据,其精度对地形分析结果的精度具有明显影响[1]。目前,利用离散高程点数据进行空间插值来构建DEM是生成地形分析数据的主要途径[2],也是ArcGIS软件地形分析的重要功能之一[3]。一直以来,空间插值方法对DEM精度的影响都得到研究人员的广泛关注。早在1988年,Hutchinson提出利用高程信息生成DEM的Hutchinson算法[4],并在其基础上开发出专业化软件ANUDEM[5];邓小波等[6]研究了不同插值方法对气候要素分析的影响,认为改进插值方法和优化插值参数能够改善插值效果。另外,汤国安等[7]的研究表明,不同的插值方法生成的DEM具有明显的差异,而这种差异对地形分析结果的影响比较显著;孟庆香等[8]研究基于GIS的黄土高原降水量和年均温空间插值,结果表明不同插值方法构建的降水量和年均温分布图差异较大;李新等[9]通过比较研究不同空间插值方法差异,结果表明要得到理想的空间插值效果,必须针对研究区的实际情况,在对实测样点数据进行统计分析和反复试验比较的基础上选择最佳的插值方法。前人的这些研究成果主要是针对宏观尺度的DEM构建及分析,而关于微观尺度DEM构建方法的研究还几乎是一个空白。

众所周知,黄土高原是我国乃至世界上最为严重的水土流失地区之一。在该区坡耕地的土壤侵蚀量占到总侵蚀量的50%～60%[10]。虽该地区造成土壤侵蚀的因素众多,但人为土地管理方式是影响土壤侵蚀的主要原因之一。这些管理往往在耕地表面留下土粒级和土块级的高低起伏,直接影响到产流、汇流和泥沙运动。因此,这种微地形特征(或称为地表糙度)的研究,对揭示水土流失成因有重要作用。鉴于此,本文在模拟该区比较常见的三种耕作措施,即等高耕作、人工掏挖和人工锄耕,并以直线坡面为对照措施,利用非接触式激光扫描仪测量这几种措施下的地表高程的基础上,结合GIS软件中数据分析工具,以均方根预测误差、误差平均值和误差图等指标对测定的地表高程数据进行探索性分析,以便确定构建微DEM的最佳插值方法,服务于黄土坡耕地土壤侵蚀机理的研究。

1 研究方法

本研究借用了以下几种前人的研究方法来建立微DEM。

1.1 反距离加权插值法

反距离加权插值法(Inverse Distance Weighting,简称IDW)[11]是利用邻近已知点的数值进行加权运算,所需的权重根据距离的幂来确定,离插值点越近的样本点赋予的权重越大。权重显著影响内插结果,其选择标准是最小平均绝对误差[12]。Husar[13]等人的研究结果表明,距离的幂越高,内插结果越具有平滑的效果。该方法的插值计算是一个均分过程,要求采样点均匀分布,并且密集程度足以满足在分析中反映局部表面变化。

1.2 局部多项式插值法

局部多项式插值法(Local Polynomial Interpolation,简称LPI)[14]是一种局部加权最小二乘拟合法,根据有限的采样数据,采用多个多项式来拟合表面,它引入了“距离权”的概念,对于未知样点的计算,考虑在局部范围内全部已知样点对其贡献,距未知样点近的点权重较大,距未知点远的点权重较小,每一个未知样点的预测值都对应一个多项式,每个多项式都处在特定重叠的邻近区域内,通过最小二乘法求解邻域内多项式组成的方程组来得到拟合表面。

1.3 径向基函数插值法

径向基函数插值法(Radial Basis Function,简称RBF)[15]属于人工神经网络方法的一种,该方法是根据有限采样数据,选择合适的径向基函数生成一个具有最小曲率,且到各样点的Z值的距离最小的曲面。该方法所拟合的表面经过所有样点数据,并可以计算出高于或低于样点Z值的预测值,适用于采样点数据集大、表面变化平缓的情况,当局部变异性大,且无法确定样点数据的准确性,或样点数据具有很大不确定性时,该方法并不适用。

1.4 克里格插值法

克里格插值法(Kriging)[16]是用协方差函数和变异函数来确定高程变量随空间距离而变化的规律,以距离为自变量的变异函数,计算相邻高程值关系权值,在有限区域内对区域化变量进行无偏最优估计的一种方法。从数学的角度来讲,克里格插值法包括析取克里格(Disjunctive Kriging)、普通克里格(Ordinary Kriging)和概率克里格(Probability Kriging)等。对不同的样本数据,按照其数据特征的不同应该选择不同的插值方法才能取得较好的插值效果,如普通克里格法前提条件是样本数据符合正态分布;若样本不服从正态分布时,一般选用析取克里格;当同一事物的两种属性存在相关关系,且一种属性不易获取时,可选用协同克里格方法,借助另一属性实现该属性的空间内插。

2 试验数据的采集与处理

2.1 试验设计

试验在中国科学院黄土高原土壤侵蚀与旱地农业国家重点实验室降雨大厅进行。试验小区规格为1.0 m×1.0 m×0.5 m(以下简称:小区),每个措施设置3个重复。小区地表处理⑴土样准备:将过筛(筛孔0.5 cm)的土样分层填装在小区中,并将土壤容重控制在1.20～1.30 g/cm3,含水量控制在10%左右;⑵耕作措施的布设:在填好土的侵蚀槽中人工构建不同的耕作措施,包括:①等高耕作:选用横坡耕作方式,垄高7～10 cm,垄距为30 cm;②人工掏挖:采用镢头掏挖地表,深度5～8 cm,间距20～25 cm;③人工锄耕:沿地表以传统方式锄耕,深度4～5 cm;④CK:平整坡面,不做任何处理(见图1)。

2.2 数据采集

本研究数据利用非接触式激光扫描仪对微地表进行扫描后获得,每个处理下数据为46行,40列,共1 840个高程值,行列方向上空间间隔均为0.02 m,其存储格式如表1。

2.3 数据处理

数据处理在ArcGIS 9.3中进行,使用探索性数据分析工具Geostatistical Analyst-Explore Data所提供的一系列图形工具探测数据分布、全局和局部异常值(过大值或过小值)、寻求全局的变化趋势、研究空间自相关,得到不同耕作措施下地表高程值统计特征,如表2。

图1 耕作措施布设

表1 高程采样数据

表2 不同耕作措施下地表相对高程统计特征表

从表2中可以看出,标准差分别为CK(0.004 8)＜等高耕作(0.012 4)＜人工锄耕(0.012 7)＜人工掏挖(0.018 3),表明采样点数据具有较高的精度。同时反映出规则采样对地表微地形起伏状况变化不明显的CK表现能力最强,其次是等高耕作,而人工锄耕和人工掏挖相对较弱;地表微地形起伏状况数据表明CK的起伏度(0.023)＜人工锄耕(0.063)＜等高耕作(0.075)＜人工掏挖(0.094),同时最大高程和最低高程数据都出现在人工掏挖措施下的采样数据中,分别为0.270和0.176,这表明人工掏挖措施微地形局部变异较大;空间变异系数的变化分别为等高耕作(0.99)＞人工掏挖(0.97)＞CK(0.87)＞人工锄耕(0.72),表明采样数据均具有较强的空间相关性;峰度和偏态数据均不为0,这表明采样数据均不符合正态分布,但等高耕作措施下的地表微地形数据更接近正态分布,其偏态和峰度分别为0.003 2和1.811。

结果表明,由于偏态和峰度以及直方图、QQplot图等表明采样点数据不符合正态分布。所以,不宜采用采用普通克里格、简单克里格和泛克里格插值法。但是,高程采样点数据足以反映地表微地形局部表面变化且具有较强的空间相关性,达到IDW插值法、局部多项式插值法、径向基函数插值方法和析取克里格插值法对空间数据的要求,因此,理论上可以采用这些方法进行插值分析。

3 结果与分析

3.1 邻域搜索半径的设置

对于高程这种区域化变量来说,邻域搜索半径的大小不仅决定着邻域采样数据对预测点的权重,而且影响插值精度和插值表面的光滑度。因此,根据采样点数据特征,设置邻域搜索合理的邻域搜索半径是空间插值的基础,也是插值精度的保证。

对于等高耕作,由于采样数据横向相邻点高程相似度较高,区域扇形设置应更加扁平且平行于等高线方向,扇区椭圆长半轴(Major semiaxis)为 0.075,短半轴(Minor semiaxis)为0.03,角度(Angle)为90,步长0.02,搜索半径内使用最小点数为横向相邻的2个点,最大为横向4个。

对于人工掏挖、人工锄耕和直线坡面,由于采样数据横向与纵向相邻点高程相似度差异不大,区域扇形设置应接近圆形。扇区椭圆长半轴为0.04,短半轴为0.04,角度为0,步长为0.02,搜索半径内使用最小点数为中心4个点,最大点数为12个。

3.2 不同插值方法比较

通过上面的分析,在了解了数据特征的基础上,本文拟选择反距离加权插值法、局部多项式插值法、析取克里格插值法和径向基函数插值法,分别对4种微地形数据进行插值分析。通过多次计算,比较误差平均值、均方根预测误差、标准误差预测图以及使用软件中的Compare功能等,得出针对不同耕作措施下4种空间插值方法的最优参数和最佳插值结果,如表3。

表3 空间插值方法比较分析

由表3可以看出,不同插值方法对于不同地表微DEM构建的差异较大。对于等高耕作来说,球面模型的析取克里格插值法优于其他三种,其绝对误差平均值为DK(0.19)＜LPI(0.32)＜IDW(1.20)＜RBF(1.34),同时其均方根预测误差也最小(2.26)。这表明析取克里格插值方法对于具有空间相关性但不符合正态分布的高程数据是最优插值方法;对于人工掏挖来说,Power=2时的局部多项式插值法明显优于其他方法,其绝对误差平均值为 LPI(0.07)＜RBF(0.13)＜DK(0.25)＜IDW(0.63),均方根预测误差为LPI(3.31)＜RBF(3.60)＜DK(6.02)＜IDW(6.21)。这表明地表微地形存在局部短程变异时,应该优先考虑局部多项式插值方法;对于人工锄耕来说,基于规则样条函数的径向基函数插值方法是最优插值方法。其绝对误差平均值为 RBF(0.05)＜LPI(0.10)＜IDW(0.28)＜DK(1.43),均方根预测误差为RBF(4.46)＜DK(5.76)＜LPI(6.49)＜IDW(7.24)。这表明基于规则样条函数的径向基函数插值方法更适合于起伏状况变化不大的微地形表面;对于CK来说,各种插值方法精度差异不大,其绝对误差平均值为RBF(0.02)＜DK(0.03)＜LPI(0.07)＜IDW(0.10),均方根预测误差为LPI(1.45)=IDW(1.45)＜DK(1.47)＜RBF(1.49)。这表明当高程采样点数据足以表现地表微地形局部变化时,对于微地形起伏状况无明显变化的区域,插值方法对插值表面精度影响不大。

图2显示了利用不同插值方法对不同地表微地形插值计算的效果。从图中可以看出不同插值方法在微DEM构建过程中精度差异。对于等高耕作和人工掏挖来说,局部多项式插值法和析取克里格插值法效果相近且都较好,而基于规则样条函数的径向基函数插值法和反距离加权插值法效果较差,这和表3中结果相同。这表明局部多项式插值法和析取克里格插值法对于微地形变化具有一定方向性以及存在局部短程变异的地表能够取得较高的精度和更加光滑的表面;对于人工锄耕来说,基于规则样条函数的径向基函数插值方法效果最佳,其次是局部多项式插值方法和析取克里格插值方法,反距离加权插值方法效果最差。这一结果与表3中反映的精度差异一致,这表明地表微地形起伏变化较为平缓区域应该优先考虑基于规则样条函数的径向基函数插值方法;对于CK来说,析取克里格插值表面最为光滑,局部多项式插值方法次之,径向基函数插值方法和反距离加权插值方法较差。虽然表3中结果表明各插值方法对CK的差异不大,但是通过图2表明插值精度较高的插值方法取得的微地形表面不一定光滑。因此,在插值过程中还要在插值精度较高的方法中选择效果最佳的插值方法。

图2 利用不同插值方法对不同地表微地形插值计算的效果

3.3 最优插值方法的选择

通过对比表3中不同耕作措施下各插值方法的误差平均值、均方根预测误差,来定量判断各插值方法的准确性。同时,对比图2中插值表面的光滑程度来综合考虑不同耕作措施条件下构建微DEM的最佳插值方法。

4 结论

优选空间插值方法和最优参数,能够更加客观、准确地反映不同耕作措施下地表微地形起伏状况的空间分布特点和变化趋势,便于构建微地形高精度DEM。因此,本文对ArcGIS 9.3中4种常用空间插值方法进行比较研究,通过对误差平均值和均方根预测误差等指标的比较,得到不同耕作措施下构建厘米级微DEM的最佳插值方法。

结果表明,反距离加权插值方法在4种耕作措施下误差平均值和均方根预测误差皆比较大,因而不适合进行微地形的构建。这主要是由于该方法仅考虑距离作为区域变量相关性的度量,而忽略了样点的空间分布格局;局部多项式插值方法在人工掏挖措施下两指标都比较小,这主要是由于该方法能够较好地保留地表微地形细节变化,邻域重叠保持了插值表面的连续性,更适合于存在局部短程变异的数据集;析取克里格插值方法在等高耕作措施下两指标都比较小,这主要是因为该方法考虑了已知样本点的空间分布及与未知样点的空间方位关系,更适合于微地形起伏状况具有一定方向性变化的数据集;基于径向基函数的规则样条函数插值法在人工掏挖措施下两指标都比较小,其主要原因是规则样条函数更适合于表面变化平缓的微地形的拟合;以上几种插值方法在直线坡面耕作措施下,两指标都最小,这主要是由于采样数据足够精确情况下,各种插值方法表现差异不大所致。

因此,等高耕作措施微地表采用基于球面模型的析取克里格插值法精度较高;人工掏挖措施微地表采用局部多项式插值法效果较好;人工锄耕措施微地表采用规则样条函数的径向基函数插值方法精度较好;对于直线坡面上述几种插值方法的精度差异不大。

[1]李志林,朱庆.数字高程模型[M].2版.武汉:武汉大学出版社,2007.

[2]Wackernagel H.Multivariate Geostatistics:An Introduction with application[M].3rd Edition.Berlin:Springer Verlag,2003:416-429.

[3]王劲峰,李连发,葛咏,等.地理信息空间分析的理论体系探讨[J].地理学报,2000,55(1):92-98.

[4]Hutchinson M F.Anew procedure for girding elevation and strearm line data with automatic removal of spurious pits[J].Journal or Hydrology,1989,106(3/4):211-232.

[5]Hutchinson M F.ANUDEM version 5.1 User Guide[C]//Centre for Resource and Environmental Studies.The Australian National University,Canberra,2004:1-22.

[6]邓小波,罗宇翔,于飞,等.西南复杂山地农业气候要素空间插值方法比较[J].中国农业气象,2008,29(4):458-462.

[7]杨昕,汤国安,刘学军,等.数字地形分析的理论方法与应用[J].地理学报,2009,64(9):1058-1070.

[8]孟庆香,刘国彬,杨勤科.基于GIS的黄土高原气象要素空间插值方法[J].水土保持研究,2010,17(1):10-14.

[9]李新,陈国栋,卢玲.空间内插方法比较[J].地球科学研究,2000,15(3):260-264.

[10]Tang Keli,Zhen Fenli,Ca Xuan.Soil erosion on the sloping farmland in the loess plateau of China[C]//Proceedings of the Fourth International Symposium on River Sedimentation.Beijing:China Ocean Press,1989.

[11]常文渊,戴新刚,陈洪武.地质统计学在气象要素场插值的实例研究[J].地球物理学报,2004,47(6):982-990.

[12]方书敏,钱正堂,李远平.甘肃省降水的空间内插方法比较[J].干旱区资源与环境,2005,19(3):47-50.

[13]Husar R B,Falke S R.Uncertainty in the spatial interpolation of PM 10 monitoring data in Southern California[EB/OL].http://capita.wustl.edu/CA PITA/CapitaReports/CaInterp/CaINTERP.HTM L,1997-03-03/1999-10-25.

[14]Fan Jianqing,Irè ne Gijbels.Local polynomial modelling and its applications[M].Chapman&Hall/CRC,2003:159-216.

[15]周志刚,陈丽红.RBF神经网络及其在数值计算中的应用[J].复杂系统与复杂性科学,2006,2(3):64-68.

[16]曾怀恩,黄声享.基于Kriging方法的空间数据插值研究[M].测绘工程,2007,16(5):5-13.