陳 凱，彭 楊，吳志毅

(華北電力大學(xué) 可再生能源學(xué)院，北京102206)

0 引 言

水庫防洪調(diào)度是減少洪災(zāi)損失的重要非工程措施。水庫防洪調(diào)度需要考慮水庫自身安全、上下游防洪目標(biāo)和眾多約束條件，是一個典型的多階段、多約束的多目標(biāo)決策問題。受流域降雨、洪水、工程狀態(tài)等客觀不確定性以及各種人為調(diào)度管理等主觀因素的影響，有可能導(dǎo)致在防洪調(diào)度實際過程中，難以避免地出現(xiàn)實際調(diào)度結(jié)果不能達(dá)到預(yù)期調(diào)度目標(biāo)，不僅影響水庫汛期的安全穩(wěn)定運(yùn)行和水庫綜合效益的發(fā)揮，也會給調(diào)度決策帶來一定的風(fēng)險。因此，有必要進(jìn)行水庫防洪調(diào)度多目標(biāo)風(fēng)險分析。

自20世紀(jì)80年代初以來，我國學(xué)者已開始水庫調(diào)度風(fēng)險問題研究，提出了概率法[1]、一次二階矩法[2]、事故樹法[3]、蒙特卡洛隨機(jī)模擬法[4]等解決方法。對于入庫徑流和洪水的不確定性帶來的風(fēng)險，常采用蒙特卡洛隨機(jī)模擬法進(jìn)行估算[4]。在水庫防洪調(diào)度風(fēng)險研究方面，李英海等[5]通過對三峽水庫入庫洪水預(yù)報誤差進(jìn)行隨機(jī)模擬，并利用灰色關(guān)聯(lián)法對防洪調(diào)度方案進(jìn)行了風(fēng)險決策。覃暉等[6]采用蒙特卡洛法對三峽水庫入庫徑流進(jìn)行隨機(jī)模擬，并通過調(diào)度確定了防洪風(fēng)險，并應(yīng)用基于相對優(yōu)勢度的風(fēng)險型多屬性決策方法對各調(diào)度方案進(jìn)行了評價。王麗萍等[7]通過對預(yù)報誤差進(jìn)行隨機(jī)模擬來修正溪洛渡-三峽梯級水庫的預(yù)報洪水，并通過建立的梯級防洪調(diào)度風(fēng)險估計模型對梯級水庫的防洪風(fēng)險進(jìn)行了分析。周研來等[8]采用Copula-MC法對萬家寨水庫上游區(qū)間洪水進(jìn)行了模擬，并與干流洪水疊加組成了萬家寨水庫的入庫洪水，通過對水庫汛限水位動態(tài)控制分析了各調(diào)洪方案的風(fēng)險。而當(dāng)水庫下游與防洪控制區(qū)之間有支流匯入時，支流來水的不確定性也會對下游防洪控制區(qū)帶來一定的風(fēng)險。閆寶偉等[9,10]考慮干支流來水不確定性，利用Copula函數(shù)建立聯(lián)合分布，分析了長江與清江洪水發(fā)生時間與發(fā)生量級的遭遇風(fēng)險。劉章軍等[11]利用Copula-MC法通過對干支流洪水進(jìn)行隨機(jī)模擬，通過水庫調(diào)洪函數(shù)計算了隔河巖水庫下游高壩洲斷面的洪水概率分布。上述學(xué)者雖考慮了干支流洪水的不確定性，但卻沒有結(jié)合水庫調(diào)度對風(fēng)險進(jìn)行進(jìn)一步分析。

本文考慮干支流洪水的不確定性對水庫防洪調(diào)度的影響，提出了一種基于Copula-MC法的水庫防洪調(diào)度多目標(biāo)風(fēng)險分析方法。將該方法應(yīng)用于向家壩水庫防洪調(diào)度風(fēng)險分析中，給出兼顧上下游防洪目標(biāo)要求的最優(yōu)調(diào)度方案，為水庫汛期防洪調(diào)度決策提供參考。

1 水庫防洪調(diào)度多目標(biāo)風(fēng)險分析

本文考慮入庫流量和區(qū)間支流流量的相關(guān)性，采用Copula-MC法對干支流洪水同時進(jìn)行模擬，在對模擬洪水進(jìn)行防洪優(yōu)化調(diào)度的基礎(chǔ)上，計算不同調(diào)度方案的防洪目標(biāo)值與風(fēng)險率，并采用基于組合權(quán)重的理想點法對調(diào)度方案進(jìn)行風(fēng)險決策，給出最優(yōu)調(diào)度方案。整個防洪調(diào)度多目標(biāo)風(fēng)險分析的過程如圖1所示。

圖1 防洪調(diào)度多目標(biāo)風(fēng)險分析流程圖Fig.1 Flowchart of multiobjective risk analysis for flood control operation

1.1 干支流洪水隨機(jī)模擬

當(dāng)水庫下游有支流匯入時，下游防護(hù)區(qū)的安全需在調(diào)度過程中同時兼顧區(qū)間來水，故本文考慮入庫洪水和區(qū)間來水之間的相關(guān)性，利用Copula-MC法對干支流洪水同時進(jìn)行模擬。

一般情況下，干支流洪水間既存在著隨機(jī)性又有相關(guān)性，可采用Copula函數(shù)描述干流洪峰X和支流洪峰Y的聯(lián)合分布[9]。由Sklar定理[12]可知，一個二維聯(lián)合分布可以由兩個獨(dú)立的邊緣分布和一個Copula函數(shù)構(gòu)成。設(shè)干流洪峰X和支流洪峰Y的邊緣分布分別為u=FX(x)和v=FY(y)，可采用P-Ⅲ型分布[13]對它們進(jìn)行描述。采用Clayton Copula函數(shù)[14]作為聯(lián)接函數(shù)，則干支流洪峰的聯(lián)合分布可表示為：

C(u,v)=(u-θ+v-θ-1)-1/θ

(1)

式中：θ為Copula參數(shù)，可根據(jù)與Kendall秩相關(guān)系數(shù)的關(guān)系求出[14]。

隨機(jī)變量X和Y之間存在相關(guān)性，由Monte Carlo產(chǎn)生的二維隨機(jī)數(shù)須保持相應(yīng)的相關(guān)性。鑒于邊緣分布u和v是服從[0,1]的均勻分布，則支流洪峰為某一定值時干流洪峰的條件分布為：

C(u|v)=[1+vθ(u-θ-1)]-(1+θ)/θ

(2)

由于條件分布C(u|v)也服從[0,1]之間的均勻分布，因此采用Copula-Monte Carlo方法對干支流洪峰進(jìn)行隨機(jī)模擬的基本思路為：首先由計算機(jī)產(chǎn)生二維隨機(jī)數(shù)m1和m2，令v=m1；然后通過求解方程m2=C(u|v)得到u；u和v分別為兩隨機(jī)變量不超過某值的概率，可通過二者的邊緣分布推求出X和Y的隨機(jī)數(shù)值。

根據(jù)上述方法可推求出干支流洪水的洪峰流量，然后選取相應(yīng)的典型洪水放大得到各自的洪水過程。

1.2 防洪優(yōu)化調(diào)度

水庫防洪任務(wù)主要有三類，即保證壩體自身的安全、減輕下游洪澇災(zāi)害和減少上游淹沒損失。本文主要考慮水庫自身安全和下游防洪控制區(qū)的安全，選取最高壩前水位最低和下游防洪控制區(qū)最大流量最小為目標(biāo)函數(shù)，建立防洪優(yōu)化調(diào)度模型，其目標(biāo)函數(shù)為

Zm=min[maxZ(t)]

(3)

qm=min{max[Qck(t)+Qqj(t)]}

(4)

式中：Z(t)表示各時段的壩前水位；Qck(t)表示水庫各時段的下泄流量；Qqj(t)表示區(qū)間洪水的各時段流量。

約束條件：主要包括水量平衡約束、最大泄流能力約束、調(diào)度期末水位約束、下泄流量約束、水位庫容約束及防洪調(diào)度規(guī)程等[15]。

在防洪調(diào)度準(zhǔn)則和約束條件的限制下，可利用權(quán)重協(xié)調(diào)法[16]將多目標(biāo)優(yōu)化問題轉(zhuǎn)化為單目標(biāo)優(yōu)化問題，進(jìn)而利用動態(tài)規(guī)劃法進(jìn)行求解，得到非劣解集。

1.3 風(fēng)險率推求

在考慮干支流來水不確定性的情況下，不同調(diào)度方案的防洪目標(biāo)值及其帶來的風(fēng)險是決策者們所關(guān)心的最重要的問題。因此有必要對每一種方案進(jìn)行風(fēng)險估計，確定不同目標(biāo)突破安全閾值的風(fēng)險率。

將上述隨機(jī)模擬得到的干支流洪水序列作為防洪優(yōu)化調(diào)度模型的輸入條件，通過求解可得到不同方案下(即非劣方案)壩前最高水位系列Zm和下游防洪區(qū)最大流量序列qm，分別統(tǒng)計壩前最高水位超出安全閾值(防洪高水位)和下游防洪區(qū)組合流量超過安全閾值(防洪控制點安全過流量)的次數(shù)，計算其越限的概率，即為防洪控制區(qū)的出險率和水庫自身的出險率。

1.4 基于組合權(quán)重的理想點法多目標(biāo)風(fēng)險決策

多目標(biāo)非劣方案的綜合評價方法有很多，但在水庫防洪運(yùn)用中，為追求防洪效益最佳，應(yīng)盡可能選取能夠反映目標(biāo)值與最優(yōu)值接近程度的評價方法，其中理想點法[17]不失為一種切實可行的方法。由于各個目標(biāo)的相對重要程度不同，可將各指標(biāo)的權(quán)重系數(shù)與評價方法相結(jié)合，綜合考慮各目標(biāo)的主客觀權(quán)重，建立起基于組合權(quán)重的理想點法，對非劣解進(jìn)行排序，選出令決策者滿意的方案。整個決策過程主要分為兩步：首先通過層次分析法和熵權(quán)法[18]分別獲得主觀權(quán)重wzj和客觀權(quán)重wkj，然后將主客觀權(quán)重進(jìn)行組合[19]，得到組合權(quán)重wj=αwzj+(1-α)wkj，其中α代表決策者對主觀經(jīng)驗和客觀數(shù)據(jù)的偏好程度；然后將根據(jù)組合權(quán)重，采用理想點法[20]對各方案進(jìn)行排序，選出最優(yōu)方案。

2 實例應(yīng)用

向家壩水庫是金沙江水電開發(fā)基地的最后一級，工程開發(fā)任務(wù)以發(fā)電為主，同時改善上、下游通航條件，結(jié)合防洪攔沙，兼顧灌溉，并對溪洛渡水庫起反調(diào)節(jié)等作用。水庫正常蓄水位380.00 m，死水位和防洪限制水位均為370.00 m，總庫容51.63 億m3，調(diào)節(jié)庫容和防洪庫容均為9.03 億m3，具有季調(diào)節(jié)性能。水庫汛期不僅可以與溪洛渡水庫聯(lián)合運(yùn)用，提高下游川江河段沿岸宜賓、瀘州、重慶等城市的防洪標(biāo)準(zhǔn)，還可以配合三峽水庫運(yùn)用為長江中下游地區(qū)發(fā)揮積極的防洪作用。水庫下游左岸有長江一級支流岷江的匯入。由于岷江是長江流域水量最大的支流，汛期暴雨洪水活動頻繁，與金沙江洪水遭遇后，會增加水庫下游李莊防洪控制點的風(fēng)險。

分別選取屏山站和高場站1951-2010年實測流量資料和歷史特大洪水資源來擬合金沙江洪峰和岷江洪峰的邊緣分布。采用年最大法取樣，分別得到這兩個站的洪峰流量序列，其邊緣分布可用PIII型分布函數(shù)表示，采用線性矩法估計參數(shù)，并采用χ2檢驗方法進(jìn)行假設(shè)檢驗，結(jié)果如表1所示。

表1 P-Ⅲ型分布參數(shù)及檢驗結(jié)果Tab.1 Parameters and hypothesis test results ofP-Ⅲ margin distributions

由表1可見，在5%的顯著性水平下，自由度為k-r-1的χ2檢驗的結(jié)果均小于臨界值(表示括號中數(shù)據(jù)為檢驗的臨界值)，落在接受域內(nèi)，因此這2個站的P-Ⅲ型分布都通過假設(shè)檢驗。經(jīng)統(tǒng)計計算可得金沙江和岷江洪水洪峰的秩相關(guān)系數(shù)為τ=0.123，由此可根據(jù)文獻(xiàn)[14]推出Clayton Copula函數(shù)的參數(shù)θ=0.28。

采用上述Copula-Monte Carlo方法，可以得到10 000組金沙江洪水和岷江區(qū)間洪水的洪峰流量，如圖2所示。選取1998年8月份洪水為典型洪水，根據(jù)同倍比原則按照洪峰進(jìn)行放大，即可得到模擬的金沙江和岷江相應(yīng)的10 000場洪水過程。

圖2 干支流洪水洪峰隨機(jī)模擬值Fig.2 Simulated flood peaks of mainstream and its tributary

將模擬得到的10 000場金沙江洪水和岷江洪水帶入防洪優(yōu)化調(diào)度模型式(3)和式(4)中，通過求解可以得到非劣解方案，每個解都可用其目標(biāo)均值來表示，如圖3所示。

圖3 向家壩防洪優(yōu)化調(diào)度非劣解集Fig.3 Pareto solutions of flood control optimal operation in Xiangjiaba reservoir

向家壩水庫防洪運(yùn)用時，下游防洪控制區(qū)李莊的安全流量是40 000 m3/s，水庫自身的安全水位為380 m。根據(jù)計算得到的9組非劣調(diào)度方案，統(tǒng)計兩目標(biāo)各自超過其安全閾值的次數(shù)，即可得到防洪控制區(qū)的出險率和水庫自身的出險率。表2給出了各非劣方案的防洪目標(biāo)均值與相應(yīng)的風(fēng)險率。

由表2可見，從方案一至方案九，防洪控制區(qū)的最大流量呈逐漸降低的趨勢，其風(fēng)險率也隨之減小，而壩前最高水位逐漸升高，其風(fēng)險率呈現(xiàn)出非減的趨勢。其中，方案一的下游防洪區(qū)的風(fēng)險最高且最大下泄流量最大，方案九的壩前水位最高，且對于水庫自身防洪具有一定的風(fēng)險。

將表2中防洪目標(biāo)的均值和風(fēng)險率作為原始評價指標(biāo)，同時將各指標(biāo)均視為成本型指標(biāo)，根據(jù)文獻(xiàn)[17]中的標(biāo)準(zhǔn)0～1變換對其進(jìn)行歸一化處理，可得到標(biāo)準(zhǔn)化指標(biāo)矩陣：

表2 各防洪調(diào)度方案目標(biāo)均值及風(fēng)險率Tab.2 Objective mean value and risk of different floodcontrol operation scheme

方案防洪區(qū)最大流量均值/(m3·s-1)風(fēng)險率/%壩前最高水位均值/(m3·s-1)風(fēng)險率/%一32921．4212．72370．260．02二32768．1412．33370．390．02三32582．6111．76370．550．03四32384．1511．25370．770．03五32098．7410．31371．090．03六31748．899．12371．620．18七31448．347．73372．430．61八31043．556．27374．092．01九30676．124．70376．074．19

采用熵權(quán)法可以得到各評價指標(biāo)的客觀權(quán)重，對于主觀權(quán)重，依據(jù)側(cè)重下游防洪區(qū)安全、側(cè)重水庫自身安全兩種偏好給出兩組主觀權(quán)重。將主客觀組合因子皆取為0.5，則可得到不同側(cè)重對象下的組合權(quán)重，具體見表3。

表3 評價指標(biāo)權(quán)重Tab.3 Weights of different evaluation index

據(jù)此采用理想點法對各組調(diào)度進(jìn)行排序，具體結(jié)果見表4。

由表4可見水庫自身的風(fēng)險率和風(fēng)險損失較低，最優(yōu)方案基本集中在中下部，表明排序結(jié)果合理。根據(jù)側(cè)重下游安全、側(cè)重水庫安全兩種決策偏好不同，最優(yōu)有所差別，當(dāng)側(cè)重水庫下游安全時，下游最大下流流量較小的方案八更優(yōu)；側(cè)重水庫自身安全時，壩前水位較低且最大下泄流量不大的方案七更優(yōu)。因此，決策者可以根據(jù)自己的主觀偏好給予各目標(biāo)不同的權(quán)重，最終決策出的方案必定是一個綜合協(xié)調(diào)了下游安全和水庫自身安全的均衡解。

表4 各調(diào)度方案優(yōu)選排序結(jié)果Tab.4 Ranking results of each operation scheme

3 結(jié) 語

本文以向家壩水庫為研究對象，以干支流洪水不確定性作為對防洪優(yōu)化調(diào)度過程中的主要風(fēng)險因子，采用Copula-Monte Carlo法對其進(jìn)行定量描述，在對模擬洪水進(jìn)行防洪優(yōu)化調(diào)度的基礎(chǔ)上，對不同調(diào)度方案的防洪目標(biāo)值與風(fēng)險進(jìn)行計算，并采用基于組合權(quán)重的理想點法對調(diào)度方案進(jìn)行風(fēng)險決策，選出最優(yōu)調(diào)度方案。結(jié)果表明：①所建模型能對水庫汛期干支流洪水不確定性給向家壩水庫給上下游防洪帶來的風(fēng)險進(jìn)行估計、評價和決策，有助于決策者制定出效益與風(fēng)險均衡的防洪調(diào)度決策方案；②水庫下游防洪控制區(qū)的過流量越小，其超過安全泄量的風(fēng)險也越小，但水庫壩前最高水位會增加，且最高壩前水位超過防洪安全水位的風(fēng)險也會增加，因此，控制下游防洪區(qū)的過流量，勢必會以增加水庫自身的防洪風(fēng)險為代價；③針對決策者的主觀傾向不同，其決策出的最優(yōu)方案也會有差異，評價時需謹(jǐn)慎處理下游防洪區(qū)安全和水庫安全兩個目標(biāo)之間的矛盾。

□

[1] 傅 湘,王麗萍,紀(jì)昌明. 洪水遭遇組合下防洪區(qū)的洪災(zāi)風(fēng)險率估算[J]. 水電能源科學(xué),1999,(4):23-26.

[2] 曾玉紅,胡敏良,梁在潮. 防洪系統(tǒng)風(fēng)險分析及其預(yù)報閾值研究[J]. 武漢大學(xué)學(xué)報(工學(xué)版),2003,(6):27-30.

[3] 程衛(wèi)帥,陳 進(jìn). 防洪體系統(tǒng)風(fēng)險評估模型研究[J]. 水科學(xué)進(jìn)展,2005,(1):114-120.

[4] 刁艷芳,王本德. 基于不同風(fēng)險源組合的水庫防洪預(yù)報調(diào)度方式風(fēng)險分析[J]. 中國科學(xué):技術(shù)科學(xué),2010,(10):1 140-1 147.

[5] 李英海,周建中,張勇傳,等. 水庫防洪優(yōu)化調(diào)度風(fēng)險決策模型及應(yīng)用[J]. 水力發(fā)電,2009,(4):19-21,37.

[6] 覃 暉,李清清,周建中. 基于相對優(yōu)勢度的水庫防洪調(diào)度多屬性風(fēng)險決策方法研究[J]. 長江科學(xué)院院報,2011,(12):58-63.

[7] 王麗萍,黃海濤,張驗科,等. 梯級水庫群聯(lián)合防洪調(diào)度風(fēng)險估計模型[J]. 中國農(nóng)村水利水電,2014,(1):69-72,76.

[8] 周研來,梅亞東. 基于Copula函數(shù)和Monte Carlo法的防洪調(diào)度風(fēng)險分析[J]. 水電能源科學(xué),2010,(8):37-39.

[9] 閆寶偉,郭生練,陳 璐,等. 長江和清江洪水遭遇風(fēng)險分析[J]. 水利學(xué)報,2010,(5):553-559.

[10] 閆寶偉,郭生練,余 維. 長江和清江洪水過程遭遇風(fēng)險分析[J]. 水力發(fā)電學(xué)報,2013,(1):50-53.

[11] 劉章君,郭生練,胡 瑤,等. 基于Copula-Monte Carlo法的水庫下游洪水概率分布研究[J]. 水力發(fā)電,2015,(8):17-22.

[12] Nelsen R B. An introduction to Copulas[M]. Springer, New York, 1999.

[13] 中華人民共和國水利部. 中華人民共和國水利行業(yè)標(biāo)準(zhǔn):水利水電工程設(shè)計洪水計算規(guī)范[M]. 北京：中國水利水電出版社, 2006.

[14] 郭生練,閆寶偉,肖 義,等. Copula函數(shù)在多變量水文分析計算中的應(yīng)用及研究進(jìn)展[J]. 水文,2008,(3):1-7.

[15] 陳森林. 水電站水庫運(yùn)行與調(diào)度[M]. 北京:中國電力出版社, 2008.

[16] 賈本有,鐘平安,陳 娟,等. 復(fù)雜防洪系統(tǒng)聯(lián)合優(yōu)化調(diào)度模型[J]. 水科學(xué)進(jìn)展,2015,(4):560-571.

[17] 劉 方. 水庫水沙聯(lián)合調(diào)度優(yōu)化方法與應(yīng)用研究[D]. 北京：華北電力大學(xué),2013.

[18] Liping Wang，Yanke Zhang, Changming Ji, et al. Risk calculation method for complex engineering system [J]. Water Science and Engineering, 2011,4(3):345-355.

[19] 王漢斌,楊 鑫. 一種基于AHP-RS的組合權(quán)重確定方法[J]. 中國安全生產(chǎn)科學(xué)技術(shù),2010,(6):155-160.

[20] 方志耕,劉思峰. 決策理論與方法[M]. 北京:科學(xué)出版社, 2009.