基于Erlang分布蓄水容量曲線的流域產(chǎn)流模型

閆寶偉,李正坤,段美壯,江慧寧,劉 昱

(1. 華中科技大學(xué)土木與水利工程學(xué)院,湖北武漢 430074；2. 中國(guó)電建集團(tuán)華東勘測(cè)設(shè)計(jì)研究院有限公司,浙江杭州 311122)

根據(jù)蓄滿產(chǎn)流原理,降水首先補(bǔ)充土壤缺水量,在流域土壤達(dá)到田間持水量以后,流域產(chǎn)流。蓄水容量曲線用于描述流域土壤缺水量空間分布不均勻的問(wèn)題,曲線的形狀決定著流域產(chǎn)流面積和產(chǎn)流量的多少,對(duì)于徑流的準(zhǔn)確模擬具有重要意義。流域蓄水容量受流域下墊面地形和土壤等多種因素影響,通常難以直接測(cè)量,蓄水容量的空間分布更難以確定,通常是經(jīng)驗(yàn)假定。在新安江模型中廣泛采用B次拋物線表示蓄水容量曲線,至于為何采用該類線型,至今尚無(wú)合理解釋。為能反映雨季和枯季蓄水容量曲線的靈活性,周買春和Jayawardena[1]提出了雙拋物線型蓄水容量曲線,以反映不同土壤水分并存狀態(tài)下的蓄水容量。流域蓄水容量曲線同時(shí)也是對(duì)流域地形地貌特征的一種反映,Gao等[2]通過(guò)對(duì)美國(guó)404個(gè)流域的洪水模擬發(fā)現(xiàn),流域地形地貌對(duì)拋物線型蓄水容量曲線的參數(shù)B有較大影響。TOPMODEL中的地形指數(shù)同樣反映了流域的缺水情況,其累計(jì)頻率曲線與蓄水容量曲線實(shí)質(zhì)是一致的,都是對(duì)流域缺水量空間不均勻性的一種描述,郭方等[3]提出了利用TOPMODEL地形指數(shù)的分布推求流域蓄水容量曲線的思路。石朋等[4]通過(guò)對(duì)流域內(nèi)柵格地形指數(shù)的統(tǒng)計(jì)分析,采用對(duì)數(shù)韋布爾分布曲線描述流域蓄水容量空間分布的不均勻。孔凡哲和宋曉猛[5]根據(jù)各子流域的平均坡度計(jì)算蓄水容量的空間分布。向小華等[6]則進(jìn)一步提出了融合土壤和地形的柵格蓄水容量計(jì)算方法,為研究流域蓄水容量的空間分布提供了一種有效工具。Wang[7]針對(duì)SCS模型中的土濕比函數(shù)沒(méi)有考慮土壤蓄水容量空間分布不均勻性的問(wèn)題,提出了一種新的分布函數(shù)來(lái)描述土壤蓄水容量的空間分布。Gao等[8]根據(jù)植被根系深度隨距離最近河道高程的空間分布規(guī)律,建立了以地形為基礎(chǔ)的土壤含水量和變?cè)串a(chǎn)流面積間的非線性關(guān)系,提出了一種新的產(chǎn)流模型,為無(wú)資料地區(qū)的水文模擬提供了一種可行的工具。

本文嘗試通過(guò)柵格蓄水容量的計(jì)算,分析多個(gè)典型流域蓄水容量的空間分布型式,選用Erlang分布曲線進(jìn)行流域產(chǎn)流的推導(dǎo),提出基于Erlang分布蓄水容量曲線的流域產(chǎn)流模型,一方面可以更加真實(shí)地反映流域的產(chǎn)流過(guò)程,另一方面,也為無(wú)資料地區(qū)的產(chǎn)流計(jì)算提供一種可行途徑。

1 流域蓄水容量的空間分布

流域蓄水容量的實(shí)質(zhì)是流域下墊面的缺水量,與流域的地形地貌有較大關(guān)系,對(duì)流域的產(chǎn)流有直接影響,由于蓄水容量難以直接測(cè)定,其空間分布更是難以直接確定?？紤]到新安江模型中的蓄水容量與TOPMODEL 中缺水量的概念相對(duì)應(yīng),而TOPMODEL可以根據(jù)地形指數(shù)和流域平均缺水量計(jì)算出各點(diǎn)的缺水量,進(jìn)而可以分析流域缺水量或蓄水容量的空間分布。蓄水容量除了與流域地形有關(guān)之外,還受到土壤類型的影響,如黏性土壤的土壤顆粒間孔隙較小,相應(yīng)的土壤缺水量也較小；而砂質(zhì)土壤的土壤顆粒間孔隙較大,相應(yīng)的土壤缺水量也較大。向小華等[6]以van Genuchten模型表述的土壤水分特征曲線為基礎(chǔ),推導(dǎo)出了融合地形和土壤特征的流域柵格單元蓄水容量的計(jì)算公式:

(1)

式中:Wi為某點(diǎn)的蓄水容量,mm；Di=Szm[max(tpi)-tpi],為某點(diǎn)的地下水埋深,mm；Szm為流域內(nèi)非飽和區(qū)的最大蓄水深,mm；tpi為地形指數(shù)；θr、θf(wàn)和θs分別為土壤凋萎含水量、田間持水量和飽和含水量；ψc為土壤水分達(dá)到田間持水量時(shí)毛管上升高度,mm；β、m為土壤水參數(shù)。各種土壤水分特征參數(shù)的取值見(jiàn)文獻(xiàn)[9],Szm則由流域平均蓄水容量試算而得,流域平均蓄水容量可通過(guò)參數(shù)率定得出,無(wú)資料地區(qū)可采用Gao等[10]和Wang-Erlandsson等[11]提出的MCT方法進(jìn)行估算。

按式(1)計(jì)算表1所示典型流域的柵格蓄水容量,得到各典型流域蓄水容量的空間分布,如圖1所示。

表1 所選典型流域?qū)傩员?/p>

圖1 典型流域蓄水容量的空間分布Fig.1 Spatial distribution of watershed storage capacity in the selected basins

2 流域蓄水容量曲線線型分析

Erlang分布是一種連續(xù)型概率分布,該分布與指數(shù)分布一樣多用來(lái)表示獨(dú)立隨機(jī)事件發(fā)生的時(shí)間間隔,常應(yīng)用于保險(xiǎn)業(yè)以及排隊(duì)論中[12]。Erlang分布作為伽馬分布的特例,計(jì)算和應(yīng)用較為便捷,且對(duì)該分布函數(shù)的積分易得出顯式表達(dá)式,可以簡(jiǎn)化后續(xù)的產(chǎn)流計(jì)算,由此選用Erlang分布作為新的蓄水容量曲線線型,其概率密度函數(shù)為

(2)

式中:Γ(·)表示伽馬函數(shù)；n、λ為參數(shù);n值為整數(shù)。若式中參數(shù)n取值為實(shí)數(shù),則式(2)即為伽馬分布的密度函數(shù)。Erlang分布是一種伽馬分布,表示獨(dú)立隨機(jī)事件發(fā)生的時(shí)間間隔,水文中廣泛應(yīng)用的Nash瞬時(shí)單位線就是伽馬分布函數(shù)。對(duì)上式積分可得其分布函數(shù):

(3)

式中:En(x)為Erlang分布函數(shù)。

根據(jù)流域蓄水容量的空間分布,進(jìn)一步統(tǒng)計(jì)蓄水容量不同量級(jí)所占的面積比例,繪制流域蓄水容量的面積分布圖,并采用正態(tài)分布、對(duì)數(shù)正態(tài)分布、指數(shù)分布、Erlang分布等多種不同型式的分布函數(shù)以及常用的拋物線型曲線進(jìn)行擬合,發(fā)現(xiàn)Erlang分布擬合效果更優(yōu),如圖2所示,圖中α表示產(chǎn)流面積比。

圖2 典型流域蓄水容量分布的擬合圖Fig.2 Fitness of distribution of the storage capacity in the selected basins

3 基于Erlang分布蓄水容量曲線的產(chǎn)流計(jì)算

參考拋物線型蓄水容量曲線定義,則基于Erlang分布的蓄水容量曲線可以表示為

(4)

式中:W′為流域中某點(diǎn)的蓄水容量；Wmm為流域中最大的點(diǎn)蓄水容量。

采用Erlang分布蓄水容量曲線計(jì)算流域產(chǎn)流,如圖3所示,則由蓄滿產(chǎn)流原理可知,流域平均蓄水容量

(5)

上式的推導(dǎo)用到了Erlang分布的如下性質(zhì)[13]:

(6)

同理可得流域初始蓄水容量

(7)

圖3 流域蓄水容量分布曲線與降雨徑流關(guān)系Fig.3 Watershed storage capacity distribution curve and its rainfall runoff relationship

式中:A為初始蓄水容量W0對(duì)應(yīng)的縱坐標(biāo)。為計(jì)算產(chǎn)流量,需首先確定A的值,在已知W0的情況下,A可以通過(guò)求解式(7)而得,而該式為A的隱式方程,很難直接求解,此處采用牛頓迭代法進(jìn)行求解。令

(8)

進(jìn)一步對(duì)上式求導(dǎo)可得

(9)

根據(jù)牛頓迭代公式可知,經(jīng)過(guò)i次迭代后,得

(10)

當(dāng)相鄰兩次迭代值之差小于允許誤差時(shí),便可求得A的數(shù)值解。根據(jù)蓄滿產(chǎn)流原理,有效降雨P(guān)e產(chǎn)生的徑流為

(11)

綜上,若Pe+A

(12)

若Pe+A≥Wmm,即全流域產(chǎn)流時(shí):

R=Pe-Wm+W0

(13)

在產(chǎn)流計(jì)算中,只要給定流域平均蓄水容量Wm和Erlang分布的參數(shù)n、λ,即可實(shí)現(xiàn)產(chǎn)流計(jì)算。在無(wú)資料地區(qū),也可根據(jù)式(1)計(jì)算出的柵格蓄水容量,統(tǒng)計(jì)出流域蓄水容量的分布曲線,采用極大似然估計(jì)等參數(shù)估計(jì)方法估算Erlang分布的參數(shù),從而可以實(shí)現(xiàn)無(wú)資料地區(qū)的產(chǎn)流計(jì)算。在資料豐富的地區(qū),可以根據(jù)降雨徑流資料直接進(jìn)行參數(shù)率定。

4 模型應(yīng)用

將新安江模型的產(chǎn)流模塊替換成上述基于Erlang分布蓄水容量曲線的流域產(chǎn)流模型,匯流模塊仍然采用原算法,其中,壤中流和地下徑流匯流采用線性水庫(kù)法,河網(wǎng)匯流采用滯后演算法,由此得到改進(jìn)的新安江模型。將改進(jìn)后的模型分別應(yīng)用于松柏、恒口、欒川和青獅潭4個(gè)典型流域的洪水模擬,以分析模型的應(yīng)用效果。

選取松柏流域1983—2012年、恒口流域1980—2017年、欒川流域1998—2010年和青獅潭流域2002—2015年的場(chǎng)次洪水用于模型的率定和檢驗(yàn)。其中,模型的參數(shù)采用SCE-UA算法進(jìn)行優(yōu)化率定,計(jì)算時(shí)段長(zhǎng)均為1 h,模型參數(shù)模擬結(jié)果如表2所示。

表2 模型參數(shù)模擬結(jié)果

選取洪峰相對(duì)誤差、峰現(xiàn)時(shí)間誤差和洪水總量相對(duì)誤差以及確定性系數(shù)作為模型精度評(píng)價(jià)指標(biāo),分別計(jì)算每場(chǎng)洪水的評(píng)價(jià)指標(biāo)值,表3給出了各個(gè)流域這些評(píng)價(jià)指標(biāo)的平均值以及與新安江模型的結(jié)果對(duì)比?？梢钥闯?松柏流域和恒口流域的模擬精度與新安江模型相當(dāng),欒川流域和青獅潭流域的模擬精度有一定提高,尤其是欒川流域的提高較為明顯,率定期確定性系數(shù)由0.73提高到0.80,其他指標(biāo)也都有不同程度的提高。圖4進(jìn)一步展示了欒川流域檢驗(yàn)期3場(chǎng)洪水的模擬結(jié)果,可以發(fā)現(xiàn),改進(jìn)后的模型無(wú)論是洪峰還是洪水過(guò)程的模擬都取得了較高的精度,相比新安江模型也都有提高。

表3 洪水模擬結(jié)果平均值統(tǒng)計(jì)及對(duì)比

圖4 欒川流域檢驗(yàn)期洪水過(guò)程模擬結(jié)果Fig.4Simulated flood hydrographs in the validation period in Luanchuan station

5 結(jié) 論

蓄水容量曲線的形狀決定著流域產(chǎn)流面積和產(chǎn)流量的多少,對(duì)于徑流的模擬有直接影響。通過(guò)對(duì)多個(gè)典型流域蓄水容量空間分布的分析,選用Erlang分布曲線進(jìn)行流域蓄水容量曲線的擬合,提出了基于Erlang分布蓄水容量曲線的流域產(chǎn)流模型,主要結(jié)論如下:

(1)基于Erlang分布推導(dǎo)了一種新的流域產(chǎn)流模型,相比新安江模型增加了1個(gè)產(chǎn)流參數(shù),模型的適應(yīng)性增強(qiáng),模擬結(jié)果更接近實(shí)際的產(chǎn)流過(guò)程,一定程度上提高了模型的模擬精度。由于增加了模型參數(shù),可能會(huì)增加其他參數(shù)如Wm的不確定性,需進(jìn)一步分析。

(2)基于Erlang分布蓄水容量曲線產(chǎn)流模型的參數(shù),可以借助流域的地形和土壤類型數(shù)據(jù)進(jìn)行估算,為無(wú)資料地區(qū)的產(chǎn)流計(jì)算提供了一種可行的途徑。