基于不同算法的剖面曲率提取研究

1、第15页（共19页）基于不同算法的剖面曲率提取研究引言地形因子是为定量表达地貌形态特征而设定的具有一定数学意义的数学参数或指标1。剖面曲率（Profile Curvature）KV 是地形因子之一，基于提取算法的坡面因子分，剖面曲率属于二阶坡面因子。在地理信息系统（GIS）中，剖面曲率的提取算法有三种：汤国安提出的SOS/SOA算法、四次表面模型算法和基于空间矢量模型的差分计算算法。数字高程模型（DEM）作为地理信息系统( GIS)空间数据库的核心数据库之一，是进行二维地形空间模拟和地形图生产的基础，DEM 能派生出各种地形因子，如坡度、坡向、平面曲率、剖面曲率、汇水面积等1。其中，剖面曲率是

2、最重要的地形因子之一，地面曲率影响着，也是制约生产力空间布局及自然地理研究的重要因子 。尽管利用GIS 平台在DEM 上提取坡度的算法已经成熟，但是由于算法会影响坡度的测算，所以需要比较不同算法对提取坡度的影响。与坡度、坡向等基本地形参数一样，剖面曲率在数学定义上是明确的和唯一的，具有确定的表达式。尽管其理论定义是明确的，然而DEM是地形曲面的微分模拟，算法设计必然存在各种各样的假设，不同假设和前提导致不同的剖面曲率计算模型和结果，这虽然对地形特征的可视化和地形分类的影响不大，但对以数值计算为主的地学分析模型的影响却是非常显著的2。剖面曲率的含义剖面曲率是地面曲率在垂直方向上的分量，是对地面坡

3、度的沿最大坡降方向的地面高程变化率的度量，是研究地形表面的重要因子之一。剖面曲率因子提取算法的基本原理为：将离散的高程数值拟为一个连续的曲面，基于微分几何的思想，模拟曲面上每一个点所处的平行于水平面的曲线，利用曲线曲率的求算方法推导得出各个曲率因子的公式，进而求算出每一点的曲率值1。剖面曲率KV 具有重要的理论和实践意义。KV 决定地表及土壤中物质移动的相对速度，KV 0表明移动加速，KV 0表明移动减速，反映了侵蚀或沉积的程度。KV 通常可以揭示出谷脊线呈突起状或阶梯状，还可以用在很多领域，如在大尺度区域中模拟土壤湿度、土层厚度、土壤侵蚀、山崩石头的散布、植被等；也可用于不同尺度的地貌研究，

4、研究区域尺度和次大陆尺度地质构造、地貌和土壤的空间相关规律；还用于识别显示断层、地质断陷线等3。笛卡尔坐标系中，x、y 表示空间位置，z = f( x，y) 表示地形曲面函数，假定函数具有单值性、连续性和二阶可导性，则剖面曲率的计算公式为：Kv=-p2r+2pqs+q2tp2+q2(1+p2+q2)3/2 （1）其中，p=zx ，是x方向高程变化率；q=zy ，是y方向高程变化率；r=2zx2 ，是x方向高程变化率的变化率；s=2zxy ，是x方向高程变化率在y方向的变化率；t=2zy2 ，是y方向高程变化率的变化率剖面曲率提取算法一般都是根据P.A.Borrough提出的窗口微分法遍历整个图

5、幅进行计算的。即是在33的DEM栅格窗口中进行的，窗口在DEM数据矩阵中连续移动后而能完成整幅图的计算工作4。Z1Z2Z3Z4Z5Z6Z7Z8Z9图1 33地形窗口(Z1Z9分别代表各个栅格的高程值)剖面曲率的计算方法GIS中剖面曲率的计算大多是基于DEM数据的，由离散的数据计算剖面曲率有不同算法，在考虑到算法的通用性和使用性的情况下，这里选择了ArcGIS采用的汤国安提出的SOS/SOA算法、四次表面模型算法和基于空间矢量模型的差分计算算法1三种常用的剖面曲率算法进行分析。汤国安提出的SOS算法剖面曲率为地面高程变化的二次导数，即坡度在垂直方向的最大变化率；Tang已经证明5，在Arc/In

6、fo等GIS软件中，地面剖面曲率的数字矩阵可以直接通过对DEM数据的求取地面坡度的坡度（Slope of Slope,SOS）而获得。SOS就是坡度的变率，是在所提取的坡度值的基础上对地面每一点再求算一次坡度，坡度是地面高程变化率的求解，因此，坡度变率表示了地球表面高程相对于水平面变化的二阶导数3。这里的“坡度”并非是坡度，而是按照坡度的算法实现了对坡度数据阵列中其变化率的量化提取。所获得的剖面曲率值在0,90内，并非是真正的曲率值，但却真实的地反映了坡度和等高线的变化率，可提取山谷、山脊等地形因子，而且因其算法简单易于实现，所以也称剖面曲率，有些教材中称之为地面变率因子3。3.2 四次表面模

7、型Zeverbergen 和 Thorne(1987)提出了一种计算算法，构造一个含有9个参数的四次方程（x、y的最高次数为2）6，见公式（2）.函数的所有参数由33窗口（如图1）内的单元值唯一确定，函数经过每个数据点，包括中心单元。系数A、B、C、I在计算曲率时未用到，可以忽略，系数DH由公式（3）（7）计算，则中心格网的剖面曲率的表达式变为表面函数系数组成的式子，见公式（8）3。Z = Ax2y2 + Bx2y + Cxy2 + Dx2 + Ey2 + Fxy + Gx + Hy + I. (2)D =（Z4 + Z6）/2-Z5/L2 . (3)E =（Z2 + Z8）/2-Z5/L2

8、. (4)F =(- Z1 + Z3 +Z7 Z9)/4L2. (5)G =(- Z4 + Z6)/2L. (6)H =(Z2 Z8)/2L . (7)剖面曲率：KV = 2(DG2 + EH2 + FGH)/(G2 + H2). (8)3.3 基于空间矢量模型的差分计算 在DEM数据的基础上，从微分几何的思想出发，求出剖面曲率公式(1)中的各个导数pt,即可求出各曲率值，差分计算常用的是三阶反距离平均权差分。对每个栅格点都确定一个33分析窗口，如图2（1），ai分别代表各个栅格的高程值，L表示栅格大小，p、q可由以下公式求出3。p=zx=c+2f+i-(a+2d+g)8L （9）p=zx=a

9、+2b+c-(g+2h+i)8L (10)abcaxbxcxaybxcxdefdxexfxdxexfxghigxhxixgxhxix （1） （2） （3）图2 提取地面曲率的分析窗口 求出所有 p、q 值，则组成 p 值栅格图（ 见图1（ 2 ） ） 和 q 值栅格图（ 见图1（ 3 ） ） 。 以此数据为基础，进行差分运算，可求出r、s、t8,见公式( 10 ) 、（ 11 ） 和（ 12 ） 。r=2zx2=cx+2fx+ix-(ax+2dx+gx)8L （11）s=2zxy=ax+2bx+cx-(gx+2hx+ix)8L （12）t=2zy2=ay+2by+cy-(gy+2hy+iy)

10、8L （13）由不同算法提取的剖面曲率由任何一种算法得出的剖面曲率，我们都要评估其精确性和制图效果，并以此来作为我们求算地形的剖面曲率的依据。因此，一个科学、合理、准确的剖面曲率分析是目前研究剖面曲率的首要任务9。在此，以西峡县的某一地区为研究区，我将依据三种剖面曲率算法所作出的图像和数据进行比较分析，得到柱状图和各统计量，为找到不同高程域的适合剖面曲率提取算法提供依据。4.1 研究区概况研究区地处伏牛山腹地，豫、鄂、陕交界地区，东临南阳盆地，西嵌秦岭，南襟荆襄平原，北负八百里伏牛山。境内地形复杂，北部是海拔高、坡度大的中低山地，南部是鹳河谷地，两侧是起伏大的低山丘岭。海拔181米，自然坡降为

11、33%。研究区内的高程分布如表1所示。由表1知，分布在120405之间的高程数值最多，为3197380个，约占总数值个数的27.9%，分布在405690之间的高程数值约占总数值个数的22.5%，分布在690991之间的高程数值约占总数值个数的21.5%，分布在9911316之间的高程数值约占总数值个数的17.8%，而分布在13162171之间的高程数值最少，为1186783个，约占总数值个数的10.4%。并且根据图3中所示可判断该研究区的地形呈阶梯性抬高，总体上是北高南低。表1 原DEM的高程值分布范围图3 原DEM重分类图图4和图5中的三维图是在ArcScene中制作的，原DEM的高程变化不

12、是很明显，仅从颜色的分类可以知道高程的分布，和地形的起伏。但用高程值的3倍显示，则可以看出此研究区的地形有一定的起伏度，但不是很大。图4 原DEM的三维显示侧视图图5 原DEM高程的3倍的三维显示侧视图由汤国安提出的SOS算法提取的剖面曲率由图6知，由SOS算法提取的剖面曲率的图像的南部，有两道狭长而倾斜的浅颜色部分，这部分的变率在0.0013778 547.339714169之间，说明这部分的剖面曲率变化比较小，地形相对平坦。图6 汤国安的SOS算法剖面曲率4.3 由四次表面模型算法提取的剖面曲率由图7知，由四次表面模型算法提取的剖面曲率的图像相比图6的相对区域，图7的这两个狭长部分的变率则

13、在-0.7680859470.646213027之间，剖面曲率的变化比图3变化更小，而且其间的其他数值范围的变率也较之增加。图7 四次表面模型算法剖面曲率4.4 由基于空间矢量模型的差分计算算法提取的剖面曲率由图8知，由基于空间矢量模型的差分计算算法提取的剖面曲率的图像中除了空数值为零外，几乎所有的变率数值都分布在-0.1462174060.000303681之间，只有零星的部分分布在0.0003036810.027690800之间。与前两种算法相比，此算法的差别较大。图8 基于空间矢量模型的差分计算算法剖面曲率试验结果及分析5.1 不同算法提取剖面曲率的统计指标分析由图9可知，由汤国安的SO

14、S算法提取的剖面曲率的数值偏正态分布，主要集中在2.335129512.303019之间，虽然数值不代表真正的曲率值，但能够表示曲率的变化率。图中的数值多处于图像的正区，表明用此算法提取的该研究区的剖面曲率偏大，且地形多为凸地形，呈阶梯性抬高地表。图9 由SOS算法得出的剖面曲率柱状图由图10知，由四次表面模型算法提取的剖面曲率的数值是正规的正态分布图，且正区、负区基本持平，大多集中在-1.0182720.793010之间，表明用此算法提取的该研究区的地形剖面凸地形和凹地形分布均匀。图10 由四次表面模型算法得出的剖面曲率柱状图由图11知，由基于空间矢量的差分计算算法提取的剖面曲率的部分数值发

15、生异常，导致数据分析失真，但从分布的数值来看，数值主要集中在-0.003119700.00030368之间，数值多分布在负区，表明用此算法提取的该研究区地形剖面凹地形多与凸地形。图11 由基于空间矢量模型的差分计算算法得出的柱状图由表2知，汤国安的SOS算法的最大值与最小值之差最大，说明其值域幅度最大；基于空间矢量模型的差分计算算法的最大值与最小值之差最小，说明其值域幅度最小。而从标准差来看，基于空间矢量模型的差分计算算法最小，准确度最高；四次表面模型算法的准确度居中；汤国安的SOS算法标准差最大，准确度最差；用标准差除以平均值，得出其各自的变异系数，SOS算法的变异系数为0.66459639

16、74616885，四次表面模型算法的变异系数为-32.43107466561735，基于空间矢量模型的差分计算算法的变异系数为-1.172042558940668。表2 三种算法的统计量比较算法最小值最大值和平均值标准差SOS0.00137785445.64225006102631668.98.9639350395.957398934四次表面模型-13.29473412.46571159-369645.7671-0.0322851681.047042694基于空间矢量模型的差分计算-0.1462174060.027690800-13322.68346-0.0015469370.00181307

17、65.2 基于三种算法提取的剖面曲率的优缺点、误差及制图效果分析对运用不同算法提取的剖面曲率数据做统计分析,发现同一研究区不同算法所获得的剖面曲率最大数值、最小数值、标准差是不同的。并且，通过对各算法所提取出的剖面曲率做重分类，更显示出其差异。由图12可知，由汤国安算法提取的剖面曲率数值分布在0.0013783.581054之间的数值比较集中，且数值最小，并与前面图3、图4、图5相对应，表明用此算法提取这一部分区域的剖面曲率是适合的，与原地形相比，在研究区南部的低平地区，制图效果明显，说明此算法在提取这种简单地形有一定的优势。图12 SOS算法提取的剖面曲率重分类由图13可知，在图12的相对应

18、区域，用此算法提取的剖面曲率数值大约在-0.9701290.242128之间，数值偏小，且中间夹杂着其它类别的数据，并与前面图3、图4、图5所示有一定的偏差，表明此算法不太适合用于该研究区的剖面曲率提取，制图效果仍有一定的明显性。图13 四次表面模型算法的剖面曲率重分类由图14可知，由基于空间矢量差分计算算法提取的剖面曲率的图像在图12、图13的对应区域，剖面曲率数值分布在-0.0010660.000988之间，中间也有其它类别的数值，比图12显示的数值稍小，表明制图效果的地形起伏不太大，但也不平坦。图14 基于空间矢量的差分计算算法的剖面曲率重分类6 结论通过在不考虑DEM来源误差的情况下，

19、比较三种算法在研究区提取剖面曲率图，并统计各剖面曲率的最大值、最小值、平均值和标准差，以及对各算法提取的剖面曲率做重分类，进行综合分析。从三种算法提取的剖面曲率重分类图像来看，曲率数值值域最大的是汤国安的SOS算法，曲率数值值域最小的是基于空间矢量的差分计算算法，但是这两种算法都能凸显出该研究区部分区域的剖面曲率的特征。由图3的对原高程图像进行的重分类可知，该研究区的大致地形是由南部的最低处向两边逐级抬升，有三个最高区域，两个稍高区域，两个低谷区域和一个平坦区域。通过比较图像特征较明显的区域，即平坦区域，认为在该研究区提取剖面曲率选用汤国安的SOS算法、四次表面模型算法较为合理，而汤国安的SO

20、S算法在计算原理与方法上比四次表面模型算法更精细、更有特征性、逻辑性更强。虽然汤国安的SOS算法所计算的是坡度的坡度，不是真正的剖面曲率，但这种算法反映了剖面的变化率，且简单易操作，更能准确的突出研究区的特征。而四次表面模型算法所提取的剖面曲率在一定程度上体现了研究区的地形特征，具有极大的参考性。因此，在该研究区进行的基于不同算法剖面曲率的统计量分析时使用这两种算法是科学的，由表2知，通过数据对比，该研究区的凹地形地区多于凸地形地区，而那部分平坦地区类似于盆地地形。本文对算法的效率问题在提取剖面曲率时没有考虑，且仅分析了3种算法的成图数据统计量，DEM的精度、DEM的比例尺、DEM的分辨率、数

21、据的组织方式等都对地形的剖面曲率提取具有重要影响，因此，对于使用不同算法提取剖面曲率的研究，有待考虑更多的因素，发现更加普遍的规律。参 考 文 献1 汤国安，李发源，刘学军等.数字高程模型 M . 北京：科学出版社，2011:142157.2 李天文，刘学军，陈正江，汤国安，李军锋等.规则格网DEM坡度坡向算法的比较分析 J . 干旱区地理，2004, 27( 3) : 398 - 403. 3 范红艳，戚鹏程，黄天勇等.关于GIS中三种地面曲率的探讨 J .南阳师范学院学报，2011, 10( 6):66-69.4 刘学军，龚健雅，周启鸣，汤国安等.基于DEM 坡度坡向算法精度的分析研究 J

22、 .测绘学报，2004, 33( 3) : 258 - 263.5 陈楠，王钦敏，汤国安,赵元伟等.6种坡度提取算法的应用范围分析以在黄土丘陵沟壑区的研究为例 J .测绘信息工程，2006, 31( 4) : 20 - 21.6 Tang Guo-an A Research on the Accuracy of Digital Elevation ModelsM Beijing: Science Press Beijing， 20007 Zeverbergen L W，Thorne C R Quantitative Analysis of Land Surface TopographyJ Ea

23、rth Surface Processes and Landforms， 1987， 12( 1) : 47 568 刘学军，王叶飞，曹志东，汤国安等.基于DEM的地形曲率计算模型误差分析 J . 测绘学报，2006, 31( 5) : 50 - 53.9 陈楠，王钦敏，汤国安等.自DEM由不同算法提取坡度的对比分析 J .测绘工程，2006，15( 1) : 6 - 9. profile curvature extraction based on different algorithm researchAbstract:Without consideration of DEM error, this article through