彭楚杰
(湖南省湘潭水文水資源勘測中心,湖南 湘潭 411100)
在中國南方地區(qū),山丘區(qū)小流域分布廣泛。山丘區(qū)中小流域洪水呈現(xiàn)洪峰高、沖擊力大的特點[1],且流域內(nèi)調(diào)蓄能力有限,在遭遇強(qiáng)暴雨襲擊時,發(fā)生山洪災(zāi)害的可能性極高。所以,為了實現(xiàn)山洪災(zāi)害預(yù)警,快速推求小流域徑流過程顯得尤為重要。
大多數(shù)山丘區(qū)中小流域水文站網(wǎng)覆蓋率低,長系列水文監(jiān)測資料匱乏,概念性模型難以推廣使用,將GIS技術(shù)與水文學(xué)原理相結(jié)合的分布式水文模型成為主要的水文應(yīng)用發(fā)展方向[2]。如利用GIS等空間技術(shù)建立暴雨洪澇災(zāi)害綜合指數(shù)評價模型,確定暴雨洪澇風(fēng)險圖[3]。美國陸軍工程兵團(tuán)研發(fā)的HEC-HMS分布式水文模型具備操作簡單和免費(fèi)開源的特點,能整合多種產(chǎn)、匯流及河道演算方法[4-6],在印度蒂斯塔流域山洪災(zāi)害研究[7]、印度南部提魯奇拉帕利市無資料地區(qū)洪水風(fēng)險評估[8]以及阿爾及利亞梅克爾流域降徑過程模擬[9]等中都得到了較好的實踐。在中國針對HEC-HMS模型的研究主要有高程數(shù)據(jù)精度影響分析[10]、參數(shù)變化影響[11]和模型比對[12]等,同時也逐步應(yīng)用于洪水預(yù)報[13]、臨界雨量計算及山洪災(zāi)害預(yù)警[14-15]等領(lǐng)域。該模型利用網(wǎng)格結(jié)合DEM數(shù)據(jù)將流域劃分為多個子流域,能充分考慮流域內(nèi)地形與氣候因素的時空差異,還原流域內(nèi)產(chǎn)匯流物理過程,在資料稀缺流域也能得以適用[16]。
研究區(qū)域選取青山河青山橋水文站以上流域(112°24′E~112°36′E,27°36′N~27°24′N),青山河為涓水一級支流,全長32.0 km,流域集水面積277 km2,河流坡降0.244%;地形三面環(huán)山,東北、西北、西南地勢均較高。研究區(qū)域內(nèi)最高點海拔743 m,最低點海拔67 m(圖1)。研究區(qū)屬中亞熱帶季風(fēng)濕潤氣候區(qū),全年降水量為1 200~1 500 mm。山丘區(qū)總體植被覆蓋較好,丘陵區(qū)主要以農(nóng)業(yè)耕作區(qū)為主。青山河流域?qū)僦衼啛釒|部常綠闊葉林亞帶,按植被區(qū)系劃分,屬華中偏東亞系。青山橋站為中小流域自動監(jiān)測站,主要監(jiān)測項目有降水量、水位和流量信息。水文站上游無閘壩,開發(fā)利用程度不高。
圖1 青山河流域地形
DEM原始數(shù)據(jù)下載自地理空間數(shù)據(jù)云平臺,下載的數(shù)據(jù)源為ASTGTM2_N27E112,分辨率30 m。采用ArcGIS軟件中的空間分析工具提取出流域的柵格數(shù)據(jù),并進(jìn)行坐標(biāo)轉(zhuǎn)換與投影。運(yùn)用ArcGIS軟件中的HEC-GeoHMS拓展模塊,對提取出的流域柵格數(shù)據(jù)進(jìn)行填洼、流向、流量計算,確定提取閾值,提取出青山河的數(shù)字水系圖[17],見圖2。在此基礎(chǔ)上,將出口以上流域進(jìn)行子流域劃分。根據(jù)水系結(jié)構(gòu)生成河網(wǎng)拓?fù)浣Y(jié)構(gòu),并根據(jù)地形數(shù)據(jù)提取流域質(zhì)心,計算最長流路;初步選取產(chǎn)流參數(shù)與河道演算方法,計算坡度、河長、等部分模型計算參數(shù),導(dǎo)入到子流域?qū)傩詧D層。
圖2 青山橋水文站以上流域水系
研究區(qū)域內(nèi)有青山橋站降水量觀測資料,同時,周邊鄰近流域有橫鋪、赤石和黃巢3處降水量站,積累了較長序列降水量資料。利用4站降水量資料采用泰森多邊形法計算面降水量,雨量站分布泰森多邊形代表面積分配見圖3。青山橋站自2016年建站以來有部分實測洪水過程流量資料,可用于模型參數(shù)率定及預(yù)報檢驗,利用HEC-DSSVue數(shù)據(jù)處理工具將計算時段流量數(shù)據(jù)導(dǎo)入模型,進(jìn)行參數(shù)試算與精度對比。
圖3 青山河流域雨量站分布及泰森多邊形劃分
模型基礎(chǔ)部分包括網(wǎng)格文件、流域模型、氣象模型、控制模塊、時序數(shù)據(jù)五大塊。模型構(gòu)建流程見圖4。流域模型中包含了各子流域的坡度、高程等地理參數(shù)。在HEC-GeoHMS中根據(jù)流域地形將研究區(qū)域劃分為57個子流域,并生成網(wǎng)格,將流域產(chǎn)匯流過程概化成水庫、匯流、下滲等樹狀水系要素,概化結(jié)果見圖5。氣象模型可以輸入雨量站點的坐標(biāo)、選擇面雨量統(tǒng)計方法,可采用雨量站時段降水?dāng)?shù)據(jù)作為輸入??刂颇P驮O(shè)定計算的起止時間及步長,也是各時序文件與其他模塊單元的輸入端口;時序數(shù)據(jù)是降水量數(shù)據(jù)與流量數(shù)據(jù)輸入各模塊的載體。在本研究中,時序文件與計算步長統(tǒng)一設(shè)定為1 h,通過網(wǎng)格劃分單元計算各子流域的產(chǎn)流量,再進(jìn)行匯流與河道演算,最后推求出流域出口斷面處的徑流過程。
圖4 模型構(gòu)建流程
圖5 流域模型概化
模型中總體分徑流量計算、地表徑流模擬、基流與河道水流模擬4個部分。模型對各部分計算均提供了多種方法選擇。
3.2.1徑流量計算
集水區(qū)內(nèi)不透水表面無水量損失并產(chǎn)生徑流的部分,即為凈雨,也就是徑流的主要來源。模型提供了數(shù)種計算損失的方法,是以水量損失的不同計算形式來區(qū)分,本次選取初損常損法。
假設(shè)最大潛在降雨損失fc在整個降雨過程中為一恒定值。若pt是某時段內(nèi)的面平均降雨深度,則時段內(nèi)的凈降雨pet可表示如下:
(1)
引入初始損失量I,可用于表征截留與填洼水量,還包括蒸發(fā)與下滲等其他損失量。設(shè)pi為初始降雨深度,最終在集水區(qū)內(nèi)所累積的凈降雨可見式(2):
(2)
其中需要輸入的參數(shù)有初損值I、損失速率與不透水率。損失速率可以從1986年Skaggs和Khaleel發(fā)表的SCS土壤分類損失速率表中查得,不透水率可以從土壤組成及土地覆蓋類型數(shù)據(jù)中進(jìn)行估算,再通過模型試算確定。
3.2.2地表徑流模擬
地表徑流過程模擬采用斯奈德單位線計算。1938年斯奈德提出了參數(shù)化單位線的方法,給出了流域特征預(yù)測單位線參數(shù)的關(guān)系[18]。輸入的參數(shù)為洪峰時滯和峰值系數(shù)。洪峰時滯是單位線峰現(xiàn)時間與集水區(qū)質(zhì)心的時間差,在斯奈德給定的標(biāo)準(zhǔn)單位線時滯方程為:
tp=5.5tr
(3)
式中tr——降雨歷時;tp——洪峰時滯。
對于斯奈德單位線中的峰值系數(shù)Cp與時滯tp有如下相關(guān)關(guān)系:
(4)
式中Up——標(biāo)準(zhǔn)單位線的峰值;A——流域面積;C——轉(zhuǎn)換常數(shù)。
在實際模型計算中,對于單位線的洪峰時滯估算,可采用式(5):
(5)
式中Ct——流域區(qū)系數(shù);L——出口到分水點的主河道長度;Lc——出口到子流域質(zhì)心的長度;C同上。
峰值系數(shù)可用以表征單位降水產(chǎn)生的水文過程斜率,推薦取值0.3~0.8;洪峰時滯可定義為流域質(zhì)心到洪峰過程之間的時間間隔,均需要采用試算的方法確定。
3.2.3基流計算模擬
非線性boussinesq基流計算法利用假定boussinesq方程用來模擬流域降水后河道水流消退的過程。具備單場所次與連續(xù)暴雨的模擬。其所設(shè)計的參數(shù)較多,其優(yōu)勢在于大部分參數(shù)可以通過實測數(shù)據(jù)進(jìn)行確定,當(dāng)然也可以進(jìn)行試算求得。參數(shù)有初始基流量、孔隙率、滲透系數(shù)、流量系數(shù)臨界比、特征地下流經(jīng)長度。其中孔隙率與滲透系數(shù)可以根據(jù)土壤構(gòu)造及覆蓋情況進(jìn)行估算,特征地下流徑長度則可以按流域邊界到河流的平均距離進(jìn)行估算。
3.2.4河道水流模擬
模型提供的馬斯京根法基于槽蓄方程和水量平衡原理,將水面簡化為一個線性非水平面,使用分段連續(xù)演算的方式,其方程系具備二階精度差分格式,其涉及參數(shù)較少,應(yīng)用廣泛。
(6)
需要設(shè)定的參數(shù)有傳播時間K與權(quán)重系數(shù)X,可以從斷面與水流特性依據(jù)公式進(jìn)行估算,也可以通過試算得出。
(7)
X是入流和出流影響的權(quán)重系數(shù),范圍通常在0.0~0.5,亦可用式(8)進(jìn)行估算:
(8)
式中 ΔL——河段長度;n——糙率;λ——轉(zhuǎn)換系數(shù);S0——河流坡降;Q0——流量。
子河段數(shù)可根據(jù)實際河段長度與時間間隔進(jìn)行確定。
模型建立后,需要率定每一個模塊的各項參數(shù),模型中提供了試算的功能。選擇單純形法進(jìn)行迭代試算,需要為各參數(shù)選定一個初始值。目標(biāo)函數(shù)采用洪峰流量誤差百分比,設(shè)置容差范圍及對應(yīng)試算時段。采用2017—2019年5場單峰洪水過程數(shù)據(jù)進(jìn)行參數(shù)試算及率定。試算參數(shù)收斂性好,模型計算穩(wěn)定,最后得到一組最佳產(chǎn)匯流參數(shù),摘取部分參數(shù)見表1、2。
表1 流域產(chǎn)流參數(shù)率定成果
表2 河道參數(shù)率定成果
在2017—2019年中選取13場降水過程(含率定場次)進(jìn)行模型驗證。計算精度評定采用水情預(yù)報規(guī)范中要求的峰量相對誤差、峰現(xiàn)時間誤差與確定性系數(shù)(DC)[19]。由于采用的是洪峰誤差百分比函數(shù)率定的參數(shù),所以峰量相對誤差值均控制較好,其中洪峰流量誤差均達(dá)到了20%以內(nèi),11場達(dá)到10%以內(nèi),占總場次的85%;峰現(xiàn)時差均小于3 h,其中不大于1 h的9場,占總場次的69%;確定性系數(shù)只有2場小于0.7;綜合評定精度等級達(dá)到甲等有2場,乙等的有9場,丙等2場。滿足洪水預(yù)報的要求。洪水過程見圖6,評定結(jié)果見表3。
a)201705230000
表3 模型驗證成果
斯奈德單位線法核心是一個標(biāo)準(zhǔn)單峰單位線,從原理上說更適用于單峰洪水。結(jié)合徑流過程線擬合情況可看出峰后無雨或少雨的單峰洪水模擬效果普遍優(yōu)于復(fù)峰洪水,丙等洪水全為復(fù)峰洪水,洪號為H201806170900F的復(fù)峰洪水,確定性系數(shù)僅為0.56。洪號為H201705230000場次洪水,洪峰流量相對誤差較大,因當(dāng)月青山橋雨量站點存在故障,所采用雨量資料為參證所得,故有可能存在面雨量計算偏小的情況,以致于模擬結(jié)果偏小。通過對復(fù)峰洪水過程線的對比,模擬的洪水漲落尤其是消退段有較大誤差,在對消退系數(shù)進(jìn)行調(diào)整后作用有限,推測模型對于地表土壤層透水性與下層蓄水的模擬準(zhǔn)確度不高,故當(dāng)下層蓄水容量大時,并不能很好地還原基流消退情況。
利用DEM數(shù)據(jù)對青山河流域構(gòu)建分布式洪水預(yù)報模型,利用5場單峰洪水進(jìn)行參數(shù)率定,試算出一組最優(yōu)參數(shù),以13場洪水進(jìn)行驗證分析。
a)模型較好地還原了下墊面地形特征,將研究區(qū)域按網(wǎng)格化的方法分割為多個小流域進(jìn)行單獨(dú)計算,流域提取、模型生成、參數(shù)率定與檢驗均可在程序里完成,可操作性強(qiáng),計算便捷。
b)初估參數(shù)經(jīng)過實測洪水資料率定修正,率定過程中參數(shù)收斂無發(fā)散,模型計算穩(wěn)定。以峰量誤差、峰現(xiàn)時差、確定性系數(shù)來檢驗?zāi)P途?,其中洪峰流量誤差均達(dá)到了20%以內(nèi),峰現(xiàn)時差均小于3 h,確定性系數(shù)只有2場小于0.7,綜合評定結(jié)果表明,模型能滿足預(yù)報要求;總體上模型對于單峰洪水的模擬優(yōu)于復(fù)峰洪水,大部分過程線擬合良好,個別場次起漲段與消落段的計算有較大偏差,需要對參數(shù)做進(jìn)一步調(diào)整,從山洪預(yù)警的目的出發(fā),退水過程對于單峰洪水的影響大于復(fù)峰洪水,在實際應(yīng)用中需特別注意。
c)模型對于基礎(chǔ)資料的要求較高,特別是土壤利用數(shù)據(jù)及降水資料,現(xiàn)行條件下能獲取到的土壤利用數(shù)據(jù)分辨率對于小流域而言精度不高,模型計算結(jié)果將產(chǎn)生誤差;而降水資料的準(zhǔn)確性則直接影響計算結(jié)果的準(zhǔn)確性。
d)模型涉及的參數(shù)較多,每一步計算過程中所需要輸入的參數(shù),即便是具有實際物理意義的參數(shù),“多參同效”的可能性依然存在,所以需要更多實測資料進(jìn)行反復(fù)校正。未來將對參數(shù)敏感性進(jìn)行進(jìn)一步分析。
e)模型假定流域為天然狀態(tài),沒有考慮流域內(nèi)蓄水工程影響,實際情況下流域內(nèi)存在多處蓄水工程,這也是模型計算誤差的來源之一。若流域內(nèi)有大型蓄、引、調(diào)水工程存在,其影響不可忽視。同時,本次研究范圍為中小流域,模擬結(jié)果較好,下一步將探求在更大尺度流域上的應(yīng)用。