GIS支持下岷江上游水文特征空间分析

曾超,赵景峰,李旭娇

(1.四川师范大学地理与资源科学学院,成都610101;2.成都理工大学遥感与GIS研究所,成都610059)

自20世纪30年代以来,岷江上游流域的水资源量呈明显减少趋势,且流域内部水土流失严重,泥石流等地质灾害频发[1]。研究该流域的水文特征可为其水文过程建模、水土流失分析提供重要参数和依据,从而更好地保障水资源的科学管理和灾害的有效防治。GIS技术拥有较强的地理空间数据和属性数据的处理和整合能力,能够快速、有效地实现水文参数的提取和水系模型的构建。因此,GIS技术成为了近年来学者从事该流域水文和灾害研究的重要手段。

舒栋才、田兵伟等利用DEM数据生成了岷江上游流域的数字水系[2-3],但对选取集水面积阈值尚无明确的方法;庞学勇、薛建辉、许申来等从生态景观的角度探讨了岷江上游植被的水文特性[4-6];邵骏、林勇、彭立等研究者以紫坪铺水文站的径流数据为基础,利用分形理论、小波、灰色系统等方法分析了岷江上游径流的丰枯及周期变化,但未考虑影响径流的气象因子[7-9];另外,针对岷江上游流域的河网密度、水系分维、纵剖面以及径流系数的研究较为少见。

本文在前人研究的基础上,首先选取河网密度随集水面积阈值变化趋于稳定的点作为阈值,提取出岷江上游水系,并计算河网密度、水系分维数、河流纵剖面等水文特征;其次结合降水量和径流量资料,计算径流系数。文中水文参数和水系网络的研究是水文模型构建的重要基础,径流系数的分析能为水土保持研究和水资源的管理提供有力支持。

1 研究区概况

岷江在都江堰以上称为岷江上游,其流域位于龙门山、岷山等一系列山脉与四川盆地之间组成青藏高原东缘地形陡变带内[10]。岷江上游流域的区域位置在北纬 30°45′-33°09′,东经 102°35′-103°56′;地处我国东部季风湿润区和青藏高原高寒区的过渡地带,流域年平均降水500～850 mm;多年平均径流量为144亿m3,丰水期6-11月占全年径流的75%,而枯水期12月一次年5月仅占全年的25%;流域内部海拔高程为740～6 190 m,平均海拔 3 440 m[2]。流域范围含盖了四川省的汶川、理县、黑水、茂县、松潘五县。它是成都平原的水源地,也是长江的重要支流。

2 资料与方法

2.1 数据资料

SRTM3-DEM数据来源于USGS(United States Geological Survey),数据的水平分辨率为90 m,垂直分辨率精度为6 m[11];水文资料采用了1982-1987年间岷江上游紫坪铺、姜射坝以及镇江关等8个水文站的径流量数据;气象资料为同时期该流域松潘、黑水、杂谷脑、渔子溪以及寿溪等21个雨量站的降水量数据[12]。

2.2 处理方法

2.2.1 分水岭提取与子流域划分

(1)研究区界限的确定和填洼处理。以手工数字化边界为基础,建立外围2 km范围的缓冲区,将缓冲区以内的区域作为研究范围。并对研究区DEM中存在的凹陷区域进行填充。

(2)水流方向与水流累积矩阵计算。流向的确定采用了D8算法。水流累计矩阵值表示流域地形每个点的水流累积量,即直接或者间接流经某单元格的所有的单元格总数[12]。DEM中每个单元格的水流累积量乘以它的面积,就得到了该单元格上的上游来水总量。

(3)子流域划分方法。以水流累计矩阵为基础,手工数字化集水面积为零的单元格,即得到整个流域以及主要子流域的分水岭。并以此为依据,结合控制主要支流的水文站的地理位置,将岷江上游划分为镇江关、黑水河、杂谷脑河、鱼子溪、寿溪5个子流域以及两个干流流域(图1)。

图1 基于DEM的岷江上游河网提取

2.2.2 水系网络生成 水系的提取结果取决于阈值的确定方法。对于阈值的确定我们采用了孔凡哲、杨邦等提出的河网密度法:即河网密度随着集水面积阈值的增大而减小,起初减小较快,当出现一个突变后趋于稳定,此时的河网密度对应了最合理的集水面积阈值[13-14]。以水流累计矩阵为基础,把大于该阈值的网格归为河道,小于它的做为汇流区。

水系生成后,本文采用斯特拉勒(Strahler)方法对岷江上游水系分级,即从河源出发的河流为1级河流;同级两条河流交汇形成的河流级比原来增加1级;不同级的两条河流交汇形成的河流级等于两者中较高者[15]。

2.2.3 流域与水系的形态参数 在岷江上游子流域划分与水系网络提取的基础上,计算全流域及各子流域的面积、河长和河网密度。基于水系分级结果,又对岷江上游河流的分叉比、河长比以及分维数进行计算。

(1)河流的分叉比(Rb)指流域内除最高级别水系外,每一级别水系的总数与比它高一级别总数的比值;河长比(RL)是指流域内除最低级别水系外,每一级别水系的平均长度与比它低一级别平均长度的比值。整个流域河网的分叉比和河长比为所有等级水系的平均值,即。

(2)采用河系定律法,即利用岷江上游河网的分叉比和河长比来计算其分维数(Db)。表达式为[15]:

2.2.4 河流纵剖面的提取与河道比降计算 从河流源头开始,以DEM最高分辨率(90 m)的间隔分别提取主河道的长度和高程值。然后对提取的数据进行滑动平均并绘制成纵剖面图。根据各流域内主河道高程差与总长度的比值得出河道平均比降。

2.2.5 径流系数的计算 径流系数α是指某一段时间的径流深度R与相应的降水深度P之间的比值[16]。径流深度与降水深度计算方法:

式中:W——控制各流域水文站的径流量;A′——流域面积;以岷江上游流域21个雨量站的年降水数据为基础,采用普通克里金(Ordinary Kriging)插值法对流域内部无观测资料的区域进行插值,得到整个流域的降水量格网数据。各流域内所有单元格对应降水量值的总和与该流域内单元格总数的比为其降水深度,即流域面上的平均降水量。计算公式为

式中:Pij——流域内部第 i行、j列的单元格内的降水量;N——流域内单元格总数。

3 结果与分析

3.1 流域基本特征

以岷山山系郎架岭作为岷江源头,紫坪铺大坝排水口为流域的出水口,计算得出流域面积为22 566 km2,干流长度为344.2 km。流域内最高处海拔为6 190 m,最低处为749 m,岷江上游干流河道(图1)比降为9.7‰。

3.2 河网的提取结果

根据集水面积阈值与河网密度的关系曲线(图2)可判断出:从0.81 km2(100个汇流网格)对应河网密度0.63开始,起初河网密度变化较剧烈,随阈值继续增大其变化逐渐变缓。当横坐标为16.2 km2(2 000个汇流网格),即河网密度为0.172时,曲线的变化(一次导数)趋于稳定,河网密度趋向稳定。于是将16.2 km2作为集水面积阈值,提取出岷江上游水系(图2)。

图2 岷江上游水系集水面积阈值与河网密度关系

3.3 河网密度与河流分级

河网密度等于干支流的总长度和流域面积之比,即表示单位面积上的河流长度。本文分别计算了镇江关、黑水河、杂谷脑河、渔子溪、寿溪5个子流域以及岷江上游全流域的主河道长度、流域面积和河网密度(表1)。

由表1可以得出:(1)岷江上游全流域的河网密度值为0.172 km/km2,其5个子流域的值为(0.17±0.01)km/km2;(2)岷江上游最北面的镇江关流域的河网密度值最大,为0.177 km/km2;最南面的寿溪流域的值最小,为 0.165 km/km2。

表1 岷江上游流域主要水文特征的计算值

据斯特拉勒(Strahler)方法,岷江上游水系可分为5个等级:1级河流有331条、平均长度为5.6 km;2级河流有76条、平均长度为12.7 km(表2)。计算得出各级水系的分叉比为 4～4.3,河长比为 2.3～2.8,其变化较小,近似为一常数。于是求得它们的平均值(全流域河网分叉比与河长比)分别为4.30和2.47。

水系的分维值可作为表征流域形态特征的定量指标。并且由于分维反映了流域形态的复杂性,以及具有尺度不变性,因此分维可以解决不同比例尺度下河流形态的推衍问题,并对理解掌握河流的发育演变规律有一定帮助[17]。何隆华和赵宏对全国14个大流域和67个小流域的分维计算结果表明:其值小于1.6时,河流发育处于侵蚀发育的幼年期;为1.6～1.89时处于壮年期[18]。本文根据式(1)计算得出岷江上游水系分维数为1.61,说明其流域地貌发育已经从幼年期过渡到壮年期阶段。

表2 岷江上游水系等级

3.4 河流纵剖面形态与河道比降变化

经过滑动平均绘出大姓沟、黑水河、杂谷脑河、鱼子溪、寿溪5条主要支流以及岷江上游干流6条河道纵剖面曲线(图3)。

图3 岷江上游干流与主要支流的河道纵剖面

根据图3和河道比降计算得出:(1)所有河流越靠近源头,河道的纵剖面凹曲度越大;(2)沿着岷江上游自北向南分布的大姓沟、黑水河、杂谷脑河、鱼子溪、寿溪这5条支流(图2)的凹曲度呈逐渐增大趋势。求得它们的平均河道比降分别为11.2‰、11.6‰、18.6‰、37.6‰、58.6‰,岷江上游干流为9.7‰,其变化与纵剖面凹曲度一致。岷江上游河谷较场地区的构造抬升是导致AB段(图3)异常的原因[19]。抬升隆起使该地区河床升高,其上游河段堆积作用明显,河道纵剖面会相对平缓;而较场至茂汶段则处于构造隆起区向非构造区过渡的陡坡地带,河流的深切侵蚀作用明显,纵剖面异常陡峭;茂汶以下河段受构造抬升作用减弱,且上游较场至茂汶段侵蚀物质在此段堆积,使其坡度恢复平缓(图3)。

图4 岷江上游多年平均降水量分布(1982-1987年)

图5 岷江上游各子流域平均径流系数(1982-1987年)

3.5 径流系数特征

3.5.1 降水量插值结果 岷江上游流域降水量插值结果表明:降水量的整个趋势是从流域的北、西、南三个方向向东面递减(图4)。位于西面的三打古和南面的紫坪铺是降水量的两个峰值区,其值均大于1 200 mm;北面的松潘和马拉墩站的降水均在650 mm以上;而流域东面的沙坝降水最低,为394 mm。全流域1982-1987年平均降水量为742 mm;降水最大值出现在1984年,为790 mm;最小出现在1985年,为705 mm。子流域中黑水河流域(图4)多年平均降水量最大,为831 mm;而干流流域(上)最小,仅为524 mm 。3.5.2 径流系数计算结果 由图5中各子流域控制点水文站的位置,分别计算了镇江关流域(镇江关水文站控制)、干流流域(上)(镇江关、沙坝、姜射坝水文站之间)、黑水河流域(沙坝水文站控制)、姜射坝以下至紫坪铺水文站之间的整个流域(包括杂谷脑河、鱼子溪、寿溪、干流流域(下)4个子流域)以及岷江上游全流域1982-1987年的年平均径流系数(图5)。得出全流域平均径流系数为0.77。从北向南分布的子流域径流系数值依次递增,其中北面的镇江关流域径流系数值最小,为0.55;南面的姜射坝水文站以下流域最大,为0.88。黑水河为0.73,与全流域的值最为接近(图5)。

4 讨论

(1)集水面积阈值的选取直接关系到从DEM中提取河网的准确度。我们利用河网密度和集水面积阈值的关系,选取河网密度变化趋于稳定时对应的16.2 km2作为阈值生成岷江上游水系。提取结果跟实际水系的匹配度高。

(2)由于雨量站的年降水量是一种空间上的连续过程,其分布在一定距离内具有一定的结构性和自相关性[20]。因此,本文选取普通克里金插值法对岷江上游流域降水量进行插值。但由于该流域内部雨量站点有限,且插值过程中未考虑地形、坡向等因素对降雨的影响,插值结果较实际降水量有一定误差。

(3)在划分的7个子流域中(图1),渔子溪和寿溪流域无流量观测资料,而对于杂谷脑河流域,缺少其控制站桑坪水文站1986年径流数据,所以将杂谷脑河、渔子溪、寿溪以及干流流域(下)合并为姜射坝水文站以下流域(图5)。

(4)通过计算岷江上游水系分维数,并结合何隆华和赵宏对全国14个大流域和67个小流域的研究成果得出:该流域地貌发育已经从幼年期过渡到壮年期阶段。水系分维数不仅能表达流域的发育状况,还可以反映水系的密布性、弯曲性等自然特性以及气候变化对流域地貌的影响[17,21]。因此,它与流域地质、地貌、气候和水文状况的关系有待进一步研究。

5 结论

利用DEM、降水和径流数据,基于GIS空间分析技术,系统定量地提取和分析了岷江上游流域的河网、子流域、分叉比、水系分维数以及河网密度、河流纵剖面、径流系数等典型水文特征。

(1)岷江上游流域1982-1987年间平均径流系数为0.77,平均降水量为742 mm,径流深度为573 mm,即该流域169 mm的降水量消耗在了蒸发、植被截留、土壤蓄水和地下水的补给等过程中。从北到南的子流域径流系数值依次递增,其中最南面的姜射坝以下流域对岷江上游的产水量贡献最大,其径流系数为0.88。

(2)流域干流的平均比降为9.7‰。沿着岷江上游自北向南分布的大姓沟、黑水河、杂谷脑河、鱼子溪、寿溪这5条支流的纵剖面凹曲度和比降均呈逐渐增大趋势,大姓沟的比降为11.2‰,寿溪为58.6‰。

(3)整个流域的河网密度值为0.172 km/km2;河流的分叉比为4.30,河长比为2.47,水系分维数为1.61。在5个子流域中,镇江关流域的河网密度值最大,寿溪流域最小,岷江干流东西两侧的子流域发育极不对称,其中面积超过500 km2以上的流域均分布在干流的西侧。

[1] 王渺林.岷江流域水资源安全问题探讨[J].四川水利,2005(4):32-34.

[2] 舒栋才,程根伟,林三益.基于DEM 的岷江上游数字流域的离散化[J].四川大学学报:工程科学版,2004,36(6):6-11.

[3] 田兵伟,范建容,王道杰,等.基于GeoProcessing方法的岷江上游流域数字水系建模[J].长江流域资源与环境,2008,17(1):106-112.

[4] 庞学勇,包维楷,张咏梅,等.岷江上游中山区低效林改造对枯落物水文作用的影响[J],水土保持学报,2005,19(4):119-122.

[5] 薛建辉,郝奇林,何长青,等.岷江上游两种亚高山林分枯落物层水文特征研究[J].水土保持学报,2009,23(3):168-172.

[6] 许申来,李秀珍,胡远满,等.岷江上游黑水河流域不同景观结构小流域径流系数的比较[J].生态学杂志,2007,26(5):712-717.

[7] 邵骏,袁鹏,李秀峰.岷江上游年径流变化的多时间尺度分析[J].东北水利水电,2007,25(6):52-54.

[8] 林勇,刘世荣,李崇巍,等.小波变换在岷江上游杂古脑流域径流时间序列分析中的应用[J].应用生态学报,2005,16(9):1645-1649.

[9] 彭立,苏春江,徐云,等.岷江上游径流枯水年时间序列的分形研究[J].水土保持研究,2007,14(4):440-442.

[10] 张会平,杨农,张岳桥,等.岷江水系流域地貌特征及构造指示意义[J].第四纪研究,2006,26(1):126-135.

[11] United States Geological Survey.Shuttle Radar Topography Mission[EB/OL].[2009-10-21].http://edc.usgs.gov/#/Find_Data/Products_and_Data_Available/SRTM

[12] 程根伟,舒栋才.水文预报的理论与数学模型[M].北京:中国水利水电出版社,2006.

[13] 孔凡哲,李莉莉.利用DEM提取河网时集水面积阈值的确定[J].水电能源科学,2005,23(4):65-67.

[14] 杨邦,任立良.集水面积阈值确定方法的比较研究[J].水电能源科学,2009,27(5):11-14.

[15] 芮孝芳.水文学原理[M].北京:中国水利水电出版社,2004:18-43.

[16] 黄锡荃,李惠明,金伯欣.水文学[M].北京:高等教育出版社,2005.

[17] 冯平,冯焱.河流形态特征的分维计算方法[J].地理学报,1997,52(4):324-330.

[18] 何隆华,赵宏.水系的分形维数及其含义[J].地理科学,1996,16(2):124-128.

[19] 赵小麟,邓起东,陈社发.岷山隆起的构造地貌学研究[J].地震地质,1994,10(4):429-439.

[20] Huang W C,Yang F T.Streamflow estimation using Kriging[J].Water Resources Research,1998,34:1599-1608.

[21] 王博,田富强,胡和平.基于分维的水系发育程度与气候特征关系[J].清华大学学报:自然科学版,2009,49(12):1948-1953.