黄土丘陵沟壑区地表过程响应单元划分及其地形特征提取

秦伟，朱清科，左长清，赵磊磊，邝高明，王昭艳，张京凤，单志杰

(1.中国水利水电科学研究院泥沙研究所，100048;2.北京林业大学 水土保持与荒漠化防治教育部重点实验室，100083:北京)

地表过程(earth surface processes)，即地球表面大气圈、生物圈与岩石圈之间因物质与能量的交互作用而形成的一系列物理与化学过程，主要包括岩土物质的侵蚀、搬运和沉积，水热养分的运移、转化和循环，以及土地利用演变及其环境效应等。地表过程是地球表面固有的自然过程，由此引发地貌格局演变，水系发育变迁，并影响生态环境，同时，地表过程与人类活动密切相关、相互影响。在全球气候变化与人类活动加剧的背景下，地表过程将更加深刻地影响人类的可持续发展。全面揭示地表过程规律，有效地向避害趋利的方向调控，以促进人与自然和谐相处，已为世界普遍关注，诸如土壤侵蚀、山地灾害、生态水文等典型地表过程也因此长期成为学界研究焦点。由于地表过程的复杂性，通常需确立一定标准，据此将研究范围划分为一系列基本单元，先就基本单位完成模拟或分析，再经累加或串联实现整个研究范围的模拟以及相关特征的揭示。这种化繁为简、化整为零的思路是解决地理学、生态学、水文学等众多学科复杂问题的重要途径。其中，所划分的基本单元，即为本文的地表过程响应单元。地表过程响应单元划分是否合理、获取是否高效很大程度上决定地表过程研究结果是否可靠、过程是否顺利。现有的划分主要包括规则单元和不规则单元2 类。规则单元是将研究区域按一定分辨率划分为均匀的方格单元，即栅格单元。这种方法简单、快速、易实现，且在地理信息系统技术广泛应用的背景下，能更好地与遥感影像等多源数据结合分析。目前采用最多，但其仅考虑了地形等环境因素影响地表过程的空间差异性，忽略了一定空间范围内下垫面影响地表过程的相似性。不规则单元是根据地形地貌特点、水文运移路径或土地利用类型等确定特征边界将研究区域划分为不规则的面状班块。这种方法在特定角度更真实反映了研究对象的固有属性，更利于后期分析研究，但相比规则单元划分，其实现难度更大，尤其针对较大空间尺度时，往往难以快捷地获取预期结果。比较有代表性的不规则相应单元包括水文响应单元(Hydrological Response Unit，HRU)［1］、水文相似单元(Hydrological Similar Unit，HSU)［2］和分组相应单元(Grouped Response Unit，GRU)［3］等。

黄土丘陵沟壑区是世界上地形最破碎复杂的地区之一，加之干旱与暴雨交叠的气候条件，造成该区土壤侵蚀、滑坡坍塌、生态水文等为代表的地表过程具有显著的特殊性和复杂性。针对该区特点，建立一种合理、高效的地表过程响应单元划分方法对推动地表过程研究具有重要的基础意义。笔者以黄丘陵沟壑区典型流域为例，拟提出针对土壤侵蚀为代表的地表过程模拟研究的响应单元划分方法及其在GIS 平台中的实现技术，以期为促进该区地表过程研究提供有益参考。

1 研究区概况

四面窑沟流域位于陕北吴起县铁边城镇、北洛河上游南岸，E 107°45'47″～107°49'24″，N 36°54'08″～36°59'58″，属北洛河水系，面积32.68 km2。流域内海拔1 369 ～1 730 m，沟壑纵横，属典型黄土高原丘陵沟壑地貌类型(图1)。土壤主要为地带性黑垆土剥蚀后广泛发育在黄土母质上的黄绵土，质地为轻壤。多年平均降雨量478.3 mm，降雨年际变化大，季节分配不均，7—9 月降雨量占年均降雨量的62.4%，且多为暴雨。受耕作等人类活动影响，流域内天然植被退化严重，现存植被主要为人工次生植被，以沙棘(Hippophae rhamnoides)、柠条(Caragana mocrophylla)等灌木，河北杨(Populus hopeiensis)、小叶杨(Populus simonii)、山杏(Prunus sibirica)、杜梨(Pyrus betulaefolia)等阔叶乔木以及紫花苜蓿(Medicago sativa)等人工牧草为主，林草覆盖率约80%。

2 流域地表过程响应单元划分

图1 研究流域区位图Fig.1 Location of the study watershed

除气候条件外，地表过程主要受下垫面影响。划分响应单元的核心在于针对下垫面影响地表过程的空间异质性，分类、合并地表过程差异显著的下垫面，形成内部地表过程特征相同或相似的下垫面单元。土地利用与地形地貌是下垫面的主要内容。其中，土地利用分类标准和划分边界比较明确，即按不同土地利用类型及其相连的界线进行划分。地形地貌所包含的具体内容更为丰富，不同地形地貌指标均能不同程度影响地表过程，且以不同地形地貌指标为依据的划分边界也不尽清晰;因此，响应单元划分的难点主要在于地形地貌特征解析及其边界提取。就黄土丘陵沟壑区的土壤侵蚀而言，坡度、坡向、坡位是影响土壤侵蚀的主要地形地貌特征。坡面作为组成该区地貌结构的基本组份，其内部的坡度和坡向虽存在空间分异，但并不显著，独立坡面在整体上具有基本一致的坡向和坡度。坡位对土壤侵蚀的影响集中表现在梁坡和沟坡间土壤侵蚀强度和机制的差异，即梁坡以薄层径流的剥蚀为主，侵蚀强度较低;沟坡以含沙股流的冲蚀为主，侵蚀强度较高。因此，以小流域为对象，参照水文研究中所提及的水文响应单元［1］，将土壤侵蚀响应单元定义为由流域沟谷线、山脊线、分水线及沟缘线交叠构成的具有相对均一坡向、坡度与坡位的独立坡面，并在Arc-GIS9.2 中，以数字高程模型为基础，实现响应单元划分及其地形特征提取。

2.1 流域Hc-DEM 的建立

数字高程模型(Digital Elevation Model，DEM)作为真实地貌的离散数字表达，是水文、土壤、地貌、生态等领域的重要基础数据，分为规则格网(如grid、raster 等)和不规则三角网(Triangulated Irregular Network，TIN)［4］。通常多基于数字化等高线生成不规则三角网模型，然后再转化为规则格网模型;然而，因内插计算误差或少数真实地形影响，由此获得的DEM 表面将存在局部凹陷洼地，阻断汇水路径从而影响后期的地表过程分析［5］。为此，采用陕西省1∶1万纸质地形图，扫描后在R2V 软件中完成数字化，进行平滑、赋值等编辑处理，形成矢量等高线;然后在ArcGIS9.2 中，利用Spatial Analysis 模块的Topo to Raster 工具，基于ANUDEM 算法生成水文关系正确的数字高程模型(Hydrologically Correct DEM，Hc-DEM)。Hc-DEM 消除了低洼等局部干扰地形，能形成通畅水文路径，可为后期分析奠定基础。

2.2 流域地表过程响应单元划分

2.2.1 沟谷线、山脊线及分水线提取 沟谷线、山脊线及分水线不仅是流域地形地貌的分界线，也是流域水文循环过程中的特征线。在ArcGIS9.2 中，提取沟谷线、山脊线及分水线均以数字高程模型为基础数据，在Hydrology 模块中完成。

1)沟谷线提取。利用Hydrology 模块的Flow Direction 工具和Flow Accumulation 工具，以Hc-DEM 为基础数据，先后获得网格单元的汇流方向和累积汇流量;再利用Raster Calculator 工具，输入不同累积汇流阈值选择生成不同级别的沟谷线。输入的累积汇流阈值越小，获得的沟谷密度越大，长度越长。流域内的沟谷系统是线状水流作用形成的不同等级和规模的沟谷总称。根据其形状、大小和发育阶段等特征，通常分为细沟、浅沟、切沟、冲沟和坳沟等5 类［6］，内部相应存在以细沟侵蚀、浅沟侵蚀、切沟侵蚀、冲沟侵蚀和河道侵蚀为主的侵蚀类型。由于细沟、浅沟多存在于独立坡面内，且难由一般精度的数字高程模型中提取，切沟、冲沟多存在于沟缘线上下，结构不稳定，不宜用于划分响应单元;因此，以沟道和河道为提取对象，通过将不同累积汇流阈值对应的提取结果与高分辨率遥感影像(SPOT-5，2.5 m×2.5 m 分辨率，下同)进行对比，最终将累积汇流阈值0.6 km2时提取的沟谷作为研究流域的输沙主沟道和永久河道。

2)山脊线与分水线提取。山脊线与分水线的提取方法与沟谷线一致。由于是汇流起点，因此，累积汇流阈值为零的栅格所组成的线即山脊线与分水线。其中，分水线是流域最外围的连续闭合山脊线。2.2.2 沟缘线提取 黄土丘陵沟壑区从分水线到沟底，地貌形态存在明显地带分异，依次包括梁峁顶部、梁峁坡、沟坡和沟底。其中，梁峁顶部和梁峁坡构成沟间地(梁坡)，沟坡和沟底构成沟谷地(沟坡)，沟间地与沟谷地之间以沟缘线为界(图2)。沟缘线上下的土壤、地形和植被存在差异，并由此导致土壤侵蚀等地表过程具有不同的特点。作为正负2大地形分区的界线，沟缘线是黄土丘陵沟壑区地表过程研究中十分重要的地貌特征线。

图2 黄土丘陵沟壑区典型坡面横剖面示意图Fig.2 Structure feature of typical slope profile in loess hilly and gully region

沟缘线提取包括利用等高线或高分辨率航片、卫片手工勾绘和基于数字化高程模型或遥感影像的编程提取。手工勾绘精度较高，但主观性大、效率低，难以满足大、中空间范围的研究需要;因此，基于特定数据源的编程自动提取应用越来越广。以Hc-DEM 为数据源，在ArcGIS9.2 中将其转为txt 文本格式作为编程计算对象。编程思路为:根据栅格单元间的汇水路径，逐个计算从分水线出发到河道的每条汇水路径上单元的坡度变化率，将坡度变化率最大的单元作为沟缘线组成点;在此基础上，排除局部破碎点和独立点，形成一条联通、完整的沟缘线;将结果与高分辨率遥感影像进行比对，确定提取结果精度。具体步骤如下。

1)沟缘线上下的地形差异主要表现为坡度陡变，沟缘线位置通常对应坡面坡度变化率最大的部位，因此，首先以3×3 的窗口为范围，采用三次加权差分法计算中心单元坡度，再以同样方法获得其与周围8 个单元的坡度变化率(单元格间位置编号见图3)，以作为沟缘线提取依据:

式中:S 为单元格坡度，(°)，S128为周边8 个单元中某单元的坡度，其他同;fx为x 方向高程变化率;fy为y 方向高程变化率;Z 为高程，m，Z128为周边8 个单元中某单元的高程，其他同;Sr0为单元格坡度变化率;g 为单元格大小，即栅格分辨率，m。

2)采用ArcGIS9.2 中Hydrology 模块的Flow Direction 工具，输入Hc-DEM，获得各单元汇流方向;然后由分水线出发，按流向至河道结束的每条汇流路径上选取坡度变化率最大的栅格，作为该路径(坡面纵断面)上的沟缘线位置点;将该点上、下的单元格分别标记为沟间地单元和沟谷地单元，并相应赋值(1、2)。

图3 DEM 3×3 局部移动窗口Fig.3 3×3 moving part window of DEM

3)坡度变化率筛选和单元类型标记后可初步形成沟缘线，但因局部高地的存在，使个别位置存在单个或几个孤立单元格构成的孤岛，影响流域沟缘线的联通性。为此，对上述孤立点进行排除处理，获得完整、联通的沟缘线以及由其划分的流域沟间地和沟谷地分区。采用ArcGIS9.2 中Conversion 模块的From Raster 工具，将栅格数据转换为矢量线、面数据，并与高分辨率遥感影像比对，判断提取结果的准确性。结果显示，以上算法提取的沟缘线与实际地貌吻合较好，与遥感影像上清楚可辨的沟缘线相比，误差不超过±10 m，即提取误差在1 个单元格内，误差较小，结果可靠(见图4)。

图4 流域沟缘线提取结果Fig.4 Extraction of shoulder line of gully

2.2.3 地表过程响应单元划分 在ArcGIS9.2 中，利用Hydrology 模块，将基于Hc-DEM 提取的流域山脊线和沟谷线按水流方向反向延长与分水线相交;叠加编程提取的沟缘线，并将各地貌特征线进行合并，形成地表过程响应单元［7］(图5)。

图5 流域地表过程响应单元划分结果Fig.5 Division of response units to earth surface processes

3 流域地表过程响应单元地形特征提取

完成地表相应单元划分后，可按以下方法提取其地形特征。

1)响应单元平均坡度提取。在ArcGIS9.2 中，直接提取的坡度是单元格与周围单元格间的最大坡降，属局部坡度，不能作为响应单元平均坡度。因此，采用3D Analysis Tool 的剖面线Interpolate Line功能，以Hc-DEM 为数据源，在响应单元内以分水线/山脊线为起点，沟缘线/沟谷线为终点，沿顺坡方向创建带有高程信息的直线，并利用Create Profile Graph 功能生成直线对应的剖面图;根据剖面图中所包含的直线长度和直线两端的高程差计算坡面坡度［8］。根据规格，每个响应单元创建5 至10 条直线，将平均值作为对应响应单元的平均坡度。

2)响应单元整体坡向提取。采用3D Analysis Tool 的Suface Analysis 功能，提取坡向，并划分为阴坡(337.5°～360°和0°～67.5°)、半阴坡(67.5°～112.5°和292.5°～337.5°)、阳坡(157.5°～247.5°)和半阳坡(112.5°～157.5°和247.5°～292.5°)4 个坡向类别，分别记为1、2、3、4。将响应单元内多数栅格单元的坡向类别作为其整体坡向，以消除局部地形影响。

3)响应单元顺坡长度提取。在响应单元内，采用Measure 功能量算分水线/山脊线与沟缘线/沟谷线间沿垂直等高线方向的距离，除以所在响应单元平均坡度的余弦，获得顺坡实际长度。根据规格，每个响应单元创建5 至10 条直线，将平均值作为对应响应单元的顺坡长度。

4)响应单元面积与周长提取。采用Conversion Tool 的To Coverage 功能将获得的响应单元图层转化为Coverage 格式，则可在属性表中自动增加各单元的周长与面积。

4 结论与讨论

与传统规则单元划分法相比，针对黄土丘陵沟壑区土壤侵蚀的地形地貌分异特点，提出的在流域内以沟谷线、山脊线、分水线和沟缘线相互交叠构成的，具有相对均一坡向、坡度与坡位的地表过程响应单元划分方法，在反映黄土丘陵沟壑区不同地形地貌分区侵蚀空间异质性的基础上，保留了相对完整的地形地貌斑块，能更好地与针对坡面的原位试验观测、基于地块的水土保持措施配置相结合，在水土保持工程、林业生态工程和山地灾害防治等领域具有较广阔的应用前景。

针对所提出的响应单元，在ArcGIS9.2 中，以Hc-DEM 为基础，分别确定的沟谷线、山脊线、分水线和沟缘线的提取方法，以及平均坡度、整体坡向、面积与周长等地形特征的测算技术，可实现地表过程响应单元的快速划分及其地形特征高效提取，从而可克服不规则单元划分通常不易实现的问题。

土壤侵蚀是黄土丘陵沟壑区的典型地表过程，土壤侵蚀模拟与预报是该学科的主要技术难点和热点。基于不规则单元的模型构建和模拟分析，能拓展土壤侵蚀预报、模拟的研究思路，有助于推动该学科领域的发展。

［1］ Beven K.Linking parameters across scales: sub-grid parameterization and scale dependent hydrological models［J］.Hydrological process，1995，9:507-525

［2］ 周大良.基于遥感信息的林冠截流模型［D］.北京:北京大学，1997:28-30

［3］ 任立良.数字流域模拟研究［J］.河海大学学报，2000，28(4):1-7

［4］ Ao T Q，Takeuchi K，Ishidaira H，et al.Development and application of a new algorithm for automated pits removal for grid DEMs［J］.Hydrological Sciences Journal，2003，48(6):985-997

［5］ 周云贵，刘瑜，邬伦.基于数字高程模型的水系提取算法［J］.地理学与国土研究，2000，16(4):77-81

［6］ 左大康.现代地理学辞典［M］.北京:商务印书馆，1990:222

［7］ 秦伟.北洛河上游土壤侵蚀特征及其对植被重建的响应［D］.北京:北京林业大学，2009:96-97

［8］ 秦伟，朱清科，赵磊磊，等.基于RS 和GIS 的黄土丘陵沟壑区浅沟侵蚀地形特征研究［J］.农业工程学报，2010，26(6):58-64