馮 平，曾 杭，李 新

(天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

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

1 非一致性洪水序列變異點(diǎn)的診斷方法

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

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

2.1 混合分布模型

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

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

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

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

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

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

2.2 參數(shù)估計(jì)

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

3 實(shí)例分析

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

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

3.1 洪水序列變異點(diǎn)的診斷

3.1.1 初步診斷

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

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

3.1.2 詳細(xì)診斷

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

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

3.1.3 綜合診斷

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

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

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

3.2.1 理論頻率計(jì)算

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

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

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

2)理論頻率計(jì)算

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

3.2.2 頻率曲線擬合

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

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

從表4可知，不同固定時(shí)段年最大洪量的3個(gè)水文變異序列均通過(guò)顯著性水平為 0.05的 K-S檢驗(yàn)，說(shuō)明樣本序列服從混合分布．

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

［1］ 梁忠民，胡義明，王 軍. 非一致性水文頻率分析的研究進(jìn)展[J]. 水科學(xué)進(jìn)展，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］ 宋松柏，李 揚(yáng)，蔡明科. 具有跳躍變異的非一致分布水文序列頻率計(jì)算方法[J]. 水利學(xué)報(bào)，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ì)算[J]. 西北農(nóng)林科技大學(xué)學(xué)報(bào)：自然科學(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ì)算原理[J]. 武漢大學(xué)學(xué)報(bào)：工學(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é)報(bào)，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］ 謝 云. 模擬退火算法的原理及實(shí)現(xiàn)[J]. 高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，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).