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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 風險因子的隨機模擬

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

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

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

因洪水預報預見期亦服從正態(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 基于蒙特卡洛模擬的風險率計算

第一，在確定性預報洪水過程Q（t）的基礎上，考慮預報誤差構建不同的入庫流量過程。

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

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

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

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

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

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

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

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

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

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

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

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

4 實例分析

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

表1 洪水預報樣本選取表

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

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

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

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

考慮洪量預報誤差、連續(xù)無雨天數(shù)以及洪水預報預見期三因子隨機組合，采用蒙特卡洛方法生成10 000組隨機數(shù)，由于19830907 號洪水滾動預報長度不夠，故舍棄該場洪水，選取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，遭遇不同頻率設計洪水時的調(diào)洪最高水位作為特征水位，特征水位取值見表3，求得對應不同頻率設計洪水的最高水位Zmax系列，安全水位取防洪高水位285.43 m，統(tǒng)計最高水位超過安全水位的概率，即為風險率，結(jié)果如表4所示。

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

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

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

5 總結(jié)

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

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