一种顾及坡度因子与极值的地表面积计算法

程传录,蒋光伟,田晓静，马新莹，王文利

(国家测绘地理信息局大地测量数据处理中心，陕西 西安 710054)

一种顾及坡度因子与极值的地表面积计算法

程传录,蒋光伟,田晓静，马新莹，王文利

(国家测绘地理信息局大地测量数据处理中心，陕西 西安 710054)

针对局部地形的非线性变化与图形的极值信息对地表面积计算结果的影响，提出了一种顾及坡度因子与图形极值的地表面积计算法，分析了地表面积计算理论与误差。最后以中国西部具有典型地貌特征的两县数据为例进行试验分析。试验结果表明，该方法在一定程度上对局部地形逼近具有一定的改进，宏观上与规格格网交叉对角划线法具有一定的相似度。文中推导的三角网面积传递误差可初步表明考虑局部地形坡度因子与极值的地表面积计算结果可靠性更高。

坡度；地表面积计算；三角网；数字高程模型；破碎图斑

DEM是连续地形变化地形表面的离散化表达，是空间信息资料和地形分析的核心数据。对DEM信息挖掘可分析地表对象的内在空间特性、相互关系、分布规律和发展趋势。图斑表面积是地理国情监测的重要对象，而影响表面积计算的结果涉及坐标系统、投影变形、DEM精度、地表粗糙度、地形参数、内插算法、曲面逼近、原始数据分布等问题[1-2]，目前地表面积计算大都基于DEM数据。相关学者对DEM作了深入研究：文献[3—6]对DEM坡度、坡向的有效尺度范围及坡度计算的适宜窗口进行了研究；王耀革研究了DEM建模与不确定性分析[1]。在图斑面积计算方面：文献[7—8]研究了椭球体表面积计算方法，杨学祥建立了地球3个半周与地表面积之间的函数关系[7]；江帆研究了基于复化公式的DEM表面积算法分析[9]；陈吉龙运用微积分的方法对林地的表面积进行计算[10]；曹振宇研究了大地坐标系中多边形面积计算，提出了顾及地形起伏的任意多边形表面积计算方法[11]；党亚民研究了不同坐标系坐标转换及其对国土资源面积量算的影响[12]；薛树强研究了一种顾及复杂地形因子的地表面积计算方法，对比了不同分辨率DEM计算地表表面积的差异，得出顾及地形因子的地表面积统计可明显提高低分辨率DEM地表表面积的统计精度[13]；秦艳红从数学角度利用多分辨率分析，从理论上给出了实现由单一比例尺数据的坡度面积生成任意多比例尺系列坡度面积的依据和途径，并证明了小波变换前后面积差平方的对数与复杂地貌复杂度及尺度之间的变化关系，对不同尺度下坡度面积的变化进行了量化估计[14]。

地形表面的非线性变化及一些不确定性因素决定了地形表面表达的复杂度，任何曲面的逼近都是对地表形态的近似表达。本文在以上基础上研究了地表不规则图斑表面积计算的理论，阐述了影响地表面积统计的误差因素，重点分析了规格格网模拟地表趋势方法的差异，提出了一种考虑地形坡度因子与图形极值信息的表面积计算方法。

1 地表面积计算理论与误差分析

1.1 地表面积计算理论

设地表曲面S可表示为方程H=h(x,y)，D为曲面S在xOy面上的投影区域，函数h(x,y)在D上具有连续偏导数hx(x,y)和hy(x,y)，则曲面S面积可表示为

(1)

由式(1)可知，局部区域的地形曲面函数h(x,y)在边界区域内进行积分可得到地形表面积。实际的地形曲面h(x,y)很难得到，较常用的方法是对规则格网的DEM表面进行整体建模与分块建模。整体建模大都用于研究地形的趋势性，不能反映局部地形的粗糙度；分块建模能较好地顾及地形粗糙度与计算复杂度，设第i块函数为hi(x,y)，则地表面积可表示为

∑dSisec(β)=

(2)

由式(2)可知，地表面积计算可采用分块DEM模拟地形曲面hi(x,y)基于积分方法计算，也可使用投影面积与坡度来计算。实际的地形曲面是未知的，采用数值模拟的方法获取局部hi(x,y)计算复杂，常采用构建地表三角网法局部地形逼近

(3)

1.2 地表面积计算误差分析

当DEM格网分辨率越高时，每个格网认为是线性变化，此时DEM曲面常采用双线性曲面进行表达。李志林采用双线性模拟DEM表面[1]，给出了地形表面精度估计公式

(4)

(5)

εmax=cdtanβ

(6)

由式(4)—式(6)可知，在用DEM模拟地形时,误差主要分为两类:一类是生成DEM时节点的传递误差；另一类是用线性曲面代表地表时的损失误差。当格网分辨率越高时，地形模拟误差越小，地表面积越精确。在计算地表面积时，不同地形的DEM传递误差与地形曲面模拟误差所占的比例虽然都不一致，但是地形曲面模拟的准确度在一定程度上能够反映地表面积计算的可靠性。

2 顾及局部坡度因子与极值的地表面积计算

2.1 规格格网地表面积计算

规格格网地表面积计算常采用同向对角划线或交叉对角划线将空间四边形转换为两空间三角形。该方法假设条件是局部地形变化在空间上是连续的、相关的，地形曲面的极值条件与对角线方向是一致的，对角连线方向贴近地面。对角划线法在一定范围内虽具有一定的补偿性，然而实际地形具有的不规则性、非线性特点，在地形变化复杂的山区尤为突出，该方法也可能与地形变化不一致，造成局部地形逼近失真。

hi(x,y)=a1x3y3+a2x2y3+a3xy3+a4y3+a5x3y2+a6x2y2+

a7xy2+a8y2+a9x3y+a10x2y+a11xy+a12y+

a13x3+a14x2+a15x+a16

(7)

2.2 边界破碎格网表面积计算

边界破碎格网表面积计算常采用坡度比例法反推地表面积，坡度计算适宜的窗口选择、坡度参数的有效范围都与地形复杂度紧密相关，一般以3×3的局部窗口计算。不规则三角网(TIN)具有多分辨率、贴近局部地形、精度高的特征，多用于三维地形仿真上[15-16]。建立不规则三角网的方法有多种，如Green-Sibson、角度判别法、Delaunay法等，其中Delaunay三角划分的理论更完善。相关研究表明，Delaunay三角网拟合性能最优越，具有从局部到整体最优的特点。因此，本文采用Delaunay方法构建不规则三角网近似表达边界破碎格网地表形态。

2.3 三角形面积误差

(8)

(9)

依据协方差传播规律可求得空间边长L的误差为

(10)

假设空间边长误差中，高程误差起主要作用，忽略x、y误差，则

(11)

(12)

根据泰勒级数与协方差传播规律可推导出空间三角形面积误差(见式(12)，限于文中篇幅要求，不再详细推导)

(13)

依据式(4)、式(11)、式(13)可初步得出顾及坡度因子与图形极值信息的地表面积统计结果精度与可靠性较高。

3 算例分析

使用我国西部具有典型地貌特征的K县、S县的边界数据与10 m分辨率DEM数据，通过与其他常规方法对比，采用4种方案计算两区域地表面积。文中分析了顾及坡度因子与极值的地表面积计算方法的可靠性与差异性。其中，K县地形复杂，地貌以盆地、山地、平原为主，县区内水域较多；S县地形较简单，以沙漠、平原为主。

方案1：采用顾及坡度因子与极值的地表面积计算法；

方案2：构建地表三角网法，规格格网对角线西北方向划分，边界破碎格网采用Delaunay构TIN；

方案3：构建地表三角网法，规格格网对角线东北方向划分，边界破碎格网采用Delaunay构TIN；

方案4：地表三角网与坡度法组合方法，规格格网依次交叉对角划分，边界破碎格网采用坡度比例法计算。

表1 4种方案表面积差异

表1给出了4种方案计算的两县地表面积差异。从统计结果可以看出：4种方案计算的结果在地形复杂区域K方案1与方案4差异较小、符合性好，方案2、方案3与方案1差异较大。在地形相对平缓的区域S，4种方案计算结果差异不大。从地表面积对比结果可得出：规格格网的表面积计算方法的差异占主要因素，破碎格网表面积方法差异影响较小。

本文提出的顾及坡度因子与图形极值信息的地表面积计算法图形结构上与交叉对角划线法相比，宏观上二者具有一定的相似度，存在尺度上细小的差异；微观上文中提出的方法在地形变化为非线性，考虑局部地形极值出现的区间，对局部地形逼近具有一定的改进，可削弱同向对角划线构建的图形二异性与不合理性。

4 结束语

地表起伏的不确定性与复杂性，影响地表面积统计结果的可靠性。本文采用实际算例进行比对试验，试验结果证明笔者提出的顾及坡度因子与极值的地表面积计算方法可以获取较可靠的地表面积统计结果。另外，文中推导的三角网地表面积的传递误差，理论上也初步得出顾及坡度因子与局部图形极值信息的地表面积计算法对地形非线性变化复杂的区域具有一定的改善作用。

[1] 王耀革.DEM建模与不确定性分析[D].郑州：信息工程大学,2009.

[2] 赵杰.数字地形模拟-地形数据获取与数字地形分析研究[D].武汉：武汉大学,2004.

[3] 刘学军,龚健雅,周启鸣,等.DEM结构特征对坡度坡向的影响分析[J].地理与地理信息科学,2004, 20(6):1-5.

[4] 刘学军,张平.DEM坡度、坡向的有效尺度范围[J].武汉大学学报(信息科学版),2008, 33(12):1254-1258.

[5] 刘学军,张平,朱莹. DEM坡度计算的适宜窗口分析[J].测绘学报,2009, 38(3):264-271.

[6] 陈艳丽,李少梅,刘岱岳.基于规则格网DEM的坡度坡向分析研究[J].测绘与空间地理信息, 2009,32(5):36-39.

[7] 杨学祥.地球表面积的计算[J].长春地质学院学报,1987, 17(3):346-352.

[8] 纪兵,边少锋.三轴椭球体表面积的计算[J].海军工程大学学报, 2008,20(4):21-24.

[9] 江帆,吕晓华,王仲兰.基于复化公式的DEM表面积算法分析[J].测绘学院学报,2005, 22(4):263-265.

[10] 陈吉龙,武伟,刘洪斌.DEM在林地表面积计算中的应用研究[J].西南农业学报,2008, 21(5):1348-1352.

[11] 曹振宇,谭明建,沈艳敏.大地坐标系中多边形面积计算方法研究[J].测绘通报,2012(S1)：570-572.

[12] 党亚民,程英燕,吴秀娟,等.不同坐标系图斑理论面积计算研究[J].测绘科学,2005, 30(6):23-25.

[13] 薛树强,党亚民,秘金钟,等.顾及非线性地形因子的地表面积计算[J].测绘学报,2015, 44(3):330-337.

[14] 秦艳红,卢林,吴纪桃.地形坡度面积的二维多分辨率分析研究[J].数学理论与应用,2005, 25(2):43-48.

[15] 蒋光伟,焦立芬,田晓静,等.基于三角剖分与方向线拼接的高精度坐标转换格网模型的建立[J].大地测量与地球动力学,2015, 35(2):222-224.

[16] 徐青.地形三维可视化技术[M].北京:测绘出版社,2000.

AMethodofTerrainAreaCalculationinConsideringofSlopeFactorsandExtremeValue

CHENG Chuanlu,JIANG Guangwei,TIAN Xiaojing,MA Xinying,WANG Wenli

(Geodetic Data Processing Center, National Administration of Surveying,Mapping and Geoinformation, Xi’an 710054, China)

In connection with the effection of the local terrain slope factors and extreme value to terrain area calculation, a new method that considering of slope factors and extreme value is proposed. The surface area theory and errors are analyzed. Finally the typical two region’s data of west China are used experiment analysis. The results indicate that the algorithm enhance the terrain fitting, on the macro, the graphic property is similar with the regular grid of cross-diagonal method. The deduce area errors formula indicated that the terrain surface area calculation is more reliable in considering the slope factors and extremevalue.

slope; terrain area calculation; triangulation network; DEM; patch

P208

A

0494-0911(2017)01-0126-04

程传录,蒋光伟,田晓静,等.一种顾及坡度因子与极值的地表面积计算法[J].测绘通报，2017(1)：126-129.

10.13474/j.cnki.11-2246.2017.0028.

2016-02-18

2015年测绘地理信息公益性行业科研专项(201512004)；国家测绘地理信息局基础测绘科技项目；现代工程测量国家测绘地理信息局重点实验室开放课题(TJES1504)

程传录(1965—),男,高级工程师,主要从事大地测量数据处理研究。E-mail: xianccl@163.com