基于Copula-Monte Carlo法的水庫防洪調(diào)度多目標風險分析

陳 凱，彭 楊，吳志毅

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

0 引 言

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

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

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

1 水庫防洪調(diào)度多目標風險分析

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

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

1.1 干支流洪水隨機模擬

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

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

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

(1)

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

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

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

(2)

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

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

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

水庫防洪任務主要有三類，即保證壩體自身的安全、減輕下游洪澇災害和減少上游淹沒損失。本文主要考慮水庫自身安全和下游防洪控制區(qū)的安全，選取最高壩前水位最低和下游防洪控制區(qū)最大流量最小為目標函數(shù)，建立防洪優(yōu)化調(dià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)度準則和約束條件的限制下，可利用權重協(xié)調(diào)法[16]將多目標優(yōu)化問題轉(zhuǎn)化為單目標優(yōu)化問題，進而利用動態(tài)規(guī)劃法進行求解，得到非劣解集。

1.3 風險率推求

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

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

1.4 基于組合權重的理想點法多目標風險決策

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

2 實例應用

向家壩水庫是金沙江水電開發(fā)基地的最后一級，工程開發(fā)任務以發(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)合運用，提高下游川江河段沿岸宜賓、瀘州、重慶等城市的防洪標準，還可以配合三峽水庫運用為長江中下游地區(qū)發(fā)揮積極的防洪作用。水庫下游左岸有長江一級支流岷江的匯入。由于岷江是長江流域水量最大的支流，汛期暴雨洪水活動頻繁，與金沙江洪水遭遇后，會增加水庫下游李莊防洪控制點的風險。

分別選取屏山站和高場站1951-2010年實測流量資料和歷史特大洪水資源來擬合金沙江洪峰和岷江洪峰的邊緣分布。采用年最大法取樣，分別得到這兩個站的洪峰流量序列，其邊緣分布可用PIII型分布函數(shù)表示，采用線性矩法估計參數(shù)，并采用χ2檢驗方法進行假設檢驗，結(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-Ⅲ型分布都通過假設檢驗。經(jīng)統(tǒng)計計算可得金沙江和岷江洪水洪峰的秩相關系數(shù)為τ=0.123，由此可根據(jù)文獻[14]推出Clayton Copula函數(shù)的參數(shù)θ=0.28。

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

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

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

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

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

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

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

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

方案防洪區(qū)最大流量均值/(m3·s-1)風險率/%壩前最高水位均值/(m3·s-1)風險率/%一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

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

表3 評價指標權重Tab.3 Weights of different evaluation index

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

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

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

3 結(jié) 語

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

□

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

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

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

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

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

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

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

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

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

[10] 閆寶偉,郭生練,余 維. 長江和清江洪水過程遭遇風險分析[J]. 水力發(fā)電學報,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è)標準:水利水電工程設計洪水計算規(guī)范[M]. 北京：中國水利水電出版社, 2006.

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

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

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

[17] 劉 方. 水庫水沙聯(lián)合調(diào)度優(yōu)化方法與應用研究[D]. 北京：華北電力大學,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的組合權重確定方法[J]. 中國安全生產(chǎn)科學技術,2010,(6):155-160.

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