(河海大學(xué) 水文水資源學(xué)院,南京 210098)

土地利用/覆被變化(LUCC)會(huì)引起流域下墊面條件發(fā)生變化,導(dǎo)致流域產(chǎn)匯流規(guī)律發(fā)生不同程度的改變,直接影響暴雨洪水特性[1].目前,國(guó)內(nèi)外主要研究LUCC 的徑流響應(yīng),關(guān)于LUCC 的設(shè)計(jì)洪水響應(yīng)研究還相對(duì)較少.吳淼等[2]對(duì)淮河上游小流域構(gòu)建SWAT 模型進(jìn)行模擬,結(jié)果表明該流域在LUCC 前后,頻率越大,流域設(shè)計(jì)洪水減小越多.雷超佳等[3]采用HEC-HMS模型對(duì)奉化江皎口水庫(kù)流域進(jìn)行研究,分析了該流域不同重現(xiàn)期的洪水過(guò)程在LUCC影響下的變化情況.魏沖等[4]以深圳市石巖流域?yàn)檠芯繉?duì)象,經(jīng)分析得出,該流域不同頻率的設(shè)計(jì)洪水受城市化影響較劇烈.

近年來(lái),沂河上游流域土地利用方式發(fā)生改變,不透水面積呈現(xiàn)增加趨勢(shì),容易導(dǎo)致洪澇災(zāi)害頻繁發(fā)生,因此,迫切需要定量研究設(shè)計(jì)洪水對(duì)LUCC 的響應(yīng),了解下墊面變化對(duì)流域防洪具有重要的意義.本文基于HEC-HMS分布式水文模型,分析沂河上游流域1985-2015年LUCC 對(duì)設(shè)計(jì)洪水的影響,通過(guò)劃分響應(yīng)分區(qū)探討不同重現(xiàn)期下設(shè)計(jì)洪水對(duì)LUCC的空間響應(yīng),為流域防洪減災(zāi)提供一定的科學(xué)依據(jù),也為流域下游水庫(kù)安全運(yùn)行提供一定的參考.

1 研究區(qū)概況

沂河發(fā)源于山東省沂源縣田莊水庫(kù),上游為發(fā)源地至沂源縣東里店,地勢(shì)自西北向東南傾斜,以低山丘陵為主,干流河長(zhǎng)63.5km,河道平均比降2.7‰.沂河上游流域夏季高溫多雨,汛期集中在6～9月.源頭田莊水庫(kù)是大(二)型山區(qū)水庫(kù),水庫(kù)流域面積424 km2,水庫(kù)總庫(kù)容1.2億m3.研究區(qū)位于田莊水庫(kù)和東里店水文站之間,如圖1 所示,流域面積為746 km2,土地利用方式主要以林地、耕地為主,土壤類型為棕壤土、褐土和潮土.研究區(qū)以東里店水文站為控制站點(diǎn),除東里店站以外,流域還有悅莊、蘆芽店、石橋和焦家上莊降水站.

圖1 研究區(qū)位置圖

2 數(shù)據(jù)與方法

2.1 基礎(chǔ)數(shù)據(jù)

遙感影像數(shù)據(jù)采用研究區(qū)1985年Landsat5 TM和2015年Landsat8遙感影像,以及Sentinel-2高分辨率遙感影像;數(shù)字高程數(shù)據(jù)(DEM)采用30m 精度的ASTER-GDEM V2高程圖,以上數(shù)據(jù)均來(lái)自于地理空間信息云平臺(tái);土壤類型數(shù)據(jù)采用中國(guó)科學(xué)院寒區(qū)旱區(qū)環(huán)境與工程研究所1∶100 萬(wàn)的土壤柵格數(shù)據(jù);氣象水文數(shù)據(jù)采用研究區(qū)降水站和水文站獲得的小時(shí)降雨徑流數(shù)據(jù).

2.2 研究方法

HEC-HMS是美國(guó)陸軍工程師團(tuán)水文工程中心研發(fā)的降雨徑流模型,是一個(gè)具有物理概念的半分布式水文模型[5],主要由5個(gè)模塊構(gòu)成,分別是氣象模塊、產(chǎn)流計(jì)算模塊、坡面匯流計(jì)算模塊、河道基流計(jì)算模塊和河網(wǎng)匯流計(jì)算模塊[6].HEC-HMS模型能夠綜合考慮研究區(qū)下墊面條件的時(shí)空變化,不同的子流域可以選用不同的產(chǎn)匯流方式[7].為了更好地分析研究區(qū)LUCC對(duì)設(shè)計(jì)洪水的影響,本文選取SCS曲線法、SCS單位線法、馬斯京根法和衰減指數(shù)法分別作為HEC-HMS模型產(chǎn)流、坡面匯流、河網(wǎng)匯流和河道基流的計(jì)算方法.

1)SCS曲線法是美國(guó)水土保持局(SCS)開發(fā)用于流域產(chǎn)流的水文方法,SCS曲線法通過(guò)土地利用/覆蓋、累計(jì)雨量和前期土壤濕度來(lái)估算流域的凈雨量[8],能夠較好地反映土地利用變化對(duì)流域產(chǎn)流的影響,計(jì)算公式如下:

式中:Pe、P分別為t時(shí)刻的累計(jì)雨量、降雨深度(mm);S為子流域最大截留降雨能力(mm);CN 為SCS曲線法的曲線數(shù),取值范圍為30～100,與土壤滲透率有關(guān).

2)SCS單位線法通過(guò)流域滯時(shí)tlag反映,與流域土地利用條件、土壤類型、水流長(zhǎng)度和坡度有關(guān)[9],能夠較好地模擬土地利用變化下的坡面匯流情況,計(jì)算公式如下:

式中:tlag為子流域的洪峰滯時(shí)(min);L為水流長(zhǎng)度(m);S為子流域最大截留降雨能力(mm);y為子流域平均坡度(%);L和y可利用HEC-GeoHMS模型插件通過(guò)DEM 來(lái)提取.

3)馬斯京根法常用于山區(qū)河道流量演算,能較好地模擬研究區(qū)河網(wǎng)匯流的情況,其參數(shù)河段傳播時(shí)間K和流量比重因子X可由以下公式估算:

式中:a為洪水波速修正系數(shù),在1.33～1.67之間;Δx為河段長(zhǎng)度(m);c′為洪水平均波速(m/s);Q0為參考流量(m3/s);B為水流頂寬(m);S0為底坡度.

4)衰減指數(shù)法描述了基流從初始流量開始以指數(shù)的形式進(jìn)行衰減,可通過(guò)實(shí)測(cè)徑流退水部分估算衰減指數(shù),t時(shí)刻基流Qt的計(jì)算公式如下:

式中:Q0為初始基流(m3/s);k′為衰減指數(shù),取值范圍在0～1之間.

3 模型構(gòu)建與率定

3.1 研究區(qū)土地利用/覆被變化分析

通過(guò)ENVI采用神經(jīng)網(wǎng)絡(luò)法對(duì)1985年和2015年遙感影像進(jìn)行監(jiān)督分類,把研究區(qū)土地利用類型劃分為林地、耕地、草地(低矮灌木林、草地等)、水體、裸地和建設(shè)用地六大類型,如圖2所示.結(jié)合沂源各類歷史專題研究圖件、地形圖等相關(guān)資料并參考地圖軟件的歷史影像,對(duì)監(jiān)督分類結(jié)果進(jìn)行人工校正.通過(guò)Sentinel-2遙感數(shù)據(jù),采用混淆矩陣進(jìn)行精度評(píng)價(jià),總體精度分別達(dá)到82%、88%,解譯精度滿足研究要求.

圖2 研究區(qū)1985年和2015年土地利用分布圖

通過(guò)ENVI分析兩期的土地利用分布圖,得到研究區(qū)1985-2015年的土地利用轉(zhuǎn)移矩陣(見表1).

表1 研究區(qū)1985-2015年土地利用類型面積轉(zhuǎn)移矩陣 (單位:km2)

從表1可知,研究區(qū)的土地利用主要以耕地、林地、草地、建設(shè)用地為主,四者所占的總面積占研究區(qū)面積的97%以上,水體與裸地所占的面積相對(duì)偏少.相較于1985年,2015年的林地、耕地面積均有所下降,分別減少了24.44km2、32.68km2,所占比例分別從38.95%、46.76%下降到35.72%、42.43%;草地、建設(shè)用地面積均有所上升,分別增加了8.87km2、50.33km2,所占比例分別從7.01%、4.47%上升到8.14%、11.14%,建設(shè)用地變化最為顯著.

在1985-2015年土地利用轉(zhuǎn)移矩陣中,林地、耕地分別約有16km2、35km2轉(zhuǎn)變?yōu)榻ㄔO(shè)用地,27 km2、23km2轉(zhuǎn)變?yōu)椴莸?這主要是由于社會(huì)經(jīng)濟(jì)的發(fā)展,研究區(qū)30年以來(lái)不斷進(jìn)行城鎮(zhèn)化建設(shè),城鎮(zhèn)用地不斷增加,一部分的林地、耕地用于旅游用地,建設(shè)旅游景點(diǎn),發(fā)展當(dāng)?shù)氐穆糜螛I(yè),同時(shí),丘陵山地以果園為主的灌叢林廣泛種植,取代原來(lái)的喬木林用地與農(nóng)作物用地,提高當(dāng)?shù)氐慕?jīng)濟(jì)水平.由于原耕地的減少以及土地利用合理規(guī)劃,林地、草地分別約有67.15 km2、27.58km2轉(zhuǎn)變?yōu)楦?

從研究區(qū)LUCC 可知,林地、耕地面積減少,草地、建設(shè)用地增多,研究區(qū)的土地利用/覆被結(jié)構(gòu)發(fā)生了較大的變化,下墊面條件的改變會(huì)引起流域產(chǎn)匯流規(guī)律的變化,對(duì)各子流域的降雨徑流過(guò)程產(chǎn)生不同的影響.

3.2 模型基礎(chǔ)數(shù)據(jù)處理

利用HEC-GeoHMS 模塊對(duì)研究區(qū)DEM 進(jìn)行流域構(gòu)建、水系提取、流域劃分等一系列的地表水文分析,以東里店水文站作為流域出口,把研究區(qū)劃分為25個(gè)子流域,13條概化河道,如圖3所示,對(duì)每個(gè)子流域和概化河道進(jìn)行物理參數(shù)、水文參數(shù)的提取,并為每個(gè)子流域和概化河道選定產(chǎn)匯流計(jì)算方法.為了反映時(shí)空分布的差異性[10],本文采用距離平方倒數(shù)加權(quán)法,通過(guò)HEC-GeoHMS計(jì)算5個(gè)降水站在各個(gè)子流域的權(quán)重系數(shù).

圖3 研究區(qū)子流域劃分圖

CN 采用SCS制定的土壤分類標(biāo)準(zhǔn),根據(jù)土壤滲透率高低分為A、B、C、D 四類,通過(guò)中國(guó)土壤科學(xué)數(shù)據(jù)庫(kù)獲取研究區(qū)土壤物理性質(zhì),把原土壤類型進(jìn)行歸并,得出研究區(qū)水分土壤組為B、C和D 三組.根據(jù)研究區(qū)的土地利用情況,確定前期土壤濕度處于AMCⅡ平均等級(jí)下各地類的CN2S值[11-13],將土壤數(shù)據(jù)與土地利用數(shù)據(jù)進(jìn)行空間疊加分析,最終通過(guò)HECGeoHMS計(jì)算研究區(qū)內(nèi)各子流域的平均CN2S值.確定各地類的不透水率Im[14-16],對(duì)各個(gè)子流域中所有土地利用的不透水率進(jìn)行加權(quán)平均,估算每個(gè)子流域的不透水率.

3.3 模型率定與驗(yàn)證

由于研究區(qū)降水集中于6～9月,植物處于生長(zhǎng)期,在前期土壤濕度AMCⅡ平均等級(jí)下,前5天總雨量為35.6～56.3mm[17].根據(jù)這一特點(diǎn),選取研究區(qū)東里店水文站1979-2012年共10場(chǎng)實(shí)測(cè)徑流過(guò)程,用于模型的率定與驗(yàn)證.由于沂河源頭為田莊水庫(kù),把田莊水庫(kù)對(duì)應(yīng)場(chǎng)次洪水的下泄流量作為本模型的入流條件.時(shí)間步長(zhǎng)的長(zhǎng)短會(huì)影響模擬值與實(shí)測(cè)值的擬合程度[18],因此,通過(guò)HEC-DSSVue的插值功能,對(duì)徑流資料存在的數(shù)據(jù)缺失值,按照15min的時(shí)間步長(zhǎng)進(jìn)行插值應(yīng)用.

選取峰值加權(quán)標(biāo)準(zhǔn)誤差方程作為擬合目標(biāo)函數(shù)[19],采用洪峰流量相對(duì)誤差ΔQ、徑流深相對(duì)誤差ΔR、峰現(xiàn)時(shí)間誤差ΔT和Nash效率系數(shù)E共4個(gè)檢驗(yàn)標(biāo)準(zhǔn),校檢HEC-HMS 在研究區(qū)的洪水模擬,1985年和2015年主要參數(shù)率定優(yōu)化結(jié)果見表2,兩者模擬的結(jié)果見表3.

表2 模型主要參數(shù)率定優(yōu)化結(jié)果

表3 研究區(qū)水文模擬結(jié)果

由表2可知,在模擬的10場(chǎng)洪水過(guò)程中,洪峰流量相對(duì)誤差和徑流深相對(duì)誤差均在20%范圍以內(nèi),峰現(xiàn)時(shí)間誤差均在3h以內(nèi),Nash效率系數(shù)E均在0.84以上.根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》(GB/T 22482—2008)的規(guī)定,滿足洪峰和徑流深相對(duì)誤差在20%以內(nèi),峰現(xiàn)時(shí)差在3h以內(nèi),確定性系數(shù)在0.7以上即為合格.因此,本文構(gòu)建的HEC-HMS模型在東里店流域洪水模擬的總合格率為100%,同時(shí),實(shí)測(cè)與模擬的洪水過(guò)程線擬合良好,如圖4所示,HEC-HMS半分布式水文模型模擬研究區(qū)降雨徑流的效果較好.

圖4 實(shí)測(cè)與模擬的洪水過(guò)程

4 LUCC下研究區(qū)設(shè)計(jì)洪水響應(yīng)特性分析

4.1 LUCC對(duì)設(shè)計(jì)洪水的影響分析

為了分析研究區(qū)LUCC 對(duì)設(shè)計(jì)洪水的影響,采用研究區(qū)內(nèi)5個(gè)站點(diǎn)1968-2016年共49年的降雨資料系列,選取不同年份各站點(diǎn)最大24h暴雨,采用P-Ⅲ型頻率曲線進(jìn)行排頻計(jì)算.

根據(jù)研究區(qū)所在的地理位置,選用山東省泰沂南山區(qū)1h雨型對(duì)設(shè)計(jì)雨量進(jìn)行逐日設(shè)計(jì)暴雨時(shí)程分配得到P=10%、P=5%、P=2%和P=1%的逐時(shí)設(shè)計(jì)暴雨.通過(guò)HEC-HMS模擬計(jì)算兩種下墊面條件下研究區(qū)不同重現(xiàn)期的設(shè)計(jì)洪水,流域出口模擬結(jié)果見表4,設(shè)計(jì)洪水過(guò)程如圖5所示,各子流域洪水變化見表5.

表4 研究區(qū)流域出口斷面模擬結(jié)果

表5 研究區(qū)子流域洪水變化率 (單位:%)

續(xù)表5 研究區(qū)子流域洪水變化率 (單位:%)

圖5 不同下墊面條件下的設(shè)計(jì)洪水過(guò)程

由模擬結(jié)果可知,2015年與1985年下墊面條件模擬的洪水相比,研究區(qū)P=10%、P=5%、P=2%、P=1%的洪峰和徑流深均有不同程度的增加.由圖5 可知,峰現(xiàn)時(shí)間分別提前了105、105、90、90 min.可以發(fā)現(xiàn),重現(xiàn)期越低,研究區(qū)洪峰和徑流深的的變化就越大,峰現(xiàn)時(shí)間越提前.

由表4可知,重現(xiàn)期由低到高,流域出口的徑流深分別增加了6.18、5.63、6.19、6.58mm,洪峰分別增加了68.08、74.43、81.08、85.34m3/s.可以發(fā)現(xiàn),研究區(qū)LUCC 在增加相同徑流深的情況下,重現(xiàn)期越高,洪峰增加得越多.

由表5可知,研究區(qū)各子流域不同重現(xiàn)期下設(shè)計(jì)洪水變化情況,石橋鎮(zhèn)北部(W670)的設(shè)計(jì)洪水變化最大,西里鎮(zhèn)北部(W1060)的變化最小.南麻街道北部(W630)、石橋鎮(zhèn)(W820、W890)、悅莊鎮(zhèn)南部(W700)、沂河沿岸(W900、W990)和燕崖鎮(zhèn)西部(W950、W960)的設(shè)計(jì)洪水變化較大,結(jié)合土地利用變化分析得知,前三者由于城鎮(zhèn)化建設(shè),城市化向外擴(kuò)張,建設(shè)用地逐漸增加,林地耕地轉(zhuǎn)變?yōu)椴煌杆玫?流域產(chǎn)流能力增加,導(dǎo)致它們的設(shè)計(jì)洪水產(chǎn)生較大變化;后兩者設(shè)計(jì)洪水變化較大的原因,主要是由于林地向耕地、草地轉(zhuǎn)變,林地面積減少,人類活動(dòng)影響增加,流域涵養(yǎng)水源能力降低.

通過(guò)比較研究區(qū)洪峰與徑流深可知,洪峰變化小于徑流深變化,在忽略基流徑流深時(shí),洪峰變化則大于直接徑流深的變化.通過(guò)比較子流域的面積、ΔCN和設(shè)計(jì)洪水變化率可知,子流域面積越大,ΔCN 增加越多,研究區(qū)設(shè)計(jì)洪水變化就越大.

4.2 LUCC下設(shè)計(jì)洪水的空間響應(yīng)分析

為了分析LUCC 下設(shè)計(jì)洪水的空間響應(yīng),提供研究區(qū)防洪減災(zāi)依據(jù),本文根據(jù)洪峰變化率和CN 變化值的關(guān)系趨勢(shì)以及變化范圍,選定ΔQ/Q1985=10%和ΔCN=1作為分區(qū)判定條件,把研究區(qū)分為3個(gè)分區(qū).當(dāng)ΔQ/Q1985＞10%且ΔCN＞1時(shí)為1區(qū),當(dāng)ΔQ/Q1985＜10%且ΔCN＞1時(shí)為2 區(qū),當(dāng)ΔQ/Q1985＜10%且ΔCN＜1時(shí)為3區(qū).不同重現(xiàn)期下的空間響應(yīng)分區(qū)如圖6所示.由圖6可知,隨著重現(xiàn)期由低到高,設(shè)計(jì)洪水對(duì)LUCC的響應(yīng)從中部地區(qū)(沂源縣主城區(qū))向周邊地區(qū)逐漸減弱,50年一遇和100年一遇的空間響應(yīng)基本一致,出現(xiàn)這種趨勢(shì)主要是由于城市化向外擴(kuò)張,周邊鄉(xiāng)鎮(zhèn)地區(qū)的林地面積大幅度減少導(dǎo)致的.設(shè)計(jì)洪水對(duì)LUCC 響應(yīng)較強(qiáng)的1區(qū),主要分布于東部地區(qū),集中于石橋鎮(zhèn)區(qū)域,建設(shè)用地明顯增加;西部地區(qū)隨著重現(xiàn)期的變高,由1區(qū)轉(zhuǎn)為2區(qū),該區(qū)域設(shè)計(jì)洪水對(duì)LUCC 的響應(yīng)逐漸減弱,敏感性降低,林地向耕地、草地轉(zhuǎn)變的同時(shí),也存在耕地向草地轉(zhuǎn)變;而響應(yīng)較弱的3區(qū),主要分布于南北部地區(qū),南部地區(qū)響應(yīng)最弱,該地區(qū)30年間土地利用變化較小,主要是林地、耕地向草地的轉(zhuǎn)變.

圖6 不同重現(xiàn)期下空間響應(yīng)分區(qū)圖

5 結(jié)論

1)1985-2015年間,研究區(qū)土地利用類型以林地、耕地、草地和建設(shè)用地為主,占總面積的97%以上.相較于1985年的土地利用,2015年時(shí)的林地和耕地面積均減少,草地和建設(shè)用地的面積均有所增加,建設(shè)用地的增幅最為顯著.

2)在1985年和2015年的不同下墊面條件模擬下,研究區(qū)內(nèi)不同重現(xiàn)期的洪峰和徑流深均有不同程度的增加,峰現(xiàn)時(shí)間提前,低重現(xiàn)期的洪水變化最顯著.在增加相同徑流深時(shí),重現(xiàn)期越高,洪峰呈現(xiàn)增加越多的趨勢(shì),且子流域面積越大,ΔCN 增加越多,設(shè)計(jì)洪水變化呈現(xiàn)越大的趨勢(shì).

3)研究區(qū)設(shè)計(jì)洪水對(duì)LUCC 響應(yīng)最強(qiáng)的區(qū)域位于東部地區(qū),西部地區(qū)次之,南部地區(qū)則最小,空間響應(yīng)趨勢(shì)隨著重現(xiàn)期變高從中間向四周減弱.對(duì)于設(shè)計(jì)洪水變化較大的區(qū)域,可采取合理的退耕還林還草政策,減少洪澇災(zāi)害的發(fā)生.