張 晨,于 昊,于若蘭,鄭云鶴,楊 蕊,高學平

(天津大學水利工程仿真與安全國家重點實驗室,天津 300350)

河湖水系連通工程作為解決中國水災害、水資源、水環(huán)境、水生態(tài)四水問題的有效手段,對區(qū)域水資源優(yōu)化配置、河湖生態(tài)環(huán)境修復以及區(qū)域綜合競爭力提升具有重要作用。關于河湖水系連通概念至今尚未形成較一致的認識,最初概念是河流與湖泊水系之間的連通狀況[1],而實際包涵的范圍很廣,既有河流自身的縱向連通、又有河流、河湖之間的橫向連通等,目前多將其引申為跨流域調(diào)水實現(xiàn)不同流域水系之間的相互聯(lián)系[2- 3]。

21世紀以來,中國河湖水系連通工程不斷興建,如南水北調(diào)、引江濟太、常州市主城區(qū)河網(wǎng)以及沂沭泗流域河湖連通等工程[4]。河湖水系連通工程作為一項復雜的系統(tǒng)工程[5],勢必伴隨一系列潛在風險,如調(diào)水與洪水遭遇引起的防洪風險、連通工程引起的污染物轉移或滯留風險等[6]。鑒于連通工程的復雜性、連通內(nèi)涵的多維性與連通伴生風險的隨機性,亟需構建能夠定量分析連通功能與風險評估的河湖水系連通伴生風險管控方法。分析河湖連通伴生風險的風險源,既存在氣候變化、地質災害等自然因素,也存在人工調(diào)蓄、河道利用整治等人為因素[7]。在對人為因素造成的河湖連通工程伴生風險管控中,隨著水力學模型計算方法的廣泛應用,部分學者基于河網(wǎng)水動力模型,通過調(diào)節(jié)控制性水工建筑物,改變不同引調(diào)水時間[8]、引水流量[9]、水量分配方案[10]等條件,對河網(wǎng)蓄水、輸水調(diào)度方案進行優(yōu)化,提出河湖水系連通伴生風險管控建議；也有學者通過引調(diào)水及閘壩的控制,探析調(diào)控干擾對河道水生態(tài)環(huán)境或行洪能力的影響特征[11- 12],進而實現(xiàn)風險管控下的防汛抗旱、水質提升、城市水景觀提升等多目標協(xié)同功能。

橡膠壩作為一種新型擋水建筑物被廣泛應用于河道治理工程中[13],其能夠根據(jù)需要調(diào)節(jié)壩高控制上游水位,兼具灌溉、防洪、擋潮等效益。目前有關橡膠壩調(diào)度對河道防洪能力影響研究表明,倘若洪水期調(diào)度運用不當,易造成連通河網(wǎng)下游洪峰疊加,存在一定的度汛安全隱患[14]；橡膠壩在汛期立壩運行,不能保證安全度汛,無法充分發(fā)揮工程景觀效益[14]；通過動態(tài)控制進行調(diào)度后,壩前水位及壩后流速皆會受到影響[15],可能造成洪水風險空間轉移。因此,在復雜的連通河網(wǎng)中,探究橡膠壩運行調(diào)度對河湖水系連通伴生風險的影響,并提出相應的風險管控措施,是當前面臨的一項科學問題。此外,現(xiàn)有的風險分析方法多為單一的統(tǒng)計學方法或基于水力學模型的情景分析方法,或不具備物理過程模擬,或缺乏對風險隨機性的探討,如何將風險隨機過程與物理過程相結合,進而揭示風險的發(fā)生機制,也是值得探討的科學問題。

本文結合Vine Copula函數(shù)與水動力模型模擬技術,構建在空間上考慮連通河網(wǎng)不同河流洪峰相關性、在物理過程上通過隨機水情條件考慮風險隨機性的河湖水系連通伴生風險分析模型,以沂沭河水系上游河湖水系連通工程為例,探討不同洪水重現(xiàn)期條件下,橡膠壩調(diào)度運行方式對沂沭河上游水系連通伴生洪水風險的影響,確定相應的風險安全域,更為直觀地反映風險因子作用范圍,以期突破河湖水系連通工程伴生風險管控技術難題,為流域實際管理提供指導。

1 研究區(qū)域及方法

1.1 研究區(qū)域與數(shù)據(jù)來源

沂沭河水系位于淮河流域東北部,山東境內(nèi)控制流域面積17 253 km2,如圖1所示。1951—1953年導沭整沂時期開挖了分沂入沭水道,并在后期建設不同程度擴大工程,從而實現(xiàn)沂沭河洪水東調(diào)工程計劃,減輕沂河劉道口以下河道行洪壓力,形成了沂沭河水系連通工程。該水系北起沂蒙山,東臨黃海,多年平均降雨量592.1 mm。沂沭河水系特點鮮明,汛期洪水峰高量大,源短流急,極易造成洪澇災害,調(diào)度管理反應期較短、防洪難度較大。此外,為優(yōu)化城市水景觀、增強水資源儲備,沂河與沭河分別建有10余座橡膠壩提升蓄水能力,一定程度上起到防汛抗旱的作用。

圖1 沂沭河水系上游連通狀況及水文站分布Fig.1 Connectivity and the distribution of hydrological station in the upstream of Yi & Shu River system

根據(jù)水文數(shù)據(jù)的可獲取性,分別選用水位、流速、流量作為連通水系伴生洪水風險的衡量指標。依據(jù)沂沭河水系水文年鑒,選取豐水年(1957年)、平水年(1985年)、枯水年(1981年)3個典型水文年的水位、流速、流量數(shù)據(jù),并通過內(nèi)插得到汛期2h間隔數(shù)據(jù)集,用于模型驗證。

1.2 河湖水系連通伴生風險分析模型

1.2.1 風險分析模型控制方程

隨機水情條件下河湖水系連通伴生風險分析模型集合Vine Copula模型與水動力模型模擬技術,采用大數(shù)據(jù)分析研究思路,以計算風險發(fā)生概率為目標,在水動力模型的基礎上,通過Copula函數(shù)構建連通河網(wǎng)不同河流洪峰的相關結構,并采用蒙特卡羅法進行隨機模擬,產(chǎn)生隨機水情條件作為水動力模型的邊界進行模擬,根據(jù)設定的風險限值,利用水位、流速、流量模擬結果計算求得河道各斷面的風險概率,進而完成河湖水系連通伴生風險分析。

二維水動力數(shù)學模型聯(lián)立連續(xù)性方程和運動方程,求解水流運動[16]。由于風險過程是一項隨機發(fā)生的可能性事件,因此,求在控制方程中增加隨機源匯項(ω),其物理意義為水系連通中引起潛在洪水風險的驅動力。改進后的連續(xù)性方程與運動方程為:

(1)

(2)

式中:u、v為沿x、y方向的流速,m/s；ω為隨機源匯項,s-1,ω=f(x,y,t)；g為重力加速度,m/s2；H為位置水頭,m；h為水深,m；D為量綱一的經(jīng)驗系數(shù),與河床、邊坡形狀相關；R為水力半徑,m；C為謝才系數(shù),m1/2/s。式(2)等式左側為加速度項,右側從左至右依次表示重力產(chǎn)生的壓力梯度、黏滯力項、河床底部摩擦及隨機源匯項。

為量化分析風險隨機源匯項的不確定性,采用Copula模型構建不同河流洪峰變量的聯(lián)合分布,并進行蒙特卡羅隨機模擬,產(chǎn)生隨機流量邊界作為源匯項。C- Vine Copulan維聯(lián)合密度函數(shù)表達式[6]:

(3)

式中:x1,x2,…,xn表示關鍵風險因子；fi(xi)為邊緣密度函數(shù)；ci,i+j|1:(i-1)為n維的Copula密度函數(shù)；F(xn|x1)為條件密度函數(shù)；θ為Copula函數(shù)中n維變量的相關關系的指標。

對式(3)兩邊進行積分得到其n維累積分布函數(shù):

(4)

累積分布函數(shù)的逆函數(shù)C-1(x1,x2,…,xn)代表隨機生成的關鍵風險因子x1,x2,…,xn序列,由此可推導出運動方程(2)中隨機源匯項(ω)的表達式。

根據(jù)質量守恒定律,采用歐拉法描述流體微團的質量變化,ω可表示為

(5)

式中:qx、qy分別為x、y方向上的單寬流量,m2/s。由于流速、水位都可以表示為流量的函數(shù),故本文風險分析模型例中以流量表述隨機源匯項。將累積分布函數(shù)的逆函數(shù)C-1(x1,x2,…,xn)代入式(5),得:

(6)

1.2.2 風險分析模型建立及驗證

為更加準確和直觀地模擬水位水量實際聯(lián)動效果,對研究區(qū)域河流渠道建立一維數(shù)學模型,對水庫湖泊建立二維數(shù)學模型,采用一維與二維數(shù)學模型耦合的方式進行模擬。沂河、沭河、分沂入沭3個一維河道每隔1 km劃分1個斷面,河段長度總計374 km,共劃分375個斷面；新沭河山東段與新沭河江蘇段2個一維河道每隔1 km左右劃分1個斷面,河段長度總計78 km,共劃分為70個斷面。對于石梁河水庫采用四邊形結構網(wǎng)格建立二維模型,網(wǎng)格尺寸為100 m×100 m,共劃分6 058個網(wǎng)格。

選取沂河葛溝站、沭河莒縣站以及新沭河大官莊站和大興鎮(zhèn)站對模型進行驗證。由于風險分析模型的隨機流量條件范圍較大,故選用豐水年、枯水年洪水過程作為率定期,平水年作為校核期。率定結果為,沂河與沭河河床及邊坡糙率0.024～0.027,新沭河糙率0.030。部分代表性誤差分析結果如圖2和圖3所示,率定期和校核期水位誤差統(tǒng)計值亦列于圖中。結果表明,平均絕對誤差(EMA)和均方根誤差(ERMS)低于0.3 m、納什效率系數(shù)(ENS)大于0.85、相對誤差(ER)小于12%,證明模型可靠。

圖2 模型率定期模擬水位與實測水位對比Fig.2 Comparison of simulated and observed water level during the calibration period

圖3 模型校核期模擬水位與實測水位對比Fig.3 Comparison of simulated and observed water level during the verification period

2 橡膠壩對沂沭河水系上游連通伴生洪水風險影響分析

2.1 風險源與模擬方案

沂河、沭河河道上分別修建了10余座橡膠壩,其作為1種阻礙河湖水系連通結構的攔蓄工程,改變了河流原有的結構連通；另一方面,由于沂河與沭河具有相似的水文特性,洪水遭遇時常發(fā)生。本文結合這2種風險源,在沂河、沭河最大50年一遇與100年一遇洪水重現(xiàn)期條件下,從流量、水位、流速三方面探究橡膠壩對沂沭河水系上游連通伴生洪水風險的影響,并通過調(diào)節(jié)橡膠壩運行高度提出其伴生綜合風險管控建議。

構建聯(lián)合分布前,需對邊緣分布進行擬合,本文對沂河與沭河上游邊界斷面洪峰流量進行分布擬合。采用極大似然法估參,并通過赤池信息準則(AIC準則)進行擬合優(yōu)度檢驗,AIC值越小代表擬合分布模型適合度越高。選取水文中常用的正態(tài)分布、logistic分布、對數(shù)正態(tài)分布、威布爾分布分別對沂河與沭河洪峰流量原始序列進行邊緣分布擬合,并將各分布對應的AIC值列于表1。

表1 各隨機變量不同分布類型對應AIC值

由表1結果可知,沂河洪峰流量最優(yōu)擬合分布為正態(tài)分布,沭河洪峰流量最優(yōu)擬合分布為威布爾分布。對新序列進行概率積分轉換使其滿足在[0,1]區(qū)間上均勻分布。根據(jù)AIC擬合優(yōu)度檢驗準則,自Gaussian Copula、Student′t Copula、Clayton Copula、Gumbel- hougaard Copula、Frank Copula、Joe Copula、BB族Copula中選取兩變量間最適宜的Pair- Copula類型,并采用極大似然法進行參數(shù)估計。結果可得兩變量間最適宜的Copula函數(shù)類型為Rotated Clayton Copula (90°)。

根據(jù)構建完成的Copula模型生成變量隨機序列,在沂河和沭河洪峰流量區(qū)間內(nèi)隨機組合出水情條件(按正態(tài)分布依據(jù)典型豐水年洪水過程生成3個月的洪水過程線)作為風險源邊界條件。根據(jù)典型風險源連通性特征,將風險模擬方案及邊界條件列于表2。流量區(qū)間下限依據(jù)歷史資料記載的洪災較小流量,如沂河1957年的洪峰流量8 120 m3/s；上限依據(jù)沂沭泗水利管理局2012年頒布的《沂沭泗河洪水調(diào)度方案》中50年一遇標準或100年一遇標準確定,每種連通方案下隨機得到200組或400組水情條件。根據(jù)各連通方案下模擬結果,提取河道內(nèi)各斷面每組最大流量、最高水位以及最大流速計算結果,采用累積分布函數(shù)(CDF)分別對流量、水位、流速結果進行處理以獲得其概率分布,再分別選取河道警戒水位、不沖流速、校核流量作為判別發(fā)生風險的限值,計算出流量、水位、流速具體的風險概率(P),如式(7)所示:

表2 模擬方案

(7)

式中:GI≤IL為風險評估指標計算結果(I)小于等于風險限值IL的模擬組數(shù)；GT為總模擬組數(shù)。

本文定義連通伴生洪水風險概率介于0～0.2之間為極低風險,0.2～0.4為低風險,0.4～0.7為中風險,0.7～1.0為高風險。

2.2 橡膠壩對沂沭河水系連通伴生洪水風險的影響

沂河、沭河河道橡膠壩壩前蓄水量約814.9萬m3,下泄所蓄水量的塌壩時間平均約6.1 h[17],在洪水來臨前預泄能夠避免其與洪峰遭遇疊加。但是汛前塌壩共同下泄蓄水,將人為形成洪峰疊加,可能增加下游防洪風險。為此,本文對該情景進行模擬,其風險概率分布如圖4(a)所示?？梢钥闯?由于洪峰的疊加,所受影響最大的是沭河中下游河段,水位與流速風險均為高風險,風險概率在0.8以上；沂河并未表現(xiàn)出相應的特點,但河道整體處于中高風險區(qū)間,臨沂市水位與大官莊水利樞紐流速是需要重點關注風險區(qū)域(概率達到1)；沂河、沭河及分沂入沭流量風險自上游至下游從0.1～0.3逐漸增加,皆處于低風險以下。由此可知,若在汛前完全塌壩下泄,沂沭河上游,尤其沭河中下游河段將產(chǎn)生極高風險。

進一步分析不塌壩對水系行洪能力的影響,對表2中方案1與方案2進行模擬計算。依據(jù)圖4(b)可知:在50年一遇的洪水重現(xiàn)期條件下,由于橡膠壩對水體的攔蓄作用,水位大幅上漲溢流使橡膠壩對泄流風險幾乎沒有影響；水位方面,若汛前未將橡膠壩前攔蓄水量提前預泄,壩前水位風險顯著提升至高風險,風險概率達到最大值1,分沂入沭及新沭河所受影響不大,處于極低風險；流速方面,沂河、沭河流速風險因橡膠壩阻水作用位于極低風險,概率在0～0.2之間,其余河道無風險。而當洪水重現(xiàn)期升至100年一遇后(圖4(c)),水系內(nèi)各河道流量、水位、流速風險整體均大幅提升至中高風險(0.5以上),橡膠壩對河湖水系連通的阻礙作用更為顯著,壩址處水位風險概率為1,流速風險概率為0。

圖4 沂沭河水系上游各情景下風險概率分布Fig.4 Risk probability distribution of each scenario in the Yi & Shu River system

綜合以上分析,當沂沭河水系面臨較大洪水時,無法預泄橡膠壩前所蓄水量將提升河道內(nèi)水位風險,而汛前塌壩下泄亦會造成人為洪峰疊加產(chǎn)生極高風險。

3 橡膠壩運行高度管控伴生洪水風險及安全域

3.1 橡膠壩運行高度對伴生洪水風險的影響規(guī)律

為有效將橡膠壩對沂沭河水系河湖連通伴生洪水風險影響維持在安全可控的范圍內(nèi),本文針對性地提出調(diào)整橡膠壩運行高度的管控措施,其中具體工況與模擬邊界條件詳見表2中方案3和方案4。使用風險齒輪圖探討橡膠壩運行高度為25%、50%、75%設計壩高3種情況的風險變化規(guī)律,圖內(nèi)涵蓋了不同壩高下水系內(nèi)河道所有斷面的風險概率值。流量在模型中作為邊界條件,故此處不討論流量風險。

如圖5所示,在50年一遇洪水條件下,沂河、沭河水位風險概率與橡膠壩高度正相關,沭河流速風險概率與橡膠壩高度負相關,每當橡膠壩高度升高25%的設計壩高時,壩前水位風險提高約70%,沭河流速風險降低約50%,沂河由于河道坡度相對較小,流速始終在極低風險區(qū)間內(nèi)。由于橡膠壩對水體的攔蓄作用,水位大幅上漲溢流使得橡膠壩對下游流量風險幾乎沒有影響,因此,壩高變化對于下游河道影響甚微,3種運行高度情況下風險概率變化幅度較小。需要重點關注的是河與河、河與湖連接處,可通過加強河道堤防防洪標準并提升其附近河段抗沖蝕能力以降低風險。

圖5 50年一遇重現(xiàn)期不同橡膠壩運行高度下沂沭河水系洪水風險齒輪圖Fig.5 Flood risk gear under different rubber dam operating heights during the 50- year flood return period in the Yi & Shu River system

在100年一遇洪水條件下,沂河水位風險變化和沭河水位、流速風險變化,與重現(xiàn)期50年一遇時在風險變化率上基本呈現(xiàn)出相同的規(guī)律,且在25%設計壩高運行高度情況下壩址處水位風險概率已經(jīng)達到最大值。而沂河的流速風險因洪量增加,其與橡膠壩運行高度負相關的趨勢更加顯著,河道上除橡膠壩壩址處流速風險外,風險概率均在0.5以上(中高風險)。由于篇幅所限,圖略。

3.2 橡膠壩運行調(diào)度管控伴生綜合洪水風險建議

前述研究發(fā)現(xiàn),在沂河、沭河修建橡膠壩的連通狀況下,2條河道水位風險和流速風險呈現(xiàn)的變化規(guī)律相反。在50年一遇洪水條件下,沂河流速風險始終位于極低風險區(qū)間(圖5),由于沭河水力坡度變化較大,相比沂河潛在流速風險高,河段部分區(qū)域水位與流速皆達到中風險。因此,在調(diào)節(jié)橡膠壩至不同運行高度情況下,針對沂河僅討論水位風險概率規(guī)避效果,再以沭河水深、流速為例,討論橡膠壩運行高度對伴生綜合洪水風險影響。

沂河水位風險概率隨水深增加而增加,且壩高越高時其正相關關系越顯著,如圖6所示,圖中各點為所研究河段不同斷面的結果。在75%壩高情況下,部分壩前河段水深較大,水位風險概率為1；當壩高降至25%與50%時,高風險河段大量減少,水位風險概率均值由管控前0.65降為0.3,即低風險。在圖6中,利用中- 低風險限值0.4與沂河警戒水深10 m圍取沂河斷面安全范圍,25%壩高時落在安全范圍內(nèi)的河道斷面相對較多,也存在極少部分斷面水位風險概率為1;壩高升至50%時,這部分河段因為更多水體被上游橡膠壩攔蓄,其水位風險概率下降,但更多河段升至中高風險。沂河水務部門提出了“梯級橡膠壩限蓄,洪水分級,峰小量大,泄量相當,跟蹤預報”等6個調(diào)度原則[18],有學者建議沂河汛期最大蓄水壩高為60%設計壩高[19]。本研究進一步根據(jù)風險分析模型計算,建議沂河橡膠壩在汛前調(diào)節(jié)至小于50%壩高,以使沂河約78%的斷面風險能夠降至低風險區(qū)間。

圖6 沂河橡膠壩運行壩高管控措施下水位風險概率隨水深變化關系Fig.6 Relationship between the water level risk probability and water depth under the control measures of rubber dam operating heights in Yi River

沭河水深、流速與綜合風險概率關系如圖7所示,根據(jù)計算結果將河道斷面用邊界線確定在有效區(qū)域范圍內(nèi),定義水位風險概率或流速風險概率其中之一大于0.4,則此斷面為風險斷面；二者皆小于0.4為安全斷面。可以看出,在安全斷面與風險斷面之間有相對明顯的界限(水位風險概率為0.4或流速風險概率為0.4),利用此界限與邊界線確定出沭河安全域。對比圖7(a)—圖7(d),當運行高度為其設計壩高的50%時,沭河安全域最大,能夠承擔的水深與流速較大(均未超過沭河最大警戒水深與不沖流速)。同時,該運行壩高下,安全斷面亦最多,其余運行高度風險斷面均超過一半以上,且能夠承擔的最大流速、水深相對較小。因此,建議沭河橡膠壩在汛前調(diào)節(jié)50%設計壩高運行,并控制河道水深小于12 m、流速小于2.23 m/s,從而有更多的安全斷面；若以25%壩高運行,需滿足河道水深小于12 m,同時流速小于1.94 m/s；75%與100%運行高度情況下,則要保證河道水深與流速分別小于11 m和1.81 m/s。

圖7 不同壩高情況下沭河河道斷面安全域Fig.7 Security domains of river sections under different rubber dam heights in Shu River

洪水重現(xiàn)期升至100年一遇后,沂河、沭河各方面整體均提升至中高風險,兩河道水位、流速風險均隨橡膠壩壩高增加呈現(xiàn)相反的變化規(guī)律,但通過風險概率計算結果得知,無論在何種運行高度下,所有河段水位與流速風險概率均在0.5以上,僅壩址處因橡膠壩阻水作用流速風險概率為0。區(qū)別于50年一遇情況,100年一遇條件橡膠壩在較低運行高度時(25%以下),壩址處水位風險概率已大幅增長至1,上文中所定義安全域無法適用于100年一遇條件。筆者又對100年一遇未修建橡膠壩情景進行計算,發(fā)現(xiàn)此情況下水系內(nèi)各河道流量、水位、流速風險概率全部為0.5左右,處于中低風險,并且無高風險河段。因此,建議在面對100年一遇洪水時,需要將橡膠壩調(diào)至25%設計壩高以下,或汛前盡早緩慢塌壩,才能降低沂沭河水系防洪壓力。

橡膠壩的風險管控方式與湖庫有較大差別,湖泊水庫調(diào)度對風險管控有一明顯的可容忍區(qū)域[20],而修建橡膠壩的河道蓄水能力相對較小,其風險彈性有限,在調(diào)控過程中易造成洪水風險空間轉移,因此,需通過調(diào)控橡膠壩壩高嚴格控制河道水位與流速,使得洪水風險達到最小。

4 結 論

本文創(chuàng)建了集合Vine Copula函數(shù)與一維二維水動力模型的隨機水情條件下河湖水系連通伴生風險分析模型,形成了不同洪水重現(xiàn)期條件下典型區(qū)域橡膠壩對河湖連通伴生洪水風險影響的分析和管控技術。主要結論如下:

(1) 利用該風險分析模型對沂沭河水系河湖連通工程人為洪峰疊加、50年一遇及100年一遇修建橡膠壩3種情景洪水風險進行了分析。若在汛前完全塌壩下泄橡膠壩所蓄水量形成人為洪峰疊加,對沭河中下游河段產(chǎn)生極高風險；若將橡膠壩維持在最大高度運行,沂河沭河橡膠壩壩址處水位風險大幅提升。

(2) 針對典型區(qū)域橡膠壩對連通伴生洪水風險的影響提出了管控建議。重現(xiàn)期50年一遇時,建議汛前調(diào)節(jié)沂河沭河橡膠壩運行高度50%以下,并控制水深與流速在一定限值(安全域)內(nèi)；而在100年一遇情況下,建議將橡膠壩調(diào)至25%設計壩高以下,或者汛前盡早緩慢塌壩下泄蓄水,以減輕防洪壓力。