張鑫厚，陳伏龍

(石河子大學(xué)水利建筑工程學(xué)院，新疆 石河子 832003)

隨著全球氣候變暖和人類活動對流域下墊面干預(yù)的加劇，水文序列失去了一致性，洪水樣本的獨立同分布假定已遭到極大挑戰(zhàn)[1]，采用傳統(tǒng)頻率分析方法得到的工程設(shè)計結(jié)果的可靠性倍受質(zhì)疑[2-3]，因此，采用非一致性洪水頻率計算方法分析水庫防洪安全顯得十分必要。水庫下游用水需求的增加對水庫汛期的防洪調(diào)度提出了更高的要求，即不僅能夠在汛期保證水庫的安全，還要充分利用洪水資源[4]，因此對水庫進行防洪風(fēng)險分析具有實際意義。目前，國內(nèi)許多專家學(xué)者對此作了研究，姜樹海[5]采用事故樹分析方法分析了漫壩失事的形成原因，得出了相應(yīng)的風(fēng)險率，并論證了合理選擇大壩防洪設(shè)計標(biāo)準(zhǔn)的原則和方法；李響等[6]將預(yù)泄能力約束法與蒙特卡洛隨機模擬相結(jié)合，同時考慮洪水預(yù)報誤差及洪水過程線形狀不確定性等因素，在控制汛期防洪風(fēng)險的條件下推求三峽水庫汛限水位動態(tài)控制域；吳成國等[7]以潘家口水庫為例，將三角模糊數(shù)理論運用于水庫防洪調(diào)度風(fēng)險分析，結(jié)果表明該方法能夠全面地估計調(diào)度過程中的各種風(fēng)險因子；李大鳴等[8]采用Copula函數(shù)構(gòu)建洪峰洪量兩變量模型反映水文不確定性，并應(yīng)用Monte Calro方法計算了桃林口水庫的防洪風(fēng)險率；黃凱[9]采用水文變異診斷系統(tǒng)對比分析了過去、現(xiàn)在條件下的極限防洪風(fēng)險率，從而進一步提高了環(huán)境影響下的預(yù)報精度；李平[10]等將貝葉斯網(wǎng)絡(luò)引入梯級水庫防洪風(fēng)險計算中，結(jié)果表明貝葉斯網(wǎng)絡(luò)方法能直觀、簡便地分析出多種風(fēng)險因素共同作用下的水庫群連潰風(fēng)險。

干旱環(huán)境、氣候變化和人類活動的共同影響，以冰雪融水為基礎(chǔ)的水資源系統(tǒng)非常脆弱，使得區(qū)域水循環(huán)系統(tǒng)的穩(wěn)定性和水資源的可再生性降低，極端水文事件發(fā)生的頻度和強度增大[11]。瑪納斯河流域以融雪洪水與暴雨洪水混合為主，其洪水總量較高，洪水過程線多寬緩等特點[12]，對該流域年最大洪量序列的研究符合洪水特點，對確保水利工程防洪度汛安全和制定防洪調(diào)度方案具有重要意義，對于肯斯瓦特水庫防洪風(fēng)險的研究結(jié)果可為瑪納斯河水庫防洪及汛期洪水資源利用提供指導(dǎo)依據(jù)。

1 研究區(qū)概況

本文以瑪納斯河上游肯斯瓦特控制流域段為研究區(qū)域，如圖1所示?？纤雇咛厮恼臼乾敿{斯河干、支流匯合后的出山口控制站，海拔約910 m，控制流域面積為4 637 km2，多年平均年徑流量12.21億m3?？纤雇咛厮畮炀哂蟹篮?、灌溉、發(fā)電等多種功能其正常蓄水位990 m，最大壩高129.4 m，總庫容1.88億m3，控制灌溉面積21.10萬hm2，屬于大(2)型工程；其設(shè)計洪水標(biāo)準(zhǔn)為500年一遇，相應(yīng)的洪峰流量為2 382 m3/s，校核洪水標(biāo)準(zhǔn)為5 000年一遇，相應(yīng)的洪峰流量為3 601 m3/s，下游防洪保護標(biāo)準(zhǔn)為50年一遇，相應(yīng)的洪峰流量為1 249 m3/s。

圖1 肯斯瓦特控制流域

2 洪量序列變異點檢驗

2.1 基于Pettitt法變異點檢驗

采用Pettitt非參數(shù)檢驗法[13]對瑪納斯河肯斯瓦特控制流域融雪洪水年最大洪量序列進行變異點檢驗，結(jié)果(圖2)顯示：年最大一日洪量序列可能在1993和1995年發(fā)生變異，年最大三日洪量序列可能在1993、1994、1995和1996年生變異，年最大七日和十五日洪量序列可能在1994、1995和1996年發(fā)生變異。取統(tǒng)計量P值最大的年份為最可能發(fā)生的變異年份，因此，年最大一日洪量序列變異點為1993年，年最大三日、七日和十五日洪量序列變異點為1995年。

圖2 Pettitt法變異點檢驗結(jié)果

圖2(續(xù))

2.2 基于可變模糊法變異點檢驗

可變模糊集檢驗法基于陳守煜[14]提出的質(zhì)變與量變定理，由LI Jianzhu等[15]提出并應(yīng)用于變異點檢驗。該方法基于序列的均值進行計算，對于序列的異常值不敏感。

對瑪納斯河肯斯瓦特水庫控制流域融雪洪水年最大洪量序列分別選用基于三角模糊數(shù)的變異點檢驗法進行變異點檢驗，其結(jié)果如圖3所示。

圖3 可變模糊集法變異點檢驗結(jié)果

從圖3可以看出：年最大洪量序列均發(fā)生質(zhì)變，即發(fā)生跳躍性變異，其中，年最大一日洪量序列發(fā)生變異年份為1993年，年最大三日、七日和十五日洪量序列發(fā)生變異年份為1995年，綜合兩種方法認為該結(jié)果能真實反映出年最大洪量序列的變異情況。

3 洪量序列一致性修正

3.1 序列分解計算

根據(jù)謝平等[16]提出的“分解-合成”理論及非一致性水文序列頻率計算方法，洪水序列由兩種或兩種以上成分線性疊加而成。根據(jù)Pettitt檢驗與可變模糊集變異點檢驗的結(jié)果，肯斯瓦特控制流域1957—2014年年最大一日洪量序列在變異點1993年可分為2個子序列，年最大三日、七日、十五日洪量序列在變異點1995年可分為2個子序列，其結(jié)果見表1。

表1 年最大洪量序列分解計算成果

根據(jù)水文時間序列的線性疊加原理，可求得年最大洪量序列的確定性成分Yt，按公式(6)可求出肯斯瓦特控制流域1957—2014年年最大洪量序列剔除跳躍成分后的隨機序列，如圖4所示。此序列為過去條件下年最大洪量序列，受自然條件變化影響較小，可以滿足洪水序列的一致性要求。

圖4 年最大洪量序列分解結(jié)果

3.2 序列合成計算

假設(shè)過去條件下年最大洪量序列服從P-Ⅲ型分布，可采用優(yōu)化適線法[17]得到過去條件下年最大洪量序列。采用蒙特卡洛法隨機模擬生成500個隨機性成分Sp，并結(jié)合t時刻(2014年)求出確定性成分Yt，利用以下合成公式

Xt,p=Yt+Sp

(1)

獲得現(xiàn)狀條件下(2014年)年最大洪量序列的樣本，通過優(yōu)化適線法對其進行P-Ⅲ型分布擬合，得到現(xiàn)狀條件下年最大洪量序列統(tǒng)計參數(shù)，結(jié)果見表2，相應(yīng)頻率的設(shè)計值見表3。

由表3可以看出：現(xiàn)狀條件下的年最大洪量序列設(shè)計值比相應(yīng)頻率下2008規(guī)劃的設(shè)計值均有所增加，變化幅度為1%～20%，其中年最大七日洪量與年最大十五日洪量變化量均超過了10%。根據(jù)鄒全等[18]的研究可知，主要是因為流域氣溫顯著上升和氣候的變化導(dǎo)致水文極端事件發(fā)生的頻率增加。因此，為保證肯斯瓦特水庫在汛期安全運行的前提下，對洪水資源進行充分利用，有必要在非一致性年最大洪量序列條件下，對水庫極限防洪風(fēng)險進行分析。

表2 現(xiàn)狀條件下年最大洪量序列統(tǒng)計參數(shù)

表3 現(xiàn)狀條件下年最大洪量設(shè)計值

4 極限防洪風(fēng)險分析

近年來，許多國內(nèi)外學(xué)者針對于水文不確定性所帶來的風(fēng)險做了大量的研究，提出了一些計算防洪風(fēng)險率的方法，主要有頻率分析法、隨機微分方程法、隨機模擬法。本文選擇頻率分析法和基于三角模糊數(shù)風(fēng)險分析法對肯斯瓦特水庫現(xiàn)狀條件下不同頻率的融雪洪水設(shè)計值分別進行水庫極限防洪風(fēng)險率計算。

4.1 基于頻率分析法風(fēng)險分析

傳統(tǒng)的頻率分析法中假設(shè)水庫年調(diào)洪最高水位與年最大洪水出現(xiàn)的頻率相同，以年調(diào)洪最高水位等于不破壞水利工程指標(biāo)水位的洪水頻率作為水庫極限風(fēng)險率。具體的計算步驟為：首先指定一個水庫極限防洪風(fēng)險控制指標(biāo)Zd，計算出不同設(shè)計頻率Pi(i=1,2,…,l)下入庫設(shè)計洪水過程，然后根據(jù)起調(diào)水位Z0并結(jié)合水庫防洪調(diào)度規(guī)則進行調(diào)洪演算，通過不斷的試算，計算出l個最高庫水位Zmi(i=1,2,…,l)，最后建立Zmi～Pi經(jīng)驗頻率曲線，并依據(jù)此頻率曲線可由Zd值反查出水庫極限防洪風(fēng)險率Pf(Zm≥Zd)。

4.1.1 調(diào)洪規(guī)則

通過河道過流能力分析，結(jié)合防洪工程的總體布局，從防洪要求出發(fā)，設(shè)定汛期限制水位為984 m、下游泄量為500 m3/s，以此作為控制條件進行調(diào)洪演算，其中調(diào)洪規(guī)則為：

當(dāng)入庫洪水小于50年一遇的標(biāo)準(zhǔn)時，水庫水位低于防洪高水位992.66 m，控制水庫的洪水下泄，最大下泄洪水為500 m3/s。

當(dāng)入庫洪水大于50年一遇標(biāo)準(zhǔn)時，水庫運行按照水庫水位分時段控制：當(dāng)入庫洪水小于下游安全泄量500 m3/s時，來多少泄多少，維持水庫汛限水位984 m不變。當(dāng)水庫水位高于984 m低于防洪高水位992.66 m，且入庫洪水大于下游安全泄量500 m3/s時，下泄量按下游安全泄量500 m3/s下泄。當(dāng)水庫水位超過防洪高水位992.66 m時，若入庫洪水流量小于泄洪建筑物下泄能力，按入庫洪水流量下泄，維持防洪高水位；若入庫洪水流量超過泄洪建筑物下泄能力時，按泄洪建筑物的泄流能力自由下泄。退水段水位逐漸下降至汛限水位后，若入庫洪水流量小于泄洪建筑物泄流能力，則維持汛限水位不變；若入庫洪水流量大于泄洪建筑物泄流能力，則按泄洪建筑物的泄流能力自由下泄。

4.1.2 防洪風(fēng)險率計算

以1996年典型融雪洪水過程為基礎(chǔ)，通過同頻率放大法計算現(xiàn)狀條件下不同頻率設(shè)計融雪洪水過程，并作為入庫融雪洪水，根據(jù)水庫水位—庫容—泄量關(guān)系曲線，結(jié)合肯斯瓦特水庫調(diào)度運行規(guī)則，經(jīng)調(diào)洪演算計算出設(shè)計洪水條件下的調(diào)洪結(jié)果，見表4。

表4 現(xiàn)狀條件下水庫調(diào)洪成果表

按表4中不同設(shè)計頻率所對應(yīng)的汛期最高庫水位建立現(xiàn)狀條件下最高庫水位與設(shè)計頻率之間的相關(guān)關(guān)系，即Zmi—Pi經(jīng)驗頻率曲線(圖5)，并按此頻率曲線由Zd值反查出水庫的極限防洪風(fēng)險率Pf(Zm≥Zd)。

圖5 Zmi—Pi經(jīng)驗頻率曲線

從表4可知，現(xiàn)狀條件下肯斯瓦特水庫0.02%設(shè)計頻率所對應(yīng)的最高庫水位Zm=994.78 m。

以水庫校核洪水位Zd=993.66 m，由圖5查得其水庫極限防洪風(fēng)險率Pf(Zm≥Zd)=0.123 7%>0.02%。這表明肯斯瓦特水庫控制流域氣溫、降雨量以及融雪洪水徑流量的變化會導(dǎo)致現(xiàn)狀條件下水庫極限防洪風(fēng)險增大。

4.2 基于三角模糊數(shù)法風(fēng)險分析

4.2.1 三角模糊數(shù)

三角模糊數(shù)[7]定義了隸屬度函數(shù)的區(qū)間數(shù)，這一較確定性實數(shù)能更好地描述和表征復(fù)雜系統(tǒng)的波動性和不確定性特征[19]。

本文設(shè)a1、a2、a3分別為一模糊變量的最小可能值(下限)、中值和最大可能值(上限)，則可建立一個三角模糊數(shù)，即

A=(a1,a2,a3)，a1≤a2≤a3。

三角模糊數(shù)不同置信度水平所對應(yīng)的區(qū)間數(shù)長度不同。設(shè)α為置信度水平，且α∈[0,1]，則可將三角模糊數(shù)A轉(zhuǎn)化為與一定置信度水平α相對應(yīng)的區(qū)間數(shù)，記

[(a2-a1)α+a1,-(a3-a2)α+a3]，

稱Aα為三角模糊數(shù)A的α-截集，它實際上是置信度水平不低于α的數(shù)據(jù)集合。

4.2.2 抽樣方法

根據(jù)歷史實測洪水特征值之間相互獨立的特性，采用自回歸模型對洪水過程進行隨機模擬[20]。本文選用一階自回歸模型AR(1)，其基本形式為

Xt=u+φ1(Xt-1-u)+ηt，

(2)

式(2)中：Xt為時間序列，φ1為一階自回歸系數(shù)，ηt為獨立隨機序列，u為序列均值。

對于洪水序列，獨立隨機序列ηt一般采用P-Ⅲ型分布，故式(2)可變?yōu)?/p>

(3)

式(3)中，r1由尤爾-沃爾克方程求得，σ為序列標(biāo)準(zhǔn)差，φt為標(biāo)準(zhǔn)化的P-Ⅲ型分布。

本文采用此一階自回歸模型進行洪水特征量的隨機模擬，在肯斯瓦特水庫風(fēng)險分析中，以1996年一場大洪水作為典型洪水，隨機模擬10萬次洪水組成入庫洪水序列。

4.2.3 極限防洪模糊風(fēng)險率計算

以5 000年一遇洪水調(diào)洪所得最高水位994.78 m作為模糊變量的上限值，防洪高水位992.66 m作為模糊變量的下限值，校核洪水位993.35 m作為模糊變量的中值，同時對置信度水平α在[0,1]中取不同值，可得與不同置信度水平α相對應(yīng)的水庫校核水位模糊變化區(qū)間結(jié)果，見表5。

對隨機生成的10萬條入庫洪水過程線設(shè)定一置信度水平α值，可得一組水庫校核水位的模糊變化區(qū)間，再通過統(tǒng)計計算得到一組水庫極限防洪的模糊風(fēng)險率區(qū)間值，結(jié)果見表6。

表5 不同置信度水平下水庫校核水位的變化區(qū)間α

表6 不同置信度水平下水庫極限防洪的模糊風(fēng)險率區(qū)間α

從表5和表6可知：

(1)當(dāng)截集水平α=0.0時，水庫風(fēng)險率區(qū)間范圍最大；隨著截集水平α的逐漸增大，水庫極限防洪風(fēng)險指標(biāo)模糊區(qū)間的范圍逐漸減小。這表明數(shù)據(jù)的模糊性逐漸減弱，提高了風(fēng)險指標(biāo)模糊區(qū)間內(nèi)數(shù)據(jù)的可信度水平。

(2)當(dāng)截集水平α=1.0時，其所對應(yīng)的區(qū)間表示風(fēng)險指標(biāo)和風(fēng)險率的相對最概然區(qū)間。這反映出水庫洪水調(diào)度過程中從風(fēng)險出現(xiàn)到風(fēng)險加劇進而致災(zāi)的一個緩變過程，而采用傳統(tǒng)的方法進行水庫極限防洪風(fēng)險分析，風(fēng)險指標(biāo)的估計值為一個點值，因此，風(fēng)險模糊區(qū)間比傳統(tǒng)的風(fēng)險率計算確定值更符合實際，可以為決策者提供更多的參考信息。

(3)兩種方法所得到的極限防洪風(fēng)險率，均比0.02%要大。其主要原因在于流域氣溫呈顯著上升趨勢，溫度升高引起以冰雪融水為基礎(chǔ)的融雪洪水徑流的季節(jié)性變化，改變了流域降水、融雪洪水徑流的時空分布和原有產(chǎn)匯流過程，使融雪洪水徑流過程產(chǎn)生變化，造成融雪洪水特征序列的變異，使得肯斯瓦特水庫控制流域年最大洪量序列呈現(xiàn)增加的趨勢。

5 結(jié)論

本文通過對肯斯瓦特水庫入庫年最大洪量序列的分析，以及采用頻率分析法和三角模糊數(shù)風(fēng)險分析法對現(xiàn)狀條件下水庫極限防洪模糊風(fēng)險率的計算分析，得出以下結(jié)論：

(1)采用Pettitt非參數(shù)檢驗法及可變模糊集變異點檢測法對年最大洪量序列進行變異點分析，在考慮P值最大的情況下，得到年最大一日洪量變異年份為1993年，年最大三日、七日和十五日洪量變異年份為1995年。

(2)采用“分解-合成”理論對跳躍變異的年最大洪量序列進行一致性修正，得到現(xiàn)狀條件下各頻率對應(yīng)的設(shè)計值，與2008年規(guī)劃的設(shè)計值相比均有所增加，變化幅度為1%～20%。

(3)采用頻率分析法和三角模糊數(shù)風(fēng)險分析法對現(xiàn)狀條件下的各頻率入庫洪水進行極限防洪風(fēng)險率的計算分析，給出了模糊風(fēng)險區(qū)間，所得結(jié)果可為水庫防洪科學(xué)調(diào)度提供依據(jù)。