基于DEM的小班坡度自动提取算法及其验证

陈 晨 陈永刚 徐文兵 梁 丹

( 1. 浙江农林大学浙江省森林生态系统碳循环与固碳减排重点实验室，浙江 杭州 311300；2. 同济大学测绘与地理信息学院，上海 200092)

地面坡度（简称坡度）是对地面倾斜程度的定量描述，是一个基本的地貌形态指标[1]，其大小直接影响着地表物质流动和能量转换的规模和强度，是工程设计和生产力空间布局的重要因子[2]，也是研究水文模型[3]、土壤侵蚀[4]、监测滑坡[5]和规划土地利用[6]的基础。小班坡度是森林资源规划设计调查（二类调查）中的一项重要因子，它是反映小班地表形态、小班立地质量、生产作业条件及水土流失的重要指标之一，也是林木采伐首先要考虑的因素之一[7]。在森林资源连续清查的小班调查中，小班坡度是重要的调查因子。因此，快速准确的提取小班坡度显得至关重要。

在以往的二类调查中，提取小班坡度的传统方法主要有两种：一是外业实地测量小班坡度。该方法主要是在实地进行小班区划后，利用坡度仪或者手持罗盘仪实地测量小班的平均坡度；二是利用地形图进行室内估算，即在小班内两处等高线的水平距离及垂直距离，根据地形图提供的坡度尺，计算小班的坡度值。利用传统方法提取的二类小班坡度虽然很直观，但主要依靠手工操作完成，导致工作效率较低，工作量大，时间周期长等问题，并且在野外实施也有一定的难度；另外，由于依靠调查人员的经验判断，导致结果具有一定的主观性及随意性，受人为因素影响大，数据不准确，难以满足二类调查的要求。研究人员也提出了一些自动化提取小班坡度的方法。陈建义等[7]针对传统方法的不足，探讨了一种基于ArcGIS软件的小班坡度自动求算方法。熊安华[8]、叶江霞等[9]均利用ArcGIS进行小班坡度的计算。这些基于数字高程模型（DEM）并结合GIS分析工具的方法能够高效方便地计算出小班的坡度，但是这些方法只计算了小班内所有栅格点坡度的平均数或者众数。这种以简单数学统计结果作为整个图斑坡度表达的方法忽略了地形因素，并不能作为小班内实际坡度的真实反映，这种缺陷在地形复杂度较高的地区尤为明显[10]。

针对林业小班坡度计算的现状，本研究提出了一种基于DEM的小班坡度自动提取算法。该算法先将小班内真实地形中存在的异常高低起伏点进行剔除，然后将小班整体作为计算对象，对剔除异常点的DEM进行回归，得到小班所在的斜平面模型，最后根据坡度公式计算该小班的坡度值。该算法能够利用DEM快速准确地提取小班的坡度，并且能够适应地形复杂度较高的区域，实现代替人工实测提取小班坡度。

1 基于DEM的小班坡度自动提取的算法

1.1 算法思路

在计算小班坡度值时，理想的小班内部地形表面是平整光滑的，然而在真实地形中，小班内的地形表面总是存在着一些高低起伏，见图1～2。在计算小班的坡度时，这些异常的高低起伏在一定程度上影响小班坡度计算的准确度。为了更加准确地计算小班的坡度，需要将这些高低起伏的地方进行剔除。因此，本研究提出了一种剔除异常高低起伏点的小班坡度自动提取算法（FSS），其算法流程见图3。采用线性回归获得小班所在的斜平面模型，从而能够综合考虑小班整体情况。

图 1 理想小班内部地形表面Fig. 1 Ideal subcompartment internal terrain surface

图 2 山坡上异常高低起伏点Fig. 2 Abnormal high and low undulations on the slope

图 3 算法流程Fig. 3 The flow of experiment

1.2 识别并剔除异常高低起伏点

首先，借助ArcGIS的综合工具，使用3×3的邻域窗口，将DEM进行均值综合。然后，使用双线性插值的方法将综合后的DEM重采样至原始DEM的分辨率。双线性插值方法基于4个最近输入单元中心的加权距离平均来确定单元新值。最后，将重采样后的DEM与原始DEM作差，由此得到原始DEM与重采样DEM之间的高程差。对DEM高程差进行分析统计，将超过2倍方差的点视为异常的高低起伏点。在本研究中，小班内的DEM栅格数量均大于30，假设服从正态分布。剔除异常高低起伏点的公式如下：

式中：X为小班内DEM的高程差，为小班内高程差的算术平均数，为小班下高程差的标准差。

1.3 建立小班斜平面模型

剔除异常高低起伏点后，利用线性回归来建立小班所在斜平面的回归模型，回归模型的方程如下：

式中：Z是斜平面模型的因变量，此处为小班内所有点的高程值集合，X1和X2是斜平面模型中的自变量，此处分别为小班内所有点的地理坐标X值和Y值集合，n是单个小班内DEM栅格点的个数，zi是第i个点的高程值，xi和yi则分别表示第i个点的地理坐标的X值和Y值。b1和b2是X1和X2的系数，ε是误差值，a 是常数项。

在算得每个小班的斜平面模型后，根据坡度计算公式得到小班的坡度。坡度计算公式如下：

式中：fx为南北方向高程变化率，fy为东西方向高程变化率。此处的fx和fy分别等于式（2）中的 b1和 b2。

2 检验数据来源及方法

2.1 数据来源

采用宁海县2007年二类调查成果数据，该数据为shapefile格式，坐标系为1980年国家大地坐标系，该数据作为本研究的实测数据。选取68个小班作为样本，这68个小班均未出现跨山脊和山谷的现象，具体小班分布见图4。DEM数据是通过对1∶50 000的等高线数据进行加工而来，同样为1980年国家大地坐标系。DEM由国家测绘部门提供。用等高线数据生成不规则三角网（TIN）；然后将TIN进行栅格化，采用自然邻域法，从而生成像元大小为10 m×10 m的DEM数据；最后，对DEM数据进行降噪，填补洼地，以提高DEM数据的精度。

图 4 实验区Fig. 4 Study area

2.2 检验方法

检验主要侧重于算法准确性和实用性的检验。准确性检验就是将FSS算法与陈建义等[7]的均值法以及常用的中心点法（将小班中心点的坡度值作为小班的坡度值）计算小班坡度并进行对比验证。由于无法获得小班的真实坡度值，本研究利用模拟斜平面进行验证，通过向斜平面上添加随机噪声来模拟现实世界中地表的起伏。若计算所得结果更加接近预设值，则认为该方法得到的小班坡度更加准确，能够适应复杂度较高的区域。为测试FSS算法的实用性，采用森林资源二类调查数据与本算法的结果进行对比。验证方法均为t检验。

为了进行算法检验，本研究首先构造了一个30°的斜平面和一组随机噪声。其次，为模拟真实地形，向斜平面中分别加入随机噪声，噪声分别为指数噪声、伽玛噪声、整数噪声、正常噪声、泊松噪声和均匀噪声（图5）。模拟斜平面和噪声的像元大小与实验区DEM数据相同。最后，在模拟斜平面上绘制1块图斑，用于模拟小班，保存成shapefile格式。

图 5 增加噪声的模拟地形Fig. 5 Simulated terrain with increased noise

3 检验结果与分析

3.1 准确性检验

使用均值法、中心点法和FSS算法分别计算模拟斜面上小班的坡度。在95%置信区间内，通过配对t检验对这3组数据进行检验。其中，原假设和备择假设是：H0: abs_A≥abs_B，H1: abs_A＜abs_B。其中，abs_A和abs_B分别表示数据A和数据B与真值差的绝对值。因为本研究需要比较FSS算法和其他2种方法计算结果的准确度，所以此处的abs_A和abs_B有以下2种情况：|FSS-真值|和|均值法-真值|、|FSS-真值|和|中心点法-真值|。

配对t检验的FSS算法与均值法和中心的检验结果见表1。由表1可知，P值均小于0.05，即拒绝原假设，应用FSS方法所得坡度值的均值更接近于真值。由表2可知，FSS方法提取的小班坡度与真值偏差的均值、标准差、最大值和最小值均比其他2种方法的结果小。FSS方法能够准确地提取小班坡度，同时面对实验预设的复杂地形，其相较于其他方法有更好的适应性。

表 1 配对t检验的结果比较Table 1 Comparison of the results of paired t-test

表 2 各方法提取小班坡度结果与真值的偏差Table 2 Each method extracts the deviation of the subcompartment slope result from the true value

3.2 实用性检验

使用FSS算法对实验区数据进行试验，得到了对应小班的坡度值。用FSS算法计算得到的值与人工实测值进行配对t检验，其中，原假设和备择假设是：H0: data_A=data_B，H1: data_A≠data_B。其中，data_A和data_B分别表示FSS算法结果和实测值数据。由表3可知，在95%的置信区间内，这2组数据的P值为0.281，P值大于0.05，接收原假设，即FSS算法结果与人工实测数据无明显差异。

表 3 实测数据与FSS算法结果配对t检验的结果比较Table 3 Comparison of the results of paired t-test between measured data and FSS algorithm results

4 结论与讨论

针对现有计算小班坡度算法的不足，本研究提出了一种基于DEM的小班坡度提取算法。该算法能够准确的提取小班坡度值，并且其结果与人工实测值无区别，能够替代人工野外实测。FSS算法比常用的算法更加合理，原因可以从以下三方面进行解释：

1）小班内的地形并不是光滑平坦的，坡面上存在着异常的高低起伏点。这些高低起伏点会影响小班坡度计算的准确性。本研究的算法能够剔除这些异常高低起伏点。

2）FSS方法利用线性回归获取了小班所在斜平面模型，综合地考虑小班整体情况，而非通过简单的统计计算得到小班的坡度值。因此，本研究所提出的算法结果比常用的算法更加准确。

3）本研究所提出的算法是一种自动化的提取方法，对于野外人工调查具有一定的参考价值。可以节省人力、物力和时间，给森林资源调查提供了便利。

实际上，小班坡度算法准确性检验是一项很具有挑战性的研究。由于无法获得小班坡度的真实值，现有条件得到的小班坡度值均为测量值或其他方法所得，这些类型的数据并不准确。在相关的基于DEM的坡度算法研究中，格网的坡度值通常使用模拟曲面来检验，以此来评估算法的合理性，例如Chen等[11]，Florinsky[12]，Zhou等[13]为了检验DEM提取坡度值的精度，在模拟曲面上进行算法比较。借鉴以上研究的思路，本研究构建了一种适合于小班坡度算法的模拟地形曲面，与之前不同的是该曲面是一个添加了噪声的斜平面，虽然它无法完全模拟真实地形，但是其坡度值可以作为一个理想真值。在异常高起伏点的剔除过程中，本研究选用了最为常用的3×3邻域窗口，并且该窗口在实验中达到了预期效果。在实际应用研究中，窗口大小的选择受地形和DEM分辨率等因素影响[14-15]。因此，在今后研究中，将对邻域窗口的选择与DEM分辨率、异常区域大小等因素进行深入探讨。