劉 磊，王 帥

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

1 引言

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

2 風(fēng)險計算方法

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

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

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

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

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

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

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

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

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

3 計算結(jié)果

3.1 調(diào)洪演算

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

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

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

圖1 三種典型入庫洪水過程線情況下的泄洪流量過程圖

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

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

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

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

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

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

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

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

圖4 ［Qo］直方圖

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

3.3 風(fēng)險計算

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

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

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

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

表3 確定性洪水過程與確定性洪水過程風(fēng)險計算結(jié)果表

4 結(jié)論

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

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

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