考慮動態(tài)匯流路徑的時變分布式單位線

易 彬,陳 璐

(1.華中科技大學土木與水利工程學院,湖北 武漢 430074；2.華中科技大學數(shù)字流域科學與技術湖北省重點實驗室,湖北 武漢 430074)

流域匯流是地面徑流、壤中流、地下徑流匯集到流域出口斷面的徑流過程。凈雨滴通過不同介質、不同路徑、不同粗糙面向流域出口斷面移動,實際匯流過程十分復雜。過去100 a間,大量學者從水動力、水文、水滴運動以及系統(tǒng)視角觀察和分析流域匯流過程[1]。1932年,Sherman[2]在研究如何用降水量推求徑流量的過程中,提出了單位線的概念,目前該方法已成為使用最為廣泛的匯流計算方法。

單位線將匯流過程概化為線性時不變系統(tǒng),實際上,匯流特性不僅受控于流域地形地貌和下墊面特征,也與凈雨時空分布等因素有關[3]?；诹饔虻孛矃?shù),Rodríguez-Iturbe等[4]提出了地貌單位線方法,但該類方法無法考慮流域地形及下墊面特征的空間分布；為了考慮下墊面空間異質性對匯流的影響,Maident等[5]提出了基于數(shù)字高程模型的分布式匯流單位線。分布式單位線計算方法將水質點在流域中的滯留時間等價為流域瞬時單位線,其核心在于速度場的計算,常用的坡面流速計算公式有SCS流速公式、曼寧公式、均勻流方程、達西公式等[6],上述流速公式考慮了植被、坡度、土地利用等因素的影響[7-8],但忽略了由時變雨強、上游入流貢獻等因素引起的徑流匯流非線性問題。

有學者通過建立坡面流速與時段凈雨強度、上游入流貢獻等因素的定量映射關系,提出了時變分布式單位線法[9-10]?？追舱艿萚11]基于改進SCS流速計算公式提出了考慮降水強度的時變分布式單位線計算方法；Bunster等[12]考慮動態(tài)的上游入流貢獻,開發(fā)了時變動態(tài)上游入流貢獻的分布式單位線,提高了洪峰預報精度。然而,上述方法均假定全流域水質點存在靜態(tài)匯流路徑,即流域上均勻分布的凈雨滴通過唯一的匯流路徑向出口流動,事實上,以蓄滿產流為主的南方濕潤地區(qū),在流域未達到蓄滿狀態(tài)之前,只有蓄滿區(qū)域產生凈雨,其匯流路徑應隨流域蓄水狀態(tài)而動態(tài)變化,當流域蓄滿后,動態(tài)匯流路徑轉變?yōu)殪o態(tài)匯流路徑。

為此,本文考慮下墊面空間分布異質性和降水強度的時變特性,根據(jù)各個時段凈雨量級大小、下墊面土壤含水量狀態(tài)推求流域時變匯流路徑,嘗試根據(jù)流域產流情況將流域靜態(tài)匯流路徑轉換為動態(tài)匯流路徑,推求同時考慮降水強度和下墊面土壤含水量時空分布特征的坡面流速計算公式,提出一種考慮動態(tài)匯流路徑的時變分布式單位線。

1 研究方法

現(xiàn)有時變分布式單位線存在2個假定：全流域蓄水量在降水結束前可以達到飽和狀態(tài)以及全流域的坡面、河網匯流路徑唯一。本文從流速計算和匯流路徑兩方面提出解決思路,一方面,根據(jù)土壤蓄水狀態(tài)提取流域動態(tài)匯流路徑,另一方面,通過考慮時變土壤含水量改進現(xiàn)有坡面流速計算公式,在此基礎上,結合改進時變坡面流速計算公式和動態(tài)匯流路徑推求考慮動態(tài)匯流路徑的時變分布式單位線。

1.1 產流計算

SCS模型由美國農業(yè)部水土保持局于1954年開發(fā),特點是結構簡單、參數(shù)較少,被廣泛用于小型集水區(qū)徑流預報[13],本文采用該模型進行時段凈雨(Rt)計算。

1.2 動態(tài)匯流路徑

本文將動態(tài)匯流路徑定義為蓄滿柵格流向流域出口過程中所經過的路徑,由于降水過程的持續(xù),流域蓄水狀態(tài)實時變化導致匯流路徑動態(tài)變化,獲取任一時刻流域匯流路徑的關鍵在于確定流域當前蓄水狀態(tài)分布情況,新安江模型常用帕累托分布表征土壤蓄水能力的空間分布不均勻性,其缺陷在于無法獲得流域蓄水能力在空間上的具體分布情況[14]。研究表明,地形指數(shù)與流域蓄水能力密切相關[15],本文采用地形指數(shù)(TI)表征流域土壤蓄水能力,根據(jù)地形指數(shù)計算柵格蓄水能力[16],表達式如下：

(1)

式中：WM,j為j柵格單元的最大蓄水容量；WM,min為全流域所有柵格中的最小蓄水容量；WM,max為全流域所有柵格中的最大蓄水容量；TI,min=min{TI,j,j=1,2,…,N}，N為流域柵格單元總數(shù)；TI,max=max{TI,j,j=1,2,…,N}；n為待率定參數(shù),反映蓄水容量不均勻分布的系數(shù),該值越大,流域蓄水容量的空間分布越不均勻。

(2)

式中：Fαt為流域蓄滿部分的面積；Pi為i時段降水量；Ri為i時段凈雨量；Wj,1為j柵格初始時段蓄水量；Aj為柵格單元面積。

采用動態(tài)匯流路徑進行流域匯流計算時,由于僅蓄滿柵格發(fā)生產流,而SCS模型認為全流域均勻產流,因此,t時段凈雨量(Rt)需重新分配至蓄滿面積,根據(jù)水量平衡原理：

(3)

(4)

1.3 基于動態(tài)匯流路徑的時變分布式單位線

分布式單位線方法已較為成熟,被廣泛用于中小流域匯流計算,具體流程詳見文獻[11],本文所提方法步驟如下：

(1) 根據(jù)流域DEMs將流域柵格化,確定柵格邊長,在此基礎上,采用D8算法確定柵格水流方向,并根據(jù)柵格水流方向提取該柵格至流域出口的匯流路徑集合Rroad,j(第j個柵格至流域出口路徑上的柵格單元集合),則全流域所有柵格匯流路徑集合可表示為

Rroad={Rroad,j|j=1,2,…,N}

(5)

由于該路徑集合是靜態(tài)的,然而實際匯流情況為只有蓄滿區(qū)域有匯流路徑,為此,結合分布式蓄水容量,推求不同蓄滿比例(αt)下蓄滿柵格的匯流路徑集合(Rroad,αt),滿足：

Rroad,αt={Rroad,j|Aj∈Fαt}

(6)

式中：Rroad,αt為流域蓄滿區(qū)域的匯流路徑集合,Rroad,αt?Rroad。隨著降水過程的持續(xù),αt不斷變化,因此,Rroad,αt構成流域不同蓄水狀態(tài)下的動態(tài)匯流路徑集合。

(2) 計算每個柵格逐時段的匯流速度(Vj,αt),獲得不同流域土壤含水量狀態(tài)下的流速分布場?，F(xiàn)有分布式單位線僅考慮降水時空分布不均勻性和下墊面空間分布異質性對流域匯流速度的影響,其假定流域在降水脈沖結束之前可達到蓄滿狀態(tài),這與實際情況不符,忽略了時變土壤含水量造成下墊面狀態(tài)的改變,已有研究表明,坡面流速與植被覆蓋、雨強、坡度、土壤含水狀態(tài)有關[7-8,17]。為此,針對蓄滿單元,研究采用考慮雨強的SCS流速計算公式[11],針對處于匯流路徑的未蓄滿柵格單元,提出考慮降水強度和土壤含水量的時變流速計算公式：

(7)

式中：Vj,αt為當流域蓄滿比例為αt時j柵格單元的匯流速度；k為流速系數(shù),依據(jù)下墊面植被、土壤類型進行確定；βj為j柵格坡度；Ic為流域參考雨強；Wj,t為j柵格t時段蓄水量。

由式(7)可知,當某一柵格屬于河道單元時,流速固定為2 m/s[6]；當該柵格處于未蓄滿區(qū)域且未處于蓄滿柵格匯流路徑上時,柵格流速為0,不參與匯流；當該柵格處于蓄滿區(qū)域時,流速跟時變雨強有關；當該柵格處于未蓄滿區(qū)域,但上游存在蓄滿區(qū)域時,需要同時考慮雨強和土壤含水量對柵格流速的影響。

其中,未蓄滿柵格的土壤含水量表達式如下：

(8)

式中：Wj,t-1為j柵格t-1時段蓄水量；Qj,t為j柵格t時段的上游入流貢獻量;Pt為t時段降水量。

(9)

(10)

式中：Lj為柵格邊長。

(11)

在此基礎上,將蓄滿柵格面積轉換為流域出口流量過程,得到任意時段的分布式單位線：

(12)

式中：Qm為第m個時段的單位線流量,m3/s；H為單位凈雨,mm。

1.4 參數(shù)率定與目標函數(shù)

研究全面考慮洪水過程低流量誤差、高流量誤差、洪量誤差等多目標因素,其中,高流量誤差強調洪峰預報誤差[18]。選擇Kling-Gupta系數(shù)(EKG)[19]、納什效率系數(shù)(ENS)[20]、納什效率系數(shù)對數(shù)形式及均方根誤差與標準差的比值(RSR)4個目標參與優(yōu)化率定,優(yōu)化算法采用SCE-UA算法[21]。其中,EKG和ENS包含的平方誤差項強調了高流量的預報精度；RSR系數(shù)通過均方根誤差與標準差的比值量化了洪量誤差；ENS通過取對數(shù)壓縮變量尺度,強調了低流量預報精度。為全面考慮4個目標,建立了涵蓋上述指標的綜合目標函數(shù),依據(jù)Brunner等[22]研究成果,采用多目標加權的方式對水文模型進行優(yōu)化率定,表達式為

O=w1(1-ENS)+w2(1-EKG)+w3lg(1-ENS)+w4RSR

(13)

式中：w1、w2、w3、w4分別為4個目標權重。本文重點關注匯流精度,洪峰準確度較為重要,因此,ENS和EKG的比重相對較高,RSR和log(1-ENS)的比重相對較小。

2 例證研究

2.1 研究區(qū)域及數(shù)據(jù)

龍虎圩水和東石河是中國南方山區(qū)小河流,位于廣東省韓江流域,分別匯入程江和柚樹河,地處山區(qū)丘陵地帶,地形復雜,屬亞熱帶氣候,氣候高溫濕熱、暴雨頻繁,但時空分布不均,具有河床坡降陡、天然落差大、洪水易漲易退等特點,降水徑流存在顯著非線性關系,容易引起山洪災害。為此,研究選擇以上2個小流域為研究對象,其中,龍虎圩水河長20.8 km,流域面積130.8 km2,河流平均比降2.92‰；東石河全長23.6 km,流域面積152.4 km2,河流平均比降3.56‰。流域DEMs分辨率為30 m×30 m,在此基礎上提取流域柵格坡度、流向等信息;流域植被覆蓋數(shù)據(jù)來自清華大學植被數(shù)據(jù)庫(data.ess.tsinghua.edu.cn)；用于參數(shù)率定和檢驗的28場降水徑流數(shù)據(jù)由廣東省水文局梅州水文分局提供,其中,20場用于參數(shù)率定,8場用于結果檢驗。

2.2 參數(shù)率定結果

選取2個流域近年各10場洪水進行參數(shù)率定。模型率定分為產流參數(shù)率定、單位線參數(shù)率定、匯流參數(shù)率定3個階段。需要率定的參數(shù)有初損系數(shù)(η)、降雨前期流域特征參數(shù)(CN)、Ic、表征蓄水容量空間分布均勻性的系數(shù)(n)、WM，min及WM，max,其余參數(shù)(流域坡度、流速系數(shù))可根據(jù)流域下墊面信息確定。

(1) 根據(jù)歷史場次洪水過程推算流域凈雨量,以實際凈雨量與SCS模型預報凈雨量差值最小為依據(jù)率定SCS模型產流參數(shù)CN和η,實際凈雨量通過流域出口斷面流量過程反算得到。

(2) 采用試算法率定分布式單位線參數(shù)Ic。由于目標流域場次洪水的降水強度多集中在5～15 mm/h,根據(jù)水量平衡原理重新分配后,蓄滿區(qū)雨強范圍多為20～60 mm/h,在參數(shù)Ic率定時,先根據(jù)流域歷史平均降水強度假定幾組離散值,本文假定20 mm、40 mm、60 mm,然后在其中優(yōu)選,具體參見文獻[11]；此外,推求分布式單位線還需獲取流速系數(shù)、流域坡度等信息,其中,流速系數(shù)依據(jù)植被覆蓋情況進行確定,不同植被類型對應的流速系數(shù)見文獻[6],流域植被類型、地形坡度通過Arcgis工具提取。

(3) 根據(jù)時段凈雨及式(12)進行匯流計算,該過程需要依據(jù)流域蓄滿比例選擇單位線,流域蓄滿比例的計算關鍵在于確定流域蓄水容量空間分布,因此,采用1.4節(jié)所提方法率定蓄水容量參數(shù)n、WM,max及WM,min,目標函數(shù)權重確定參考文獻[22],w1—w4分別為0.5、0.25、0.15、0.1。根據(jù)率定后的流域CN值計算流域潛在蓄水能力(S),以該值作為流域平均蓄水容量,進一步試算參數(shù)n和WM,max,研究發(fā)現(xiàn),在流域平均蓄水容量不變的前提下,參數(shù)WM,min不敏感,變化范圍較小,本文取值5 mm,參數(shù)n、WM,max對流域蓄水容量的空間分布有顯著影響,參數(shù)n取值越大,蓄水容量空間分布越不均勻,WM,max取值也越大。本文參數(shù)率定結果如表1。

表1 龍虎圩流域和東石流域參數(shù)取值方案Table 1Parameter schemes of the Longhuwei River basin and Dongshi River basin

表2 龍虎圩流域和東石流域場次洪水率定結果及方案選擇Table 2Calibration accuracy and schemes selection results in the Longhuwei River basin and Dongshi River basin

由表2可知,龍虎圩流域和東石流域率定期洪峰相對誤差在±20%以內,峰現(xiàn)時間誤差均在±6 h之間,2個流域的ENS均值分別為0.86和0.82,場次洪水合格率超80%。

2.3 動態(tài)匯流路徑的時變分布式單位線推求

根據(jù)2.2節(jié)確定的單位線參數(shù)可以推求考慮動態(tài)匯流路徑的時變分布式單位線和現(xiàn)有時變分布式單位線,為方便工程實際應用,結合相關研究[10-11,23],本文采用忽略上游入流貢獻的式(8)進行未蓄滿柵格土壤含水量計算；同時,為減少單位線條數(shù),將降水強度和流域蓄滿比例離散化,凈雨量和流域蓄滿比例因子離散區(qū)間見表3和表4。

表3 每個時段的雨強對應的離散值Table 3Discrete value Is of excess rainfall intensity It/40 corresponding to each time interval

表4 每個時段的土壤含水量狀態(tài)對應離散土壤含水狀態(tài)Table 4Discrete value αs of soil moisture content αt corresponding to each time interval

以龍虎圩流域為例,獲得蓄水容量空間分布后,根據(jù)其蓄滿比例獲得具體蓄滿區(qū)域,不同蓄滿比例下(以αt取0.25、0.50、0.75為例),流域蓄滿區(qū)域分布如圖1所示,結合表3、表4和流域蓄滿面積可以推求考慮動態(tài)匯流路徑的時變分布式單位線。以推求蓄滿比例為0.25時的單位線為例,即圖1(a)中綠色區(qū)域所產生的單位凈雨在流域出口形成的流量過程線(圖2中綠色的單位線),相應地,可以推求不同蓄滿比例下的匯流單位線,現(xiàn)有時變分布式單位線結果如圖3所示。

圖1 龍虎圩流域不同蓄滿比例對應的產流面積空間分布Fig.1 Runoff areas corresponding to different full storage ratio in the Longhuwei River basin

由圖2可知,在同一雨強下,考慮動態(tài)匯流路徑的時變分布式單位線在流域蓄滿比例較小時單位線峰值也較小,單位線上漲歷時略有提前趨勢,但不顯著；在同一蓄滿比例下,隨著降水強度的增加,流域匯流速度增加,匯流時間減少,單位線峰現(xiàn)時間提前,總歷時減小。本文所提考慮動態(tài)匯流路徑的時變分布式單位線的物理意義為單位時段、單位凈雨均勻降落在蓄滿區(qū)域,在流域出口形成的流量過程線；圖3為現(xiàn)有時變分布式單位線,其物理意義為單位時段、單位凈雨均勻降落在全流域,在流域出口形成的流量過程線。兩者在降水及下墊面時空分布、匯流速度及匯流路徑等方面均存在顯著差異。

圖2 龍虎圩流域考慮動態(tài)匯流路徑的時變分布式單位線Fig.2 Time-varying distributed unit hydrographs considering dynamic flow routing paths for the Longhuwei River basin

圖3 龍虎圩流域現(xiàn)有時變分布式單位線Fig.3 Current time-varying distributed unit hydrographs for the Longhuwei River basin

2.4 與現(xiàn)有分布式單位線對比分析

根據(jù)SCS模型產流計算結果,分別采用本文所提單位線和現(xiàn)有分布式單位線進行模型檢驗,選取龍虎圩流域和東石流域各4場洪水進行預報,當流域前期較為干燥時選用參數(shù)方案1,否則采用參數(shù)方案2,預報結果如表5。

表5 本文方法和現(xiàn)有方法預報結果對比Table 5Comparisons of the results predicted by the proposed and current methods

由表5可知,采用考慮動態(tài)匯流路徑的時變分布式單位線進行匯流計算后,龍虎圩流域和東石流域場次洪水平均納什效率系數(shù)在0.85以上,洪峰相對誤差均值約為8%,平均峰現(xiàn)時間為1.65 h,預報精度達到工程應用乙級及以上精度,且優(yōu)于現(xiàn)有時變分布式單位線預報結果,其中龍虎圩流域預報效果優(yōu)于東石流域。以龍虎圩流域20030517號、20120527號場次洪水為例展開分析,采用本文所提動態(tài)匯流單位線進行徑流預報時,這2場洪水預報精度有較大提升,2場洪水預報徑流過程如圖4所示；龍虎圩流域20060601號、東石流域20161021號、20190609號場次洪水存在精度提升不明顯或各指標沒有同步提升現(xiàn)象,原因如下：龍虎圩流域20060601號洪水前期土壤含水量達106 mm,流域接近蓄滿狀態(tài),動態(tài)流徑的時變分布式單位線和現(xiàn)有方法推求的單位線差距較??；東石流域20161021號、20190609號洪水產流參數(shù)不具有普適性。

圖4 場次洪水預報流量與實測流量對比Fig.4 Comparisons of the predicted and observed flow for flood events 20030517 and 20120528

由于采用相同的產流參數(shù),采用2種單位線進行預報的產流量是一樣的,因此,預報洪水過程誤差主要由單位線的差異造成。由圖4可知,采用現(xiàn)有時變分布式單位線進行匯流計算時,20030517號洪水和20120527號洪水預報誤差均高于所提方法,其中,20030517號洪水洪峰相對誤差超過14%,ENS為0.88；20120527號洪水峰現(xiàn)時間誤差為3 h,ENS為0.87。采用所提單位線進行預報后,2場洪水過程的ENS分別提高至0.94和0.93。結果表明,考慮動態(tài)匯流路徑的時變分布式單位線可以更好地表征流域的匯流過程。

2.5 流域匯流時間規(guī)律解析

分布式單位線的核心為柵格匯流時間的推求,為解析動態(tài)匯流路徑對流域匯流時間的影響機制,下面深入分析流域匯流時間的動態(tài)變化規(guī)律。采用動態(tài)匯流路徑后,有一部分蓄滿柵格由于流經未蓄滿區(qū)域,導致匯流時間發(fā)生變化,為比較流域不同蓄滿比例下土壤含水量對柵格匯流速度的影響,離散凈雨強度設定為最大值Is=2.0,通過提取流域最先達到蓄滿狀態(tài)的25%流域面積(F0%—F25%)的柵格匯流時間,分析該部分面積上柵格在不同流域蓄滿比例下的匯流時間變化,并與采用現(xiàn)有時變流速計算公式的匯流時間進行比較,以龍虎圩流域為例,流域最先蓄滿的25%面積上的柵格匯流時間分布如圖5。

圖5 龍虎圩流域25%蓄滿柵格匯流時間分布Fig.5 Flow routing time distribution of the 25% full storage grid in the Longhuwei River basin

圖5展示了圖1(a)中綠色區(qū)域的柵格匯流時間,αt=1.00的橙色區(qū)域代表現(xiàn)有時變流速計算公式的結果,橫坐標代表柵格,按照蓄水容量從小到大排序,以圖5F18%處對應的柵格為例,其在蓄滿狀態(tài)為0.25、0.50、0.75、1.00時對應的匯流時間分別為9 h、7 h、6 h、5 h,可以看出隨著流域蓄滿比例的增加,蓄滿柵格的匯流時間逐漸減少,這是因為隨著降水的持續(xù),處于蓄滿路徑的未蓄滿柵格單元也逐漸趨于飽和,未蓄滿柵格單元逐漸減少,由式(7)可知,柵格蓄滿后流速達到最大值,因此,隨著流域蓄滿比例的增加,柵格匯流時間逐漸縮短,直至流域達到全蓄滿狀態(tài),此時動態(tài)匯流路徑轉化為靜態(tài)匯流路徑,流域匯流時間趨于穩(wěn)定。

圖6為龍虎圩流域不同蓄滿比例下考慮動態(tài)匯流路徑時計算的匯流時間和靜態(tài)匯流路徑計算的匯流時間箱線圖,結果表明,現(xiàn)有時變流速公式會高估坡面匯流速度,匯流時間不隨蓄滿比例而改變,考慮土壤含水量對流域匯流速度的影響后,隨著流域逐漸達到全部蓄滿狀態(tài),即αt趨于1.00時,2種方法匯流時間趨于一致,匯流時間分布更加集中；此外,隨著流域蓄滿面積的增加,土壤含水量對匯流時間的影響有減小趨勢,所提匯流速度計算方程對前期較為干旱的情況影響較大。

圖6 本文所提方法和現(xiàn)有方法不同蓄滿比例下流域匯流時間箱線圖Fig.6 Box plots of flow routing time under different storage ratio of the proposed and current methods

3 結 論

(1) 針對現(xiàn)有模型中的蓄水容量曲線只能定性描述流域蓄水容量空間不均勻分布的缺陷,引入地形指數(shù)刻畫流域柵格蓄水容量的空間不均勻分布,定量表征每一個柵格蓄水容量,結合產流計算結果,進而從柵格尺度劃分流域蓄滿柵格和未蓄滿柵格,可為柵格流速場計算提供更高精度輸入場。

(2) 提出了考慮時變降水強度和土壤含水量的匯流速度計算公式,通過劃分4種柵格狀態(tài)(蓄滿柵格、未蓄滿但處于匯流路徑柵格、未蓄滿且不處于匯流路徑柵格、河道柵格),更加精確地刻畫了流速場的空間分布規(guī)律。

(3) 建立了考慮動態(tài)匯流路徑的時變分布式單位線,針對傳統(tǒng)單位線假定全流域均勻產流的不合理假定,本文以南方山區(qū)龍虎圩流域和東石流域為研究對象,基于蓄滿區(qū)域更易產流的理論基礎,推求了流域不同蓄滿比例下的動態(tài)匯流路徑,完善了單位線推求的理論基礎。

(4) 結合本文所提改進時變流速計算公式和動態(tài)匯流路徑理論,計算了不同蓄滿比例下柵格匯流時間分布場,具體而言,通過改進流速公式提高了匯流時間計算的精度,通過劃分蓄滿區(qū)域和未蓄滿區(qū)域,將匯流方法的匯流路徑從靜態(tài)匯流路徑發(fā)展為動態(tài)匯流路徑,應用結果表明,所提方法應用精度較高,提高了分布式單位線匯流時間計算的準確性。