魏 沖，董曉華,2，成 潔，時 磊，劉 冀，喻 丹，楊芝辰，譚雪松

(1.三峽大學(xué)水利與環(huán)境學(xué)院，湖北 宜昌 443002；2.水資源安全保障湖北省協(xié)同創(chuàng)新中心，武漢 430072；3.深圳市水務(wù)規(guī)劃設(shè)計院有限公司，深圳 518001；4.宜昌晟泰水電實業(yè)有限責任公司，湖北 宜昌 444299)

0 前 言

隨著流域內(nèi)城市的飛速發(fā)展，該流域的土地利用會有所改變。人類在流域內(nèi)進行諸如流域河道整治、水庫塘壩以及跨流域調(diào)水等水利工程建設(shè)的各種活動，必然會對流域下墊面產(chǎn)生一定的影響[1]，使土地利用/覆被狀況發(fā)生變化，從而改變流域蒸發(fā)截留、填洼下滲等產(chǎn)匯流條件[2-4]，導(dǎo)致一系列的水文生態(tài)效應(yīng)。研究流域內(nèi)土地利用變化對水文過程的影響對流域內(nèi)城市的發(fā)展有重要指導(dǎo)意義。

國內(nèi)外學(xué)者應(yīng)用各種形式的方法對土地利用變化的響應(yīng)進行研究，而應(yīng)用最廣的還是能考慮流域自然條件的分布式水文模型[5]，它可以體現(xiàn)出流域的各種物理條件對水文過程的影響，并可以結(jié)合遙感與GIS，靈活設(shè)置不同的土地利用情景，從事更精細的研究[6]。SWAT模型開發(fā)于20世紀90年代[7]，是近年來在土地利用變化的水文響應(yīng)研究領(lǐng)域應(yīng)用比較成熟的模型[8]。王強[9]等人在太湖上游西苕溪流域構(gòu)建SWAT模型，應(yīng)用GWR模型在空間上定量評估了土地利用變化對徑流過程的影響；董曉華[10]等人在淮河上游構(gòu)建SWAT模型，通過設(shè)定不同的土地利用模式來模擬未來土地利用情景下的水文響應(yīng)；吳淼[11]等對潢川上游流域構(gòu)建SWAT模型，模擬流域不同土地利用下最大幾日設(shè)計暴雨洪水過程，研究結(jié)果表明：該流域土地利用變化對設(shè)計洪水過程的影響比較顯著。但是洪水過程與徑流過程相比，具有時間短、強度大和產(chǎn)流快等特性，需要更為精細的模型進行模擬。雷超桂[12]等人在奉化江皎口水庫流域構(gòu)建HEC-HMS水文模型，分析不同頻率的暴雨洪水對下墊面變化的響應(yīng)情況，指出暴雨洪水受土地利用變化的影響程度不一，對設(shè)計頻率高的洪水影響更明顯；毛慧慧[13]等人基于新安江-海河模型，分析了永定河上游土地利用變化對設(shè)計洪水的影響；馮平[14]等人對灤河流域下游構(gòu)建流域水文模型，通過分析土地利用變化對洪水過程的影響，發(fā)現(xiàn)水土保持工程對洪水過程的峰值及總量都有一定削減作用。

城市化會影響該地區(qū)的產(chǎn)匯流過程，從而加大其受到洪澇災(zāi)害的風(fēng)險[15]，了解下墊面變化對流域設(shè)計洪水的影響對城市防洪以及城市發(fā)展有重大意義。深圳市石巖水庫建于1960年3月，同年設(shè)立石巖水庫雨量站，距今已經(jīng)有近60年的歷史。近年來深圳市發(fā)展迅猛，城市化過程非常劇烈，水庫修建之初使用的設(shè)計暴雨及設(shè)計洪水是否適用于現(xiàn)在，亟待研究。本文對石巖流域1988-2016年間四期的土地利用遙感圖進行解譯，構(gòu)建改進后的SWAT模型，模擬該流域不同設(shè)計頻率的洪水過程。以此對不同時期的土地利用情景下不同頻率設(shè)計暴雨對應(yīng)的設(shè)計洪水進行模擬，研究其變化情況，為城市防洪管理及水庫管理提供理論支持。

1 研究區(qū)概況

石巖水庫位于深圳市石巖鎮(zhèn)境內(nèi)，地理位置為東經(jīng)113°54′，北緯22°41′，流域面積約44 km2(圖1)，庫容為3 120 萬m3，水面面積約2.98 km2。庫區(qū)位于亞熱帶常綠林海洋性季風(fēng)氣候區(qū)，氣候特征為均溫高，日照時間長，無霜期時間長，多雨，且雨量在年內(nèi)不均。根據(jù)資料統(tǒng)計，該地年均氣溫22.4 ℃，最高溫38.7 ℃，最低溫0.2 ℃；年均日照時數(shù)2 120 h，年均蒸發(fā)量1 322 mm，流域年均降雨量為1 628.3 mm，平均每年4-9月占其年降雨量的85%左右。石巖水庫有6條支流，石巖河為其干流。

圖1 石巖流域概況圖Fig.1 Overview diagram of Shiyan River basin

2 數(shù) 據(jù)

2.1 土地利用數(shù)據(jù)

本研究所使用的土地利用圖均來自對研究區(qū)域遙感影像數(shù)據(jù)進行解譯。在土地利用/覆被分類的相關(guān)研究中，通常使用監(jiān)督分類法和非監(jiān)督分類法。本研究選用貝葉斯監(jiān)督分類[16]，一種應(yīng)用極廣的監(jiān)督分類方法。本文使用的遙感影像資料源于地理空間數(shù)據(jù)云提供的Landsat4-5和Landsat8的衛(wèi)星數(shù)字產(chǎn)品，為減小云量和植被季節(jié)性變化給解譯帶來的誤差，綜合考慮這兩方面的因素，選取成像時間為1988年11月24日、2002年11月7日，2010年12月23日和2016年2月7日的影像數(shù)據(jù)進行解譯，此四期遙感影像的云量少，且植被處在相近生長周期。數(shù)據(jù)經(jīng)過衛(wèi)星系統(tǒng)自動校正后，選擇2、3、4三個波段組合進行影像合成，再進行監(jiān)督分類，得到1988、2002、2010和2016年的土地利用圖。結(jié)果見圖2。

圖2 石巖流域土地利用圖Fig.2 Landuse map of Shiyan River basin

2.2 水文氣象數(shù)據(jù)

本研究選用氣候預(yù)報系統(tǒng)CFSR(Climate Forecast System Reanalysis)提供的研究區(qū)1979-2010年共32 a的降雨、溫度、太陽輻射、相對濕度等氣象數(shù)據(jù)，計算天氣發(fā)生器所需的參數(shù)。DEM數(shù)據(jù)從地理空間數(shù)據(jù)云獲取，土壤數(shù)據(jù)來源于世界和諧土壤數(shù)據(jù)庫(Harmonized World Soil Database-HWSD)1∶100萬土壤數(shù)據(jù)(聯(lián)合國糧農(nóng)組織(Food and Agriculture Organization-FAO)和維也納國際應(yīng)用系統(tǒng)研究所(International Institute for Applied System Analysis-IIASA)所構(gòu)建，土地利用數(shù)據(jù)選用通過遙感解譯所得到的1988、2002、2010和2016年土地利用數(shù)據(jù)。

基于石巖水庫雨量站1963-2010年的歷史實測暴雨資料，采用P-III型頻率曲線法進行適線分析，得到設(shè)計頻率P=0.05%、0.1%、1%、2%、5%下的1、6、24、72 h年最大設(shè)計暴雨，結(jié)果見表1；再用推理公式求相應(yīng)頻率的設(shè)計洪水過程，用于率定模型參數(shù)。

表1 設(shè)計暴雨計算結(jié)果 mm

3 SWAT模型的構(gòu)建

SWAT(Soil and Water Assessment Tool)模型是由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心的研究團隊于1994年研發(fā)的分布式水文模型[17]。該模型以日尺度為模擬步長，主要適用于大尺度流域的長序列徑流模擬。

流域水文過程主要分為兩個階段：一是水文循環(huán)陸相階段；二是水文循環(huán)匯流階段。SWAT模型通過不同的土地利用、土壤類型及坡度組合將流域離散成若干個水文響應(yīng)單元(HRU)，以此反映不同作物與土壤的蒸散發(fā)差異。模型分別計算每個HRU內(nèi)的產(chǎn)流量，然后將每個子流域內(nèi)所有HRU的產(chǎn)流量疊加，得到各個子流域的總產(chǎn)流量。再進行子流域坡面匯流，計算其主河道的水量，最后通過特征河長法或馬斯京根法來進行匯流階段的運算。

由于洪水事件的歷時相對較短，需要更加精細的模擬步長來捕捉流量的瞬時變化，故以日為模擬步長的SWAT模型無法很好的應(yīng)用于洪水模擬。本文對SWAT模型進行修改，建立研究區(qū)小時尺度的SWAT模型。

3.1 SWAT模型的改進

模型有兩種方法來估算地表徑流：SCS曲線法以及Green&Ampt下滲法，都是以日為默認的時間步長進行演算，故只能進行日尺度長序列的水文模擬。SCS曲線法是20世紀50年代在美國鄉(xiāng)村廣泛使用的經(jīng)驗?zāi)Ｐ停鳪reen&Ampt下滲法是在假設(shè)土壤剖面均質(zhì)、前期水分在剖面中均勻分布的前提下連續(xù)模擬下滲過程的一個物理模型，具有一定的物理意義。

本文將Green&Ampt下滲法[18]作為計算日以下時間步長地表徑流的方法，對模型進行改進。其計算公式見式(1)。

(1)

式中：Finf,t是給定時間步長下的累積下滲量，mm；Ke是有效滲透系數(shù)，mm/h；ψwf是濕潤峰處的基質(zhì)勢，mm；Δθv是濕潤峰處土壤體積含水量的變化量，mm/mm。

當降雨強度小于下滲率時，時段內(nèi)的降雨全部滲入土壤，地表徑流量為0；當降雨強度大于下滲率時，按土壤的下滲能力下滲，其余部分形成地表徑流。對各個時段下滲雨量積分，得到累積下滲量，并計算地表徑流量。

3.1.1 地表徑流滯后模塊的修改

由Green&Ampt下滲法計算的時段地面凈雨通過一個滯后方程[式(2)]，將有部分凈雨被滯留在HRU中，在后一時段再匯入河道中。該方程通過滯后系數(shù)surlag和匯流時間tconc計算滯留的地面凈雨。而SWAT2005中的滯后方程基于日時間步長，并不適用于日以下時間步長。

(2)

(3)

式中：Δt是模擬步長，h；其余符號意義同前。

如圖3可知，在匯流時間一定的情況下，當滯后系數(shù)取值從1～10，隨著模擬步長Δt的增大，tconc/Δt逐漸減小，而匯流因子逐漸增大，匯入主河道的地表徑流量也逐漸增大，符合“隨著時間間隔的增加，將會有更多的地表徑流匯入主河道”這一規(guī)律。

圖3 匯流因子曲線圖Fig.3 Confluence factor graph

3.1.2 河網(wǎng)匯流模塊的修改

SWAT模型主河道匯流一般采用兩種方法：特征河長法和馬斯京根法。為方便計算，本文選用特征河長法，該方法將河段復(fù)雜的調(diào)蓄作用簡化成線性關(guān)系(式(4))。

Vout,2=SC(Vin+Vstored,1)

(4)

式中：Vout,2表示時段末的出流量，m3；Vin表示時段內(nèi)的入流量，m3；Vstored,1表示時段初的蓄水量，m3；SC表示蓄水系數(shù)。

蓄水系數(shù)與模擬時間步長之間存在式(5)的關(guān)系，且在河道演算過程中存在的諸如滲漏損失、蒸發(fā)損失及河岸調(diào)蓄等過程的模擬，模型均以日時間步長進行模擬。修改模型代碼中的模擬步長，使其適用于任意Δt步長的模擬。

(5)

式中：TT表示水流的運動時間，h；Δt表示模擬時間步長，h。

修改前的模型，在時間空間上各存在三個循環(huán)層，即時間上的Δt→日→年和空間上的HRU→子流域→流域。修改后，取消原代碼中時間上的年循環(huán)部分，將其替換成連續(xù)的場次循環(huán)。每場洪水的初始條件，包括初始土壤含水量、初始深層以及淺層含水層儲量、初始河道儲量，均是從模型的日模結(jié)果中提取?？臻g上，仍以HRU為最小計算單位，分別計算各個HRU的產(chǎn)水量。其中，冠層截留、地表徑流、坡面匯流和河道流量演算以Δt為時間步長模擬，但是融雪、基流、側(cè)流、蒸發(fā)以及含水層儲水是以日為時間步長計算，然后平均到每個時間步長。隨后匯總至子流域循環(huán)層進行坡面匯流，最后進入流域循環(huán)層進行河道流量演算。計算完每個Δt時間層上的循環(huán)后，再以同樣的方法進行日循環(huán)和年循環(huán)，直至模型模擬期結(jié)束。

3.2 SWAT模型參數(shù)敏感性分析及自動率定方法

SWAT模型中常見的對徑流模擬有影響的參數(shù)有26個，只有部分參數(shù)對模型的不確定性具有較大程度的影響，需對參數(shù)進行不確定性分析[19]，選出敏感性較強的參數(shù)來進行率定，以降低模型的不確定性。SWAT模型采用結(jié)合了LH抽樣法的穩(wěn)定性和OAT算法的精確性的LH-OAT算法[20]進行敏感性分析，可準確獲取對模型結(jié)果影響較大的重要輸入?yún)?shù)，提高模型適用范圍。模型選用SCE-UA算法[21]率定模型參數(shù)，SCE-UA算法通過搜尋所有選定參數(shù)的可行范圍，盡可能獲取全局最優(yōu)參數(shù)組。

本研究使用2010年土地利用和五種不同頻率的設(shè)計暴雨及設(shè)計洪水作為模型輸入，分析模型所有徑流參數(shù)的敏感性，并對敏感性較強的參數(shù)進行率定。將2010年9月期間發(fā)生的實測小時降雨過程輸入模型，模擬得到的小時徑流過程聚合成日徑流，與實測日徑流數(shù)據(jù)進行對比，驗證模型精度。(實測日徑流數(shù)據(jù)由水庫水位及取用水記錄，根據(jù)水量平衡反推得到。)

敏感性排名前十的參數(shù)及其自動率定結(jié)果見表2，參數(shù)率定結(jié)果見表3，參數(shù)驗證結(jié)果見表4和圖4。

表2 敏感性參數(shù)及最優(yōu)值Tab.2 Sensitive parameters and the optimal value

表3 模型率定結(jié)果Tab.3 Model calibration result

表4 模型驗證結(jié)果Tab.4 Model validation result

圖4 模型驗證洪水過程Fig.4 The flood process of model validation

用于驗證的三場實測洪水過程，第一場的模擬結(jié)果和實測結(jié)果相差較大，模擬的洪水過程有兩天存在較大洪量，而實測洪水過程則只有一天存在較大洪量。由于實測洪水過程由水庫水位數(shù)據(jù)反推得到，可能在實際過程中存在相應(yīng)的泄洪措施，使得實測的徑流過程在有降雨的情況下沒有洪峰，模擬前由模型日模結(jié)果計算的初始條件可能存在些許偏差，影響到第一場暴雨的模擬結(jié)果，使結(jié)果產(chǎn)生較大誤差；而第二場和第三場實測洪水過程的模擬精度很高，均合乎《水文情報預(yù)報規(guī)范》(GB/T 22482-2008)規(guī)定的要求。綜合考慮，該模型模擬結(jié)果可以接受，模型可用于實際應(yīng)用中。

4 結(jié)果及分析

4.1 土地利用變化情況

利用ArcGIS，統(tǒng)計流域各個時期各類土地利用的面積和比例，結(jié)果見表5。

由表5可知，林草地始終是石巖流域的主要土地利用類型，1988-2010年均超過流域一半的面積，但總體呈遞減趨勢，而城鎮(zhèn)用地則呈持續(xù)遞增趨勢。1988-2016年間，林草地減少39%，耕地減少80%，城鎮(zhèn)用地增加3 380%，未利用地增加93%。其中，1988-2002年間的土地利用變化最為劇烈，林草地減少27%，耕地減少61%，而城鎮(zhèn)用地增加1 200%，未利用地增加440%；2002-2010年，林草地和耕地分別持續(xù)減少約10%和88%，城鎮(zhèn)用地增加約114%，而未利用地減少30%；2010-2016年則相對平穩(wěn)，只有城鎮(zhèn)用地和未利用地發(fā)生了較大程度的變化。

表5 流域各期土地利用狀況Tab.5 Basin periods of landuse status

耕地和林草地的持續(xù)減少，以及城鎮(zhèn)用地的不斷增加，標志著該流域城市化進程正迅猛發(fā)展，下墊面條件也隨之發(fā)生劇烈變化。下墊面的改變直接導(dǎo)致洪水特性發(fā)生改變，適用于以前的設(shè)計洪水在下墊面發(fā)生改變后也會產(chǎn)生相應(yīng)的變化。

4.2 不同土地利用的水文響應(yīng)模擬結(jié)果

將五場不同頻率的設(shè)計暴雨分別輸入到不同土地利用條件的SWAT模型中,分別模擬四種不同時期土地利用下設(shè)計暴雨產(chǎn)生的洪水過程,并將之與對應(yīng)的設(shè)計洪水過程進行對比，分析比較不同時期土地利用變化對不同頻率的暴雨洪水的影響規(guī)律。模擬結(jié)果如圖5。

圖5 不同土地利用下設(shè)計洪水模擬結(jié)果Fig.5 The design flood simulation result under different landuse

由圖5可知，同一頻率的暴雨在四種不同時期土地利用情景下產(chǎn)生的洪水過程有明顯差異。隨時間的推進，同頻率暴雨下產(chǎn)生的洪水總量與洪峰流量隨之增加，且洪水起漲時間及峰現(xiàn)時間均呈逐漸提前的趨勢；2010年與2016年情景下洪水過程線趨勢相似，洪峰及洪量基本都接近設(shè)計洪水過程，洪水起漲階段歷時短且峰陡，退水階段歷時長且相對平緩；隨著設(shè)計暴雨重現(xiàn)期的增大，2010年與2016年情景下洪水模擬過程逐漸趨近于實際設(shè)計洪水過程。

通過對不同年份的土地利用下暴雨洪水進行模擬，分析其結(jié)果，見表6和表7。

表6 不同土地利用情景下暴雨洪水模擬洪峰及變化量 萬m3

表7 不同土地利用情景下暴雨洪水模擬洪量及變化量 萬m3

根據(jù)表6和表7可得，同等頻率的暴雨在1988-2016年四期土地利用情境下產(chǎn)生的洪水峰值呈逐次增大趨勢。2002年模擬結(jié)果洪峰相對1988年模擬結(jié)果變化量與2016年模擬結(jié)果平均值相對于2010年模擬結(jié)果相近，分別為45.6和36 m3/s；而2010年模擬結(jié)果的洪峰相對于2002年則有巨大改變，其峰值平均相差205.6 m3/s。洪水重現(xiàn)期越短，土地利用變化下其洪峰和洪量的變化量越大。洪量與洪峰具有相同趨勢，同等頻率下暴雨在1988-2016年四期土地利用情境下產(chǎn)生的洪水洪量也呈逐步上升趨勢。其中2002年模擬結(jié)果平均值相對1988年來說洪量僅增加了23.8 萬m3，而2002-2010年則增加了634.2 萬m3，2016年相對2010年也增加了84 萬m3。

結(jié)合土地利用變化分析，1988-2002年與2002-2010年，林草地和耕地均大量減少，城鎮(zhèn)用地均增加6 km2以上，但前者新增了6.6 km2的未利用地，后者則減少了2.4 km2；相比之下，后者的徑流變化量遠大于前者，其洪峰平均變化量為205.6 m3/s，洪量平均變化量為634.2 萬m3；以此分析除卻林草地和耕地有增大流域暴雨下滲量外，未利用地也對增大流域下滲量有積極作用。而2010-2016年間，林草地減少的速度變緩，耕地面積新增0.9 km2，城鎮(zhèn)面積增加3.5 km2，同時未利用地減少2.8 km2，其間洪峰平均變化量為36 m3/s，洪量平均變化量為84 萬m3，與分析結(jié)果相同。

5 結(jié) 語

本文分析了深圳市石巖流域1988-2016年間的土地利用變化情況，構(gòu)建改進后的SWAT模型，模擬不同時期土地利用下設(shè)計暴雨的洪水過程，并分析其變化情況。主要結(jié)論如下：

(1)改進后的SWAT模型在石巖流域的洪水模擬中具有良好的適用性。其模擬精度均達到《水文情報預(yù)報規(guī)范》(GB/T 22482-2008)規(guī)定的要求，可用于該流域的洪水模擬研究。

(2)1988-2016年間，石巖流域下墊面發(fā)生了比較劇烈的變化。流域內(nèi)林草地和耕地大量減少，城鎮(zhèn)用地大量增加，且變化主要集中發(fā)生在1988-2010年間。

(3)隨著石巖流域不同時期下墊面的變化，流域設(shè)計洪水過程也發(fā)生了劇烈的變化。流域內(nèi)大量的林草地和耕地向城鎮(zhèn)用地發(fā)展，改變了流域下墊面的滲透系數(shù)，進而影響其產(chǎn)匯流過程。在下墊面變化比較劇烈的1988-2010年間，各個頻率的設(shè)計暴雨所產(chǎn)生的設(shè)計洪水洪峰均增加了一倍左右，洪量平均增加了77%，且設(shè)計頻率越低的設(shè)計洪水所受到的影響大于設(shè)計頻率較高的設(shè)計洪水，增加了其洪水風(fēng)險。