基于DEM的地形曲率计算模型误差分析

基于DEM的地形曲率计算模型误差 基于DEM的地形曲率计算模型误差分析 第31卷第5期 2006年9月 测绘科学 ScienceofSurveyingandMapping Vo1.31No.5 Sep 基于DEM的地形曲率计算模型误差分析 刘学军??,王叶飞?,曹志东?,汤国安? (?南京师范大学地理科学学院,南京210097;?长沙理工大学公路工程学院,长沙410076) 【摘要】地形曲率是地形表面几何形态和地学建模的基本变量之一.本文首先对曲率计算模型进行了归纳,然 后通过数据独立的DEM误差分析方法和实际DEM的分析验证,对目前九种曲率计算的三类曲率计算模型进行了 量化分析比较.研究结果表明,当DEM数据精度比较高时,高次曲面(四次曲面)能给出较高精度的曲率计算结 果.而当DEM误差较大时,低次曲面(二次曲面)由于具有误差的平滑作用而能产生较高精度的曲率值. 【关键词】数字高程模型;地形曲率;算法,精度 【中图分类号】TP751【文献标识码】A【文章编号】1009—2307(2006)05—0050—04 1引言 地形曲率(TerrainCurvature)是表达地形曲面结构的主 要参数之一,也是地表过程模拟,土壤侵蚀模型,土地利 用分类等环境模型的基本变量.在地貌,地理,土壤,资 源与环境等领域有着重要的应用.在地理信息系统(GIS) 中,地形曲率计算在是数字高程模型(DEM)上通过一定的 曲率计算模型来实现. 与坡度,坡向等基本地形参数一样,各种类型的地形 曲率在数学定义上是明确的和唯一的,具有确定的表达 式….然而DEM是地形的模拟模型,不同观点在对地形的 表达上有不同的认识,从而导致不同的DEM语义模型. 例如从DEM上提取坡度坡向的算法目前就有十多种,而流 域网络提取算法则多达二十余种…,地形曲率计算也不例 外,正如Moore,Burrough,Zhou等人所指出的,在 基于DEM的地形参数提取和GIS地学建模中,模型的选择 是不可忽视的因素,直接影响到地形参数的计算精度和地 学建模的可靠度. 选择合理的地形参数计算模型必须对各种计算方法的 原理,特征,误差分布有深入的了解.目前对坡度坡向, 流域网络提取等的模型误差分析研究比较深入,如 Moore_6j, Flofinsky_7J,刘学军.等人的研究,而对地形 曲率计算研究相对较少.Evans,Sharry【12]等较为详细地 研究了各种地形曲率的定义,命名和分类体系,Florinsky?7 对常用的四种曲率模型和误差分布,通过误差传播定律进 行了分析比较;Scidmit则通过实验的方式对三种曲率计 算模型进行了研究.然而如文献[1]中所指出的,由于地形 本身的复杂性,对地形参数计算模型的误差分析应该在一 个客观,公正,量化的环境中进行.Florinsky的研究虽然 给出了曲率误差和DEM误差的关系式,但对模型的分析比 较却是建立在各个参数计算的精度分析上,忽略了模型参 数的相关性,同时他所弓l用的曲率计算模型也有错误. 作者简介:刘学军(1965.),男,汉族, 陕西合阳人,博士,教授,博士生导 师,主要研究方向为数字地形分析, GIS空间分析,空间数据不确定性等. E?mail:liuxuejun@njnu.edu.ca 收稿日期:2005—03—20 基金项目:国家自然科学基金 (40571120,40271089);南京师范大学 高层次人才科研启动基金项目 Scidmit_2通过实验环境来研究各种曲率模型的语义信息特 征,但没有量化地对曲率计算模型进行分析评价.基于上 述原因,本文以文献[1]中提出的数据独立DEM误差分析 方法为依据.对各种曲率计算模型的计算精度,对DEM误 差的敏感度进行了分析比较,以对各种曲率计算模型量化 地进行评价,并探讨曲率计算模型的适用范围. 2地形曲率计算模型 地形曲率在不同的层面上反映着地貌形态和结构特征. 为了实际应用的需要,人们定义了各种各样的曲率(参见 Sharry的研究).考虑到实际应用的需要和曲率之间的相 关性,本文选择其中的九种曲率作为分析研究对象,如表1 所示. 表1 名称符号 平均曲率eln 总高斯曲率cg 总曲率ctot 水平曲率cc 曲率定义(据SharryC整理) 单位计算公式 m—I+ _2( q2 p )r-2sp +1) q+川/ m一(—s)/(pq+1) m一r2+25+' m,一(rq一2spq+tp)/(pq) 剖面曲率cm,( p^+p2sp+q?'q+21)/) 流线路径曲率cfm一1[(p一 ( q p 2 2 )s + - g p ) q 3 ( ,2 r—f)]/ 经线曲率 交叉断面曲率 切曲率 一 (rp+2spq+tq)/(p+q) 一 (rq一2spq+tp)/(pq) 一 (rq一2spq+tp)/ (pq)/(p++1) 地形曲面函数为z=,(,y),P,q,r,s,t,分别为,, l—ll" 从表1的曲率定义不难看出,地形曲率计算的关键是 地形曲面在给定点处的各阶导数的估算.由于DEM是地形 曲面的离散化表达,而且地形曲面一般也不知道,因此在 DEM上对各阶导数的估算一般是通过局部窗口中(3×3窗 口,如图1)的曲面拟合方式进行.另外,为计算上的方 便,坐标原点常常设在3×3局部窗口中的中心点上,中心 点周围各点坐标与格网编号如图1所示.目前常用的曲面 拟合函数是多项式为: mn ,Y)=o'(1) 通过局部窗口中的各个高程点(,Y,)(i=1,2, 第5期刘学军等基于DEM的地形曲率计算模型误差分析5l … ,9),以最小二乘法求解多项式中的系数,从而求取各 个偏导数,然后按表1式中的公式计算中心点各种曲率. 由于局部窗口常常为3×3,已知点的高程为分布在中心格 网点周围的8个点,因而多项式的阶数不可能太高,一般 不超过4阶.在该局部窗口中的常用的多项式函数为: 二次曲面:z=觚+6y+cxy++ey+(简称E模型) (2) 限制二次曲面:=ax+6+cxy+++厂 且曲面通过中心点(简称S模型)(3) 四次曲面:z=axY+bxY+cxy++ ey+y++i(简称z模型)(4) 由(2),(3),(4)式,3×3局部窗口中的中心点处的 各阶导数估算公式如表2所示. 表23x3局部窗口导数计算模型 导 数 算法 E模型S模型Z模型 P(Z3%+Z9一z一:4一.1)/(6g)(z4一.6)/(2g) q(zj+2+3一7一8一z9)/(6g)(2一8)/(2g) [zj+z3+z4+z6+[zj+3+3(z4+z6) rz7+z9—2(2++7+z9—2(2+(z4+6—2z5)/gZ z5+z8)]/(3g)3z5+z8)]/(5g) s(3+Z7一zj+9)/(4g) [zj+z2+z3+Z7+[z1+z3+z7+z9+ tZ8+z9—2(z4+5+2(2+8)一2(Z4+(2+8—2z5)/g z6)]/(3g)3z5+z6)]/(5g) 实际上,不管采用何种曲面函数,其本质都是对中心格 网点处的导数的估计.因此,各种曲面函数模型下的导数计 算公式,都可统一到数值微分的框架下.下面以方向的导 数计算为例进行分析. 对于函数_厂(),存间距为g的等距节点—g,,g上 的函数值分别为f(—g),厂(),厂(g),则中问点处的导数 值为: _厂()=[+g)一—g)]/(2g)1, f()=[g)一2厂()+_厂(—g)]/gJ 上两式可改写成: _厂()={[+g)一)]/g+[)一—g)]/g{/21 厂()={g)-f()]/g一)一—g)]/g}/gJ (6) 令n=+g),b=),c=(—g),则(6)可表示成: _厂()=[(n—b)+(b—c)/g]/21,7, _厂()=[(n—b)/g一(b—c)/g]/gJ y方向的导数l石]理可导.式(7)为格网DEM上3×3 局部窗LJ的导数计算的通用公式,无论是何种计算方法,其 不同之处在于对a,b,c的估计的不同,亦即权值不同以及所 取点的位置的不同(周围位于坐标轴上的点或周围所有点). 对于z方法,仅取位于坐标轴线上的点并且等权,即a=z,b =z,c=z;对于二次曲面(方法),是用中心点5处周同位 于(或Y)方向的所有高程值的算术平均值作为n,b,c值 的,即:c=(zl++z7)/3,b=(z2+z5+三8)/3,r上=(z3+z6 +z9)/3;对于限制二次曲面(S方法),一阶偏导数估计与二 次曲面的计算方法类似,二阶偏导数给予位于坐标轴方向I 的点更大的权重,此处权值为3,即:c=(z,+3z+,)/5,b = (2+3z5+z8)/5,0=(z3+3z6+z9)/5. 仿上述过程,还可得到更多的曲率计算模型.本文重点 分析上述三种曲率模型. 3数据独立的误差分析方法与结果 图2是本文的实验数据分析,共包括三个环节,即 DEM数据无误差,施加随机噪音的DEM以及实际DEM分析 验证. I(-g~:I)(-g:,} 4l-g.g)S(-g,z1 g 7(-g~)8(-g~g.z) 图13×3局部窗口g 为格网分辨率 … 嘏率计锋梭精度I曲珥蔓引榄精,筻的RMSE捧}芋 RMSE排序分析l分析J}}=棚戈程J}杖点同)分析 图2数据独立DEM解译 算法误差分析流程 数据独立的DEM解译算法误差分析是文献…中提出的 用来解决各种地形参数模型精度的方法.其核心是建立一 系列和实际地形相似的模拟数学曲面,并按给定的分辨率离 散化建立相应曲面的模拟DEM.该DEM的特征是各个格网 点的高程数据没有误差,同时由于曲面的解析式表达式已 知,可获取地形参数的真值,从而通过真值和计算值的比较, 量化的分析各种曲率模型的精度. 本文的数学曲面采用高斯合成曲面,曲面表达式为(图 3): z= [1一()】e-(~-)2-(一[0.2(x--)一 (三)一(上)1e一()一(一Ce-((吉)+)一()(8)m凡J 式中A,B,C为地势起伏参数,m,n为范围控制参数.本 文格网分辨率为lOm,曲率计算模型精度采用中误差指标 (RootMeanSquareError,RMSE).在该曲面的模拟DEM上 的各种曲率计算模型统计分析结果如图4所示. 萋 图4高斯曲面模拟DEM曲 率模型中误差统计分布图 (图中数值经过对数变换) 数据独立的DEM虽然高程数据没有误差,能定量的描 述曲率计算模型精度,但不能反映曲率计算模型对DEM数 据误差的适应程度.为此,对曲面模拟DEM施加先验特征 已知的随机噪音,并按RMSE统计了DEM中误差.在含有 误差的DEM上,由于曲面表达式已知,也可定量的表达各 种曲率的计算精度.图5是在施加了不同程度随机噪音的 DEM上统计的各种曲率中误差. 由于在实际DEM的地形曲面表达式并不知道,无法获 取曲率的真值,但可通过散点图来反映各种曲率模型之间的 相关程度.散点图用不的算法构成纵横轴,通过点的聚散 程度来表述曲率模型的相关程度,如两算法相似,其散点应 集中在直线两侧,图形的发散与聚合反映了算法的相似程 度.图6是存施加了误差的模拟DEM上,不同曲率计算模 型的散点分布图. 图7是在四个实际DEM上各种曲率计算模型的散点 图.实验样区地貌类型为黄土高原,分别位于陕西北部的甘 泉和延安地区(图8),其中甘泉地区平均坡度30.9.,高差 179.5m,属黄土残塬地貌,地形支离破碎,延安样区平均坡度 30.5,高差193m,属黄土梁,地形相对简单. 52测绘科学第31卷 一一 Ji重{I!fI{I;&ji3一一}3!{l{I!fI{I国. 061l2427.6I1319506l_124276l1519.506llI2427.6l15195 一羞 7.6115 i19 j 5 :.7.6 l 19506242595062427.6595 . 06 , l242枣595lllllIIll1I . 一 曩I_0.60,4 一 l{I.5一一l. 一一 幽一}置,.一一矗. 蟋 图5加入随机误差的高斯曲面模拟DEM上各种 曲率计算模型精度统计(横轴为 DEM误差,纵轴为曲率模型误差) 21m巾误差cr Ef{l嘲 3qm巾澧蔗f—f 警 图6模拟DEM施加误差散点图(横纵坐标均为相应算法 的曲率计算值,以cf曲率为例,其他曲率有相似规律) H泉cH釉t-H'泉", 掣 E模 延安时 静 摹 E挺 诞安 静 t 擎 礁 砸密rI }:税楼 图7样区散点图(横纵坐标均为相应算法的 曲率计算值.以cc曲率为例,其他曲率有相似规律) 4分析与讨论 图8试验样区地貌 从图4可看出,当DEM数据没有误差时,无论何种曲率, 高次曲面(四次曲面,z模型)都给m比低次曲面(二次曲面,S 模型,E模型)精度为高的计算结果.实际上这一点不难理 解,由于此时DEM数据本身没有误差(仅含有DEM地形模拟 误差,即DEM格网逼近曲面所产生的误差,这不影响分析 计算结果).由数值分析理论知,任何复杂的曲面都可用高阶 多项式去逼近,多项式次数越高,越能精确描述曲面结构,同 时在DEM局部窗口中,四次曲面完全通过每一个格网点.因 此四次曲面比二次曲面更能准确刻画曲面形态,因而其曲率 计算的精度也就高些.s模型和E模型采用二次曲面,区别 在于是否通过中心点,因此两种算法的精度比较接近,两种算 法都不通过局部窗口的9个点,从而对局部窗口的数据产生 平滑.由于在数学曲面上离散生成的DEM格网点无误差,因 此这两个模型的在地形结构表达上不如四次曲面准确,基于 二次曲面的各种曲率计算误差也就相对较大. 当DEM含有误差时(图5),除等高线曲率CC,流水路径曲 率c厂外,随着DEM误差的增加,其余七种曲率的中误差随 DEM误差的增加而增加,并且E模型的曲率计算精度最高,S 模型次之,而z模型的误差最大,这与DEM数据无误差时的 分布基本相反.其本质还在于局部窗口的曲面拟合上.二次 曲面仅有六个未知数(系数),而局部窗口有几个已知点,已知 数据个数多于未知数据个数,在最小二乘原理下,二次曲面对 数据误差有一定的抑制作用,因而能给出较高精度的曲率计 算值.而四次曲面完全通过局部窗FI的各个格网点,具有唯 一 解,不存在多余观测数据,故而曲率计算精度较低. 在图5中,对于E,s,z三种模型而言,等高线曲率CC和 流水路径曲率c厂的误差变化与其他曲率误差变化并不一致, 等高线曲率CC误差随机波动,而流水路径曲率则完全呈相 反变化趋势.这说明这两种曲率的计算不但与曲面形式有 关,还与其他的参数有关.参看表1.并考虑到坡度计算 式…:tgS=+,,s为坡度,有: 等高线曲率CC=ct/sinS(9) 流水路径曲率cf=2cpXct/sin(2S)(10) 显然CC和c厂与坡度相关,坡度的存在增加了CC和c厂对 DEM误差的敏感性.由文献1知,同一点上的坡度,z模型 的计算值总是比S模型和E模型大,并考虑到S并非一个常 数和正弦函数的变化趋势,因此CC曲率的计算误差变化在 不同模型之间并非恒定,而z模型则给m比E和s模型误差 小的c厂值. 通过上面的分析,有如下的结论,当DEM数据误差比较 小的时候(<cm),高次曲面不但能够很好地描述局部地形 特征,而且也能给出比较准确的曲率估计值;而当DEM含 有误差时,二次曲面虽然在地形曲面形态刻画上不如四次曲 面,但却对DEM数据误差有一定的抑制作用,而这一点对基 于DEM的地形参数提取是至关重要的,因此二次曲面常常 给m比四次曲面精度为高的地形参数计算结果,这一点在一 些日前的研究结论中也不难发现'j.本文在实际DEM (图7)卡"模拟DEM上(图6)的结果也验证了这一结果,同为 二次曲面的s模型和E模型,各种曲率的计算结果比较接 近,并且S模型的精度要差一些,而基于四次曲面的z模型, 无论是与E模型还是与S模型,在散点图上和误差分布图上 相差都比较大,充分说明在基于DEM的地形参数提取上,低 次曲面比高次曲面具有较大的优越性. 5结论 本文首先将地形曲率计算模型统一到数值分析的框架 下,使得各个曲率计算模型有机的联系在一起.然后利用数 据独立DEM误差分析方法,定量地对各种曲率模型进行了 误差分析,同时通过实际DEM进行了各种曲率模型的相关 性分析.通过本文的研究,有如下的结论:?高次曲面虽然 能够较为准确的刻画复杂的局部地形曲面特征,但对DEM 误差比较敏感,并不适合作为DEM地形参数提取的曲面模 型;?在各种曲率模型巾,当DEM数据精度较高时,可采用 基于四次曲面的z模型进行曲率计算,而当DEM误差较大 时,以基于二次曲面的E模型为佳.考虑到目前各种精度级 . 第5期刘学军等基于DEM的地形曲率计算模型误差分析53 (上接第49页)法计算了南中国海及其邻近海域的垂线偏 差子午分量和卯酉分量,避免了按照Sandwell方法计算测 高卫星地面轨迹交叉点,并与EGM96模型的大地水准面梯 度,垂线偏差分量进行了比较,说明该方法的有效性,能 够使利用卫星测高资料反演高精度高分辨率的重力场成为 可能.由于该方法需要高分辨率的卫星测高资料,如果联 合多种卫星测高资料效果将会更好. 5(J}302f1i0【ll(12{)3fJ4【lSO 图3ENVISAT测高卫星计算的海域 垂线偏差南北分量图 2)按照间接平差计算,叼时存在法方程病态问,本文 采用分块矩阵求,'义逆的7Y~-q艮好地解决了病态法方程问题. 3)将瞬时海面高的每一周每一弧段分别进行一次差分 处理,作为星测高数据的初始观测值,由于同一弧段相 邻两点的采样几乎发生在同一时刻,故海面高差分值JL乎 消除了所有的测高系统偏差与环境改正误差.当联合不同 测高数据时,可以不考虑框架统一与海平面时变基准问题, 从而使多种卫星测高联合处理变得非常简单.与测高海面 处理不同的是,在记录瞬时海面高差分值时,同时需要记 录差分方向数据. 4)海面高差分是对海面梯度的直接采样,相对于Sand— well测高垂线偏差法反演重力场参数而言,由海面高差分得 到的海面梯度观测精度比交叉点方法间接求得的海面梯度精 度高,且空间密度也将得到大幅度增加.因此,利用海面高 差分反演海洋重力场参数比直接利用测高海面高反演有利. 参考文献 李建成,王正涛,胡建国.联合多种卫星测高数据 分析全球和中国海海平面变化[J].武汉测绘科技 大学,2000,25(4):343—347. 李建成,宁津生,陈俊勇,等.中国海域大地水准 面和重力异常的确定[J].测绘,2003,32 (2):ll4一ll9. HwangCW.InverseVening—Meineszformulaandde— flectiongeoidformula:Applicationsofthepredictionsof gravityandgeoidovertheSouthChinaSea[J].Jour- nalofGeodesy,1998,72(4):304—312. 姜卫平,李建成,王正涛.联合多种测高数据确定 全球平均海面WHU2000[J].科学通报,2002,47 (15):l187一l191. 李建成,姜卫平,章磊.联合多种测高数据建立高 分辨率中困海平均高模型[J]武汉大学,信 息科学版,2001,26(1):40—45. 彭富清,夏哲仁.卫星测高中的垂线偏差法[J]. 海洋测绘,2004,24(2):5—8. 鲍李峰,陆洋.由测高大地水准面密合平面拟合法 确定垂线偏差[J].自然科学进展,2003,13(3): 269—274. 1J1J]J]J]