大中流域尺度分布式土壤侵蚀学坡长提取的数据分片方法

曹夏雨，杨勤科，2†，兰 敏，王春梅，2

(1.西北大学城市与环境学院,710069,西安；2.中国科学院水利部水土保持研究所,712100,陕西杨凌；3.中煤西安设计工程有限公司,710054,西安)

土壤侵蚀预报研究中,坡长是指地表从径流源点开始,沿径流线到达到坡度减小直至有沉积出现地方之间或者到一个明显的沟道(自然的沟道或人工渠道)之间的水平距离[1]。针对流域内每个点(DEM栅格)提取的坡长称为流域分布式土壤侵蚀学坡长(简称流域坡长)[2]。

流域坡长的计算,是基于USLE(或RUSLE、CSLE)进行流域和区域尺度土壤侵蚀评价制图的最关键问题之一[3-5]。其计算受到诸多因素影响,主要包括DEM分辨率、流向算法和基础数据的范围等[6]。关于数据分片方面,研究者强调须以完整的流域为单元,以便能在完整的径流路径上计算坡长,以避免完整的坡面被数据边界截断,使数据边界处的部分坡长被遗失[7-8]。事实上,这种认识是针对小流域或大中流域较粗分辨率(数据量较小)的坡长提取而言的；而对于较高分辨率、较大面积(数千到上万个标准图幅组成)的情况下,坡长提取的数据如何分片(数据单元),却未见讨论。笔者对流域坡长提取的数据单元划分方法做出探讨,以便避免或减少边际效应前提下,在较大区域范围内快速提取坡长,满足区域土壤侵蚀评价与制图的需要。

1 研究区概况

本研究分别在地形比较平缓的东北漫岗丘陵区和地形较陡的黄土丘陵区展开研究。东北漫岗丘陵位于大小兴安岭和长白山的山前台地,由低平漫川和波状起伏的丘陵组成,本研究样区选在黑龙江克山和拜泉一带。黄土丘陵区位于黄河中游流域,是我国乃至世界上土壤侵蚀最严重的地区,本研究样区选在位于典型黄土丘陵区的陕西安塞。据25 m分辨率坡度统计,东北样区和黄土样区域的平均坡度分别是2.0°和21.5°。2个研究区地势特征见图1。

图1 东北漫岗丘陵样区(左)和西北黄土丘陵样区(右)地形Fig.1 Landform of the Northeast China undulated hill site(Left,NE site for short below)and the Northwest China loess hill site(Right,NW site for short below)

本研究以2个研究区各选覆盖9个1∶25万标准图幅的DEM为基础。所用DEM数据系根据1∶5万DLG(包括等高线、高程点和水系3个专题层),在ANUDEM支持下建立的水文地貌关系正确的DEM(hydrological correctly DEM,Hc-DEM)[9-10],分辨率为25 m[11]。

2 研究方法

本研究基于坡长提取的一般原理,利用数字地形分析方法,结合对水系、流域和坡长统计特征的分析来完成。

2.1 流域坡长概念模型

流域坡长的计算,从局部高点(位于分水地带)开始,顺径流路径向下坡方向对每个栅格单元的长度不断累加,直到发生泥沙沉积的地方或明显的沟道(图2)[7,12]。从典型小流域坡长提取结果看,分水线往下,坡长的值逐渐增加(图3,由白到黑坡长增加)；所以流域坡长提取的基本数据单元是DEM上可以辨认的、具有水文地貌意义的最小流域—最小有效流域(WME),只要数据以WME的边界或相应级别的流水线为边界,即可提取完整的坡长。

图2 坡长提取技术流程Fig.2 Flowchart of extracting slope length

图3 坡长与分水线Fig.3 Slope length and watershed boundary

2.2 基于WME的坡长提取数据单元

据坡长概念模型,较大范围矩形数据块划分的WME,其最外围流域边界(或流水线)的连线(BNDWME)将接近于矩形(图4中粗黑色折线)。BNDWME之内为一组完整小流域,之外为一个不完整小流域带,也是不完整坡长带。假定图4中内部黑色矩形为某标准图幅的图廓,则只要对标准图幅向外缓冲宽度大于不完整坡长带宽度(Dwr)后,即可以规则图幅(包括缓冲带)为数据单元计算得到完整的坡长。Dwr通过下面3种方法确定。

图4 规则数据块和不规则流域边界Fig.4 Rectangular data and incomplete watershed zone

1)直接量算:即在一个工作区选择典型样区,提取最小有效流域,勾绘近似矩形的BNDWME,量取其最大内接和最小外接矩形间的最大距离(Dw/m)。

2)基于沟壑密度推求:在DEM上提取与WME适应的河流,统计河流总长度,用式(1)推算平均坡长[13]。

3)基于坡长实际值推求:以流域为单元提取的坡长,对坡长做统计分布分析,用累计频率99.99%作为最大值(LENmax/m)。

最终,规则数据向外扩充(缓冲)的数据宽度(Dwr/m)为上述各值的最大值。

2.3 坡长提取与对比

以流域为单元和以图幅为单元(边界缓冲后)提取坡长,并将结果进行对比。主要步骤如下。

流域划分:经过填洼、流向计算、累计流量计算等步骤,通过在ARC/INFO workstation下编程,选择合适汇水面积阈值,提取最小有效流域,用以确定不完整坡长带宽度；然后将样区9个1∶25万标准图幅的DEM拼接为一个数据单元后,将其划分为20个左右的中等流域,用以作为流域坡长提取的基本单元。

坡长提取:分别以流域和以标准图幅(以Dwr为宽度向外缓冲)为单元,利用笔者开发的坡度坡长提取工具软件提取坡长[2,14]。

对比分析:对流域为单元和图幅为单元的坡长,从提取结果的图形、统计特征和工作效率等方面进行对比分析。

3 结果与分析

3.1 以流域为单元的坡长

3.1.1 流域划分与坡长提取 东北和黄土样区1∶25万标准分幅、25 m分辨率DEM划分的较大流域单元如图5示。为了提取正中间矩形图框(1∶25万标准图幅)部分的坡长,将涉及的中等流域(图5中的11—14,21—24)边界,经过缓冲后切割出4个流域的数据块,然后逐一提取坡长,拼接后从中间切割出工作区的坡长(图6,均为1∶25万标准图幅的左上角)。以流域为单元的坡长提取工作流程如图7。有6个基本步骤。

图5 东北漫岗丘陵样区(左)和西北黄土丘陵样区(右)地貌中等流域单元(面积分别为11.09万km2和13.42万km2)Fig.5 Medium watershed units in the NE site(Left)and NW site(Right)(Areas are 11.09×104km2and 13.42×104km2respectively)

3.1.2 坡长特征分析 以下从空间格局和统计分布(图6、表1)2个方面,对样区坡长做简单分析。

图6 东北漫岗丘陵样区(上)西北黄土丘陵样区(下)坡长Fig.6 Slope length of NE site(Up)and NW site(Down)

图7 以流域为单元的坡长提取工作流程Fig.7 Working flowchart for extracting slope length by a watershed as unit

1)东北漫岗丘陵区坡长:坡长平均值479 m,中值320 m,最大值3 570 m。空间上坡度比较平缓的漫岗地坡度最长,最大值出现在漫岗与漫川的转折部位；缓坡丘陵坡度较短。

2)黄土丘陵区的坡长:样区平均坡长86.1 m,中值68.4 m,最大值1 223.9 m。空间分布是,从丘陵顶部向下增加,到坡面向沟道或川地交界处理最长达到最大。

3.2 以图幅为单元的坡长

3.2.1 坡长不完整带的宽度 不完整流域带宽度(Dw):根据汇流面积—河流密度关系分析,并结合对河流—DEM表面关系的观察,东北漫岗丘陵区和黄土丘陵区在25 m分辨率DEM上汇水面积阈值分别为1.0和0.5 km2。根据前述方法量算,在东北漫岗丘陵区和黄土丘陵区,Dw的值分别为4 791.6 m和2 776.5 m。

表1 2个研究区坡长基本统计特征Tab.1 Statistical characteristics of two sites'slope lengths

坡长实际值的统计分布:利用前述以流域为单元提取的坡长,用流域边界切割剔除流域周遍不完整坡长部分,用累计频率99.99%作为最大值(LENmax),结果在黄土和东北样区坡长最大值分别为491.4和2 124.5 m(表3)。

表2 基于沟壑密度的坡长统计Tab.2 Slope length statistics based on gully density

表3 实际坡长统计表Tab.3 Statistics of actual slope length

坡长数据缓冲宽度(Dwr):据以上统计,黄土样区和东北样区不完整坡长带宽度分别是2 776.5和4 791.6 m(表4),实际应用中,可设置为3和5 km。

表4 矩形图幅周边不完整流域带参数统计Tab.4 The statistic parameters of the incomplete watershed around the rectangular map sheet

3.2.2 坡长提取流程 对标准图幅向外缓冲3 000(黄土样区)和5 000 m(东北样区),然后在LS_Tool系统中提取流域坡长,工作流程见图8。与图7相比,工作步骤由6步减少到3步。尤其重要的是,减少了流域划分、按流域切割数据和提取坡长后的再拼接,从而使工作效率大为提高。

图8 标准图幅坡长提取流程Fig.8 Flowchart of extracting slope length from standard map sheet

3.3 2种方法比较

把相同样区2种方法提取的坡长做差值运算,结果全部等于0,表明坡长的值完全一样。由于表面一致,因而使其统计直方图(比较上部位两条线为黄土样区坡长累计频率)也完全一致(图9)。有此可见,以标准图幅为单元提取坡长是可行的。

图9 流域为单元和标准图幅为单元坡长直方图Fig.9 Slope length histogram while a watershed as an unit and a standard map sheet as an unit

据我们工作记录,主要工作步骤的消耗时统计见表5。以流域为单元的坡长提取总耗时46 min,以标准图幅为的单元提取坡长总耗时8 min,为常规方案的17.4%；所以,以标准图幅为单元提取坡长,效率远大于以中等流域为单元的工作方式。

表5 流域坡长和规则数据块坡长计算时间比较Tab.5 Comparison of calculation times between two calculating methods(Watershed slope length and rectangular data) min_

4 结论与讨论

流域坡长提取的数据单元,是一个长期被忽视的问题,本研究结论及有关问题归纳如下。面向大区域、多各图幅和海量数据的坡长提取数据单元研究,可得出如下初步结论。

1)根据流域分布式土壤侵蚀学坡长的含义和提取算法,本文提出流域坡长概念模型—流域坡长提取的基本数据单元是DEM上可以辨认的、具有水文地貌意义的最小流域。据此,对于一个由多个图幅构成的较大工作区,只要对将每个图幅范围向外扩展一定距离(Dwr),即可直接用规则图幅为单元逐幅提取坡长。

2)基于对不完整流域带宽度(Dw)、沟壑密度和坡长均值、和坡长实际值的统计(LENmax),在地形较平缓东北样区和地形较陡黄土样区,标准图幅向外缓冲的宽度分别为2 776.5和4 791.6 m。

3)基于中等流域和标准图幅的坡长提取,其结果完全一致。工作步骤分别为6步和3步,总耗时间分别为45和8 min。所以,基于标准图幅提取大范围坡长,工作步骤大为简化,效率大为提高。

笔者论述了以标准图幅或规则矩形块为数据单元提取流域坡长的基本原理,也提出了标准图幅为单元提取坡长应该缓冲的宽度推算方法；然而以下问题有待深入研究:

1)坡长提取运算的数据量上限:经试验或者理论推导,推算出在给定的硬件条件下利用LS_Tool工具提取坡长时,可顺利运行数据的DEM最大数据量(矩阵行列数),以便更有效划分数据单元。

2)不完整坡长带的推算:上述方法虽然可推算出不完整坡长带的宽度,但是原则上要根据工作区域的地形特征做出率定和优化。今后的工作中,应在足够样本基础上,选择合适参数,建立比较通用的模型以便能推算出不完整坡长带宽度值。

5 参考文献

[1] SMITH D D,WISCHMEIER W H.Factors affecting sheet and rill erosion[J].Trans.Am.Geophys.Union,1957.38(6):889.

[2] 杨勤科,郭伟玲,张宏鸣,等,基于DEM的流域坡度坡长因子计算方法研究初报[J].水土保持通报,2010.30(2):203.YANG Qinke,GUO Weiling,ZHANG Hongming,et al.Method of extracting LS factor at watershed scale based on DEM[J].Bulletin of Soil and Water Conservation,2010,30(2):203.

[3] MOORE I.D,BURCH G J.Physical basis of the lengthslope factor in the universal soil loss equation[J].Soil Science Society of America Journal,1986.50(5):1294.

[4] WILSON J P.Estimating the topographic factor in the universal soil loss equation for watersheds[J].Journal of Soil and Water Conservation,1986.41(3):179.

[5] MOORE I D,WILSON J P.Length-slope factors for the revised universal soil loss equation:Simplified method of estimation[J].Journal of Soil and Water Conservation,1992.47(5):423.

[6] 王程,陈正江,杨勤科,等.流域分布式坡长不确定性的初步分析[J].水土保持研究,2012,19(2):15.WANG Cheng,CHEN Zhengjiang,YANG Qinke,et al.Analysis on uncertainty of DEM derived watershed Dia[J].Research of Soil and Water Conservation,2012,19(2):15.

[7] HICKEY R,A SMITH,AND P JANKOWSKI.Slope length calculations from a DEM within ARC/INFO GRID[J].Computers,Environment and Urban Systems,1994.,18(5):365.

[8] GRUBER S,PECKHAM S.Land-surface parameters and objects in hydrology[M]//HENGL T,REUTER H I.Geomorphometry:Concepts,Software,Applications,Series Developments in Soil Science Vol.33.Amsterdam:Elsevier,2009:207.

[9] YANG Qinke,McVICAR T R,Van NIEL T G,et al.Improving a digital elevation model by reducing source data errors and optimising interpolation algorithm parameters:an example in the Loess Plateau,China[J].International Journal of Applied Earth Observation and Geoinformation),2007,9(3):235.

[10] 杨勤科,师维娟,McVicar T R,等.水文地貌关系正确的DEM建立方法[J].中国水土保持科学,2007,5(4):1.YANG Qinke,SHI Weijuan,McVICAR T R,et al.On constructing methods of hydrologically correct DEMs[J].Science of Soil and Water Conservation,2007,5(4):1.

[11] 国家测绘局.基础地理信息数字产品1∶10 000、1∶50 000数字高程模型:CH/T 1008—2001[S].北京:中国标准出版社,2002:3.National Administration of Surveying.CH/T 1008- 2001 Digital products of fundamental geographic information 1∶10 000,1∶50 000 digital elevation models:CH/T 1008-2001[S].Beijing:Standards Press of China,2002:3.

[12] HICKEY R.Slope angle and slope length solutions for GIS[J].Cartography,2000.29(1):1.

[13] 陆中臣,贾绍凤,黄克新,等著.流域地貌系统[M].大连:大连出版社,1991:32.LU Zhongchen,JIA Shaofeng,HUANG Kexin,et al.Watershed landform system[M].Dalian:Dalian Publishing House,1991:32.

[14] 张宏鸣,杨勤科,刘晴蕊,等,基于GIS的区域坡度坡长因子提取算法[J].计算机工程,2010,36(9):246.ZHANG Hongming,YANG Qinke,LIU Qingrui,et al,Regional slope length and slope steepness factor extraction algorithm based on GIS,[J].Computer Engineering,2010,36(9):246.