混合分布在非一致性洪水頻率分析的應(yīng)用

馮 平，曾 杭，李 新

(天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072)

近年來，受氣候變化尤其是人類活動的影響，許多流域內(nèi)的下墊面條件發(fā)生了顯著變化．下墊面變化的發(fā)生，會改變流域原有降雨產(chǎn)匯流過程，水文序列會發(fā)生顯著變化，導(dǎo)致水文序列不滿足一致性要求，則采用傳統(tǒng)頻率計算方法得到的結(jié)果將會受到質(zhì)疑[1]．基于非平穩(wěn)極值系列的直接水文頻率分析方法有3種：①基于混合分布的非一致性水文頻率分析方法；②基于時變矩的非一致性水文頻率分析方法[2]；③基于條件概率分布的非一致性水文頻率分析方法[3]．筆者采用混合分布模型對非一致性水文頻率分析進行探討．Singh等[4]首次提出將混合分布模型應(yīng)用到非一致性洪水頻率計算，得到較好的效果，但參數(shù)估計是當時存在的最大問題；2002 年，Alila等[5]應(yīng)用此方法對 Gila流域進行了研究，發(fā)現(xiàn)混合分布模型比傳統(tǒng)的單分布模型更好．對于混合分布的參數(shù)估計問題，目前國內(nèi)外學(xué)者采用了很多方法，如Singh[6]采用非線性優(yōu)化算法、Rossi等[7]采用極大似然算法、Leytham[8]采用極大似然的 EM 算法以及Fiorentino等[9]采用最大熵準則法(POME)．非線性優(yōu)化算法是利用樣本系列的前五階中心距求解參數(shù)，對于小樣本而言高階矩會產(chǎn)生較大誤差，后3種估計方法計算過程較為復(fù)雜．成靜清等[10]采用模擬退火算法進行混合分布模型參數(shù)估計，應(yīng)用于非一致性年徑流序列的頻率計算，取得較好效果；謝平等[11]對非一致性年徑流序列的頻率計算也取得了一些研究進展．然而非一致性洪水序列頻率計算在我國缺乏系統(tǒng)成熟的方法，為此，筆者根據(jù)水文變異點的綜合診斷結(jié)果，假設(shè)非一致性洪水序列變異點前后分別服從2個皮爾遜Ⅲ型分布，全序列服從由這2個分布組成的混合分布，以大清河流域南支上游龍門水庫的入庫年最大流量序列和固定時段(1日、3日、6日)洪量序列為例，分析變異后計算得到的設(shè)計洪水值與不考慮變異情況下的設(shè)計洪水值、原設(shè)計洪水值(1985年審定通過)之間的變化及變化原因，以期為龍門水庫的設(shè)計洪水修訂提供依據(jù)，也為大清河流域南支上游其他大型水庫的設(shè)計洪水修訂提供技術(shù)方法．

1 非一致性洪水序列變異點的診斷方法

水文時間序列變異點綜合診斷方法首先是對洪水序列進行初步診斷，以檢驗序列是否存在確定性成分[12]，即采用過程線法、累積曲線斜率差異幅度分析法和Hurst系數(shù)法進行初步檢驗．判斷序列是否存在變異，若判斷結(jié)果不存在變異，則轉(zhuǎn)入成因調(diào)查分析，對結(jié)果進行確認；如存在變異，轉(zhuǎn)入詳細診斷部分．詳細診斷主要采用 Lee-Heghinian檢驗、有序聚類檢驗、最優(yōu)信息二分割檢驗、R/S檢驗和滑動 F檢驗等方法．最后進行綜合診斷，結(jié)合詳細診斷的可能變異點和物理成因分析，確定最可能變異點．

2 非一致性洪水序列的混合分布與參數(shù)估計

2.1 混合分布模型

早在 1972年 Singh等[4]首次提出將混合分布模型應(yīng)用到非一致性洪水頻率計算中，得到較好的效果．混合分布是由k個一致性的分布(子分布)混合而成的，即

式中：Fk( x)為k個子分布的累積分布函數(shù)；αk為各個分布的權(quán)重，且滿足α1+α2+ …+αk=1.0．

根據(jù) SL44—2006《水利水電工程設(shè)計洪水計算規(guī)范》規(guī)定，我國水文頻率計算一般采用 P-Ⅲ型分布，故混合分布模型中的子分布可假設(shè)服從 P-Ⅲ型分布．為降低參數(shù)估計的復(fù)雜性，一般令一個混合分布由 2個子分布組成[5]．因此，對于樣本容量為n的非一致性水文時間序列X，若其變異點為τ，假設(shè)變異點之前的序列為X1，服從概率密度函數(shù)為f1( x)的分布，其樣本長度為n1=τ；變異點之后的序列為X2，服從概率密度函數(shù)為 f2(x)的分布，其樣本長度為整體序列X服從混合分布f( x)，概率密度函數(shù)為

式中α為權(quán)重系數(shù)．

子分布的概率密度函數(shù) f1( x)和 f2(x)均服從 P-Ⅲ型分布，其表達式為

式中：αi、βi和a0i分別為fi( x)函數(shù)分布的形狀、尺度和位置參數(shù)(i=1,2)，可以由統(tǒng)計參數(shù)均值、變差系數(shù)Cvi和偏態(tài)系數(shù)Csi來表示，具體關(guān)系為因此，混合分布f( x)中有和Cs27個需要估計的參數(shù)．

2.2 參數(shù)估計

成靜清等[10]在非一致性年徑流序列頻率參數(shù)計算中，將混合分布的參數(shù)估計問題看成是組合優(yōu)化問題，采用全局優(yōu)化算法——模擬退火算法[13]進行參數(shù)估計．模擬退火算法描述簡單、使用靈活、運行效率高且受初始條件限制較少，本研究也選用此方法進行洪水序列頻率分布參數(shù)估計的計算．

3 實例分析

大清河流域南支上游龍門水庫位于河北省保定市滿城縣境內(nèi)的漕河干流上，控制流域面積470,km2，總庫容 1.27×108,m2，是一座以防洪為主結(jié)合灌溉等綜合利用的大(Ⅱ)水利樞紐工程．水庫始建于 1958年，1959—1960年擴建為大型水庫，1974—1977年進行水庫加固擴建工程，2002年開始實施水庫除險加固工程，2006年水庫防洪標準達到100年設(shè)計，2000年校核．

龍門水庫的洪水資料有洪峰流量(1956—2005年)、年最大 1日洪量(1951—2005年)、年最大 3日洪量(1951—2005年)和年最大 6日洪量(1951—2005年)4個序列．

3.1 洪水序列變異點的診斷

3.1.1 初步診斷

分別采用累積曲線斜率差異幅度分析法和 Hurst系數(shù)法對龍門水庫的 4個洪水序列進行診斷．累積曲線斜率差異幅度分析法，最大差異幅度大于 50%，判斷為變異，反之無變異；Hurst系數(shù)法，Hurst系數(shù)值對應(yīng)變異程度分級表，判斷變異程度；2種方法均變異者即判斷為變異序列．由表 1可知，1日洪量、3日洪量和6日洪量序列存在水文變異，洪峰流量序列不存在水文變異．

表1 龍門水庫洪水序列水文變異初步診斷結(jié)果Tab.1 Initial diagnosis results of hydrological alteration of Longmen reservoir flood series

3.1.2 詳細診斷

由初步診斷結(jié)果，采用 Lee-Heghinian檢驗、有序聚類檢驗、最優(yōu)信息二分割檢驗、R/S檢驗和滑動F檢驗等方法，進行變異點的詳細診斷，結(jié)果見表2．

表2 不同方法下洪水序列水文變異詳細診斷Tab.2 Detailed diagnosis of hydrological alteration of flood series by different methods年

3.1.3 綜合診斷

綜合診斷主要是在詳細診斷可能變異點的基礎(chǔ)上，結(jié)合流域水文自然地理環(huán)境和人類活動變化影響，進行成因分析．

由表2可知，龍門水庫年最大1日洪量序列可能變異點為1964年和1973年，年最大3日洪量序列的可能變異點為1964年、1958年和1972年，年最大6日洪量序列的可能變異點有 1964年、1966年、1970年和1954年．龍門水庫始建于1958年，其固定時段年最大洪量實測水文序列為 1951—2005年，可能變異點1958年和1954年位于水文序列首端，變異點不可靠；龍門水庫于 1963年遭遇了建庫以來最大的一次洪水，稱“63.8”洪水，此場洪水是汛期特大暴雨導(dǎo)致的特大洪水，可能由氣候變化引起，而由于氣候變化原因復(fù)雜，需要長系列水文氣象數(shù)據(jù)加以分析確認，因此從下墊面變化導(dǎo)致洪水變異的角度來說，1964年和 1966年變異點不作為最可能變異點；1974—1977年龍門水庫進行加固擴建工程，且于 20世紀60年代末至70年代，大清河流域?qū)ι嫌未笾行退畮?包括龍門水庫)進行了續(xù)建、擴建，基本形成以水庫、蓄滯洪區(qū)、骨干河道為主的防洪工程體系，流域下墊面條件的改變十分顯著，故診斷出的可能變異點1970年、1972年和1973年是合理可靠的．綜上物理成因分析，龍門水庫年最大 1日、3日和 6日洪量序列的變異點分別為1973年、1972年和1970年．

3.2 設(shè)計洪水成果比較分析

3.2.1 理論頻率計算

1)曲線參數(shù)估計

對變異水文序列，由序列的變異點，以頻率離差絕對值和最小準則(ABS)為目標函數(shù)，采用模擬退火算法對 2個 P-Ⅲ型混合分布的參數(shù)進行估計，結(jié)果詳見表3．

表3 變異序列的參數(shù)估計Tab.3 Parameter estimation of non-stationary series

2)理論頻率計算

將計算得到的混合分布參數(shù)代入式(2)可得整個序列的理論頻率，其計算函數(shù)為

3.2.2 頻率曲線擬合

根據(jù)經(jīng)驗頻率和理論頻率，運用 K-S檢驗法進行擬合檢驗．在 K-S檢驗法中表示隨機樣本的累積頻率函數(shù)，本文中為變異序列的經(jīng)驗頻率．表 4給出了水文序列理論分布頻率與經(jīng)驗分布頻率的 KS檢驗．

表4 龍門水庫3個非一致水文序列K-S檢驗結(jié)果Tab.4 K-S test results of three non-stationary hydrological series of Longmen reservoir

從表4可知，不同固定時段年最大洪量的3個水文變異序列均通過顯著性水平為 0.05的 K-S檢驗，說明樣本序列服從混合分布．

不同時段(1日、3日和6日)年最大洪量水文序列的變異情況理論分布頻率與經(jīng)驗分布頻率擬合曲線和不考慮變異情況的 P-Ⅲ型分布頻率曲線詳見圖1～圖3．圖1～圖3分別為年最大1日洪量、年最大3日洪量和年最大 6日洪量序列的頻率曲線．由圖1～圖 3可見，在離差絕對值和最小準則的情況下，相比不考慮變異情況下的傳統(tǒng) P-Ⅲ型分布頻率曲線，3組變異序列混合分布曲線與經(jīng)驗點據(jù)擬合明顯較好，尤其是經(jīng)驗點距頻率較大部分．這說明2個P-Ⅲ型分布混合模型應(yīng)用于非一致性洪水序列是可靠的，特別是存在水文序列變異的情況下，可以與經(jīng)驗點據(jù)整體擬合更好．

圖1 龍門水庫年最大1日洪量序列的頻率曲線Fig.1 Frequency curves of annual maximum 1-day flood volume series of Longmen reservoir

圖2 龍門水庫年最大3日洪量序列的頻率曲線Fig.2 Frequency curves of annual maximum 3-day flood volume series of Longmen reservoir

圖3 龍門水庫年最大6日洪量序列的頻率曲線Fig.3 Frequency curves of annual maximum 6-day flood volume series of Longmen reservoir

3.2.3 設(shè)計洪水成果比較分析

對于考慮變異情況下龍門水庫設(shè)計洪水的合理性，將通過與 1985年水利部規(guī)劃總院審定的龍門水庫設(shè)計洪水值，及不考慮變異情況下的 P-Ⅲ型分布給出的設(shè)計洪水值進行對比分析．表 5為龍門水庫相應(yīng)的不同設(shè)計洪水值及變化情況下的對比．

由表 5可見，診斷無變異的洪峰流量序列(1956—2005年)，考慮歷史特大洪水情況下，通過P-Ⅲ型分布頻率適線，得到的不同標準設(shè)計洪水值相較1985年審定的設(shè)計洪水值均小，變小的比例在7.7%～30.0%；而固定時段(1日、3日和6日)年最大洪量序列均診斷出變異，且3個變異序列的變異點均在1970年附近，運用混合分布模型計算出不同設(shè)計標準的設(shè)計洪水值相比不考慮變異情況下的設(shè)計洪水值和 1985年審定的設(shè)計洪水值，也均有不同程度的減小，相較不考慮變異情況下的設(shè)計洪水值，年最大1日洪量值變化比例范圍為 1.0%～20.8%；年最大 3日洪量變化比例范圍在1.2%～23.7%；年最大6日洪量變化比例范圍在 1.8%～22.1%，變異的程度均比較明顯．從統(tǒng)計角度分析，現(xiàn)狀洪水序列情況下得到的設(shè)計洪峰值和設(shè)計洪量值較1985年審定的設(shè)計洪峰值和設(shè)計洪量值小，前者可能由于樣本長度增長所致，1985年審定的設(shè)計值需完善修訂；得到的非一致性序列設(shè)計洪量值較不考慮變異情況的設(shè)計洪量值小，表明流域受下墊面變化和人類活動(水保措施和水利工程的興建等)的影響，導(dǎo)致入庫洪水減少的現(xiàn)象比較明顯．這在一定程度上說明，在下墊面變化和人類活動影響比較大的大清河流域，進行設(shè)計洪水的校核與修訂是非常必要的．

表5 龍門水庫設(shè)計洪水成果比較Tab.5 Comparison of design flood values of Longmen reservoir

4 結(jié) 論

以大清河流域南支上游龍門水庫為對象，著重研究了混合分布模型解決非一致性洪水序列頻率計算的問題，并進一步分析非一致性水文序列后的設(shè)計洪水值與未變異前的設(shè)計洪水值之間的變化，得出以下主要結(jié)論：

(1) 采用變異點綜合診斷方法和物理成因分析，可以較準確地判斷洪水序列是否變異以及得到變異序列的變異點，診斷龍門水庫固定時段(1日、3日、6日)年最大洪量序列存在變異，且 3個序列變異點年份相近．

(2) 采用模擬退火算法，以離差絕對值和最小準則為目標函數(shù)，可以比較可靠地進行參數(shù)估計，避免了人工適線的誤差，實例證明得到的混合分布曲線與經(jīng)驗點據(jù)擬合較好，比較適合存在變異情況下的水文序列頻率分析．

(3) 基于混合分布的非一致性洪水頻率分析計算，不需要進行洪水資料的還原，完全從統(tǒng)計角度完成洪水頻率分析計算，理論基礎(chǔ)比較明確．

(4) 對于龍門水庫的非一致性洪水序列，與不考慮變異情況的設(shè)計洪水值相比，不同設(shè)計標準的設(shè)計洪水值相比不考慮變異情況下的設(shè)計洪水值，也均有不同程度的減小，變化比例范圍均在 1%～25%，變異的程度比較明顯．因此，在下墊面變化和人類活動影響比較大的大清河流域，進行設(shè)計洪水的校核與修訂是非常必要的．

［1］ 梁忠民，胡義明，王 軍. 非一致性水文頻率分析的研究進展[J]. 水科學(xué)進展，2011，22(6)：864-871.Liang Zhongmin，Hu Yiming，Wang Jun. Advances in hydrological frequency analysis of non-stationary time series[J].Advances in Water Science，2011，22(6)：864-871(in Chinese).

［2］ Strupczewski W G，Singh V P，F(xiàn)eluch W. Nonstationary approach to at-site flood frequency modeling(I)：Maximum likelihood estimation[J].Journal of Hydrology，2001，248(1/2/3/4)：123-142.

［3］ 宋松柏，李 揚，蔡明科. 具有跳躍變異的非一致分布水文序列頻率計算方法[J]. 水利學(xué)報，2012，43(6)：734-739，748.Song Songbai，Li Yang，Cai Mingke. Methods of frequency analysis for hydrologic data with jump up components[J].Journal of Hydraulic Engineering，2012，43(6)：734-739，748(in Chinese).

［4］ Singh K P，Sinclair R A. Two-distribution method for flood frequency analysis[J].Journal of Hydraulics Division，1972，98(1)：29-44.

［5］ Alila Y，Mtiraoui A. Implications of heterogeneous flood-frequency distributions on traditional streamdischarge prediction techniques[J].Hydrological Processes，2002，16(5)：1065-1084.

［6］ Singh K P. A versatile flood frequency methodology[J].Water International，1987，12(3)：139-145.

［7］ Rossi F，F(xiàn)iorentino M，Versace P. Two-component extreme value distribution for flood frequency analysis[J].Water Resources Research，1984，20(7)：847-856.

［8］ Leytham K M. Maximum likelihood estimates for the parameters of mixed distributions[J].Water Resources Research，1984，20(7)：896-902.

［9］ Fiorentino M，Arora K，Singh V P. The twocomponent extreme value distribution for flood frequency analysis：Derivation of a new estimation method[J].Stochastic Hydrology and Hydraulics，1987(1)：199-208.

［10］ 成靜清，宋松柏. 基于混合分布非一致性年徑流序列頻率參數(shù)的計算[J]. 西北農(nóng)林科技大學(xué)學(xué)報：自然科學(xué)版，2010，38(2)：229-234.Cheng Jingqing，Song Songbai. Calculation of hydrological frequency parameters of inconsistent annual runoff series based on mixed distribution[J].Journal of Northwest A&F University：Natural Science Edition，2010，38(2)：229-234(in Chinese).

［11］ 謝 平，陳廣才，夏 軍. 變化環(huán)境下非一致性年徑流序列的水文頻率計算原理[J]. 武漢大學(xué)學(xué)報：工學(xué)版，2005，38(6)：6-9，15.Xie Ping，Chen Guangcai，Xia Jun. Hydrological frequency calculation principle of inconsistent annual runoff series under changing environments[J].Engineering Journal of Wuhan University，2005，38(6)：6-9，15(in Chinese).

［12］ 謝 平，陳廣才，雷紅富，等. 水文變異診斷系統(tǒng)[J]. 水力發(fā)電學(xué)報，2010，29(1)：85-91.Xie Ping，Chen Guangcai，Lei Hongfu，et al. Hydrological alteration diagnosis system[J].Journal of Hydroelectric Engineering，2010，29(1)：85-91(in Chinese).

［13］ 謝 云. 模擬退火算法的原理及實現(xiàn)[J]. 高等學(xué)校計算數(shù)學(xué)學(xué)報，1999(3)：212-218.Xie Yun. Principle and realization of the simulated annealing algorithm[J].Numerical Mathematics A Journal of Chinese University，1999(3)：212-218(in Chinese).