一种确定地形起伏度最佳分析区域的方法

摘要

本发明涉及一种确定地形起伏度最佳分析区域的方法，用以解决最佳分析区域确定过程中，由于逐次计算的复杂性和人工判别的随意性等问题。本发明的方法充分利用分析区域边长和最大地形起伏度值的对数函数关系，运用逐步回归分析方法分析对数函数因子和实验区域平均高程、高差、平均坡度、平均坡度变率之间的相关关系，建立了最大地形起伏度值和实验区域平均高程、高差、平均坡度、平均坡度变率的函数关系，实现了地形起伏度最佳分析区域范围的确定。该方法不仅可以用于地形起伏度最佳分析区域的确定，而且可以用于其他宏观地形因子最佳分析区域的确定。

说明书

技术领域

本发明涉及一种确定地形起伏度最佳分析区域的方法。

背景技术

地形起伏度指分析区域内最高点和最低点之差，反映宏观区域内地形的起伏特征，是定量描述地形形态、划分地貌类型的重要指标，计算方法如式（1）所示：

R＝H_max-H_min                                               (1)

其中，R为分析区域内的地形起伏度值，H_max、H_min为分析区域内的最高和最低高程值。

地形起伏度提取的关键在于确定分析区域内的最高和最低高程值。但是随着分析区域范围的变化，最高和最低高程值随之发生变化，导致地形起伏度值也发生变化。因此，确定最佳分析区域是地形起伏度提取算法的核心步骤和决定地形起伏度提取结果有效性的关键（参见朱红春，陈楠，刘海英，汤国安.自1：10000比例尺DEM提取地形起伏度——以陕北黄土高原的实验为例[J].测绘科学,2005,30(4)：86-88.）。

按照地貌发育的基本理论，最佳分析区域指存在使高差达到相对稳定的区域范围。对于相同区域的DEM而言，随着分析区域面积的增大，区域内最高点与最低点的高差增大。通常开始时高差以较快的速度增加，然后增速逐渐变缓。当分析区域面积达到某一阈值后，高差的变化基本稳定不变。高差变化率由快变缓的点即为拐点，拐点对应的面积就是最佳分析区域（参见涂汉明，刘振东.中国地势起伏度研究[J].测绘学报，1991,20(4)：311-319.）。

涂汉明等（参见涂汉明，刘振东.中国地势起伏度最佳统计单元的求证[J].湖北大学学报(自然科学版)，1990,12(3)：266-271.）认为中国存在2km²、6km²、16km²、20km²、22km²五种不同规模的地形起伏度最佳统计分析区域，且各有其适用范围；并指出21km²的分析区域具有全国普适性。刘新华等（参见刘新华，杨勤科，汤国安.中国地形起伏度的提取及在水土流失定量评价中的应用[J].水土保持通报,2001,21(1)：57-62．）使用1∶100万DEM数据提取地形起伏度研究水土流失状况时，认为考虑到样区均衡性因素，地形起伏度的最佳分析区域为25km²。郎玲玲等（参见郎玲玲，程维明，朱启疆，龙恩.多尺度DEM提取地势起伏度的对比分析——以福建低山丘陵区为例[J].地球信息科学,2007,9(6)：1-6.）使用福建地区1∶25万DEM和1∶10万DEM比较多尺度DEM提取的地形起伏度差异时，认为该地区的地形起伏度最佳分析区域为4.41km²和0.4km²。唐飞等（参见唐飞，陈曦，程维明.基于DEM的准噶尔盆地及其西北山区地势起伏度研究[J]．干旱区地理,2006,29(3)：388-392.）通过对准噶尔盆地及其西北山区地形起伏度的研究，得出该地区地形起伏度的最佳分析区域为4km²。刘爱利（参见刘爱利.基于1：100万DEM的我国地形地貌特征研究[D].西安：西北大学,2004.）基于1∶100万DEM数据研究我国的地形地貌特征，认为我国地形起伏度的最佳分析区域为42.25km²。张磊（参见张磊.基于地形起伏度的地貌形态划分研究——以京津冀地区为例[D].河北：河北师范大学,2009.）在使用1：25万DEM数据研究京津冀地区地形起伏度时，认为该地区的最佳分析区域为9.61km²。王雷等（参见王雷，朱杰勇，周雁.基于1∶25万DEM昆明地区地貌形态特征分析[J]．昆明理工大学学报(理工版)，2007,32(1)：6-14.）使用1∶25万地形图数据研究昆明地区地貌形态时，认为该地区的最佳分析区域为16km²。王玲等（参见王玲，同小娟.基于变点分析的地形起伏度研究[J].地理与地理信息科学,2007,23(6)：65-67.）使用1∶25万DEM数据研究新疆地区地形起伏度时，认为该地区的椎间分析区域为2.56km²。

最佳分析区域就是高差变率由快变缓的拐点对应的区域面积。拐点的判断通常包括人工判别法、最大高差法、最大高差-面积比法以及预测法等。

人工判别法是绘制“面积-最大地形起伏度”或“边长-最大地形起伏度”曲线，人工判断拐点的位置，这种方法主观性较强。如图1所示，假设高差变率变缓的阈值为30度，则B点可以为拐点，对应最佳分析区域边长为1000米。

最大高差法认为相邻分析尺度之间地形起伏度变化最大时，拐点出现。其计算公式为：

ΔH_i＝R_i-R_i-1                                              （2）

其中，R_i为第i个尺度下的地形起伏度值，ΔH_i为第i个尺度下相邻地形起伏度值之差。取R_i最大对应的第i点即为拐点。如图2所示，最大的高差点出现在1500米左右，所以最佳分析区域边长为1500米。

最大高差-面积比法是计算相邻尺度地形起伏度之差与邻域面积差的比值，得I_i，从所有I_i中选择最大值，其所对应的第i点即为拐点，其计算公式如式（3）所示。

$>I_i=\frac{R_{i+1}-R_i}{((S_{i+1}-S_i)/S_{i+1})\times 100}---\left(3\right)$>

其中，R_i为第i个尺度下的地形起伏度值，S_i为第i个尺度下的邻域面积值。如图3所示，最大的I值出现在1500米处，所以最佳分析区域边长为1500米。

预测法是根据实验区域高差，建立地形起伏度和分析区域面积的对数函数关系，从而实现最佳分析区域的计算（参见张锦明，游雄.地形起伏度最佳分析区域研究[J]．测绘科学技术学报,2011,28(5)：369-373.）。

上述分析表明：

①不同学者的研究成果表明，不同实验区域对应不同的地形起伏度最佳分析区域，不同尺度的DEM数据对应不同的地形起伏度最佳分析区域，分析区域面积（或分析区域边长）和最大地形起伏度之间存在对数函数关系。

②使用人工判别法得到的最佳分析区域存在较大的人为因素，导致与最大高差法、最大高差-面积比法存在很大的出入。

③使用参数较少，导致最佳分析区域的计算并不准确。

发明内容

本发明的目的是提供一种确定地形起伏度最佳分析区域的方法，用以解决最佳分析区域确定过程中，由于逐次计算的复杂性和人工判别的随意性等造成的不能准确确定的问题。

为实现上述目的，本发明的方案是：一种确定地形起伏度最佳分析区域的方法，步骤如下：

1）对于一个实验区域，根据$>\left(\begin{array}{c}a=25.38-0.01*\bar{H}+0.16*\mathrm{\Delta H}+16.30*\bar{S}-41.36*S^{\prime }\\ b=-0.23-0.06*\bar{H}+0.89*\mathrm{\Delta H}+61.82*\bar{S}-220.23*S^{\prime }\end{array}\right),$>得到因子a、b；为平均高程，ΔH为区域高差，为平均坡度，S′为平均坡度变率；

2）基于对数函数关系y＝a*lnx-b，y为最大地形起伏度值，x为分析区域边长，对x进行间隔取样，并分别求取y＝a*lnx-b的斜率k，k＝a/x；

3）计算相邻斜率差Δk，当相邻斜率差Δk小于设定阈值时，分析区域边长x为最佳分析区域边长，所述斜率差Δk的阈值为0.001。

本发明实际上提出了一种基于微观地形特征因子的最佳分析区域预测模型，不同于现有技术，不使用人工判据法，避免了人为因素的干扰，开拓性地使用了大量微观地形特征因子，并借助分析区域面积（或分析区域边长）和最大地形起伏度之间的对数函数关系确定最佳分析区域，区域确定更加准确。

本发明的方法充分利用分析区域边长和最大地形起伏度值的对数函数关系，运用逐步回归分析方法分析对数函数因子和实验区域平均高程、高差、平均坡度、平均坡度变率之间的相关关系，建立了最大地形起伏度值和实验区域平均高程、高差、平均坡度、平均坡度变率的函数关系，实现了地形起伏度最佳分析区域范围的确定。本发明不仅可以用于地形起伏度最佳分析区域的确定，而且可以用于其他宏观地形因子最佳分析区域的确定。

附图说明

图1为地形起伏度拐点确定的人工判别法示意图；

图2为地形起伏度拐点确定的最大高差判别法示意图；

图3为地形起伏度拐点确定的最大高差-面积比法示意图；

图4为实验区域分布图；

图5(a)、(b)、(c)为陕西实验区域示意图；

图6为地形起伏度计算的边界效应处理；

图7为陕西地区各尺度DEM计算得到的地形起伏度值；

图8为陕西地区各尺度DEM计算得到的地形起伏度值(局部放大图)；

图9为河南地区DEM计算得到的地形起伏度值；

图10为黑龙江地区DEM提取得到的地形起伏度值；

图11为内蒙西部地区DEM计算得到的地形起伏度值；

图12为对数函数中参数a和高差的散点图；

图13为对数函数中参数b和高差的散点图；

图14为对数函数中参数a拟合值的置信区间图；

图15为对数函数中参数b拟合值的置信区间图。

具体实施方式

下面结合附图对本发明做进一步详细的说明。

为了具体说明本发明，下面介绍基于微观地形特征因子的最佳分析区域预测模型的获得过程。

一、实验设计

（一）实验数据来源与特征

从全国范围内随机选取了26个地区的30km×30km、60km×60km、90km×90km三个不同范围的ASTER GDEM作为实验数据（如图4所示）（各实验区域的面积总和为21.06万km²，约占全国总面积的2%），并且建立格网尺寸分别为30m、60m和300m的多尺度DEM数据（分别如图5（a）、（b）、（c）所示），最后计算各实验区域的微观地形特征因子，包括最低高程、最高高程、区域高差、平均坡度、平均坡度变率、平均平面曲率和平均剖面曲率。

（二）地形起伏度计算方法

地形起伏度的计算通常运用邻域分析原理实现。给定DEM中任一格网点，开僻一定大小的分析区域，在该区域内搜索最高和最低高程值，得到该格网点的地形起伏度值；然后分别提取DEM中每一格网点的地形起伏度值，最终得到实验区域的最大地形起伏度值。

分析区域形状可以分为正方形邻域、圆形邻域、环形邻域或扇形邻域等。本文使用正方形邻域计算地形起伏度，并且分别计算分析区域尺度从“3×3”、“5×5”、……、“199×199”的系列最大地形起伏度值。

同时，在计算某一分析区域尺度的地形起伏度时考虑边界效应，删除不能构成完整分析区域的格网点。以图6为例，用“3×3”的分析区域尺度计算格网点A的地形起伏度值，由于A位于DEM边缘，1、2、3、7、8为无值格网点。此时格网点A的地形起伏度值存在误差。因此，对于类似含有边界效应的格网点，本文不计入最终的地形起伏度计算。

（三）实验流程

为了研究和地形起伏度最佳分析区域相关的问题，设计如下实验步骤：

①随机选择实验区域，并计算各实验区域的微观地形特征因子；

②对实验区域进行分析区域尺度从“3×3”、“5×5”、……、“199×199”的系列最大地形起伏度值的计算；

③研究相同实验区域、多尺度DEM数据提取的地形起伏度值之间的关系；

④研究地形起伏度值和实验区域微观地形特征因子之间的关系，建立并验证最佳分析区域预测模型。

二、实验分析

（一）多尺度DEM和地形起伏度的关系

研究多尺度DEM和地形起伏度之间的关系，即在同一实验区域中，从不同尺度DEM中提取的地形起伏度值是否存在差异？差异多大？

实验表明相同分析区域对应类似地形起伏度值。以陕西地区90km×90km实验区域为例（其他实验区域具有类似的结果），当分析区域为9km2时，从30mDEM提取的地形起伏度值为691m（表1），从60mDEM提取的地形起伏度值为689m（表2），从300mDEM提取的地形起伏度值为648m（表3）。

表1陕西地区30m DEM提取的地形起伏度值

表2陕西地区60m DEM提取的地形起伏度值

表3陕西地区300m DEM提取的地形起伏度值

建立“分析区域边长─最大地形起伏度”散点图（如图7、8所示），可以发现相同分析区域提取的最大地形起伏度值相似，随着DEM尺度差异的增加，相同分析区域计算得到的最大地形起伏度值差异增加。

并结合多尺度DEM的微观地形特征因子，可以发现对于同一实验区域而言，随着DEM尺度的增大，最大高程逐渐减小、最小高程逐渐增大、高程的分布趋于均化，表现在试验区域高差逐渐减小，平均坡度逐渐减小，地形起伏度逐渐减小。因此可以初步判定地形起伏度和最大高程、区域高差、平均坡度等微观地形特征因子存在某种程度的相关性。

（二）微观地形特征因子和地形起伏度的关系

已有地形起伏度研究成果表明，分析区域面积（或分析区域边长）和最大地形起伏度值之间存在对数函数关系（参见刘爱利.基于1：100万DEM的我国地形地貌特征研究[D].西安:西北大学,2004；张磊.基于地形起伏度的地貌形态划分研究——以京津冀地区为例[D].河北:河北师范大学,2009；王雷,朱杰勇,周雁.基于1:25万DEM昆明地区地貌形态特征分析[J]．昆明理工大学学报(理工版),2007,32(1):6-14；王玲,同小娟.基于变点分析的地形起伏度研究[J].地理与地理信息科学,2007,23(6):65-67.），形如：

y＝a*lnx-b                                               （4）

其中，x分析区域面积或边长，y为实验区域最大地形起伏度值。

根据以上研究，本发明进一步提出，地形起伏度最佳分析区域计算的最佳方法是基于对数函数关系通过相邻斜率差的计算得到最佳分析区域。

前面提及地形起伏度和最大高程、区域高差、平均坡度可能存在具有某种程度的相关性，那么式（4）中的未知参数a、b能否由最大高程、区域高差和平均坡度表示？

首先对各实验区域的最大地形起伏度值和分析区域边长做对数函数的曲线拟合，其拟合参数结果如表4所示（由于不同尺度DEM的拟合参数类似，因此本文仅列出尺度为30m的DEM的实验结果）。

表4实验区域最大地形起伏度与分析区域边长对数拟合参数

在拟合过程中发现两个问题：

一是地形起伏度值在某一小分析区域内达到稳定，而后随着分析区域的增大，地形起伏度值迅速增加，最终导致拟合效果极差。如60km×60km的河南实验区域（如图9所示），从图中可以发现，当分析区域边长在240m∽1000m时，最大地形起伏度值达到相对稳定；但是随着分析区域继续增大，最大地形起伏度值同样继续增大，表明240m∽1000m的最佳分析区域不符合山体完整性原则。如果拟合过程中删除0m∽1000m的数据，则曲线拟合效果变好（R²=0.9955）。

二是由于实验区域的平均坡度较小（小于3度）时，提取得到的地形起伏度值呈现明显的阶梯状或在很小分析区域范围内即达到稳定，并且随着分析区域范围的增加，地形起伏度值变化不大，同样导致拟合效果极差。如60km×60km的黑龙江实验区域（如图10所示）和60km×60km的内蒙西部实验区域（如图11所示）。

通常认为“地形起伏度是用于划分丘陵、小起伏山地、中起伏山地和大起伏山地的依据”，因此本文在实验过程中删除坡度小于3度的实验区域（实验表明坡度小于3度的实验区域的“对数函数关系”拟合效果都较差）。

在剔除或修正了拟合过程中的两种特殊情况之后，运用相关分析法分析未知参数a、b和地形特征因子的相关关系（表5），可以发现:

a和区域高差存在强相关性，相关系数达到0.9389（单独建立a和区域高差的散点图，如图12所示）；和最高高程、平均高程、平均坡度、平均坡度变率存在一定的相关性。b和区域高差存在强相关性，相关系数达到0.9136（单独建立b和区域高差的散点图，如图13所示）；和最高高程、平均高程、平均坡度、平均坡度变率存在一定的相关性。

5拟合参数a、b和微观地形特征因子相关系数表

从可决系数角度来看，使用区域高差拟合未知参数a、b不能达到理想效果。考虑到最低高程、最高高程、平均高程、平均坡度、平均坡度变率对未知参数a、b均有一定程度的影响。于是运用逐步回归分析法建立a、b和区域高差、最低高程、最高高程、平均高程、平均坡度、平均坡度变率的多元线性回归函数。可以发现区域高差、平均高程、平均坡度、平均坡度变率可以参与对a、b值的拟合，即当置信水平为0.05时，它们对a、b值具有显著性影响。其表达式（5）所示：

$>\left(\begin{array}{cc}a=25.38-0.01\stackrel{}{*\bar{H}+0.16*\mathrm{\Delta H}+}& \\ 16.30*\bar{S}-41.36*S^{\prime }& (R^2=0.9267)\\ b=-0.23-0.06*\bar{H}+0.89*\mathrm{\Delta H}+& \\ 61.82*\bar{S}-220.23*S^{\prime }& (R^2=0.8618)\end{array}\right)---\left(5\right)$>

其中，为平均高程，ΔH为区域高差，为平均坡度，S′为平均坡度变率。

最后利用式（5）拟合所有实验区域的参数a、b，并建立置信水平为0.05时的拟合值置信区间图（如图14、15所示）。可以发现：对于参数a来说，仅仅存在90km×90km的新疆北部实验区域地区不能满足要求，对于参数b来说，存在30km×30km、60km×60km的新疆南部，30km×30km的青海北部三个实验区域不能满足要求。

从总体上看，式（5）较单独使用区域高差拟合未知参数a、b略有提高；从拟合残差图看，95%以上的拟合值都能达到置信水平的要求。同时由于在实际计算过程中，我们需要计算拟合曲线的斜率，但是斜率仅仅和a发生关系。因此，可以将式（4）和（5）确定为基于微观地形特征因子的最佳分析区域预测模型。

具体实施例

以下以全国90m SRTM DEM数据为例，并结合附图对本发明的具体实施方式作进一步详细说明。

第一步：数据准备

（1.1）选择全国90m SRTM DEM数据作为实验数据，并分别计算区域高差、平均高程、平均坡度、平均坡度变率等微观地形因子值（表6）；

表6全国90m SRTM DEM数据的地形特征因子值

（1.2）根据实验区域的相关微观地形因子，根据（4）、（5）式，得到a=1195.9550、b=6390.4439。

第二步：计算斜率差，并确定最佳分析区域

（2.1）对x进行间隔取样，分别计算y＝a*lnx-b的斜率，即k＝a/x；

（2.2）计算相邻斜率之差；

（2.3）当斜率之差小于0.001时，x≈11000m，即本地区的最佳分析区域边长为11km。

我们可以用人工判别法验证a、b的计算值

对实验区域进行分析区域尺度从“3×3”、“5×5”、......、“199×199”的系列最大地形起伏度值的计算；

绘制“分析区域边长-最大地形起伏度值”的散点图，并拟合对数函数曲线，拟合得到a=1167（在0.05的置信水平时置信区间为[1070,1263]）、b=6437（在0.05的置信水平时置信区间为[5458,7415]）；

实验验证表明由式（3）拟合得到的参数a、b都位于0.05的置信区间内。

本发明的优点在于：

（1）本发明充分利用分析区域边长和最大地形起伏度值的对数函数关系，运用逐步回归分析方法分析对数函数因子和实验区域平均高程、高差、平均坡度、平均坡度变率之间的相关关系，建立了最大地形起伏度值和实验区域平均高程、高差、平均坡度、平均坡度变率的函数关系，实现了地形起伏度最佳分析区域范围的确定，确定的最佳分析区域具有很高的可信度。

（2）本发明不仅可以用于地形起伏度最佳分析区域的确定，而且可以用于其他宏观地形因子最佳分析区域的确定。