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

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

1 引言

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

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

2 汛限水位動(dòng)態(tài)控制風(fēng)險(xiǎn)定義及風(fēng)險(xiǎn)源分析

2.1 汛限水位動(dòng)態(tài)控制風(fēng)險(xiǎn)定義

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

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

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

2.2 汛限水位動(dòng)態(tài)控制風(fēng)險(xiǎn)源分析

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

一是入庫(kù)流量過(guò)程預(yù)報(bào)結(jié)果偏小，實(shí)際來(lái)水比預(yù)報(bào)來(lái)水大；二是下一場(chǎng)洪水發(fā)生的時(shí)間提前，水庫(kù)喪失了部分或全部興利預(yù)泄時(shí)間；三是下一場(chǎng)洪水入庫(kù)時(shí)間比預(yù)計(jì)的快，防洪預(yù)泄水量比預(yù)計(jì)的減少。

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

3 基于蒙特卡洛模擬的風(fēng)險(xiǎn)分析方法

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

3.1 風(fēng)險(xiǎn)因子的隨機(jī)模擬

3.1.1 預(yù)報(bào)入庫(kù)水量的隨機(jī)數(shù)生成

式中：εt為服從水量預(yù)報(bào)相對(duì)誤差系列分布的正態(tài)隨機(jī)數(shù)；ηt為標(biāo)準(zhǔn)化正態(tài)隨機(jī)數(shù)；為水量預(yù)報(bào)相對(duì)誤差系列均值；σ為水量預(yù)報(bào)相對(duì)誤差系列均方差。

3.1.2 洪水預(yù)報(bào)預(yù)見(jiàn)期的隨機(jī)數(shù)生成

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

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

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

作如下變換

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

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

3.2 基于蒙特卡洛模擬的風(fēng)險(xiǎn)率計(jì)算

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

計(jì)算洪水退水階段及下一場(chǎng)洪水起漲階段的洪量：

式中：t1、tm分別為預(yù)報(bào)洪水過(guò)程的起點(diǎn)與終點(diǎn)。

考慮洪量預(yù)報(bào)相對(duì)誤差，計(jì)算修正的入庫(kù)洪量W'：

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

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

第三，將考慮預(yù)報(bào)誤差的流量過(guò)程系列、洪水預(yù)報(bào)預(yù)見(jiàn)期系列以及連續(xù)無(wú)雨天數(shù)系列隨機(jī)組合，得到組合隨機(jī)數(shù)系列。

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

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

式中：V0為t2時(shí)刻的水庫(kù)蓄量。

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

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

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

4 實(shí)例分析

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

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

選取河口村水庫(kù)1980-2003年共8場(chǎng)洪水資料作為分析對(duì)象，選擇預(yù)見(jiàn)期T=4 h，分析預(yù)見(jiàn)期內(nèi)實(shí)測(cè)洪量與預(yù)報(bào)洪量的誤差。

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

表2 預(yù)見(jiàn)期4 h預(yù)報(bào)洪量相對(duì)誤差分布統(tǒng)計(jì)表

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

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

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

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

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

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

表4 基于蒙特卡洛模擬法的風(fēng)險(xiǎn)率計(jì)算表

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

5 總結(jié)

闡述了汛限水位動(dòng)態(tài)控制風(fēng)險(xiǎn)分析的基本流程，包括對(duì)洪量預(yù)報(bào)誤差、洪水預(yù)報(bào)預(yù)見(jiàn)期及連續(xù)無(wú)雨天數(shù)等風(fēng)險(xiǎn)因子的隨機(jī)數(shù)模擬和基于蒙特卡洛模擬的汛限水位動(dòng)態(tài)控制風(fēng)險(xiǎn)率計(jì)算。

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