劉 磊，王 帥

（1.河南省水利移民事務(wù)中心，河南 鄭州，450000；2.河南省焦作水文水資源勘測(cè)局，河南 焦作，454000）

1 引言

文章旨在考慮影響水庫(kù)泄洪風(fēng)險(xiǎn)主要不確定性變量概率分布基礎(chǔ)上，采用蒙特卡羅法模擬相對(duì)均值和標(biāo)準(zhǔn)差對(duì)風(fēng)險(xiǎn)的擾動(dòng)程度，結(jié)果表明均值對(duì)水庫(kù)泄洪風(fēng)險(xiǎn)的擾動(dòng)更敏感，洪水過(guò)程的敏感性也高于其他不確定性因素，洪水過(guò)程不確定性是引起水庫(kù)泄洪風(fēng)險(xiǎn)的主要因素。

2 風(fēng)險(xiǎn)計(jì)算方法

如何將洪水過(guò)程有效地進(jìn)行量化，王長(zhǎng)新等學(xué)者在其研究成果中提出了洪峰削減系數(shù)η的概念，即

式中，Qomax為經(jīng)水庫(kù)調(diào)蓄后的最大泄流量；Qimax為入庫(kù)洪水的洪峰流量。為獲得洪峰削減系數(shù)的實(shí)測(cè)數(shù)據(jù)樣本，作者基于洪水過(guò)程線實(shí)測(cè)數(shù)據(jù)經(jīng)調(diào)洪演算得到Qomax值，然后兩流量之比為洪峰削減系數(shù)，如式（1）所示。n條洪水過(guò)程線獲得n個(gè)洪峰削減系數(shù)值，即容量為n的洪峰削減系數(shù)η的樣本。

為了更好地分析洪水過(guò)程不確定性導(dǎo)致的水庫(kù)泄洪風(fēng)險(xiǎn)，有必要得出洪峰削減系數(shù)規(guī)律，采用Kolmogorov-Smirnov（KS）檢驗(yàn)法對(duì)數(shù)值結(jié)果進(jìn)行檢驗(yàn)。K-S檢驗(yàn)法定義為：觀測(cè)樣本理論分布F（x）與實(shí)測(cè)Fn（x）之間差值為指標(biāo)，從而判斷樣本實(shí)測(cè)數(shù)據(jù)服從某種分布。具體為：求出Dn=max|F（x）-Fn（x）|，當(dāng)實(shí)際觀測(cè)值Dn＞Dnα則拒絕此類(lèi)假設(shè)，否則接受假設(shè)。K-S檢驗(yàn)是根據(jù)數(shù)據(jù)本身統(tǒng)計(jì)規(guī)律進(jìn)行擬合的一種方法，是一種非參數(shù)檢驗(yàn)方法，常用于在樣本量較小時(shí)。

若式（3）成立，則可以認(rèn)為在顯著水平α上擬采用的分布是合理的，否則認(rèn)為樣本不服從擬采用分布。式中Dn擬服從隨機(jī)變量其分布依賴(lài)率，Dnα為顯著水平α上的臨界值。文章在α=1%α=5%和α=10%三個(gè)顯著性水平上計(jì)算實(shí)測(cè)數(shù)據(jù)是否服從某種分布規(guī)律，通常為極值-I型分布、正態(tài)分布、對(duì)數(shù)正態(tài)分布及P-III型分布。

在獲得不確定性變量的分布規(guī)律后，采用蒙特卡洛（Monte-Carlo）法計(jì)算系統(tǒng)的風(fēng)險(xiǎn)，Monte-Carlo 計(jì)算風(fēng)險(xiǎn)特別是失效概率很小，需要大量的數(shù)據(jù)進(jìn)行計(jì)算，模擬次數(shù)甚至可達(dá)到幾百萬(wàn)乃至更多，隨著計(jì)算機(jī)科學(xué)的不斷發(fā)展，這些問(wèn)題很好地被解決，也推進(jìn)了蒙特卡洛隨機(jī)模擬的應(yīng)用。

Monte-Carlo 法是一種很靈活的方法，其結(jié)果的準(zhǔn)確性得到了驗(yàn)證，且其能夠廣泛應(yīng)用于各種問(wèn)題，特別是針對(duì)非線性問(wèn)題或復(fù)雜系統(tǒng)問(wèn)題，Monte-Carlo 法具有很大的優(yōu)勢(shì)。針對(duì)本文研究?jī)?nèi)容，建立極限狀態(tài)方程：

Z＞0 時(shí)系統(tǒng)處于可靠狀態(tài)；當(dāng)Z=0 時(shí)系統(tǒng)處于極限狀態(tài)；當(dāng)Z＜0 時(shí)系統(tǒng)處于失效狀態(tài)。產(chǎn)生服從分布的隨機(jī)數(shù)代入式（4），統(tǒng)計(jì)其結(jié)果小于0的次數(shù)，計(jì)算次數(shù)與結(jié)果小于0的次數(shù)即為失效概率。

Monte-Carlo模擬風(fēng)險(xiǎn)值時(shí)，需要大量的計(jì)算，隨計(jì)算次增加，計(jì)算結(jié)果的精度更高。抽樣次數(shù)n和失效概率Pf之間存在如下關(guān)系：

式中，n為計(jì)算次數(shù)；Pf為系統(tǒng)風(fēng)險(xiǎn)值。

3 計(jì)算結(jié)果

3.1 調(diào)洪演算

式中，Qi為入庫(kù)流量；Q0為水庫(kù)調(diào)洪后出庫(kù)流量；下標(biāo)1為t時(shí)刻步變量，在t+△時(shí)刻為已知量；下標(biāo)2為在t+△t時(shí)刻變量。

根據(jù)式（6），計(jì)算了相同洪峰流量的三種典型入庫(kù)洪水過(guò)程線情況下的泄洪流量過(guò)程，如圖1所示。

如圖1（a）時(shí)所示洪水過(guò)程線，洪峰入庫(kù)時(shí)，庫(kù)區(qū)具有較高的裕度（水庫(kù)水位較低），更多的流量被蓄留，導(dǎo)致調(diào)洪后的出庫(kù)流量更小。具體地，調(diào)洪后Qomax=9601 m3/s，洪峰削減系數(shù)η=0.69。洪峰后移，水庫(kù)水位在洪峰來(lái)臨時(shí)蓄水越多，庫(kù)區(qū)裕度減小。當(dāng)洪峰靠中間時(shí)，Qomax=10 995 m3/s，η=0.79。由此可見(jiàn)，洪峰削減程度隨洪峰靠后而減弱。當(dāng)如圖洪水過(guò)程如圖1（c）所示時(shí)，Qomax=13 150 m3/s，η=0.94，此時(shí)洪峰入庫(kù)時(shí)，水庫(kù)水位高，入庫(kù)流量直接泄走，導(dǎo)致水庫(kù)削減洪峰能力小。比較圖1（a），1（b），1（c）可知，洪水過(guò)程具有不確定性，且洪峰削減系數(shù)能夠很好地量化洪水過(guò)程。為綜合考慮水庫(kù)泄洪安全，有必要考慮洪水過(guò)程不確定性計(jì)算泄洪風(fēng)險(xiǎn)。

圖1 三種典型入庫(kù)洪水過(guò)程線情況下的泄洪流量過(guò)程圖

3.2 不確定變量分布規(guī)律

為研究洪峰削減系數(shù)分布規(guī)律，擬采用K-S檢驗(yàn)法檢驗(yàn)?zāi)乘畮?kù)1970-2009 年的數(shù)據(jù)。每年入庫(kù)洪峰流量與調(diào)洪后最大泄洪量之比數(shù)據(jù)經(jīng)過(guò)擬合后，函數(shù)值與實(shí)測(cè)值之間的均方根誤差表示為：

均方根誤差值表示擬合函數(shù)的可靠程度，1970-2009年均方根誤差如圖2 所示，由圖2 可知，2002 年的洪峰削減系數(shù)數(shù)據(jù)最為可靠。為了使結(jié)果具有更強(qiáng)的適用性，采用某地水庫(kù)2002 年份監(jiān)測(cè)數(shù)據(jù)進(jìn)行研究為宜，2002 年實(shí)測(cè)洪峰削減系數(shù)和入庫(kù)洪峰流量直方圖如圖3 所示，并采用K-S 法進(jìn)行檢驗(yàn)，檢驗(yàn)結(jié)果如表1所示。

圖2 洪峰削減系數(shù)均方根誤差圖

由表1 可知，在三個(gè)顯著水平α= 0.10、0.05、0.01 上，洪峰削減系數(shù)和入庫(kù)洪峰流量服從正態(tài)分布，其均值分別為0.78、16 m3/s，標(biāo)準(zhǔn)差分別為0.12、223 m3/s。根據(jù)趙國(guó)藩研究結(jié)果，入庫(kù)洪峰流量假設(shè)服從正態(tài)分布的隨機(jī)變量，其相對(duì)均值為1，相對(duì)方差為0.05。而允許最大泄洪流量公式為：

圖3（a）（b） η和Qimax直方圖

表1 η和QimaxK-S檢驗(yàn)結(jié)果表

式中，c為泄洪孔流量系數(shù)；A為泄洪口面積；H0水庫(kù)正常蓄水為；B為水庫(kù)水深對(duì)應(yīng)面積。由于施工工藝的發(fā)展和規(guī)范，認(rèn)為B、H0、c和A為不確定性小，為方便計(jì)算假設(shè)其均為確定性變量。式（9）中僅考慮V2-V1不確定變量，V2-V1由實(shí)測(cè)累計(jì)入庫(kù)流量和出庫(kù)流量所得，計(jì)算結(jié)果直方統(tǒng)計(jì)結(jié)果如圖4所示。通過(guò)K-S檢驗(yàn)，［Qo］分布服從極值-I型分布，其均值為1 766 m3/s、標(biāo)準(zhǔn)差為53 m3/s。這種分布特點(diǎn)是泄洪容許流量大概率是偏小，泄洪時(shí)產(chǎn)生更大的風(fēng)險(xiǎn)。據(jù)此，導(dǎo)致泄洪風(fēng)險(xiǎn)的主要因素的分布均已獲得，如表2所示。

圖4 ［Qo］直方圖

表2 不確定性參數(shù)分布規(guī)律表

3.3 風(fēng)險(xiǎn)計(jì)算

敏感性分析可確定單一不確定性變量參數(shù)變化對(duì)水庫(kù)泄洪失效概率的影響程度，主要不確定性因素的相對(duì)均值μ和相對(duì)方差σ對(duì)風(fēng)險(xiǎn)的擾動(dòng)。假設(shè)相對(duì)標(biāo)準(zhǔn)差為0.95-1.05，采用Monte-Carlo法模擬得出主要不確定因素與風(fēng)險(xiǎn)關(guān)系，如圖5所示。

圖5 泄洪不確定性因素敏感性圖

由圖5 可知，來(lái)流量Qi、泄洪洞容許泄洪量［Qo］及洪峰削減系數(shù)η不確定性對(duì)失效概率的影響均不可忽略。其中洪水過(guò)程為最敏感因素，其原因?yàn)椋涸谙嗤榉搴拖嗤塾?jì)流量，當(dāng)洪峰出現(xiàn)較早，導(dǎo)致水庫(kù)水位急劇上升，產(chǎn)生更大的泄洪失效概率；當(dāng)洪峰出現(xiàn)較晚，洪峰之前累計(jì)經(jīng)過(guò)調(diào)洪獲得水庫(kù)更大的裕度。當(dāng)洪峰μ/μ0=0.95 時(shí)，考慮洪水過(guò)程不確定性的泄洪失效概率Pf=0.005 3；當(dāng)洪峰μ/μ0=1.05 時(shí)，考慮洪水過(guò)程不確定性的泄洪失效概率Pf=0.071?？芍?μ0=0.95-1.05范圍內(nèi)，泄洪失效概率增大了13倍。說(shuō)明洪水過(guò)程不確定性是最敏感的，但水庫(kù)泄洪風(fēng)險(xiǎn)對(duì)均值擾動(dòng)更敏感，而入庫(kù)洪峰流量和最大泄洪量的敏感性低于洪水過(guò)程。在水文資料不足情況下，可假設(shè)兩者服從正態(tài)分布對(duì)水庫(kù)泄洪風(fēng)險(xiǎn)進(jìn)行分析。根據(jù)圖5 所示曲線趨勢(shì)，當(dāng)相對(duì)方差為1.1時(shí)，洪水過(guò)程不確定性引起的水庫(kù)泄洪風(fēng)險(xiǎn)為其他兩者的56倍、63倍，洪水過(guò)程不確定是引起泄洪風(fēng)險(xiǎn)分析的主要因素。

為了闡明洪水過(guò)程不確定性重要性，計(jì)算了考慮η為確定性變量和不確定性變量?jī)煞N情況下的泄洪風(fēng)險(xiǎn)值相對(duì)差值為19%，如表3 所示。相對(duì)確定性分析風(fēng)險(xiǎn)值，采用不確定性分析風(fēng)險(xiǎn)值偏低，這是由于確定性分析基于假設(shè)洪峰靠后而得到風(fēng)險(xiǎn)值，并沒(méi)有考慮這一情況發(fā)生的概率，說(shuō)明采用確定性分析可能導(dǎo)致計(jì)算結(jié)果偏保守。蒙特卡羅分析失效概率，相對(duì)誤差隨模擬次數(shù)增加而減小，且在模擬次數(shù)為104次時(shí)，誤差迅速減小，模擬次數(shù)大于104時(shí)，相對(duì)誤差均小于0.5%。綜合考慮計(jì)算資源和結(jié)果精度，計(jì)算次數(shù)取104為宜。

表3 確定性洪水過(guò)程與確定性洪水過(guò)程風(fēng)險(xiǎn)計(jì)算結(jié)果表

4 結(jié)論

采用洪峰削減系數(shù)很好的量化洪水過(guò)程，并獲得入庫(kù)洪峰流量和容許泄洪量采用了蒙特卡羅法模擬了考慮洪水過(guò)程不確定性水庫(kù)泄洪風(fēng)險(xiǎn)，得到以下結(jié)論：洪峰削減系數(shù)能夠很好地量化洪水過(guò)程的，洪峰削減系數(shù)服從正態(tài)分布，其均值為0.78，方差為0.12。洪水過(guò)程對(duì)風(fēng)險(xiǎn)擾動(dòng)最大，可認(rèn)為洪水過(guò)程不確定是導(dǎo)致水庫(kù)泄洪風(fēng)險(xiǎn)的最主要因素。

不考慮洪水過(guò)程不確定性，蒙特卡羅模擬失效概率偏大，闡明考慮洪水過(guò)程為確定性的研究結(jié)果得到更大的泄洪風(fēng)險(xiǎn)值，說(shuō)明確定性分析結(jié)果偏保守。

蒙特卡羅分析失效概率，相對(duì)誤差隨模擬次數(shù)增加而減小，且在模擬次數(shù)大于104次時(shí)，誤差迅速減小且均小于0.50%。綜合考慮計(jì)算資源和結(jié)果精度，計(jì)算次數(shù)取104為宜。