基于GIS的黄土高塬沟壑区砚瓦川流域地形特征提取与分析*

方磊,刘文兆,2,李怀有

(1.西北农林科技大学资源环境学院,陕西杨陵712100;2.中国科学院水利部水土保持研究所,陕西杨陵712100;3.黄河水利委员会西峰水土保持科学试验站,甘肃 庆阳745000)

地形形态分布特征对地表径流和土壤侵蚀、植被分布和土地利用、水土保持措施布设等具有重要影响[1],它直接影响着地表物质的迁移与能量的转换,标志着地表基本环境格局的形成与发展。

基于数字高程模型(Digital Elevation M odel,简称DEM)的数字地形分析技术经过40多年的发展,其诸多基础理论问题得到了深入的研究[2],而利用GIS中DEM的自动提取技术已经成为地形因子的主要提取手段。同时DEM的应用也从简单的地形要素提取扩展到更广的范围,如辅助建立水文模型、土壤侵蚀模型和洪水预报模型等。目前国内学者对基于DEM的数字地形分析研究主要是集中在地形因子提取与定量分析以及DEM信息容量等方面。如汤国安、赵牡丹等[3-4]系统研究了坡度分级方法并对陕北不同地貌区进行坡度提取、分级,最终对不同地面坡谱进行了统计分析,刘学军等[5-9]将DEM误差分为DEM模型误差以及空间要素提取算法误差,并在考虑DEM误差自相关的前提下对基于DEM 的空间因素(包括坡度、坡向、河网等)的不同提取算法进行了系统的误差分析,陈楠等对不同尺度DEM做了信息容量比较,定量化研究了信息损失量的规律性[10-11]。但以上诸多研究多集中在陕北丘陵沟壑区,对高塬沟壑区研究较少,尤其是对位于陇东,有天下第一塬号称的董志塬相关研究更是鲜有报道。高塬沟壑区地貌复杂,塬高、沟深、坡陡,地势起伏较大,是水土流失严重地区之一,其中水蚀和重力侵蚀强烈,所以如何正确快速掌握地形特征规律,为当地相关部门尤其是水保部门提供决策依据是亟待解决的问题。

基于以上情况,本研究选择典型高塬沟壑区——砚瓦川流域为研究区域,在众多空间因子中选取了3个最能突出反映该地区地形特征的指标因子,指标选取基本原则就是这些地形因子能从不同角度反映该区域地形特征,独立性强,最终确定坡度、坡长和剖面曲率三个地形因子作为该区域地形特征的主要指标,因为坡度大小直接影响着地表物质的流动及能量转换的规模与强度,坡长能够影响坡面径流的流速和流量,进而影响土壤侵蚀强度,而剖面曲率能够反映出坡度的变化率。最后结合GIS,以数字地形分析技术为手段,建立相应地形因子的数字地面模型(Digital Terrain M odel简称DTM),并统计分析各因子分布特点与规律,望给当地经济建设提供支持,为后续相关研究提供参考。

1 材料与方法

1.1 研究区域概况

砚瓦川流域系泾河水系,位于甘肃省庆阳市西峰区境内 ,107°37′-107°55′E,35°31′-35°44′N,流域中轴线呈NW-ES走向,总面积约368 km2,该区属半湿润季风气候区,年均温为8.1℃,年均降水量为579 mm,最多年降水量为852.8 mm(1964年),最少年降水量336.7 mm(1986年),7-9月约占全年降水量的60%左右,并多以暴雨形式出现,大量降水以地表经流的形式汇入江河,发生暴洪的机率大。流域以黑垆土为主,占流域面积的70%,主要分布在塬面,土层厚约250 cm,其上松下粘,耕性良好,适耕期长。流域地面海拔高度为 940～1 409m,平均高程 1 243 m,平均沟壑密度2.17 km/km2,沟谷是侵蚀最剧烈的部位,重力侵蚀频繁发生,对区域经济发展、居住安全以及下游水利工程正常运行带来不利影响,水土保持任务艰巨。

1.2 数据基础及预处理

目前制作DEM的主要方法有两种,即通过遥感立体像得到DEM,或是基于数字地形图通过内插的方法得到DEM,后者是迄今为止最为成熟和经济实用的DEM建立方法,本文综合考虑该区域的数据可获性,选择后者建立本研究区域的DEM,所使用到的基础数据为覆盖研究区域的4幅1∶5万国家标准地形图,基本参数为:54北京坐标系,56黄海高程系,等高距为20 m,使用投影为高斯-克吕格投影。

一般来说,地形图扫描后都有一定变形,必须进行预处理,也就是扫描纠正,并对纠正后的栅格图进行检查,才能得到矢量化准确的数字等高线。本研究采用A rcGIS 9.2提供的Georeferencing工具进行图像纠正,先在每幅地形图上选取12个均匀分布的控制点,为方便计算,本研究选取的都是格网交点,然后建立链接表,并逐个对控制点进行误差检查,合格后对图像进行重采样,最后得到经过精确纠正的扫描图。接下来使用Geoway 3.5对该扫描图进行等高线、高程点、河流、湖泊以及流域边界线的矢量化,通过要素拼接最终得到砚瓦川流域完整的数字地形图。

1.3 建立水文增强型DEM

DEM是一组有序数据阵列对地表高度空间变化的定量表现。是一定范围内规则格网点的平面坐标(X,Y)及其高程(Z)的数据集,它主要是描述区域地貌形态的空间分布,可通过等高线或相似立体模型进行数据采集(包括采样和量测),然后进行数据内插而形成的。根据DEM对水文特征和地貌关系的表现能力,可将DEM分为2类,即水文增强型DEM(hydrologically enhanced DEM s/drainage enforced DEM s)[12],也称为水文地貌关系正确DEM(hydrologically correct DEM s,H c-DEM)[13],以及非水文增强型DEM。本研究采用澳大利亚国立大学Hutchinson教授开发的ANUDEM软件[13-14]来建立水文增强型DEM(H c-DEM),考虑到本研究区域为塬区,相对于丘陵沟壑区地形平坦,用ANUDEM软件插值时设置的三个关键参数[15]分别为:分辨率10m,第二糙度系数0.5,迭代次数35,最终得到水文增强的DEM:H c-DEM,图1是基于H c-DEM建立的本研究区域的三维示意图。

1.4 地形指标计算

坡度是指地球表面任一点通过该点的切平面与水平地面的夹角,表示地表面在该点的倾斜程度,是一个重要的地形因子。坡度大小直接影响着地表物质的流动及能量转换的规模与强度。坡度的计算如公式(1),其实现过程是在3×3的DEM 栅格窗口中进行的,窗口在DEM数据矩阵中连续移动后完成整幅图的计算工作。

式中:fx,fy——南北和东西方向上的偏导数。目前对 fx,fy的计算有二阶差分、三阶不带权差分(Sharp nack算法)、三阶反距离平方权差分(Horn算法)、三阶反距离权差分、Frame差分、简单差分等6种常见算法,对于这6种算法的分析比较,国内学者已经做了大量的研究[5,16-17],最终,本研究确定使用三阶反距离平方权差分(Horn算法),该算法已经集成到A rcGIS中,利用该算法提取的坡度图如附图4。在实际应用中,坡度有两种表达方法,一种是坡度即水平面与地面表面之间的夹角;另一种是坡度百分比即高程增加量与水平增加量之比的百分数,本研究采用第一种方法来表达坡度。剖面曲率是对地面坡度沿最大坡降方向地面高程变化率的度量,通过在A rcGIS下对坡度求“slope”来实现提取(附图5)。

坡长通常是指在地面上沿水流方向到其流向起点间的最大地面距离在水平面上的投影长度。坡长因子是水土保持上的重要因子之一,坡长通过影响坡面径流的流速和流量,影响水流挟沙力,进而影响土壤侵蚀强度,因而是定量计算土壤流失的重要指标,也是土壤侵蚀模型中的必要参数。在不考虑其他影响因素时,坡长越大,径流速度就愈大,汇集的径流也愈大,水力侵蚀也愈强[18]。研究坡长对产流产沙过程对建立土壤侵蚀预报模型、制订有效的水土保持体系有重要的理论价值。本研究对坡长的提取是基于Hickey和Van Remortel等提出的基于DEM提取区域坡度坡长因子的算法,H ickey和Van Rem ortel等编制了相应AM L语言程序[19-21],后又改进为利用C++语言编写的程序,使运行效率大为提高[22]。本研究正是通过运行改进后的程序来实现坡长的提取(附图6),其流程图如图2。

图2 砚瓦川流域坡长提取流程图

2 结果与分析

2.1 流域坡度特征分析

将坡度图属性表导出,利用Visual Foxp ro进行频率统计得到图6,可知流域的坡度分布主要集中在1°～5°,其在累积频率曲线上表现为曲线增势较陡,当坡度＞5°时,累积频率曲线的增势趋于缓和,且5°以上的各坡度频率分布均小于3%。流域16%的面积坡度在 1°以下,坡度在9°以下的约占流域面积的50%,＞30°的地区不到流域总面积的20%。流域坡度分布在 0～78.5°,平均坡度为14.9°,坡度标准差为 14.6,变异系数为 0.98,说明研究区的地面较为平坦,沟壑发育还不完全。

图3 砚瓦川流域坡度频率曲线图

根据文献资料记录,在黄土高原丘陵沟壑区县南沟流域[23],基于5m分辨率DEM提取的坡度最小值为0°,最大值为 72.9°,平均坡度为28.9°;榆林地区典型样区[24]10 m分辨率DEM提取的坡度最小值为0°,最大值为 59.8°,平均坡度为20.3°,坡度标准差为10.48,变异系数为0.52;延川样区的平均坡度为 32.27°,坡度标准差为 11.94,变异系数为0.37。比较结果知,相对于丘陵沟壑区,本研究区的坡度分布范围较广,平均坡度小,说明区域内平缓坡面积比较大,而坡度标准差和变异系数均大于丘陵沟壑区,反映了研究区域的地形起伏波动程度要大于丘陵沟壑区,具体统计数据见表1。

表1 研究区及丘陵沟壑区若干样区地形特征值统计表[23-24]

2.2 流域剖面曲率的特征分析

导出流域的剖面曲率图属性表,并在V isual Foxp ro下进行频率统计(图4)。本研究区的剖面曲率分布主要集中在1～3,其中曲率为1的地区面积最广,约占流域的29.2%,其在累积频率曲线上的表现为曲线增势很陡,当曲率＞1时,累积频率曲线的增势趋于缓和,4以上各曲率的频率分布均小于3%;由累积频率曲线知,流域50%的区域坡面曲率小于7,＞24的地区不到流域总面积的20%。流域剖面曲率分布在0～75.4,平均剖面曲率为12,曲率标准差为13.6,变异系数为1.13。

图4 砚瓦川流域剖面曲率频率曲线图

根据文献资料记录[24],在黄土高原丘陵沟壑区,榆林地区典型样区10 m分辨率DEM提取的剖面曲率最小值为0,最大值为67.8,平均剖面曲率为21.6,曲率标准差为13.2,变异系数为0.61。比较结果知,相对于丘陵沟壑区,本研究区的剖面曲率分布稍广,而平均剖面曲率不到丘陵沟壑区一半,说明陡坡区域相比丘陵沟壑区较少,曲率标准差非常接近,但变异系数是丘陵沟壑区的一倍,表明高塬沟壑区的坡度变化率的离散程度比较大,地貌起伏的波动程度要大于丘陵沟壑区。

2.3 流域坡长特征分析

通过导出坡长图的属性表,在Visual Foxp ro下进行频率统计(图5),本研究区的坡长分布主要集中在10～150 m,其在累积频率曲线上表现为曲线增势较陡,其中坡长为10 m的坡分布最广,占流域的11.21%,150m以上各坡长的频率分布均小于2%;由累积频率曲线知,流域50%的区域坡长小于80m,当坡长大于200m时,累积频率曲线的趋势变缓,坡长大于240 m的坡只占流域总面积的20%。坡长较大值均分布在海拔1 300 m以上的高塬区,主要在流域的边界处,沟壑发育不完全;坡长较小值主要分布在沟谷区,这些地区地表破碎。流域坡长分布在0～9 583.1m,平均坡长为181.8m,标准差为350.8,变异系数为1.93,说明研究区的地面较为平坦,沟壑发育还不完全。

图5 砚瓦川流域坡长频率曲线图

根据黄土高原丘陵沟壑区文献资料记录,在县南沟流域[23],基于5 m分辨率DEM提取的坡长最小值为2.5 m,最大值为485.2m,平均坡长为33.7m;基于25 m分辨率DEM提取的坡长最小值为12.8 m,最大值为443.8 m,平均坡长为54.2 m;在延河流域10m分辨率的典型样区[24],坡长最小值为5m,最大值为1 642.25m,平均坡长为64.16m,坡长标准差为65.30,变异系数为1.02;榆林样区基于10m分辨率DEM提取的坡长,最小值为5 m,最大值为571.3 m,平均坡长为37.05 m,坡长标准差为38.2,变异系数为1.03。比较结果知,相对于丘陵沟壑区,本研究区的坡长分布范围非常广,平均坡长是丘陵沟壑区的3～6倍,坡长标准差也是丘陵沟壑区的5～9倍,变异系数约是丘陵沟壑区的2倍,反映了本研究区地表在水平方向上的平均转折变化水平较小,坡长平均值较大,说明研究区的连续平坦区域面积较大,这也与本区域塬面占主导的实际地形特征相一致。

3 结论

(1)流域的坡度分布主要集中在1°～5°,且5°以上的各坡度的频率分布均小于3%。

(2)剖面曲率分布主要集中在1～3,其中曲率1约占流域的29.2%,4以上各曲率的频率分布均小于3%。

(3)坡长分布主要集中在10～150 m,10 m的坡长坡占流域的11.21%,150 m以上各坡长的频率分布均小于2%。

(4)相对丘陵沟壑区,高塬沟壑区的坡度分布范围较广,平均坡度小,区域内平缓坡面积比较大,但坡度变化率的离散程度比较大,地貌起伏的波动程度要大于丘陵沟壑区。

这些结果的获取有助于区域水保措施部门因地制宜,制定出合理的水保措施,为流域规划方案的制订提供理论依据和技术指导,本研究仅考虑了坡度、剖面曲率与坡长因子,并未考虑其他地形因子如坡向、高程极差、高程曲线积分等,需要今后进一步进行研究。

[1] 朱显谟.中国黄土高原土地资源[M].西安:陕西科学技术出版社,1986.

[2] 李志林,朱庆.数字高程模型[M].武汉:武汉大学出版社,2001.

[3] 汤国安,宋佳.基于DEM坡度图制图中坡度分级方法的比较研究[J].水土保持学报,2006,20(2):157-160.

[4] 赵牡丹,汤国安,陈正江,等.黄土丘陵沟壑区不同坡度分级系统及地面坡谱对比[J].水土保持通报,2002,22(4):33-36.

[5] 刘学军,任志峰,王彦芳,等.基于DEM的任意方向坡度计算方法[J].地域研究与开发,2009,28(4):139-141.

[6] 刘学军,晋蓓,王彦芳.DEM流径算法的相似性分析[J].地理研究,2008,27(6):1347-1357.

[7] 刘学军,卞璐,卢华兴,等.顾及DEM误差自相关的坡度计算模型精度分析[J].测绘学报,2008,37(2):200-206.

[8] 刘学军,王永君,龚健雅,等.DEM流域网络提取算法的误差特性分析[J].测绘学报,2007,36(2):224-230.

[9] 刘学军,龚健雅,周启鸣,等.基于 DEM 坡度坡向算法精度的分析研究[J].测绘学报,2004,33(3):258-263.

[10] 陈楠,王钦敏.基于地形因子信息量的数字高程模型分辨率的选择:以黄土高原的研究为例[J].武汉大学学报:信息科学版,2009,34(6):692-695.

[11] 陈楠,林宗坚,汤国安,等.数字高程模型的空间信息不确定性分析[J].测绘通报,2005(11):14-17.

[12] 杨勤科,M cVicar T R,VanN iel T G,等.ANUDEM和TIN两种建立DEM方法的对比研究[J].水土保持通报,2006,26(6):84-88.

[13] Hutchinson M F.A new p rocedure for gridding elevation and stream line data with automatic removal of spurious pits[J].Journal of Hyd rology,1989,106(3/4):211-232.

[14] 周买春,黎子浩.数值地形图的生成及其水文地貌特征评价[J].水利学报,2002(2):71-74.

[15] 杨勤科,M cVicar T R,VanNiel T G,等.用 ANUDEM建立水文地貌关系正确DEM的方法研究[J].测绘科学,2006,31(6):155-157.

[16] 陈楠,王钦敏,汤国安,等.6种坡度提取算法的应用范围分析:以在黄土丘陵沟壑区的研究为例[J].测绘信息与工程,2006,31(4):20-22.

[17] 刘学军,王叶飞,曹志东,等.基于 DEM 的坡度坡向误差空间分布特征研究[J].测绘通报,2004(12):11-13.

[18] 唐克丽,陈永宗.黄土高原地区土壤侵蚀区域特征及综合治理途径[M].北京:中国科学技术出版社,1990.

[19] Van Remortel R,H am ilton M,Hickey R.Estimating the LS fac tor for RUSLE through iterative slope length p rocessing of digital elevation data[J].Cartography,2001,30(1):27-35.

[20] H ickey R.Slope angle and slope length solutions for GIS[J].Cartography,2000,29(1):1-8.

[21] H ickey R,Sm ith A,Jankow ski P.Slope Length Calcu lations from a DEM within ARC/INFO GRID[J].Computers,Environment and U rban Systems,1994,18(5):365-380.

[22] Van Remortel H M,H ickey R.Calcu lations LS Factor Using DEM dataaccording to RUSLE-based C riteria.www.on linegeographer.com/slope/rusleV4pc.am l2004.

[23] 谢红霞.延河流域土壤侵蚀时空变化及水土保持环境效应评价研究[D].西安:陕西师范大学,2008.

[24] 张婷.基于DEM的陕北黄土高原多地形因子空间关联特征研究[D].西安:西北大学,2005.