基于地貌视角的喀斯特流域水系分维估算方法适应性分析

安 全，贺中华，赵翠薇，梁 虹，焦树林，杨朝晖

(1.贵州师范大学地理与环境科学学院，贵阳 550001；2.贵州师范大学国家喀斯特石漠化防治工程技术研究中心，贵阳 550001；3.贵州省山地资源与环境遥感应用重点实验室，贵阳 550001；4.贵州省应急管理厅，贵阳 550001)

0 引言

自然界许多系统、要素的空间组合结构具有某种意义上的分形特点。诸多河流发育的支流与干流在空间结构上呈现出一种典型的分形特征。通过传统计算水系分维方法工作量大、费时耗力，利用地理信息系统(geographic information system，GIS)技术并结合1∶5万地形图水系与30 m空间分辨率ASTER-GDEM数字高程数据估算流域水系分维，探讨流域地貌发育与水文水系特征的关系是当前研究的热点问题。分形是自然界诸多客体表现局部特征与整体的相似性客观几何规律，在现代数学的发展中逐步演化成一支新的具有独特理论体系与认识的世界观与方法论。在地学领域，地貌发育与河流水系的演化规律中本质上都潜藏着分形的自相似性特点。分形理论在地学界的快速发展始于20世纪90年代，早在1977年Mandelbrot首先将分形理论引入水文学，成为水系分维理论的创始人。进入20世纪90年代后，国内外对水系分维的研究迅速展开。何隆华等[1]利用计盒法对全国14个大流域与67个小流域进行了计算，并得出水系分维值1.6是流域地貌发育阶段的临界值。进入21世纪后，基于数字高程模型(digital elevation model,DEM)提取水系估算水系分维的研究得到快速发展[2-5]。往后随着科学技术的快速发展，特别是GIS与DEM的发展与建立使得基于ArcGIS系列技术与DEM模型的流域水系分维研究进入大发展时期[6-15]。如王倩[16]等利用GIS技术，对秦淮河流域水系分维展开了研究，王林等[17]基于ArcGIS8.3和ArcVIEW3.2进行了DEM模型的流域水系分维估算分析。

目前，水系分维的研究主要是区域流域水系分维与地貌发育关系，人工渠系与自然水系分维对比，水系分维与洪涝、径流过程，水系分维与岩性识别，水系分维与居名点分布特征，水系分维的生态意义等方面的研究[18-20]。而针对中国西南喀斯特筑坝流域地区与地貌发育对水文特征影响方面的研究报道还较少。鉴于此，本文以黔中筑坝工程区龙场桥流域为研究对象，以1∶5万地形图水系、30 m空间分辨率的ASTER-GDEM和1∶10万贵州省综合地貌图为数据源，利用基于ArcGIS10.2的Horton-Strahler理论、水系栅格法、渔网法估算黔中筑坝工程区龙场桥流域水系分维，探讨筑坝区地貌发育对流域水文特征的影响，以期为喀斯特地区的生态建设，经济的可持续发展，筑坝带的洪涝风险评估以及流域综合管理等方面提供一定参考。

1 研究区概况及数据源

1.1 研究区概况

研究区是以龙场桥水文控制断面提取的流域，位于黔中水利工程区上游地区(图1)。黔中水利枢纽工程是贵州省首个大型跨地区、跨流域长距离水利调水工程，也是黔中地区生存和发展的生命线工程[23]。工程区内是长江流域和珠江流域两大流域交错地带，研究区流域面积4 074.15 km2，地形西高东低，海拔在1 138～2 800 m之间，流域的西部地貌类型为典型的喀斯特岩溶高原[27]，中东部峡谷、峰丛洼地相间分布，气候属于亚热带季风气候，夏季高温多雨，冬季温和湿润，四季分明，多年平均降雨量在1 100～1 400 mm之间，最热月平均气温22 ℃以上，最冷月0～15 ℃左右。植被为亚热带常绿阔叶林、灌木林、草地和裸岩相间分布。土壤以黄壤为主，兼有非地带性石灰土分布。流域内发育的乌江水系南源三岔河干流全长325.6 km，流域集水面积7 264 km2，是黔中地区大型水利工程集中地。

(a)龙场桥流域位置 (b)龙场桥流域DEM

图1 研究区概况图

Fig.1Researchareaoverview

1.2 数据源及其预处理

本文数据源包括1∶5万地形图水系(以下简称地形图水系)和30 m空间分辨率的ASTER-GDEM等。在对数据源的坐标进行一致性处理的基础上，使用30 m空间分辨率 ASTER-GDEM数据进行水系提取，利用ArcGIS10.2软件的Hydrology水文分析工具集对从中国科学院计算机网络信息中心的地理空间数据云(http://www.gscloud.cn)获取的ASTER-GDEM数字高程数据进行填洼、流向、汇流累积量计算处理，利用ArcGIS10.2软件的Con命令设定汇流阈值1 500提取研究区河网栅格水系(以下简称ASTER-GDEM水系)。

2 研究方法

利用ASTER-GDEM水系、结合0.6 m空间分辨率Google Earth影像、30 m空间分辨率Landsat 8 OLI遥感影像，对提取的筑坝水库区平行状水系依据Strahler理论进行了分级。依次将提取水系栅格的边长设定为500～10 000 m，步长为500 m；然后利用ArcGIS软件的SQL查询功能分别统计不同边长水系栅格所对应的栅格数目，再对统计的栅格边长、数目与Strahler分级的矢量水系导入相关统计软件，分别求栅格法与Horton定理的水系分维[28-29]。

2.1 Horton定理

Horton-Strahler计算原理为利用ArcGIS 10.2 提取1∶5万地形图水系，应用Strahler分级原理对提取的水系进行河流分级；同理，以ASTER-GDEM为数据源，利用ArcGIS 10.2的con函数对30 m空间分辨率DEM分别设定汇流阈值100～10 000自动提取水系，并与1∶5万地形图尺度下提取的水系不断做叠加试验，通过汇流阈值与河网密度的拟合曲线变点分析，确定用以提取研究区河网水系的最佳汇流阈值为1 500，再根据Strahler河网分级理论对矢量河网进行分级；最后统计地形图水系和ASTER-GDEM通过最佳阈值1 500提取水系的每一级河流分级的河道数及其河道长度，计算2个数据源提取水系的的分叉比与河长比，对分叉比与河长比取对数，比值即为所求分维。计算公式如下：

Rb=Ni-1/Ni，i=2,3,…,n，

(1)

(2)

(3)

式中：Rb为分叉比；RL为河长比；Ni为各级河道数；i为河道等级；DH为Horton-Strahler水系分维值。

2.2 栅格法

栅格法又称覆盖法、网格法、计盒维数法，其原理是：取边长为r的正方形网格与水系图求取交集，得到水系覆盖的网格数目为N(r)，当r不断变化时，与之相对应会得到一系列的N(r)值，二者关系为

N(r)∝r-D。

(4)

左右两边求取对数，以点(lgr,lgN(r))为坐标作双对数图，采用最小二乘法可拟合出一条直线，即

lgN(r)=-Dlgr+b，

(5)

式中：r为正方形网格的边长；N(r)为对应边长网格与水系图求取交集所覆盖的网格数目；b为待定系数；D为双对数曲线的斜率值，即所求的水系分维[25]。

3 结果与分析

3.1 黔中筑坝流域水系分维提取分析

利用ArcGIS10.2软件对提取水系进行Strahler二次分级，通过Horton-Strahler法求得地形图水系与研究区ASTER-GDEM水系的水系分维，再通过水系栅格法、渔网法对提取地形图水系与ASTER-GDEM水系，分别依次将栅格边长设定为500～10 000 m，步长为500 m，分别统计不同栅格边长的栅格个数，并求取栅格边长与不同边长栅格数的对数，利用Origing9.1统计软件作栅格边长与不同边长栅格数双对数拟合线，求其斜率即为分维值，统计结果如表1—4所示。

表1 基于地形图水系的Horton-Strahler法参数Tab.1 Horton-Strahler method parameters based on water system of topographic map

表2 基于ASTER-GDEM水系的Horton-Strahler法参数Tab.2 Horton-Strahler method parameters based on water system of ASTER-GDEM

表1和表2的统计结果表明，河流数与河流平均长度随着河流级数的增大而减少。河流的分叉比随着河流级数的增大而增大，与河流数、河流平均长度成反比；河长比随河流级数越大而减少。总体而言，河流分叉比、河长比除了受到河流数及其河流长度的影响，与河流的分级存在一定的相关性。

表3和表4是地形图水系与ASTER-GDEM水系栅格参数和渔网法参数。统计结果表明，栅格水系法与渔网法统计的栅格数目随着试验设定的栅格边长的增加而减少，栅格边长与栅格数目之间存在明显的线性关系。

表3 地形图水系与ASTER-GDEM水系栅格参数Tab.3 Water system of topographic map and ASTER-GDEM grid parameters

表4 地形图水系与ASTER-GDEM渔网法参数Tab.4 Water system of topographic map and ASTER-GDEM fishing net method parameters

表5是结合表1至表4统计数据，利用分维计算公式根据表1和表2分别求出1∶5万地形图提取水系和ASTER-GDEM提取水系的Horton-Strahler分维；利用Origing9.1统计软件根据表3和表4分别作1∶5万水系和ASTER-GDEM的栅格水系法及渔网水系法试验得到栅格边长与栅格数的双对数拟合曲线。结果见表5。

表5 不同方法提取地形图水系与ASTER-GDEM分维值对比表Tab.5 Comparison of fractal values from water system of topographic map and ASTER- GDEM by different methods

从表5中可以看出，同一地区不同方法、不同数据源估算的水系分维相差较大。其中Horton-Strahler法估算地形图水系差别显著，根据Horton-Strahler法估算地形图水系与ASTER-GDEM水系的分维值分别为1.69和0.66，二者相差1.03。根据河数定律与河长定律的分叉比、河长比理论取值范围可知，针对二元三维地表地下结构，溶洞、裂隙发育的典型喀斯特地貌而言，由于诸多干流沿途遇到溶洞，大型溶蚀管道、裂隙时会转入地下流一段距离等特殊现象，使得在地形图上提取水系时无法判断地下的水流状况导致提取水系误差增大，根据地形图水系计算的河流的分叉比、河长比的比值范围未在理论值范围，以至于所计算出的水系分维值与水系栅格法、渔网法估算的分维值大相径庭，不能稳定的判断研究区地貌发育状况。与此同时，表1和表2根据Horton-Strahler法提取的地形图水系与ASTER-GDEM水系的河段数差别尤为显著，总体上相差1 174段，河流分级的河段数差别主要体现在1～3级。这种显著差异可能是由于复杂的喀斯特地貌组合结构所致。

然而，栅格水系法与渔网法估算的水系分维无论是地形图水系还是ASTER-GDEM水系，估算的分维值均比较接近。就表5数据表明，通过ASTER-GDEM水系利用2种方法估算的分维略比地形图水系大。另外，从水系提取的R2值来看，1∶5万地形图尺度下栅格水系法与渔网法提取水系分维的R2值分别为0.996 6和0.996 4，而基于ASTER-GDEM采用栅格水系法与渔网法提取水系分维的R2值分别为0.994 1和0.993 4，这表明前者提取水系的无标度区间在0.996 4～0.996 6之间，而后者为0.993 4～0.994 1之间。

3.2 黔中筑坝流域水系分维与地貌发育关系分析

图2为研究区根据1∶10万贵州省综合地貌图矢量化的研究区地貌发育现状图和30 m空间分辨率ASTER-GDEM数字高程数据利用ArcGIS10.2制作的研究区山体阴影图。表6为研究区各类型地貌类型及其面积比例。

(a)研究区地貌类型 (b)研究区山体阴影

图2 研究区各流域地貌类型与山体阴影图

Fig.2Landformtypesandmountainshadowmapsofvariouswatershedsinthestudyarea

表6 不同流域地貌发育参数及其面积比例Tab.6 Geomorphological development parameters and area ratio of different watersheds

①字母F，K1，K2分别表示地貌发育的侵蚀-剥蚀类型、溶蚀为主型、溶蚀-侵蚀类型。

从图2(a)中可以看出，整个研究区地貌发育类型以中山谷地、峰丛谷地、峰丛洼地为主，其中龙场桥流域、阳长流域少部分地区出现深切中山型地貌。从图2(b)山体阴影图可以看出，研究区地貌发育破碎程度大，河谷多呈“V”型，越往研究区的东部即阳长、龙场桥一带，地势起伏越大。地貌发育组合以深切中山、中山谷地、峰丛洼地复型组合为主。从表6中可以看出，向阳流域地貌发育组合以峰丛谷地、峰林溶原(盆地)、浅切中山为主，所占流域面积比例分别为23.31 %，20.83%和38.04%，研究区比例分别为4.76%，4.25%和7.77%；阳长流域地貌发育组合以K化中山谷地、峰丛谷地、峰丛洼地为主，所占流域面积比例分别为26.75%，47.43%和11.62%，研究区比例分别为10.48%，18.59%和4.56%；龙场桥流域地貌发育以K化中山谷地、峰丛谷地、峰丛洼地为主，所占流域面积比例分别为43.23%，11.95%和41.14%，研究区比例分别为17.46%，4.83%和16.61%。

3.3 不同流域地貌发育成因分析

从表6中可以看出，向阳流域、阳长流域和龙场桥流域3个流域的地貌发育成因主要为侵蚀-剥蚀类型，溶蚀为主类型、溶蚀-侵蚀类型，向阳流域由于侵蚀-剥蚀成因下发育的地貌类型深切中山面积占整个流域的38.04%，这说明整个流域地貌的发育主要以侵蚀-剥蚀为主，其次为溶蚀、溶蚀-侵蚀型。阳长流域由于溶蚀、溶蚀-侵蚀成因下发育的地貌类型峰丛谷地面积占整个流域的47.43%，其地貌发育主要以溶蚀、溶蚀-侵蚀为主。龙场桥流域由于溶蚀、溶蚀-侵蚀成因下发育的地貌类型K化中山谷地、峰丛洼地面积占整个流域的34.07%，故其地貌发育主要以溶蚀、溶蚀-侵蚀为主。从地貌发育的角度看皆属于地貌发育幼年期的中晚期阶段。

综合水系估算分维值来看，根据根据何隆华的计算的水系分维与地貌发育阶段的关系可知，当D或DH≤1.6时，流域地貌发育属于侵蚀发育阶段的幼年期；当D>1.6或DH≤1.89时，流域地貌发育属于侵蚀发育阶段的壮年期；当D或DH>1.89时，流域地貌发育属于侵蚀发育阶段的老年期。Horton-Strahler法、水系栅格法、渔网法提取除了Horton-Strahler法估算水系分维与研究区实际出入较大外，水系栅格法、渔网法估算的水系分维与研究区地貌发育较为吻合。结合表5和表6数据来看，综合分析典型喀斯特复杂的地貌结构与水系分维关系，渔网法估算的水系分维与研究区实际地貌现状最为吻合。根据渔网估算的水系分维可知，研究区利用渔网法估算的1∶5万地形图提取水系分维值为1.54，通过ASTER-GDEM提取水系估算的分维值约为1.60，这说明研究区正处于地貌发育阶段的幼年晚期、壮年期早期。此结果与研究区实际发育地貌吻合。

4 结论与展望

本文以黔中筑坝工程区龙场桥流域为研究对象，以1∶5万地形图水系、30 m空间分辨率的ASTER-GDEM和1∶10万贵州省综合地貌图为数据源，利用基于ArcGIS10.2的Horton-Strahler理论、水系栅格法、渔网法估算流域水系分维，探讨筑坝区地貌发育对流域水文特征的影响，结果表明：

1)喀斯特地区复杂地貌组合结构下不同方法、不同数据源估算的水系分维相差较大。Horton-Strahler法、水系栅格法、渔网法估算地形图水系分维值分别为1.69，1,53和1.54；估算ASTER-GDEM水系的分维值分别为0.66，1.59和1.60。其中Horton-Strahler法估算分维值差别显著，差值达到1.03。

2)综合分析Horton-Strahler理论、水系栅格法、渔网法估算喀斯特筑坝区不同数据源水系分维与实际地貌发育的关系可知，渔网法估算的水系分维与研究区实际地貌现状最为吻合。根据渔网估算的水系分维可知，研究区利用渔网法估算的地形图水系分维值为1.54，ASTER-GDEM水系估算的分维值约为1.60，这说明研究区正处于地貌发育阶段的幼年晚期、壮年期早期。此结果与研究区实际发育地貌吻合。此外，3种方法估算喀斯特筑坝流域的水系分维精度排序为：渔网法>水系栅格法>Horton-Strahler法。