基于隨機微分方程的水庫防洪調(diào)度風(fēng)險分析

陳 娟，鐘平安，2，徐 斌

(1.河海大學(xué)水文水資源學(xué)院，江蘇南京 210098;2.河海大學(xué)水資源高效利用與工程安全國家工程研究中心，江蘇南京 210098)

水庫防洪調(diào)度是流域洪水管理的重要內(nèi)容，其各環(huán)節(jié)中都存在諸多不確定性［1-2］，包括水雨情信息采集中由于設(shè)備故障、通訊不暢、誤碼和量程不足等原因?qū)е碌男畔o法獲取或無法及時傳達、信息錯誤，實時洪水預(yù)報中水文氣象條件、模型結(jié)構(gòu)、模型參數(shù)等導(dǎo)致的預(yù)報誤差，調(diào)洪演算中的水庫泄流和庫容曲線等水力不確定性等。水庫防洪調(diào)度中的諸多不確定性因素導(dǎo)致了水庫防洪決策的不確定性［3-4］，因此對水庫防洪調(diào)度中的不確定性因素和決策風(fēng)險進行分析研究，具有學(xué)術(shù)意義和實踐價值［5-7］。

目前，水庫防洪調(diào)度風(fēng)險分析主要有兩種基本思路:數(shù)字試驗法［8-10］和解析法。(a)數(shù)字試驗法。首先，根據(jù)試驗數(shù)據(jù)或經(jīng)驗判斷確定防洪調(diào)度主要不確定性因素的概率分布，并采用隨機模擬技術(shù)隨機生成符合相應(yīng)分布的各主要不確定性因素的數(shù)值;然后，用模擬生成的洪水過程進行水庫調(diào)洪演算，得到模擬的防洪目標(biāo)因子(水位或泄量)樣本;最后，運用統(tǒng)計學(xué)方法分析得到在指定條件下防洪目標(biāo)破壞的頻率。數(shù)字試驗法可以解決復(fù)雜問題的求解，但計算工作量大，計算結(jié)果的可靠性、精度與樣本容量有密切關(guān)系，不能保證數(shù)值解與真值之間的一致性。(b)解析法。解析法是建立風(fēng)險因子的概率密度函數(shù)，通過積分或者微分計算得到系統(tǒng)的風(fēng)險。姜樹海等［11-12］、王軍［13］進行了解析法在水庫防洪調(diào)度風(fēng)險分析方面的研究，認(rèn)為水庫蓄洪量的隨機變化導(dǎo)致了水庫水位的隨機消長，在假定水庫蓄水量具有Wiener過程特性的前提下，推導(dǎo)了調(diào)洪演算隨機方程，求解了調(diào)洪過程中庫水位的概率密度分布。

調(diào)洪過程中不同時刻水庫蓄水量的隨機變化，是由入庫洪水的水文條件、泄流的水力條件以及水位庫容關(guān)系等不確定性因素導(dǎo)致的，因此筆者從這些不確定性因素出發(fā)，建立調(diào)洪演算隨機微分方程，推導(dǎo)水庫水位分布特征值與這些不確定性因素分布特征值的直接聯(lián)系，并以此為依據(jù)進行水庫防洪調(diào)度的風(fēng)險分析。

1 水庫調(diào)洪演算隨機微分方程

確定性調(diào)洪演算微分方程如下:

式中:V(t)——t時刻水庫蓄水量;Q(t)——t時刻入庫流量;q(t)——t時刻出庫流量。

在整個調(diào)洪過程中，入庫洪水的水文條件、泄流的水力條件、水位庫容關(guān)系的邊界條件等不確定性因素導(dǎo)致了不同時刻水庫蓄水量的隨機變化，這一隨機變化又將導(dǎo)致水庫水位的隨機消長，給水庫水位控制造成風(fēng)險。因此，采用隨機過程的概念研究水庫水位，有利于完整、全面地描述水位隨機變化過程。本文考慮的主要不確定性因素及其數(shù)學(xué)描述如下:

a.入庫洪水水文條件的不確定性。實時洪水預(yù)報中水文氣象信息、模型結(jié)構(gòu)、模型參數(shù)等不確定性導(dǎo)致了預(yù)報誤差，使得預(yù)報入庫洪水過程成為連續(xù)的隨機過程。如果預(yù)報入庫洪水過程為該隨機過程的均值線，則考慮不確定性的入庫流量可表述為

圖1 入庫流量預(yù)報誤差分布示意圖Fig.1 Schematic diagram of forecast error distribution of reservoir inflow

一般情況下，預(yù)報流量過程中各時段的ξ(t)是不同的(圖1)，直觀上認(rèn)為，離作業(yè)預(yù)報時間越近，預(yù)報精度越高，預(yù)報相對誤差的均方差越小，誤差分布密度函數(shù)越尖瘦，預(yù)報結(jié)果接近均值的概率越大;反之，離作業(yè)預(yù)報時間越遠(預(yù)見期越長)，相對誤差的均方差越大，誤差分布密度函數(shù)越矮胖，預(yù)報結(jié)果在均值左右較大區(qū)間變動［14］。

b.泄流水力條件的不確定性。在泄洪建筑物規(guī)模確定的情況下，受水庫水位和流量系數(shù)等水力參數(shù)的不確定性影響，q(t)亦可表達為一隨機過程［11-12］:

水庫防洪調(diào)度一般在較高水位上進行，在高水位時C一般較穩(wěn)定，因此q(t)亦可簡化表達為

即泄流的不確定性主要由水位的不確定性引起。

c.水位庫容關(guān)系的不確定性。由于量測、計算和繪圖的誤差，水庫淤積和庫岸坍塌等不確定性，導(dǎo)致了水位庫容關(guān)系的不確定性。因此，與庫水位H相應(yīng)的水庫蓄水量亦可表達為一隨機過程［11-12］:

根據(jù)以上分析與假定，建立如下調(diào)洪演算隨機微分方程:

式中(t)為t時刻考慮不確定性的水庫蓄水量。

2 基于隨機微分方程的水庫防洪調(diào)度風(fēng)險分析

2.1 各時刻水庫水位的均值與方差推導(dǎo)

將式(2)代入式(6)，并進行離散得

當(dāng)Δt較小時，水位變幅也較小，式(4)和式(5)可以分段擬合為如下線性關(guān)系:

式中:μt，νt——t時刻出庫泄流過程的線性擬合參數(shù);αt，βt——t時刻水位庫容關(guān)系的線性擬合參數(shù)。

如圖2所示，第一象限為水位庫容曲線，第二象限是根據(jù)以水位為判別指標(biāo)的分級防洪調(diào)度規(guī)則建立的水位與出庫流量的關(guān)系線，最上部的曲線為自由泄流時的泄流能力曲線。對于任意H(t)，可以以H(t)為中心，在泄流曲線和庫容曲線取兩點(在泄流曲線上取點時應(yīng)取同一段分段函數(shù)上的兩點)，線性擬合式(8)和式(9)，進而確定 μt，νt和 αt，βt。

將式(8)、式(9)帶入式(7)，整理得

圖2 泄流過程和庫容曲線分段線性擬合示意圖Fig.2 Schematic diagram of piecewise linear fitting for discharge process and storage capacity curve

任一時刻水庫水位服從概率分布，對式(10)求其均值過程，整理得

將求得的水位和水位均值帶入方差公式整理得

一般情況下，經(jīng)過水庫調(diào)蓄，水庫調(diào)洪水位的變幅比入庫流量的變幅要小得多，因此可以假定各時刻入庫預(yù)報流量的誤差分布與水位的誤差分布相互獨立，即cov(H(t)，ξ(t))=0，式(12)整理得

在入流過程、泄流過程、庫容曲線以及ξ(t)分布已知的情況下，可以由公式(11)和(13)計算任意時刻水位的均值和方差，進而進行調(diào)洪過程的風(fēng)險分析。

2.2 水庫調(diào)洪風(fēng)險定義與計算

按照所建立的調(diào)洪演算隨機微分方程和推導(dǎo)的水位特征值計算公式，可以得到任意時刻水庫水位的均值和方差，假定各時刻水庫水位分布可由均值和方差確定(如正態(tài)分布)，則可以得到任意時刻的水位分布(圖3)。以水庫水位為荷載，定義控制水位或決策者心理承受水位HR(水庫實時防洪調(diào)度時，HR可根據(jù)校核洪水位和后續(xù)降雨預(yù)留防洪庫容確定)［14］。定義t時刻的水庫調(diào)洪風(fēng)險率如下:

圖3 調(diào)洪風(fēng)險示意圖Fig.3 Schematic diagram of risks of flood routing

式中Pt，R為t時段的水庫調(diào)洪風(fēng)險率。

設(shè)整個洪水過程的時段數(shù)為T，令事件H(t)＞HR為At(t=0，1，2，…，T)，定義整個洪水過程的總風(fēng)險率為

假定各時段的水位分布相互獨立，即事件A0，A1，…，AT相互獨立，因此A0，A1，…，AT也相互獨立，則

3 算 例

大伙房水庫是一座以防洪灌溉為主，兼顧發(fā)電、養(yǎng)殖、城市供水的綜合利用水庫。水庫總庫容為22.68億m3，汛限水位為126.40 m，設(shè)計洪水位為136.63 m，校核洪水位為139.32 m。

根據(jù)洪水預(yù)報模型對大伙房水庫歷史洪水進行模擬預(yù)報，統(tǒng)計預(yù)報入庫洪峰相對誤差符合正態(tài)分布N(0，0.23772)。圖4為某次預(yù)報入庫洪水過程，預(yù)報洪峰流量為10800m3/s，時段長Δt=1h，時段數(shù)T=130，洪峰所在時序m=55。

圖4 預(yù)報入庫洪水過程Fig.4 Process of forecasted reservoir inflow

由于缺乏不同時刻的誤差分布信息，本文基于預(yù)報相對誤差隨時間線性增加的假定，按式(17)確定各時刻的誤差分布:

式中:uξ(t)——t時刻預(yù)報入庫流量相對誤差的均值;σξ(t)——t時刻預(yù)報入庫流量相對誤差的均方差。

由于式(11)和式(13)采用的是誤差的絕對值，采用式(18)轉(zhuǎn)換:

根據(jù)式(11)、式(13)逐時段遞推計算得到各時段末的水位均值和方差，具體步驟如下:(a)t=0，起調(diào)水位為H(t)，E(H(t))=H(t)，D(H(t))=0;(b)以H(t)為中心(圖2)，率定式(8)和式(9)中的參數(shù) μt，νt和 αt，βt;(c)由式(18)計算E(ξ(t))和D(ξ(t));(d)將E(ξ(t))，E(H(t))，αt，μt，νt代入式(11)，得E(H(t+1));(e)將D(ξ(t))，D(H(t))，αt，μt，νt代入式(13)，得D(H(t+1));(f)如果t≥T，計算結(jié)束，否則t=t+1，轉(zhuǎn)到步驟(b)。

表1 各時刻流量和水位分布特征值Table 1 Characteristic values ofreservoir inflow and water level at each moment

根據(jù)以上步驟計算各時刻水庫水位均值和均方差，結(jié)果見表1。

由表1可知，隨著預(yù)見期的增大，流量相對誤差的均方差越大，表明洪水預(yù)報結(jié)果穩(wěn)定性越差，因而隨著預(yù)見期的延長，水庫調(diào)度中應(yīng)用洪水預(yù)報結(jié)果的風(fēng)險增大。同時，調(diào)洪水位相對誤差的均方差比洪水預(yù)報流量相對誤差的均方差小得多，可見由于水庫的調(diào)蓄作用，預(yù)報流量誤差造成的風(fēng)險被衰減了。

調(diào)度期內(nèi)取不同的預(yù)報后續(xù)降雨值，計算后續(xù)降雨預(yù)留防洪庫容，根據(jù)校核洪水位和預(yù)留防洪庫容計算相應(yīng)的控制水位，利用式(14)和式(16)計算整個洪水過程的調(diào)洪總風(fēng)險率。由得到的控制水位值和相應(yīng)的調(diào)洪風(fēng)險率繪制PR～HR曲線，見圖5。

由圖5可見，隨著后續(xù)降雨的增多，水庫控制水位越低，整個洪水過程的調(diào)洪總風(fēng)險率越高。已知不同的預(yù)報后續(xù)降雨值，查PR～HR曲線可求得任意控制水位下的調(diào)洪風(fēng)險率，如后續(xù)降雨為80.6 mm，計算得控制水位為HR=135.60 m，查圖5得到風(fēng)險率為5.02%，表明整個洪水過程水位超過控制水位的總風(fēng)險率為5.02%。圖5在調(diào)度會商過程中可以為決策者提供必要的風(fēng)險信息。

圖5 不同控制水位時的調(diào)洪總風(fēng)險率Fig.5 Total risKof flood routing at different control water levels

4 結(jié) 語

水庫防洪調(diào)度是風(fēng)險調(diào)度，不確定性因素帶來了防洪決策風(fēng)險，筆者就水庫防洪調(diào)度中的不確定性因素對水庫安全影響開展研究，取得如下主要研究成果:

a.提出了水庫入庫洪水、出庫泄流以及水位庫容關(guān)系等不確定性因素的數(shù)學(xué)描述方法。

b.基于調(diào)洪演算隨機微分方程，推導(dǎo)了各種不確定性因素綜合影響下各時刻水庫水位分布均值和方差計算公式。

c.假定各時刻水庫水位分布可由均值和方差確定(如正態(tài)分布)，從而得到任意時刻的水位分布，并基于此假定提出了調(diào)洪過程中各時刻的風(fēng)險描述方法，以及整個洪水過程的總風(fēng)險率定量計算方法。

d.以大伙房水庫為背景開展了實例研究，計算結(jié)果符合水庫防洪調(diào)度的實際運行狀況，可用于定量評估調(diào)洪過程中任意時刻的風(fēng)險和總風(fēng)險，在調(diào)度會商過程中為決策者提供必要的風(fēng)險信息。

［1］黃強，劉招，閆正龍，等.洪水預(yù)報信息用于水庫防洪預(yù)報調(diào)度的風(fēng)險分析［J］.西北農(nóng)林科技大學(xué)學(xué)報:自然科學(xué)版，2008(6):200-204.(HUANG Qiang，LIU Zhao，YAN Zhenglong，et al.RisKanalysis for reservoir flood forecast operation based on flood forecast information［J］.Journal of Northwest A ＆ F University:Natural Science Edition，2008(6):200-204.(in Chinese))

［2］ ZHANG Yanping，WANG Guoli，PENG Yong，et al.RisKanalysis of dynamic control of reservoir limited water level by considering flood forecast error［J］.Science China Technological Sciences，2011，54(7):1888-1893.

［3］WU S，YANG J，TUNG Y.RisKanalysis for flood-control structure under consideration of uncertainties in design flood［J］.Natural Hazards，2012，58(1):117-140.

［4］LI Qiong，ZHOU Jianzhong，LIU Donghan，et al.Research on flood risKanalysis and evaluation method based on variable fuzzy sets and information diffusion［J］.Safety Science，2012，50(5):1275-1283.

［5］ZOU Qiang，ZHOU Jianzhong，ZHOU Chao，et al.The practical research on flood risKanalysis based on IIOSM and fuzzyα-cut technique［J］.Applied Mathematical Modeling，2012，36(7):3271-3282.

［6］鐘平安，曾京.水庫實時防洪調(diào)度風(fēng)險分析研究［J］.水力發(fā)電，2008(2):8-9.(ZHONG Ping’an，ZENG Jing.Research on risKanalysis of reservoir real-time flood control operation［J］.Water Power，2008(2):8-9.(in Chinese))

［7］LI Xiang，GUO Shenglian，LIU Pan，et al.Dynamic control of flood limited water level for reservoir operation by considering inflow uncertainty［J］.Journal of Hydrology，2010，391(1/2):124-132.

［8］王本德，張靜.基于隨機模擬的分類預(yù)報調(diào)度方式風(fēng)險分析［J］.水力發(fā)電學(xué)報，2009(2):8-13.(WANG Bende，ZHANG Jing.RisKanalysis on classified forecast dispatching mode based on stochastic simulative method［J］.Journal of Hydroelectric Engineering，2009(2):8-13.(in Chinese))

［9］周惠成，董四輝，鄧成林，等.基于隨機水文過程的防洪調(diào)度風(fēng)險分析［J］.水利學(xué)報，2006，37(2):227-232.(ZHOU Huicheng，DONG Sihui，DENG Chenglin，et al.RisKanalysis on flood control operation of reservoir based on stochastic hydrological process［J］.Journal of Hydraulic Engineering，2006，37(2):227-232.(in Chinese))

［10］朱元甡.基于風(fēng)險分析的防洪研究［J］.河海大學(xué)學(xué)報:自然科學(xué)版，2001，29(4):1-8.(ZHU Yuanshen.Flood defense based on risKanalysis［J］.Journal of Hohai University:Natural Sciences，2001，29(4):1-8.(in Chinese))

［11］JIANG Shuhai.Application of stochastic differential equations in risKassessment for flood releases［J］.Hydrological Sciences Journal，1998，43(3):349-360.

［12］姜樹海，范子武.水庫防洪預(yù)報調(diào)度的風(fēng)險分析［J］.水利學(xué)報，2004，35(11):102-107.(JIANG Shuhai，F(xiàn)AN Ziwu.RisKanalysis for flood control operation of reservoir［J］.Journal of Hydraulic Engineering，2004，35(11):102-107.(in Chinese))

［13］王軍.基于貝葉斯理論的洪水實時預(yù)報調(diào)度研究［D］.南京:河海大學(xué)，2010

［14］鐘平安.流域?qū)崟r防洪調(diào)度關(guān)鍵技術(shù)研究與應(yīng)用［D］.南京:河海大學(xué)，2006.