基于多源DEM数据的都江堰灌区地形因子提取及对比分析

肖文全,曹依帆,秦 涛,罗 尚,翟 星,冯丹禹,赵 辉

(四川水发勘测设计研究有限公司,成都 610072)

0 引言

20世纪50年代,数字高程模型(DEM)的概念被提出,它被定义为一种以数字数组形式显示地面标高的实体模型[1]。近年来,随着计算机技术、地理信息技术(GIS)和遥感技术(RS)的迅速发展,DEM数据的种类逐渐增多,分辨率也逐步提高。基于DEM数据提取的地形因子可以有效地描述地貌的形态特征,从而推动地貌形态分类的研究进展[2]。因此,DEM数据在众多基础地理信息数据中占有极为重要的地位。针对DEM的研究也成为当下的热点,高艺伟等[3]基于ASTER GDEM2数据提取陕南地区地形起伏度,并采用均值变点法确定最佳邻域分析单元;韩萍等[4]基于多源DEM数据,提取青岛市部分地区小流域水文特征并进行对比分析;胡最等[5]依据1∶25万DEM数据,开展湖南省地貌形态特征分类研究。虽然各专家学者已对各式DEM展开不同探讨,但利用多源DEM数据进行地形因子提取和对比分析的研究还略显不足,且已有研究表明,不同分辨率的DEM数据所表达的地形信息的容量和精度等均存在显著差异[6]。鉴于此,针对同一研究区,利用不同分辨率DEM数据提取地形因子并进行定量分析具有一定研究价值。本文以都江堰灌区作为研究区,首先,基于三种不同分辨率的DEM数据,提取高程、坡度、坡向、地形起伏度、地表粗糙度、地表切割深度等六个地形因子;接着,对提取结果进行定性和定量分析;最后,利用皮尔逊相关性分析筛选出研究区最佳地形因子。研究结果可为都江堰灌区的地貌类型精细化划分、灾害风险评估及土壤侵蚀等提供依据和参考。

1 研究区概况与数据来源

1.1 研究区概况

都江堰灌区位于四川盆地西部,介于东经103°29′～105°24′、北纬29°24′～31°29′之间。从地形地貌条件来看,灌区地貌类型以丘陵为主,丘陵面积占总面积的比例高达71.1%[7];海拔高度介于210～1 004 m,以龙泉山脉为分界,具有典型的西北高、东南低特征。从气候水文方面来看,灌区属亚热带湿润气候类型,全年气候适宜,年均降雨量900～1 240 mm,降雨量受季节影响较大,呈现春旱、夏洪、秋涝、冬干的特点[8]。目前,都江堰灌区实际农田灌溉面积达75.37万hm2,涉及成都、绵阳等7市40县(市、区)。研究区概况如图1所示。

图1 研究区概况

1.2 数据来源

本文利用的三种DEM数据分别为ASTR GDEM、SRTM DEM及ALOS DEM。其中,ALOS DEM空间分辨率为12.5 m,来源于NASA地球数据中心(https：//search.asf.alaska.edu/#/);ASTR GDEM空间分辨率为30 m,SRTM DEM空间分辨率为90 m,均来源于地理空间数据云(http：//www.gscloud.cn/)。基础数据的准确性是提取地形因子的关键,基于GIS软件平台,对三种原始DEM数据进行镶嵌、裁剪及投影变换等预处理,以获取研究区同一坐标系下准确的基础数据。

2 研究方法

2.1 归一化

不同的地形因子具有不同的原始值,并且对量纲和量纲单位也具有不同的要求,为了使各地形因子具有可比性,需要采取适当的手段对各地形因子进行处理,从而将各地形因子的值限定在一定的范围内[9]。因此,为了便于开展地形因子间的相关性分析,本文对各地形因子的原始值进行归一化处理,使处理结果全部介于[0,1]。

(1)

式中,X为归一化值;x为各地形因子原始值;xmin为各地形因子的最小值;xmax为各地形因子的最大值。

2.2 皮尔逊相关性分析

皮尔逊相关系数又叫皮尔逊积矩相关系数,它是一种可以准确量算两个变量之间的相关程度的经典统计学方法,被广泛地应用于国民消费结构[10]及灾害敏感性评估[11]等方面。以变量x和变量y为例,那么x与y之间的相关系数计算公式为：

(2)

式中,R表示变量x和变量y之间的相关系数;n表示变量x和变量y的个数。相关系数与具体的相关程度如表1所示。

表1 皮尔逊相关系数与相关程度

3 地形因子提取与分级

3.1 高程

高程是地貌形态的最基本要素,DEM所提供的原始数据就是高程数据[12]。对三种高程数据按分辨率从高到低排序,获得研究区的高程值分别介于210～1 004 m、192～1 014 m及254～1 042 m,考虑前人对高程的分级标准并结合研究区实际高程分布情况,将高程值分为<400 m、400～600 m、600～800 m、800～1 000 m及>1 000 m等五个等级,如图2所示。

(a)12.5 m分辨率 (b)30 m分辨率 (c)90 m分辨率

3.2 坡度

坡度定量反映了地表的陡缓程度,是常见的微观地形因子之一[13]。利用Arc GIS 10.2软件的表面分析工具,从三种不同DEM数据中提取坡度因子,按分辨率从高到低排序,最终提取结果为0～69.36°、0～71.28°、0～47.97°。依据曾佩枫等[1]对坡度的等级划分,本文将坡度分为<3°、3°～5°、5°～10°、10°～15°、15°～20°、20°～25°、25°～30°及>30°等八个等级,如图3所示。

(a)12.5 m分辨率 (b)30 m分辨率 (c)90 m分辨率

3.3 坡向

坡向是地表任意一点切平面的法线矢量在水平面的投影与过该点正北方向的夹角,可通俗地理解为自高而低的方向[1]。与坡度的提取及排序方式相似,将研究区坡向划分为平坦-1、北(0～22.5°及337.5°～360°)、东北(22.5°～67.5°)、东(67.5°～112.5°)、东南(112.5°～157.5°)、南(157.5°～202.5°)、西南(202.5°～247.5°)、西(247.5°～292.5°)及东北(292.5°～337.5°)等共九个类别,如图4所示。

(a)12.5 m分辨率 (b)30 m分辨率 (c)90 m分辨率

3.4 地形起伏度

地形起伏度(Rf)是指某一特定区域内最高点高程与最低点高程的差值,能够较为直观地反映地表形态[14],是地貌类型划分研究中最常用的因子之一,其具体计算公式为：

Rf=Hmax-Hmin

(3)

最佳分析窗口的确定是影响地形起伏度准确性的最关键因素,基于前人[2,15-16]研究成果,对于ALOSDEM、ASTR GDEM及SRTM DEM,分别采取13×13窗口,15×15窗口,17×17窗口为最佳分析窗口,其面积分别为0.026 4 km2、0.202 5 km2及2.34 km2。利用Arc GIS 10.2软件平台的焦点统计工具,从三种DEM数据中提取地形起伏度因子,按分辨率从大到小排序,最终结果为0～216 m、0～363 m、0～509 m。依据以往研究成果的数字格局地貌研究分类方法,将三种提取结果划分为五类,分别为<30 m(平原)、30～50 m(台地)、50～200 m(丘陵)、200～500 m(小起伏山地)及>500 m(中起伏度山地),如图5所示。

(a)12.5 m分辨率 (b)30 m分辨率 (c)90 m分辨率

3.5 地表粗糙度

地表粗糙度(R)是指某一固定区域内地球表面积与其投影面积的比值,是反映地表形态的一个重要宏观因子[2]。在提取的坡度数据基础上,利用Arc GIS 10.2软件的栅格计算器工具,获取研究区地表粗糙度数据。结果分别为1～2.87、1～3.12及1～1.49。参照贾丽娜[17]的分级标准,将三种分辨率的地表粗糙度数据分为1～1.03、1.03～1.08、1.08～1.14、1.14～1.21、1.21～1.29、1.29～1.39、1.39～1.57、1.57～1.93及>1.93等共九个等级,如图6所示。

(a)12.5 m分辨率 (b)30 m分辨率 (c)90 m分辨率

3.6 地表切割深度

地表切割深度(D)定义为某一特定区域内平均高程与最小高程的差值[18]。在提取地表切割深度的过程中,最佳分析窗口的选择与提取地形起伏度的窗口相同,提取结果分别为0～112.23、0～204.94及0～285.27。参照梁宏艳等[2]的划分标准,将三种提取结果分为三个等级,分别为<30 m(微切割山地)、30～100 m(浅切割山地)及100～500 m(中等切割山地),如图7所示。

(a)12.5 m分辨率 (b)30 m分辨率 (c)90 m分辨率

4 结果分析

4.1 地形因子对比分析

为了更直观地分析三种不同分辨率地形因子的等级划分结果,对各地形因子不同等级的栅格数量、面积、栅格数量占比及面积占比进行统计。不同分辨率下各等级高程数据的栅格数量占比如表2所示。

表2 研究区不同分辨率下各等级栅格数量占比(高程)

由表2可知,对于30 m和90 m分辨率的高程因子来说,各等级所占比例接近,且400～600 m等级的占比最高,达70%以上。而对于12.5 m分辨率的高程因子,400～600 m所占比例却不足50%,更多的区域被划分为<400 m等级。

不同分辨率下各等级坡度数据的栅格数量占比如表3所示。

表3 研究区不同分辨率下各等级栅格数量占比(坡度)

由表3可知,三种分辨率的坡度数据在栅格数量占比上未表现出明显规律。除<3°分级的差异较大外,其余各等级所占比例相差不大。

不同分辨率下各等级坡向数据的栅格数量占比如表4所示。

表4 研究区不同分辨率下各等级栅格数量占比(坡向)

由表4可知,三种分辨率的坡向数据,除在平面所占比例相差较大外,其余八个分级所占比例接近,且均在10%上下波动。

不同分辨率下各等级地形起伏度数据的栅格数量占比如表5所示。

表5 研究区不同分辨率下各等级栅格数量占比(地形起伏度)

由表5可知,30 m和90 m分辨率的地形起伏度数据在对平原和丘陵的划分结果相似,而12.5 m和90 m分辨率的地形起伏度数据对台地的划分结果相似。此外,90 m分辨率的地形起伏度数据中,有极少数区域被划分为中起伏度山地。

不同分辨率下各等级地表粗糙度数据的栅格数量占比如表6所示。

表6 研究区不同分辨率下各等级栅格数量占比(地表粗糙度)

由表6可知,三种分辨率的地表粗糙度数据中,各等级栅格数量所占比例差别不大,90%以上的栅格属于1～1.08之间。

不同分辨率下各等级地表切割深度数据的栅格数量占比如表7所示。

表7 研究区不同分辨率下各等级栅格数量占比(地表切割深度)

由表7可知,对于30 m分辨率和90 m分辨率的地表切割深度划分结果,各级占比差异不大。其中,70%以上的区域被划分为微切割山地,20%以上的区域被划分为浅切割山地,仅有极少区域被划分为中等切割山地。而对于12.5 m分辨率的地形切割深度分级结果来说,仅有微切割山地和浅切割山地两种类型,且微切割山地的占比高达94.78%。

4.2 最佳地形因子确定

为了确定研究区的最佳地形因子,利用皮尔逊相关系数对各因子进行相关性分析,具体步骤为：①利用Arc GIS 10.2软件的随机选点功能,生成2 000个独立的随机点;②提取各地形因子的原始值至随机点;③利用公式(1),将提取的原始值进行归一化处理;④基于SPSS 25软件,对归一化过后的数据进行皮尔逊相关分析。各因子的相关性分析结果如表8所示。

表8 地形因子相关性分析结果

由表8的分析结果并结合表1可知,坡向与其他地形因子间的相关系数均小于0.2,即极低相关性,表明坡向与其他地形因子相比,对研究区地形特征差异表达不明显,应予以剔除。除此之外,切割深度与起伏度、切割深度与坡度的相关系数均大于0.8,即极强相关,表明三者表达的地形信息有所重复。考虑到起伏度和坡度已被大量研究确定为地形划分的主要因子,因此,选择保留起伏度和坡度而删除切割深度。最终,研究区最佳地貌类型划分的地形因子为高程、坡度、地形起伏度及地表粗糙度。

5 结论与讨论

以都江堰灌区为研究区,首先,基于三种不同分辨率的DEM数据,提取高程、坡度、坡向、地形起伏度、地表粗糙度及地表切割深度等六个地形因子;接着,对提取结果进行等级划分和统计;然后,利用皮尔逊相关性分析确定研究区地貌分类的最佳地形因子。研究结果如下：

(1)由不同分辨率提取的研究区地形因子结果存在明显差异。对高程而言,90%以上的栅格分布于<600 m等级;对坡度而言,除<3°等级相差较大外,其余等级差距较小;对坡向而言,平面等级的占比相差在10%以内,其余等级差距不大;对地形起伏度而言,12.5 m分辨率DEM的提取结果与其他两种提取结果差异较大,所有的区域都被划分为平原、台地及丘陵。呈现出分辨率越低,提取结果种类越多的趋势;对地表粗糙度而言,各等级栅格数量所占比例接近;对地表切割深度而言,12.5 m分辨率DEM的提取结果仅划分出微地形切割和浅地形切割两种,且微地形切割的比例高达95%。

(2)由皮尔逊相关性分析结果可知,都江堰灌区的最佳地形因子为高程、坡度、地形起伏度及地表粗糙度。