黄土高原历史时期沟谷侵蚀量计算方法探讨

解哲辉，崔建新，常 宏

（1.中国科学院地球环境研究所，西安 710075；2.陕西师范大学 西北历史环境与经济社会发展研究中心，西安 710062；3.中国科学院大学，北京 100049)

黄土高原历史时期沟谷侵蚀量计算方法探讨

解哲辉1,3，崔建新2，常 宏1

（1.中国科学院地球环境研究所，西安 710075；2.陕西师范大学 西北历史环境与经济社会发展研究中心，西安 710062；3.中国科学院大学，北京 100049)

历史时期沟谷侵蚀量的计算多采用传统实地测量或计算方法，而现代地理信息技术的发展为沟谷侵蚀量计算提供了更为简便、快捷、高效的方式。等高线图形概括法即为一种将传统制图学与现代地理信息技术相结合的方法。本文以神木县东山旧城冲沟为切入点，采用等高线图形概括方法，基于30 m分辨率ASTER DEM，利用Arcgis10平台计算历史时期沟谷侵蚀量，并将其与基于野外测量计算侵蚀量进行对比。结果对比发现，基于野外测量的沟谷侵蚀量约占基于DEM侵蚀量的89.11%，两者之间存在一定差距，但是可被接受的。进一步分析显示，等高线简化前后，DEM平均坡度、等高线长度和表面积均发生变化，平均坡度、等高线长度变化率分别为2.78%和4.3%，表面积变化不是特别明显。三维分析显示，简化后等高线更加平滑，平均坡度趋于平缓，坡度的分布更为均匀。总体而言，等高线图形概括方法在沟谷侵蚀量计算方面具有较高的可靠性，对更大时空尺度内的沟谷侵蚀量计算具有借鉴意义。

等高线图形概括；数据挖掘；沟谷侵蚀量；DEM

沟谷侵蚀是一种侵蚀强度非常高的线状侵蚀，并且断面积足够大，对地表形态形成永久改造（Torri and Borselli，2003）。沟谷侵蚀是由于地表径流对土壤的不断冲刷而形成的，沟谷的不断发展，改变了原有陆地下垫面特征，而下垫面则又会通过影响径流而对沟谷侵蚀产生作用，如此形成恶性循环。沟谷侵蚀引起的泥沙堆积会阻塞河道，不仅对中下游地区的地形、水文条件，而且会对社会经济发展和人类生产生活活动产生深远影响。沟谷侵蚀量是衡量沟谷侵蚀程度的一个重要指标。开展沟谷侵蚀量研究，不但能够加深对沟谷侵蚀程度的直观认识，也可以直接为研究区的生产生活活动提供参考依据，还可以为沟谷侵蚀的防治提供数据支持。因此如何计算黄土高原小流域沟谷侵蚀量，成为分析沟谷形成演化过程中重要的基础性课题。

沟谷侵蚀是土壤侵蚀的一种，目前大部分的已有研究侧重于研究黄土高原各个因子与土壤侵蚀的相互关系，包括降水条件（Xu，2005；Wei et al，2007）、土壤特性（赵晓光和石辉，2003；贾志军等，2006）、地形特征（陈正发等，2010）、土地利用方式（Chen et al，2001；Wei et al，2007 ）以及植被条件（Chenet al，2007）等对侵蚀的影响。近十年来，部分学者对历史时期黄土高原的土壤侵蚀及沟谷侵蚀进行了相关研究（桑广书等，2003；姚文波，2010），尽管如此，但是关于较长历史时期沟谷侵蚀的研究依然相对较少,尤其是对历史时期所形成的沟谷其侵蚀量计算方面，大都采用传统横截面计算法或空腔体积法。传统横截面计算方法误差小，精度高，但是野外工作周期长，费时费力，并且受自然条件限制因素很多。Torri and Borselli（2003）试图通过方程对沟谷侵蚀量进行计算，其研究结果认为沟谷形态阈值不能通过简单的方法进行解释，因此也不能通过简单的方程进行估计。近年来，随着地理信息系统的发展成熟，其被越来越多的应用于土壤侵蚀研究。沟谷侵蚀量的计算方法也有了很多新进展。王辉等（2008）利用黄土高原两期航片和地形图数据，借助两期DEM的高程差图层计算沟谷侵蚀量与沉积物的产生量。但是其得出的沟谷侵蚀量与前人研究结果不一致，他将其归因为沟谷侵蚀因素之间的相互影响和研究区之间的尺度问题。刘万青等（2011）等利用等高线图形概括结合数据挖掘技术，对黄土高原马家沟流域进行沟谷侵蚀量估算，为以后不同空间尺度上的侵蚀量计算提供了一定的参考。胡文生等（2008）利用数字摄影测量方法估算半干旱区小流域沟谷侵蚀产沙量，并将实测值与估算结果进行对比评价。

本文立足于前人已有研究成果，选取陕北黄土丘陵沟壑区1座军事古城堡——东山旧城内的沟谷为研究对象,对城内沟谷侵蚀量进行估算。采用基于DEM的等高线图形概括方法以及基于实测数据的计算方法来估算侵蚀量，并将两种估算方法所得结果进行对比分析。该研究旨在将历史地理学与现代地理信息技术进行结合，为小流域历史时期沟谷侵蚀量估算提供参考,指出需要注意的影响因素，并对未来研究进行了展望。

1 研究区特征

本文研究区位于黄土丘陵沟壑区。该区黄土层深厚，具有典型代表性，土壤侵蚀方式独特，过程迅速，面蚀、沟蚀、潜蚀的共同作用使得该地区成为我国甚至全球土壤侵蚀最为严重的地区之一。该区属于中温带半干旱内陆季风气候，冬季漫长寒冷，夏季降雨集中，年际变率大。降水集中在7、8、9三个月，暴雨频率高发。东山旧城位于38°49′ N，110°30′ E，神木县向东约1.5公里的东山上，海拔901～1337 m，地势西高东低。城堡面积约为0.12 km²，城墙西侧发育有冲沟，城内也分布有细沟及雨水冲刷形成的坑洼地。城内地面不平，植被覆盖类型为草本及乔木。

如图1所示，城墙轮廓内共有三条沟。沟1为石沟，沟壁为沉积岩，具有水平沉积层理，深度为50 m左右。城堡西南城墙已失落，且与深沟相邻，虚线代表推测的西南城墙，并不确定。由于西南城墙的不确定性，无法断定沟1是存在于城堡修筑之前还是城废弃之后。因此沟1的侵蚀量不在此次计算范围内，只对沟2和沟3的侵蚀量进行对比分析。

图1 研究区地理位置（修改自裴新富(1991)）Fig.1 The location of the study area (after Pei (1991))

2 数据获取

本文采用两种方法来估算侵蚀量。第一种为实测方法，第二种为基于DEM数据的等高线数据概括方法。而数据源的获取是本项研究的基础。

2.1 实测数据获取

2012年9月，对东山旧城城墙范围进行了第一次信息采集，利用手持GPS采集城墙点位信息，及对城内冲沟进行线要素采集。2013年8月，对东山旧城冲沟进行第二次信息采集，沿沟缘线平均每5 m记录GPS点。两次信息采集结果为：1）确定了东山旧城的城墙范围。2）主要沟谷参数采集。由于沟谷地形的特殊性及现实条件的制约，沟谷宽度的测定主要通过采集高精度亚米级GPS数据，然后在地理信息系统软件中测量得到。沟谷深度的获取方式有两种：自然条件较好，可以下到沟底的部位，利用沟上与沟底GPS高程相减得到；不可以进入沟底的则采用钢丝测距绳直接测得。沟壁坡度也相应分为两种：地形制约严重的部位，利用沟谷深度与坡长之比的正弦得到；一般部位，则使用地质罗盘仪直接测得。坡长的获取方式也根据现实条件分为两种：能够进入沟底的采用GPS获取；不能够进入沟底的直接使用钢丝测距绳测量。钢丝测距1 m最大允许误差：±4.0 mm，任意5 m最大允许误差：±15 mm。数字采集所用工具为手持GPS（GEO-XT-2005和GEO-XT 2008），水平精度为0.8±0.2 m，垂直精度1.2±0.3 m。

2.2 DEM数据获取

30 m空间分辨率ASTER DEM数据通过中国科学院数据平台获得。由于其较易获得，再加上东山旧城地区历史时期沟谷在DEM中具有较强的可识别性，因此选其作为本文的数据源。

3 计算过程

3.1 根据实测数据计算的沟谷侵蚀量

根据野外调查数据，采取横截面法进行计算,沟头部分采用三棱锥法进行计算。沟道部分计算方法是横截面面积乘以横截面间距，然后求和计算得出。计算公式如下：

式中：V为沟谷侵蚀量； S为横截面面积；D为横截面间距。横截面间距是相邻两个横截面之间的距离，由于实际测量情况受到沟谷地形及其他限制因素，间距长度存在差异，沟谷横截面平均间距为20 m左右。

若沟谷为“U”形谷，横截面为梯形，其面积则用梯形面积公式计算得出：

其中：S为横截面面积；L1为截面梯形的上底；L2为截面梯形的下底；H为代表沟谷截面深度，计算时取两侧深度的均值。

若沟谷为“V”形谷，横截面为三角形，面积则用三角形面积公式计算得出：

其中：S为横截面积；L1、L2为沟谷两侧沟壁长度；A为沟壁两侧之间的夹角，等于180°减去两侧沟壁倾角之和。

3.2 基于DEM的侵蚀量

3.2.1 等高线图形概括法和数据挖掘方法

等高线图形概括法（刘万青等，2011）是将传统地图学与现代计算机技术相结合的方法，其原理即将沟谷部位通过对等高线进行简化来达到“虚拟填充”，从而从地理学上将沟谷“恢复”到发育之前的状态。等高线简化的依据是沟缘线，使得等高线沿着沟谷两侧等高线的延伸方向自然伸展。等高线简化应当遵循三个原则（刘万青等，2011）：（1）由低级到高级，由低海拔到高海拔；（2）应对构成沟谷的所有等高线进行操作，且简化后空间分布协调；（3）应该根据沟谷两侧等高线自然延伸方向和曲率进行简化。数据挖掘法主要是地形表面分析法和地学统计分析法。地形表面分析法主要用于生成等高线；地学统计分析法用于沟谷侵蚀量的统计运算。

3.2.2 计算过程

以研究区30 m空间分辨率ASTER数字高程模型（DEM）为数据源。在Arcgis10平台下，对DEM进行操作，主要分为以下几个部分：（1）根据两次野外GPS点位，连接并确立东山旧城原始城堡界限。（2）根据GPS点位勾绘堡内冲沟，结合第一次野外信息对沟缘线进行修改更正。（3）基于30 m分辨率ASTER DEM，以城堡为中心矩形裁剪出小块DEM，然后基于裁剪的DEM提取等高线，等高距为5 m，生成TIN。投影采用UTM 49 N。（4）基于Arcgis3D分析模块进行沟谷侵蚀量的计算，图2即为具体的处理流程。等高线简化后，地图表现为沟谷小弯曲的消失，进而可以通过计算机方法对目标沟谷进行侵蚀量计算。数字计算结果由“Arcgis10软件→3D Analyst工具→功能性表面→表面体积”工具计算得出。利用等高线简化前后生成的TIN体积差计算沟谷侵蚀量。

4 结果与分析

4.1 沟谷侵蚀量结果

由图3可以看出，等高线简化后，原来的沟谷区域已经被进行了“填充”而成为台地，“填充”后等高线沿其自然曲率向两侧伸展。由表1分析得到，基于野外测量的沟谷侵蚀量约占基于DEM侵蚀量的89.11%。若以野外测量结果为参照，则基于DEM的计算结果的可信度约为87.78%。从统计学上而言，基于DEM的等高线图形概括方法对历史时期沟谷侵蚀量的计算结果是比较可靠的。两者之间相差20839.41 m3，其约占实测结果的12.22%。可以认为，两种计算结果之间依然存在一些差距。

基于Arcgis10水文分析计算，根据已有研究（陈涛等，2008；张康聪，2010），将汇流阈值设为150个像元时，得出东山旧城地区沟壑密度为2.74 km·km−2，根据土壤侵蚀分类分级标准SL190-2007（2008），东山旧城所在小流域为中度侵蚀。然后根据等高线图形概括方法所得侵蚀量计算得出东山旧城地区沟谷侵蚀模数为3784.89 t·(km2·a)−1，对照土壤侵蚀标准水力侵蚀分级亦为中度侵蚀。最后，根据野外测量数据所得结果计算得出的东山旧城地区的侵蚀模数为3372.86 t·(km2·a)−1，对照土壤侵蚀标准水利侵蚀分级同样为中度侵蚀。这从侧面反映出应用等高线图形概括方法对历史时期沟谷侵蚀量的计算结果是可信的。

图2 数据处理流程图Fig.2 The flow chart of data processing

图3 东山旧城地区等高线简化前(a)后(a)对比Fig.3 The contour before (a) and after (b) simplification of Dongshan castle

4.2 简化前后地表形态分析

在对等高线进行虚拟填充后，由于沟谷的消失，沟壁表面积的减少，必然会导致简化后表面积小于简化前；由于利用等高线简化，形成了对地表形态的虚拟改造，简化后的沟谷体积相较于简化前必然会增加。由表2可以分析得到，等高线简化前后，DEM的平均坡度和等高线长度均发生了不同程度变化。其中DEM平均坡度前后相差0.36°，标准差原来为8.04，等高线简化后为7.67。相较于简化前，坡度降低了2.78%，标准差降低了4.6%。平均坡度趋于平缓，坡度的分布更为均匀。等高线长度则前后相差897.78 m，相比原来缩短了4.3%，等高线形态更加平滑。从总体来看，简化达到了虚拟“填充”的效果和标准。由于该区地势北高南低，且对于海拔最高和最低地形部位并没有进行等高线简化操作，因此简化前后，该区最高海拔和最低海拔没有发生变化，相对应的最小坡度和最大坡度也没有发生变化。从图4中可以很直观看出等高线简化前后沟谷形态的变化。经“填充”处理之后，原来的沟谷部位已经被浅色的部分所填充，浅色部分体积即为基于DEM计算得出的沟谷侵蚀量。

表1 实验结果和实测结果对比Table 1 Results by experiment and measurement

表2 简化前后对比Table 2 Contrast before and after simplification

图4 沟谷侵蚀结果三维立体图Fig.4 The 3-D visual picture of Dongshan gully erosion

4.3 原因分析

已有研究成果表明（汤国安等，2006），DEM分辨率会对侵蚀结果造成影响。同一地区,由于土壤侵蚀的计算方法不同，或者DEM尺度因素，造成其最终得到的侵蚀速率存在差异。Ren et al（2011）研究发现，小比例尺得到的土壤侵蚀面积要大于大比例尺DEM得到的结果。因此，在利用等高线图形概括对沟谷进行虚拟填充过程中，会使得沟谷地区与实测相比被“夸大”，从而造成基于野外测量的沟谷侵蚀量小于基于DEM计算侵蚀量。除此之外,由于东山旧城地区沟谷系统的复杂性,在对沟谷进行实测的过程中,容易受到路线及侧视的影响，不免存在一定误差,这些因素都会影响野外考察数据的结果。

本文数据源为30 m分辨率DEM数据。虽然沟谷形态在等高线中能大体识别，但是沟谷的具体信息还是表现的比较粗糙。从图3（a）中可以很清晰地识别出沟1、沟2和沟3。对比图1中沟谷形状与图3中等高线可以看出，沟3西北角的支沟部分存在部分差异，这可能是由于DEM栅格较大，因此细节部分不能得到很好反映。

5 结论

（1）基于DEM的等高线图形概括方法，其计算侵蚀量与基于野外测量的计算侵蚀量之间差异较小；在判断沟谷发育程度方面，其侵蚀模数与野外测量所得侵蚀模数也保持一致。因此确定该方法在沟谷发育研究中具有较好的可靠性和现实可行性，可以为长时段的沟谷侵蚀研究提供较为可信赖的数据。然而其计算结果与野外测量结果还存在一定误差，这可能与DEM空间分辨率及野外测量误差有关，还需进行深入研究。

（2）该项研究以古城内的沟谷发育作为研究对象，从某种程度上可以为较大时空尺度的沟谷侵蚀研究提供方法上的借鉴。长时期以来，历史时期沟谷侵蚀尤其是较大时空尺度的沟谷侵蚀研究，大都依靠人工野外测量或地形图，费时费力。应用等高线图形概括方法可以有效解决这一问题，因此具有很好的实际应用意义。

陈正发, 郭宏忠, 史东梅, 等. 2010. 地形因子对紫色土坡耕地土壤侵蚀作用的试验研究[J]. 水土保持学报, 24(5): 83–87. [Chen Z F, Guo H Z, Shi D M, et al. 2010. Experimental study on effect of topographical factors on purple cultivated slopes land soil erosion [J]. Journal of Soil and Water Conservation, 24(5): 83–87.]

胡文生, 蔡强国, 陈 浩. 2008. 利用数字摄影测量方法估算半干旱区小流域沟谷侵蚀产沙[J]. 地球信息科学, 10(4): 533–538. [Hu W S, Cai Q G, Chen H.2008. Digital photogrammetric techniques as a tool for estimating sediment production of gully erosion in a semiarid catchment of the Loess Plateau, China [J]. Geo-Information Science, 10(4): 533–538.]

贾志军, 蔡强国, 卜崇峰. 2006. 晋西黄土丘陵沟壑区土壤侵蚀能力的试验研究[J]. 水土保持研究, 13(6): 1–3. [Jia Z J, Cai Q G, Bu C F. 2006. Experimental study on soil anti-erodibility in the hilly and gully region of Loess Plateau Western Shanxi [J]. Research of Soil and Water Conservation, 13(6):1–3.]

刘万青, 张超超, 王旭红. 2011. 等高线图形概括在黄土高原沟谷侵蚀量估算中的应用[J]. 水土保持通报, 31(1): 128–131. [Liu W Q, Zhang C C, Wang X H. 2011. Estimate total gully erosion in Loess Plateau using simplified contour map [J]．Bulletin of Soil and Water Conservation, 31(1): 128–131.]

裴新富. 1991.关于黄土高原范围问题[J]. 中国水土保持, 12: 37 – 42. [Pei X F. 1991. On the scope of Loess Plateau [J]. Soil and Water Conservation in China, 12: 37 – 42.]

桑广书, 甘枝茂, 岳大鹏. 2003.元代以来黄土塬区沟谷发育与土壤侵蚀[J]. 干旱区地理, 26(4): 355–360. [Sang G S, Gan Z M, Yue D P. 2003. Gully development and soil erosion in the Loess Yuan region since Yuan Dynasty [J]．Arid Land Geography, 26(4): 355–360.]

汤国安, 刘学军, 房 亮, 等. 2006. DEM及数字地形分析中尺度问题研究综述[J]. 武汉大学学报(信息科学版), 31(12): 1059 – 1066. [Tang G A, Liu X J, Fang L, et al. 2006. A review on the scale issue in DEMs and digital terrain analysis [J]. Geomatics and Information Science of Wuhan University, 31(12): 1059 – 1066.]

王 辉, 王天明, 杨明博, 等. 2008.基于航片的黄土高原丘陵沟壑区沟谷侵蚀定量监测[J]. 应用生态学报, 19(1): 127–132. [Wang H, Wang T M, Yang M B, et al. 2008. Quantitative monitoring of gully erosion in hilly-gully area of Loess Plateau based on aerial Images [J]. Chinese Journal of Applied Ecology, 19(1): 127–132.]

姚文波, 侯甬坚, 高松凡. 2010.唐以来方家沟流域地貌的演变与复原[J]. 干旱区地理, 33(4): 509–517. [Yao W B, Hou Y J, Gao S F. 2010. Geomorphologic changes and restoration of Fangjiagou Drainage Area from early Tang Dynasty [J]. Arid Land Geography, 33(4): 509–517.]

张康聪. 2010. 地理信息系统导论（第五版）[M]. 陈健飞,张筱林,译. 北京: 科学出版社, 313. [Chang Kang-T. 2010. Introduction to geographic information systems (5th Edition) [M]. Chen J F, Zhang X L, translated. Beijing: Science Press, 313.]

赵晓光, 石 辉. 2003. 水蚀作用下土壤抗蚀能力的表征[J]. 干旱区地理, 26(1): 12–16. [Zhao X G, Shi H. 2003. Prescription of soil anti-erosion capability under water erosion [J]. Arid Land Geography, 26(1): 12–16.]

中华人民共和国水利部. 2008. 土壤侵蚀分类分级标准：SL 190-2007 [S]. 北京: 中国水利水电出版社, 8. [The Ministry of Water Resources of the People’s Republic of China. 2008. Standards for classification and gradation of soil erosion: SL 190-2007 [S]. Beijing: China Water & Power Press, 8.]

Chen L D, Huang Z L, Gong J, et al. 2007. The effect of land cover / vegetation on soil water dynamic in the hilly area of the Loess Plateau China [J]. Catena, 70(2): 200–208.

Chen L D, Wang J, Fu B J, et al. 2001. Land use change in a small catchment of northern Loess Plateau China [J]. Agriculture Ecosystem and Environment, 86(2): 163–172.

Ren S M, Yin L, Sun B. 2011. Research on sensitivity for soil erosion evaluation from DEM and remote sensing data source of different map scales and image resolutions [C]. 2011 3rd International Conference on Environmental Science and Information Application Technology (ESIAT 2011)．Procedia Environmental Sciences IO, 1753 – 1760.

Torri D, Borselli L. 2003. Equation for high-rate gully erosion [J]. Catena, 50(2): 449–467．

Wei W, Chen L D, Fu B J, et al. 2007. The effect of land uses and rainfall regimes on runoff and soil erosion in the semiarid loess hilly area, China [J]. Journal of Hydrology, 335(3 – 4): 247–258．

Xu J X. 2005. Precipitation-vegetation coupling and its influence on erosion on the Loess Plateau, China [J]. Catena, 64(1): 103–116.

The discussion of the computing method to historical gully erosion module in Loess Plateau

XIE Zhe-hui1,3, CUI Jian-xin2, CHANG Hong1
(1. Institute of Earth Environment, Chinese Academy of Sciences, Xi'an 710075, China; 2. Center for Historical Environment and Socio-Economic Development in Northwest China of Shaanxi Normal University, Xi'an 710062, China; 3. University of Chinese Academy of Sciences, Beijing 100049, China)

Usually, in terms of historical gully erosion study, historical geographers often prefer to field measurement or room observation based on terrain map in order to obtain the erosion module. So it seems outdated, especially with the rapid development of the modern geography information technology. However, simplified contour maps which combine traditional mapping with modern geography information system together may be regarded as a new solution to the gully erosion study. The main purpose of this study is aiming to calculate gully erosion volume more rapidly, conveniently and efficiently. The study area in this paper sites in the Shenmu County, 1.5 km away from the city, called Dongshan Castle which was abandoned in the Song dynasty in 1444 AD. Based on ASTER DEM with 30 m resolution combined with simplif ed contour maps and data mining method that are run on Arcgis10 platform, the gully erosion module is calculated and then is made contrast with the result obtained from field measurement. Consequently the results show as follows: (1) The gully erosion volume for DEM based and f eld measurement are 191431.26 m3and 170591.85 m3. The later resulttakes about as 89.11% as the former. It demonstrates the reliability of the DEM based calculation. However, apparently there is indeed error existed and some subtle factors should be taken into consideration when the method is applied into estimating gully modules in small watershed. (2) Before and after simplif cation, average slope, contour length are 12.92°, 20670.11 m and 12.56°, 19772.33 m respectively. As for the average slope, it shrinks 0.36° about 2.78% of the value before simplified. The standard deviation of average slope reduces from 8.04 to 7.67 with the rate as much as 4.6%. The contour length decreases as much as 897.78 m, which takes 4.3% approximately. Given the gully shape, the rate of Surface area is not particularly signif cant. All of the changes suggest that after simplif ed, the contour turns smoother and average slope is leveling off and gradient distribution is more uniform. (3) Given the rationality and feasibility of the simplif ed contour map method, it may provide a new solution to calculate gully erosion module and is probable to supply reliable data for the historical gully erosion study.

contour simplif ed method; data mining; gully erosion module; DEM

S157

A

1674-9901(2014)01-0016-07

10.7515/JEE201401003

2014-01-11

中央高校基本业务费项目（10SZYB10）；中国科学院重点部署项目（KZZD-EW-04-02）

崔建新，E-mail: cuijx@snnu.edu.cn