新地貌指标及山地地貌自动区划方法研究

　　摘要：地理栅格统计分类及聚类分析难以实现山地地貌的自动区划制图，实践中一般使用半定量或手工综合的方法绘制地貌图，而山地地貌区划研究的关键是地貌空间特征指标的构造及空间栅格的综合归并. 基于此，提出了地形隆起高度与切割深度、皱褶多度、地表破碎度等新地貌特征指标，并以传统统计ISO聚类结果为基础，构造了一种基于空间“核颗粒”和“纹理指标”的二次聚合方法. 研究表明：①新提出的地貌指标能以具体数值表达山体向上隆起以及向下切割的程度，度量水平和垂直方向上的破碎度;②山地地貌是多种聚类分析及空间纹理再聚合的结果，利用“核颗粒”和“纹理指标”可划分出符合认知体系的地貌类型，并将中国陆地全境划分为一级6类、二级33类的地貌区，与以往的制图结果具有一致性. 二次聚合划分地貌类型的方法，实现了地貌的计算机自动区划及制图，适用于具有纹理特征的地学类别的空间划分研究.

　　关键词：山地地貌;新地貌指标;自动区划

　　地貌区划及制图是地球科学系统研究与应用中一项基础性工作，对国土空间规划、资源开发利用、自然灾害管理、生态环境保护、生态恢复治理等都具有重要意义. 中国地貌学家开展了多种地貌区划方法研究，同时编制了各种地貌区划图[1-2]，如省级地貌分区[3]及部分地区的地貌分区图等[4].

　　李炳元等[5]利用半定量与定性方法将中国划分为两级地貌区，程维明等[6]基于2013年提出的地貌区划结果构建了全国五级地貌区划方案，但各种区划方案结果都有所差别. 其输出结果不确定性主要受分类原则、分类指标、分类方法以及数据分辨率的影响，由于缺乏定量化数据标准和模型化的计算机分类方法，多使用定性或半定量化划分地貌，造成了各种分类方法无法相容和统一.在地学分类指标方面有指标多样、等级较少，且较多研究趋向于利用海拔高度和地势起伏度进行划分[5].

　　然而中国地貌复杂多样，区域差异明显，仅利用海拔高度和地势起伏度数据进行地貌分类不能完全反映地表形态的整体性，必须寻找地貌特征明显的指标量化其地貌形态，将其细分为更易于管理与应用的组成部分[7]. 在分类方法方面，随着数字高程模型数据[8]、ArcGIS软件技术和遥感技术的发展，为地貌自动划分提供了更为广泛的数据和技术支持[9-12].

　　现有的基于聚类、规则知识及典型样点等地貌自动划分方法，都受限于语义定义模糊、多尺度分割困难、空间信息难以量化的影响[13].聚类是利用统计学原理将地貌属性相似的地理栅格归并为空间斑块，再将其定义为地形的某个部位[14-16]. 事实上，在山地地貌区域，空间变异性强，聚类生成的斑块大小不一且极其琐碎，因而不能保证地貌的完整性[17]，也无法产出符合认知概念和规则知识的地貌类型，难以进行地学解释[18].

　　其根本原因是地貌区划是定义在空间等模式变量上的函数，而不是像素上的测度. 现有研究中面向对象的影像方法克服了栅格分类的破碎图斑和“椒盐现象”[19]，而分割尺度的确定与面向对象的空间聚类是关键[20]. 其次，由于卫星影像存在包含地表覆盖物信息，难以有效融合地貌特征和尺度差异等问题而制约着地貌分类的进一步研究[21-22].

　　数字高程模型(Digital Elevation Models，DEM)的优势在于可“过滤”地表覆盖物信息而只表达地形自身信息，而“纹理”特征又是区分不同地貌形态的重要依据[23-24]，成为宏观尺度上地形特征量化与识别研究的突破点. 在数据分辨率方面，数字高程模型(DEM)分辨率越高，地形模拟越成功，但地貌形态具有分形几何特征，随分辨率的增加，数据量呈几何级数递增.选择符合应用要求又能兼顾计算机容量与处理能力的分辨率数据，制约着DEM应用范围及地貌区划研究[25].

　　基于以上分析，本文提出了能够表达地貌特征的新地学指标，并结合传统的地貌指标构建分类体系. 为克服一次聚类产生的琐碎图斑及山地地貌完整性被分裂的问题，引入了“核颗粒”和“纹理指标”特征概念，并提出了二次琐碎图斑再聚合的方法，进行中国全陆域范围的地貌区划研究，该方法对分异性较强的地理空间分类制图具有参考应用价值.

　　1 研究区域及数据

　　1.1 研究区概况

　　中国位于亚洲东部、太平洋西岸，地势西高东低呈三级阶梯分布. 第一级阶梯内有号称“世界屋脊”的青藏高原，其间耸立着多条著名的高大山系，其中喜马拉雅山脉是地球上最新隆起的年轻山系之一. 第二阶梯包含若干高原和盆地，从青藏高原北缘、东缘到大兴安岭、太行山、巫山、雪峰山. 第三级阶梯地势低平主要以高原、盆地和平原为主[26]. 中国地势起伏较大，分布着许多高大的山脉及被这些山脉所围绕隔开的大型地貌单元，由西向东逐级下降[27]. 其中，山地和高原是构成中国地貌基本轮廓的主体，山地、丘陵和高原约占中国总面积的2/3[28].

　　1.2 数据源及预处理

　　下载覆盖中国区域的SRTM3数据集，栅格分辨率为90m，运用ArcGIS软件底层函数对数据文件拼接运算，将1°×1°的分块地图拼接成中国的整体地图.由于接缝边缘存在像素灰阶不一致，构成肉眼可见的异常缝，采用滤波进行消除.对于某些超过正常海拔值的区域，用邻域灰阶进行填充或条带过滤处理.像元降尺度处理：数字高程模型是描绘地形地貌的基础数据，经过提取地形信息和再演算，可实现各种地表形态的数字化模拟表达.随着DEM分辨率的提高，可分辨和提取流域的尺度不断提高，同时数据量也呈几何级数增长.本研究以全国陆域范围为对象，考虑到研究区及大尺度地貌划分对数据的需求，分别利用90m、500m、1000m分辨率的DEM数据提取地形指标，经测算比对分析1000m尺度栅格数据能满足研究应用要求，地貌碎块的表达也较好，故以其为地理指标的栅格计算尺度.

　　2 原理与方法

　　2.1 地貌新变量

　　地貌要素特征的量化表达一直是难点，传统的海拔、坡度、坡向等不足以表达宏观地貌特征，要实现地貌类的划分需寻找新的特征指标.

　　2.1.1 地形隆起高度与切割深度

　　对于发育于大倾斜面上的山地地貌，在水平和垂直方向上具有复杂的向上隆起和向下切割的变化特征. 传统的海拔难以表达地貌隆起、沟谷切割或盆地塌陷的状态. 基于此，提出以每个流域的平均海拔为样本点构建地表趋势面，以此为起算点测算地形隆起高度与切割深度的变化特征[29]. 该指标以山体基部为起算点，度量山体隆起或切割下沉程度.

　　2.1.2 皱褶多度

　　褶皱是指岩石因受力发生一系列的波状弯曲而显示的变形，使原来近于平直的面表现出向上隆起或塌陷的曲面形态. 皱褶多度与地貌类型密切相关，平原地貌几乎无皱褶，山地地貌破碎，皱褶明显. 王艳霞等[30]基于滑动窗口方法，利用地表表面积与地物投影的平面面积比值提取褶皱多度(公式(1))，但该种方法采用面积线性对比反映研究区的褶皱变化程度，而中国山体坡度多介于20°～50°中等大小的山体发育，坡度差异小，故投影面积区分也不显著，使得褶皱量化对比不明显. 基于此，本文以具体坡度数值对褶皱多度指标进行改进(公式(2))，以DEM为数据源，利用ArcGIS软件构造褶皱地貌单元模型.

　　2.2 地貌类的栅格定义与二次聚合

　　利用ArcGIS软件提供的ISO聚类方法进行多指标统计聚类，将地貌划分为碎片单元. 因各碎块面积大小不一，有大面积的片区，也有小面积的椒盐，通过消除将小面积逐步归并到大片区中，被归并的小面积碎片称之为“核颗粒”. 如果碎片面积大小相近，多类碎片可利用重复模式结构表达地貌类，其集合体即为“纹理指标”. 依据“核颗粒”将空间相邻且相互关系密切的小碎片联通为大片区“纹理”，可划分出地貌类别。

　　3 结果与分析

　　3.1 地貌分类指标分析

　　基于数字地貌信息提取的地形因子是一个定量化指标，能客观地反映地貌形态特征与发生发育演化的规律，揭示各指标与地形地貌的因果联系. 因此，地貌指标的选择成为地貌分类的关键. 不同的地形指标从不同维度刻画地貌特征，其间存在相关关系，为消除指标相关性及信息冗余的影响，借助ArcGIS软件栅格代数运算工具，对9个指标进行中国全陆域范围的相关性统计分析。

　　相关性在0.6以上的指标，但在进行单个地貌指标分析时发现，地表破碎度、地形起伏度虽相关性较高，但两者从不同方面刻画地貌形态特征，故都保留运用于最佳地貌区划指标体系的筛选构建中. 经对各个指标组合对比分析，本文最终选择地形隆起高度与切割深度、皱褶多度、地表破碎度、海拔高度、地形起伏度共5个指标研究样区地貌区划.

　　3.1.1 地形隆起高度与切割深度

　　地表形态指标是地貌区划的主要依据之一. 地形隆起高度与切割深度能够度量山体向上隆起、沟谷向下切割的程度，体现地貌在高度上的真正垂直变化，刻画大倾斜面上大地域间构造状况及整体变化趋势.

　　地貌隆起最大高度为2755 m，切割最大深度为−2051 m. 从青藏高原到周边区域过渡带，隆起数值较大，其中以横断山区隆起最为剧烈、隆起面积最大. 其次，青藏高原北面和柴达木盆地西面、南面的昆仑山，北面的天山地形变化呈窄条带发育，海拔下降大、沟谷深切，隆起与向下切割形态显著. 塔里木盆地、准噶尔盆地、四川盆地、东北平原、华北平原和长江中下游平原等地为地貌切割平缓的典型代表单元，隆起数值较小.

　　3.1.2 皱褶多度

　　地貌形态特征主要受内外营力作用影响反映研究区地貌差异性，如花岗岩长期受外营力作用侵蚀，在中国北方地区常呈现高大险峻的山地地貌如华山、泰山等. 研究定量化的物质组成分异指标对地貌区划意义重大，褶皱多度能够衡量特定区域内地表整体形态，度量地貌坡面组合差异.中国地貌单元内褶皱度最低值集中分布于平原和中低山. 褶皱多度较明显区域地形起伏剧烈、坡度大，主要分布于昆仑山脉、天山山脉、横断山区、祁连山区、秦岭、台湾山脉等高山区.

　　3.1.3 地表破碎度

　　地表破碎度能从水平垂直方向上揭示研究区地表组合完整程度，刻画空间上的分布规律. 从外营力方面表征地貌分带现象，成为地貌区划单元转折部位划分的依据.中国地表破碎度呈现较高空间异质性，如东部的东北平原、华北平原、长江中下游平原等大部分地区;中西部的准噶尔盆地、塔里木盆地、柴达木盆地、四川盆地等典型盆地地貌区，地表完整度良好，破碎程度小;中部的秦岭、太行山、云贵高原等地地表相对崎岖，存在一定程度的破碎;西部的天山、阿尔金山、昆仑山、祁连山、横断山等区域多高山峡谷分布，高差变化大，地表破碎显著.

　　4 讨论与结论

　　4.1 讨论

　　(1)山地地貌类是栅格地理单元组成的“核颗粒”和“纹理指标”的空间再聚合集合体. 基于统计聚类结果为基础，利用斑块大小、纹理特征、斑块之间的邻近及可联通性等拓扑关系，将碎块图斑迭代聚合为片区可实现地貌类型的计算机自动区划及制图. 同时，自动区划和制图受综合自然地理要素及地貌分类指标体系的影响，需考虑不同地貌区边界.

　　在进行山地地貌区划时，应注重整个区域垂直地形特征并结合沟谷发育和大地构造的差别，如秦岭北坡山地边界的划分，以分水岭为界从东向西根据垂直地形差异关系划分地貌区. 在以低山为主，且有个别中山分布区的界线确定时，可根据地形起伏状况把突出的中山或与周边区别较大区域单独划分出来. 在以平原为主，有低山及丘陵地貌区，可根据实际地貌制图的平面结构及各指标的“纹理”特征进行划分，如鲁东低山丘陵周边各指标在该区域表现平缓，低山与周围其余地貌相比较突出，结合其所处地理位置故将其划分为鲁东低山丘陵区.

　　(2)由于地貌指标选取与分类方法不同，使其部分区划边界和属性判断存在差异. 而又缺少可对比的矢量数据，只能通过可视分析进行比对，其结果与现实会存在些许偏差. 如在类别上，程维明等[5-6]提出江河源丘状山原、羌塘高原湖盆区2个二级地貌区，但因该区受地形影响各指标表现出高度相似性，聚类后可形成大面积片区，因此将其合并划分为羌塘高原区;横断山和江河上游高山峡谷区多高大山体发育，其形态相似聚合后由多类面积相似的碎块组成，不存在优先归并顺序故将其合并划分;喀喇昆仑山极高山与昆仑极高山各指标表现特征相似而合并为一个地貌类。

　　宁镇山脉(江苏)面积较小，地形起伏变化和缓且与周围地貌区分不显著，故将其归并到华北平原地貌区. 在分界线上，程维明等[5-6]以伏牛山中部为界划分秦岭−大巴亚高山区，本文依据定量化的指标表达为保证山体完整性仍以伏牛山为界. 除此之外，中国南部地区地形破碎，聚类后多面积近似多类碎块镶嵌形成的模式结构，设置算法时具有不确定性，故在一级、二级地貌区及界线划分上与已有研究结果存在差异. 但与中国地形数据对比分析发现，本文的地貌分类结果与原始数据及地貌的科学认知体系相吻合，具有一定的现实意义.

　　地形测绘论文：多源数据补充和优化数字化大比例尺地形图测绘

　　4.2 结论

　　(1)本文构建了地形隆起高度与切割深度、褶皱多度、地表破碎度、地形起伏度、海拔高度组合的中国地貌计算机自动划分的指标体系. 在新地貌指标特征表达方面，地形隆起高度与切割深度指标更能表征地貌的垂直变化，实现宏观大地域间的地貌构造差异;地表破碎度指标用水平沟壑密度、垂直方向的隆起与切割高度状况以及窗口地形变化特征共同描述地貌破碎程度，比以往的指标更为全面、客观;褶皱多度指标从外营力方向表达了地貌起伏变化特征，对多种地理科学现象的研究都可提供借鉴作用.

　　(2)以聚类分析输出的地貌类碎片斑块为基本单元，构建基于“核颗粒”和“纹理指标”逐步迭代的二次空间聚合算法，可实现中国地貌的计算机自动区划和制图研究. 最终，发现联通面积Sk为9×103 km2时，二级地貌区开始形成;联通面积Sk为3×104 km2时，大部分二级地貌区形成;联通面积Sk为1.2×105 km2时，一级地貌区形成. 据此，将中国地貌划分为6个一级地貌区和33个二级地貌区. 二次聚合解决了山地地貌的计算机分类问题，该方法也适合于地理分异性较强的土壤制图、植被分类制图及遥感图像分类研究等.

　　参考文献：

　　尤联元, 杨景春. 中国地貌[M]. 北京: 科学出版社,2013.You L Y, Yang J C. Geomorphology in China [M].Beijing: Science Press, 2013.

　　[1]郭子良, 崔国发. 中国地貌区划系统—以自然保护区体系建设为目标[J]. 生态学报, 2013, 33(19): 6264-6276.Guo Z L, Cui G F. Geomorphologic regionalization ofChina aimed at construction of nature reserve system[J].Acta Ecologica Sinica, 2013, 33(19): 6264-6276.

　　[2]王丽娜, 丁文广, 许丹阳. 基于DEM的甘肃省地貌形态特征分类[J]. 水土保持通报, 2019, 39(1): 268-263.DOI: 10.13961/j.cnki.Stbctb.2019.01.041.Wang L N, Ding W G, Xu D Y. DEM-based geomorphologicalfeatures of Gansu Province[J]. Bulletin Soiland Water Conservation, 2019, 39(1): 268-263.

　　[3]常直杨, 孙伟红, 王建, 等. 青藏高原及其邻近地区地貌类型划分[J]. 山地学报, 2017, 35(1): 1-8. DOI: 10.16089/j.cnki.1008-2786.000189.Chang Z Y, Sun W H, Wang J, et al.

　　Objectorientedmethod based on classification of geomorphic type inthe Tibet Plateau and adjacent regions[J]. Mountain Research,2017, 35(1): 1-8.

　　作者：肖丽琼1，周汝良1**，陈永刚1，蔡志勇2，黄 彦2