沂河上游土地利用变化的设计洪水响应研究

(河海大学 水文水资源学院,南京 210098)

土地利用/覆被变化(LUCC)会引起流域下垫面条件发生变化,导致流域产汇流规律发生不同程度的改变,直接影响暴雨洪水特性[1].目前,国内外主要研究LUCC 的径流响应,关于LUCC 的设计洪水响应研究还相对较少.吴淼等[2]对淮河上游小流域构建SWAT 模型进行模拟,结果表明该流域在LUCC 前后,频率越大,流域设计洪水减小越多.雷超佳等[3]采用HEC-HMS模型对奉化江皎口水库流域进行研究,分析了该流域不同重现期的洪水过程在LUCC影响下的变化情况.魏冲等[4]以深圳市石岩流域为研究对象,经分析得出,该流域不同频率的设计洪水受城市化影响较剧烈.

近年来,沂河上游流域土地利用方式发生改变,不透水面积呈现增加趋势,容易导致洪涝灾害频繁发生,因此,迫切需要定量研究设计洪水对LUCC 的响应,了解下垫面变化对流域防洪具有重要的意义.本文基于HEC-HMS分布式水文模型,分析沂河上游流域1985-2015年LUCC 对设计洪水的影响,通过划分响应分区探讨不同重现期下设计洪水对LUCC的空间响应,为流域防洪减灾提供一定的科学依据,也为流域下游水库安全运行提供一定的参考.

1 研究区概况

沂河发源于山东省沂源县田庄水库,上游为发源地至沂源县东里店,地势自西北向东南倾斜,以低山丘陵为主,干流河长63.5km,河道平均比降2.7‰.沂河上游流域夏季高温多雨,汛期集中在6～9月.源头田庄水库是大(二)型山区水库,水库流域面积424 km2,水库总库容1.2亿m3.研究区位于田庄水库和东里店水文站之间,如图1 所示,流域面积为746 km2,土地利用方式主要以林地、耕地为主,土壤类型为棕壤土、褐土和潮土.研究区以东里店水文站为控制站点,除东里店站以外,流域还有悦庄、芦芽店、石桥和焦家上庄降水站.

图1 研究区位置图

2 数据与方法

2.1 基础数据

遥感影像数据采用研究区1985年Landsat5 TM和2015年Landsat8遥感影像,以及Sentinel-2高分辨率遥感影像;数字高程数据(DEM)采用30m 精度的ASTER-GDEM V2高程图,以上数据均来自于地理空间信息云平台;土壤类型数据采用中国科学院寒区旱区环境与工程研究所1∶100 万的土壤栅格数据;气象水文数据采用研究区降水站和水文站获得的小时降雨径流数据.

2.2 研究方法

HEC-HMS是美国陆军工程师团水文工程中心研发的降雨径流模型,是一个具有物理概念的半分布式水文模型[5],主要由5个模块构成,分别是气象模块、产流计算模块、坡面汇流计算模块、河道基流计算模块和河网汇流计算模块[6].HEC-HMS模型能够综合考虑研究区下垫面条件的时空变化,不同的子流域可以选用不同的产汇流方式[7].为了更好地分析研究区LUCC对设计洪水的影响,本文选取SCS曲线法、SCS单位线法、马斯京根法和衰减指数法分别作为HEC-HMS模型产流、坡面汇流、河网汇流和河道基流的计算方法.

1)SCS曲线法是美国水土保持局(SCS)开发用于流域产流的水文方法,SCS曲线法通过土地利用/覆盖、累计雨量和前期土壤湿度来估算流域的净雨量[8],能够较好地反映土地利用变化对流域产流的影响,计算公式如下:

式中:Pe、P分别为t时刻的累计雨量、降雨深度(mm);S为子流域最大截留降雨能力(mm);CN 为SCS曲线法的曲线数,取值范围为30～100,与土壤渗透率有关.

2)SCS单位线法通过流域滞时tlag反映,与流域土地利用条件、土壤类型、水流长度和坡度有关[9],能够较好地模拟土地利用变化下的坡面汇流情况,计算公式如下:

式中:tlag为子流域的洪峰滞时(min);L为水流长度(m);S为子流域最大截留降雨能力(mm);y为子流域平均坡度(%);L和y可利用HEC-GeoHMS模型插件通过DEM 来提取.

3)马斯京根法常用于山区河道流量演算,能较好地模拟研究区河网汇流的情况,其参数河段传播时间K和流量比重因子X可由以下公式估算:

式中:a为洪水波速修正系数,在1.33～1.67之间;Δx为河段长度(m);c′为洪水平均波速(m/s);Q0为参考流量(m3/s);B为水流顶宽(m);S0为底坡度.

4)衰减指数法描述了基流从初始流量开始以指数的形式进行衰减,可通过实测径流退水部分估算衰减指数,t时刻基流Qt的计算公式如下:

式中:Q0为初始基流(m3/s);k′为衰减指数,取值范围在0～1之间.

3 模型构建与率定

3.1 研究区土地利用/覆被变化分析

通过ENVI采用神经网络法对1985年和2015年遥感影像进行监督分类,把研究区土地利用类型划分为林地、耕地、草地(低矮灌木林、草地等)、水体、裸地和建设用地六大类型,如图2所示.结合沂源各类历史专题研究图件、地形图等相关资料并参考地图软件的历史影像,对监督分类结果进行人工校正.通过Sentinel-2遥感数据,采用混淆矩阵进行精度评价,总体精度分别达到82%、88%,解译精度满足研究要求.

图2 研究区1985年和2015年土地利用分布图

通过ENVI分析两期的土地利用分布图,得到研究区1985-2015年的土地利用转移矩阵(见表1).

表1 研究区1985-2015年土地利用类型面积转移矩阵 (单位:km2)

从表1可知,研究区的土地利用主要以耕地、林地、草地、建设用地为主,四者所占的总面积占研究区面积的97%以上,水体与裸地所占的面积相对偏少.相较于1985年,2015年的林地、耕地面积均有所下降,分别减少了24.44km2、32.68km2,所占比例分别从38.95%、46.76%下降到35.72%、42.43%;草地、建设用地面积均有所上升,分别增加了8.87km2、50.33km2,所占比例分别从7.01%、4.47%上升到8.14%、11.14%,建设用地变化最为显著.

在1985-2015年土地利用转移矩阵中,林地、耕地分别约有16km2、35km2转变为建设用地,27 km2、23km2转变为草地,这主要是由于社会经济的发展,研究区30年以来不断进行城镇化建设,城镇用地不断增加,一部分的林地、耕地用于旅游用地,建设旅游景点,发展当地的旅游业,同时,丘陵山地以果园为主的灌丛林广泛种植,取代原来的乔木林用地与农作物用地,提高当地的经济水平.由于原耕地的减少以及土地利用合理规划,林地、草地分别约有67.15 km2、27.58km2转变为耕地.

从研究区LUCC 可知,林地、耕地面积减少,草地、建设用地增多,研究区的土地利用/覆被结构发生了较大的变化,下垫面条件的改变会引起流域产汇流规律的变化,对各子流域的降雨径流过程产生不同的影响.

3.2 模型基础数据处理

利用HEC-GeoHMS 模块对研究区DEM 进行流域构建、水系提取、流域划分等一系列的地表水文分析,以东里店水文站作为流域出口,把研究区划分为25个子流域,13条概化河道,如图3所示,对每个子流域和概化河道进行物理参数、水文参数的提取,并为每个子流域和概化河道选定产汇流计算方法.为了反映时空分布的差异性[10],本文采用距离平方倒数加权法,通过HEC-GeoHMS计算5个降水站在各个子流域的权重系数.

图3 研究区子流域划分图

CN 采用SCS制定的土壤分类标准,根据土壤渗透率高低分为A、B、C、D 四类,通过中国土壤科学数据库获取研究区土壤物理性质,把原土壤类型进行归并,得出研究区水分土壤组为B、C和D 三组.根据研究区的土地利用情况,确定前期土壤湿度处于AMCⅡ平均等级下各地类的CN2S值[11-13],将土壤数据与土地利用数据进行空间叠加分析,最终通过HECGeoHMS计算研究区内各子流域的平均CN2S值.确定各地类的不透水率Im[14-16],对各个子流域中所有土地利用的不透水率进行加权平均,估算每个子流域的不透水率.

3.3 模型率定与验证

由于研究区降水集中于6～9月,植物处于生长期,在前期土壤湿度AMCⅡ平均等级下,前5天总雨量为35.6～56.3mm[17].根据这一特点,选取研究区东里店水文站1979-2012年共10场实测径流过程,用于模型的率定与验证.由于沂河源头为田庄水库,把田庄水库对应场次洪水的下泄流量作为本模型的入流条件.时间步长的长短会影响模拟值与实测值的拟合程度[18],因此,通过HEC-DSSVue的插值功能,对径流资料存在的数据缺失值,按照15min的时间步长进行插值应用.

选取峰值加权标准误差方程作为拟合目标函数[19],采用洪峰流量相对误差ΔQ、径流深相对误差ΔR、峰现时间误差ΔT和Nash效率系数E共4个检验标准,校检HEC-HMS 在研究区的洪水模拟,1985年和2015年主要参数率定优化结果见表2,两者模拟的结果见表3.

表2 模型主要参数率定优化结果

表3 研究区水文模拟结果

由表2可知,在模拟的10场洪水过程中,洪峰流量相对误差和径流深相对误差均在20%范围以内,峰现时间误差均在3h以内,Nash效率系数E均在0.84以上.根据《水文情报预报规范》(GB/T 22482—2008)的规定,满足洪峰和径流深相对误差在20%以内,峰现时差在3h以内,确定性系数在0.7以上即为合格.因此,本文构建的HEC-HMS模型在东里店流域洪水模拟的总合格率为100%,同时,实测与模拟的洪水过程线拟合良好,如图4所示,HEC-HMS半分布式水文模型模拟研究区降雨径流的效果较好.

图4 实测与模拟的洪水过程

4 LUCC下研究区设计洪水响应特性分析

4.1 LUCC对设计洪水的影响分析

为了分析研究区LUCC 对设计洪水的影响,采用研究区内5个站点1968-2016年共49年的降雨资料系列,选取不同年份各站点最大24h暴雨,采用P-Ⅲ型频率曲线进行排频计算.

根据研究区所在的地理位置,选用山东省泰沂南山区1h雨型对设计雨量进行逐日设计暴雨时程分配得到P=10%、P=5%、P=2%和P=1%的逐时设计暴雨.通过HEC-HMS模拟计算两种下垫面条件下研究区不同重现期的设计洪水,流域出口模拟结果见表4,设计洪水过程如图5所示,各子流域洪水变化见表5.

表4 研究区流域出口断面模拟结果

表5 研究区子流域洪水变化率 (单位:%)

续表5 研究区子流域洪水变化率 (单位:%)

图5 不同下垫面条件下的设计洪水过程

由模拟结果可知,2015年与1985年下垫面条件模拟的洪水相比,研究区P=10%、P=5%、P=2%、P=1%的洪峰和径流深均有不同程度的增加.由图5 可知,峰现时间分别提前了105、105、90、90 min.可以发现,重现期越低,研究区洪峰和径流深的的变化就越大,峰现时间越提前.

由表4可知,重现期由低到高,流域出口的径流深分别增加了6.18、5.63、6.19、6.58mm,洪峰分别增加了68.08、74.43、81.08、85.34m3/s.可以发现,研究区LUCC 在增加相同径流深的情况下,重现期越高,洪峰增加得越多.

由表5可知,研究区各子流域不同重现期下设计洪水变化情况,石桥镇北部(W670)的设计洪水变化最大,西里镇北部(W1060)的变化最小.南麻街道北部(W630)、石桥镇(W820、W890)、悦庄镇南部(W700)、沂河沿岸(W900、W990)和燕崖镇西部(W950、W960)的设计洪水变化较大,结合土地利用变化分析得知,前三者由于城镇化建设,城市化向外扩张,建设用地逐渐增加,林地耕地转变为不透水用地,流域产流能力增加,导致它们的设计洪水产生较大变化;后两者设计洪水变化较大的原因,主要是由于林地向耕地、草地转变,林地面积减少,人类活动影响增加,流域涵养水源能力降低.

通过比较研究区洪峰与径流深可知,洪峰变化小于径流深变化,在忽略基流径流深时,洪峰变化则大于直接径流深的变化.通过比较子流域的面积、ΔCN和设计洪水变化率可知,子流域面积越大,ΔCN 增加越多,研究区设计洪水变化就越大.

4.2 LUCC下设计洪水的空间响应分析

为了分析LUCC 下设计洪水的空间响应,提供研究区防洪减灾依据,本文根据洪峰变化率和CN 变化值的关系趋势以及变化范围,选定ΔQ/Q1985=10%和ΔCN=1作为分区判定条件,把研究区分为3个分区.当ΔQ/Q1985＞10%且ΔCN＞1时为1区,当ΔQ/Q1985＜10%且ΔCN＞1时为2 区,当ΔQ/Q1985＜10%且ΔCN＜1时为3区.不同重现期下的空间响应分区如图6所示.由图6可知,随着重现期由低到高,设计洪水对LUCC的响应从中部地区(沂源县主城区)向周边地区逐渐减弱,50年一遇和100年一遇的空间响应基本一致,出现这种趋势主要是由于城市化向外扩张,周边乡镇地区的林地面积大幅度减少导致的.设计洪水对LUCC 响应较强的1区,主要分布于东部地区,集中于石桥镇区域,建设用地明显增加;西部地区随着重现期的变高,由1区转为2区,该区域设计洪水对LUCC 的响应逐渐减弱,敏感性降低,林地向耕地、草地转变的同时,也存在耕地向草地转变;而响应较弱的3区,主要分布于南北部地区,南部地区响应最弱,该地区30年间土地利用变化较小,主要是林地、耕地向草地的转变.

图6 不同重现期下空间响应分区图

5 结论

1)1985-2015年间,研究区土地利用类型以林地、耕地、草地和建设用地为主,占总面积的97%以上.相较于1985年的土地利用,2015年时的林地和耕地面积均减少,草地和建设用地的面积均有所增加,建设用地的增幅最为显著.

2)在1985年和2015年的不同下垫面条件模拟下,研究区内不同重现期的洪峰和径流深均有不同程度的增加,峰现时间提前,低重现期的洪水变化最显著.在增加相同径流深时,重现期越高,洪峰呈现增加越多的趋势,且子流域面积越大,ΔCN 增加越多,设计洪水变化呈现越大的趋势.

3)研究区设计洪水对LUCC 响应最强的区域位于东部地区,西部地区次之,南部地区则最小,空间响应趋势随着重现期变高从中间向四周减弱.对于设计洪水变化较大的区域,可采取合理的退耕还林还草政策,减少洪涝灾害的发生.