基于蒙特卡洛法的水庫汛限水位動態(tài)控制風險分析

建劍波 曾智珍 盧金閣 孫志偉 何芳嬋

（1河南省河口村水庫管理局；2河海大學水利水電工程學院；3南京水利科學研究院；4河南省水利科學研究院）

1 引言

中國現(xiàn)階段絕大部分水庫多以靜態(tài)的方式對汛限水位進行控制，該方式雖然能夠保證水庫的防洪安全，但在汛前或汛期常常要進行大量棄水以降低至汛限水位，從而造成洪水資源大量浪費，不利于水庫的興利調(diào)度。中國是一個水資源缺乏和水資源時空分布不均的國家，近年來隨著經(jīng)濟的快速發(fā)展和城鎮(zhèn)化建設(shè)，水資源的供需矛盾越來越突出。為充分利用水資源，滿足可持續(xù)發(fā)展的需求，盡可能利用現(xiàn)有的降雨預(yù)報和洪水預(yù)報技術(shù)，優(yōu)化水庫調(diào)度運行，實施水庫汛限水位的動態(tài)控制，對提高水庫水資源的利用率，發(fā)揮水庫綜合效益，以及水資源可持續(xù)發(fā)展戰(zhàn)略都具有重要的意義。

文章以基于預(yù)報預(yù)泄的水庫洪水動態(tài)水位控制為研究對象，結(jié)合河口村水庫實例分析，對水文氣象預(yù)報等不確定性信息進行風險估計，為決策者提供理論依據(jù)和決策信息。

2 汛限水位動態(tài)控制風險定義及風險源分析

2.1 汛限水位動態(tài)控制風險定義

水庫汛限水位動態(tài)控制風險指的是由于水庫調(diào)度過程中自然、人為等不確定因素的影響，水庫系統(tǒng)在實施汛限水位動態(tài)控制前后的風險變化，主要從風險率變化和后果變化兩方面得以反映。

已建成的水庫，工程規(guī)模和功能是確定的。假設(shè)面臨調(diào)度期內(nèi)調(diào)洪最高水位Z為荷載，水庫運行階段允許最高水位Zm為承載能力，定義“Z＞Zm”為水庫洪水動態(tài)水位控制導致的水庫安全風險事件，則風險率可表達為：

式中字母含義同下。不考慮調(diào)度規(guī)則的變化，當起調(diào)水位Z0超過原設(shè)計汛限水位Z汛時，水庫最高水位Zmax可由式（2）～（4）確定：

2.2 汛限水位動態(tài)控制風險源分析

基于預(yù)報預(yù)泄的水庫汛期水位動態(tài)控制的實質(zhì)，是利用洪水尾水實現(xiàn)超原設(shè)計汛限水位蓄水，關(guān)鍵問題是分析實時增蓄水量的大小。增蓄水量在上一場洪水退水段形成，其消落有興利預(yù)泄與防洪預(yù)泄兩種方式。如果增蓄后無雨期很長，則增蓄的水量可以通過興利預(yù)泄消落，在興利預(yù)泄來不及時，可通過防洪預(yù)泄消落。根據(jù)基于預(yù)報預(yù)泄的水庫洪水動態(tài)水位風險控制方法的原理可知，導致水庫水位不能回落到原設(shè)計汛限水位的原因主要有以下方面：

一是入庫流量過程預(yù)報結(jié)果偏小，實際來水比預(yù)報來水大；二是下一場洪水發(fā)生的時間提前，水庫喪失了部分或全部興利預(yù)泄時間；三是下一場洪水入庫時間比預(yù)計的快，防洪預(yù)泄水量比預(yù)計的減少。

可見預(yù)報入庫水量的誤差、預(yù)泄采用的連續(xù)無雨天數(shù)，以及預(yù)泄采用的洪水預(yù)報預(yù)見期是導致水庫水位可能高于的三個不確定性（風險）因子，這三個因子既可能單獨，也可能組合發(fā)生作用。

3 基于蒙特卡洛模擬的風險分析方法

汛限水位動態(tài)控制的各風險變量之間存在著比較復雜的影響機制，而使用蒙特卡洛模擬方法（MC法）來解決復雜系統(tǒng)的風險分析是一種行之有效的方法。MC法在水庫洪水動態(tài)水位控制風險分析中按照以下步驟進行：①根據(jù)實驗數(shù)據(jù)或者經(jīng)驗判斷確定水庫洪水動態(tài)水位控制中主要不確定性因素的概率分布。②隨機生成符合相應(yīng)分布的各主要不確定性因素的數(shù)值。③將不確定性因素疊加到水庫洪水動態(tài)水位控制的過程中，調(diào)洪演算得到防洪目標（調(diào)洪最高水位）值。④重復步驟②～③，可得到目標樣本，運用統(tǒng)計學方法分析得到在指定條件下防洪目標破壞的頻率，作為風險計算的基本數(shù)據(jù)。

3.1 風險因子的隨機模擬

3.1.1 預(yù)報入庫水量的隨機數(shù)生成

式中：εt為服從水量預(yù)報相對誤差系列分布的正態(tài)隨機數(shù)；ηt為標準化正態(tài)隨機數(shù)；為水量預(yù)報相對誤差系列均值；σ為水量預(yù)報相對誤差系列均方差。

3.1.2 洪水預(yù)報預(yù)見期的隨機數(shù)生成

因洪水預(yù)報預(yù)見期亦服從正態(tài)分布，故其隨機數(shù)生成方法同3.1.1。

3.1.3 連續(xù)無雨天數(shù)的隨機數(shù)生成

連續(xù)無雨天數(shù)服從負指數(shù)分布Exp（λ，μ），即

作如下變換

當y 為［0，1］均勻分布隨機數(shù)時，d 即為服從式（7）的負指數(shù)分布。

因此，要生成符合負指數(shù)分布的隨機數(shù)，首先需要生成符合U［0，1］均勻分布的隨機數(shù)系列{yt}，將其代入式（7）中即可得到服從負指數(shù)分布Exp（λ，μ）的隨機數(shù)系列{dt}。

3.2 基于蒙特卡洛模擬的風險率計算

第一，在確定性預(yù)報洪水過程Q（t）的基礎(chǔ)上，考慮預(yù)報誤差構(gòu)建不同的入庫流量過程。

計算洪水退水階段及下一場洪水起漲階段的洪量：

式中：t1、tm分別為預(yù)報洪水過程的起點與終點。

考慮洪量預(yù)報相對誤差，計算修正的入庫洪量W'：

利用同倍比法放大，計算考慮預(yù)報誤差后修正的流量過程：

第二，用蒙特卡洛模擬法生成服從正態(tài)分布的洪水預(yù)報預(yù)見期系列{εt}，以及服從負指數(shù)分布的連續(xù)無雨天數(shù)系列{dt}。

第三，將考慮預(yù)報誤差的流量過程系列、洪水預(yù)報預(yù)見期系列以及連續(xù)無雨天數(shù)系列隨機組合，得到組合隨機數(shù)系列。

第四，將組合隨機數(shù)中的流量過程、洪水預(yù)報預(yù)見期以及連續(xù)無雨天數(shù)代入式（12）～（13）中。

計算得到從t0到t2水庫的蓄量變化△W'。已知t0時刻水庫完成洪水動態(tài)水位控制的預(yù)蓄過程且水位達到Zms，根據(jù)式（14）～（15）計算得到下一場洪水起漲時刻即t2時刻的水位Z0。然后，根據(jù)式（2）～（4），得到最高水位Zmax。

式中：V0為t2時刻的水庫蓄量。

第五，考慮不同的組合隨機數(shù)，重復步驟（4），得到最高水位Zmax系列。

第六，假設(shè)決策者所能接受的安全水位指標為Zd，則相應(yīng)的風險率為：

式中：m 為隨機試驗總次數(shù)；mf為調(diào)洪最高水位大于安全水位（防洪控制最高水位）的次數(shù)。

4 實例分析

單獨利用預(yù)蓄預(yù)泄確定河口村水庫汛期水位動態(tài)控制的上限值為246.00 m。單獨利用實時預(yù)報調(diào)度技術(shù)確定常遇洪水（20 a一遇以下）的汛期水位控制上限值為250.59 m，綜合考慮后將河口村水庫的主汛期水位動態(tài)控制的上限值取為246.00 m。

表1 洪水預(yù)報樣本選取表

選取河口村水庫1980-2003年共8場洪水資料作為分析對象，選擇預(yù)見期T=4 h，分析預(yù)見期內(nèi)實測洪量與預(yù)報洪量的誤差。

8場洪水共獲得有效樣本434個，經(jīng)統(tǒng)計4 h預(yù)見期洪量相對誤差分布見表2。

表2 預(yù)見期4 h預(yù)報洪量相對誤差分布統(tǒng)計表

假定洪量預(yù)報誤差符合正態(tài)分布，求得4 h預(yù)見期洪量預(yù)報誤差服從正態(tài)分布N（65.71，75.232）。同時，參考河口村水文資料可知，洪水預(yù)報預(yù)見期系列服從正態(tài)分布N（9.80，3.212），連續(xù)無雨天數(shù)系列（P≤3 mm）服從負指數(shù)分布Exp（0.19，1）。

考慮洪量預(yù)報誤差、連續(xù)無雨天數(shù)以及洪水預(yù)報預(yù)見期三因子隨機組合，采用蒙特卡洛方法生成10 000組隨機數(shù)，由于19830907 號洪水滾動預(yù)報長度不夠，故舍棄該場洪水，選取19800728、19920811、20010726 三場洪水進行風險分析，根據(jù)基于蒙特卡洛模擬的風險分析原理，分別得到下一場洪水起漲時刻的水位Z0系列，偏安全考慮，選取每種隨機模擬情景下三場洪水中的最高水位作為最終的結(jié)果，結(jié)果如圖1所示。

圖1 Z0系列水位分布直方圖與累積頻率曲線圖

在10 000次隨機模擬中，下一場洪水起漲時刻的水位Z0超過汛限水位238.00 m的次數(shù)為237次，概率為2.37%，圖1為Z0超過238.00 m的分布情況。

按照起調(diào)水位為238.00 m，遭遇不同頻率設(shè)計洪水時的調(diào)洪最高水位作為特征水位，特征水位取值見表3，求得對應(yīng)不同頻率設(shè)計洪水的最高水位Zmax系列，安全水位取防洪高水位285.43 m，統(tǒng)計最高水位超過安全水位的概率，即為風險率，結(jié)果如表4所示。

表3 前汛期起調(diào)水位238.00 m對應(yīng)的水庫各頻率洪水特征值（初步設(shè)計）表

表4 基于蒙特卡洛模擬法的風險率計算表

從表4可以看出，在主汛期水位動態(tài)控制值246.00 m的方案下，若考慮預(yù)報入庫水量的誤差ε1、預(yù)泄采用的連續(xù)無雨天數(shù)d，以及預(yù)泄采用的洪水預(yù)報預(yù)見期 這三個不確定性因子，該方案未增加500 a一遇及以下量級的洪水的風險，稍增加2 000 a一遇洪水的風險，但增加的概率也僅為1.19×10-3%。

5 總結(jié)

闡述了汛限水位動態(tài)控制風險分析的基本流程，包括對洪量預(yù)報誤差、洪水預(yù)報預(yù)見期及連續(xù)無雨天數(shù)等風險因子的隨機數(shù)模擬和基于蒙特卡洛模擬的汛限水位動態(tài)控制風險率計算。

以河口村水庫為例，采用蒙特卡洛模擬方法分析了水庫汛期水位動態(tài)控制的風險，得到在汛限水位抬高至246.00 m情況下，水庫在下一場洪水起漲時刻超原汛限水位的概率僅為2.37%，該方案不增加500 a一遇及以下量級的洪水風險，稍增加2 000 a 一遇洪水的風險，汛期水位動態(tài)控制產(chǎn)生的風險極低。