考慮非平穩(wěn)性特征的雅礱江流域洪水量化分析

高 潔，楊 龍

(1.水電水利規(guī)劃設(shè)計總院，北京 100120；2.南京大學(xué)地理與海洋科學(xué)學(xué)院，江蘇 南京 210023)

0 引 言

設(shè)計洪水是確定水電工程規(guī)模和制定防洪標(biāo)準(zhǔn)的重要基礎(chǔ)。基于頻率分析的設(shè)計洪水計算在分布線型、參數(shù)估計等方面積累了豐富的研究成果[1]。氣候變化和人類活動對水文時間序列的平穩(wěn)性特征產(chǎn)生影響，傳統(tǒng)水文設(shè)計方法面臨變化環(huán)境帶來的挑戰(zhàn)[2]。郭生練等[1]將非一致性洪水頻率分析、設(shè)計洪水不確定性研究、基于水文物理機(jī)制的洪水頻率分析等內(nèi)容列為當(dāng)前設(shè)計洪水研究前沿與熱點。魯帆等[3]在水文氣象領(lǐng)域非平穩(wěn)時間序列極值模型研究綜述中介紹，傳統(tǒng)水文頻率分析通常采用皮爾遜III型(P-III)，對數(shù)皮爾遜III型，廣義極值分布(GEV)，廣義帕累托分布(GPD)等理論線型。我國工程界設(shè)計洪水標(biāo)準(zhǔn)主要采用P-III分布。以年最大為樣本的GEV分布也在水文極值統(tǒng)計中逐漸被采用。關(guān)于GEV分布參數(shù)，GEV-CDN模型[4]認(rèn)為形狀參數(shù)決定了分布的尾部形態(tài)，具有較大的不確定性，可將形狀參數(shù)取值設(shè)定在[-0.5，0.5]。GEV模型應(yīng)用于洞庭湖區(qū)、岷江流域年最大日降水系列[5- 6]，發(fā)現(xiàn)各地的形狀參數(shù)差異較大，形狀參數(shù)為0的極值I型Gumbel分布、形狀參數(shù)<0的極值III型Weibull分布，形狀參數(shù)>0的極值II型Frechet分布均有適用。通過對全國降雨資料統(tǒng)計，Yang等[7]發(fā)現(xiàn)極端降雨在長江流域呈現(xiàn)形狀參數(shù)<0的薄尾分布；在環(huán)渤海地區(qū)形狀參數(shù)>0，呈現(xiàn)厚尾分布。長江中下游控制站大通水文站[8]1950年至2012年歷年最大和非汛期最大流量序列符合形狀參數(shù)<0的Weibull分布，新疆天山開都河[9]1957年至2011年歷年、春季和夏季最大日流量序列服從形狀參數(shù)≥0的Gumbel分布和Frechet分布。通過對珠江流域北江水系武江上游坪石水文站和下游犁市(二)水文站實測年最大徑流系列分別采用薄尾P-III分布、Weibull分布、厚尾廣義邏輯分布(GLO)、三參數(shù)混合尾GEV、Gumbel、兩參數(shù)對數(shù)正態(tài)分布(LN2)[10]擬合發(fā)現(xiàn)，厚尾GLO和三參數(shù)混合尾GEV分布在擬合時可使高水端上翹，提高了高水端擬合優(yōu)度。

非平穩(wěn)時間序列的概率分布均值、方差或協(xié)方差隨時間變化，重現(xiàn)期和風(fēng)險度有別于傳統(tǒng)水文頻率計算結(jié)果，需要對變化序列進(jìn)行還原/還現(xiàn)，或針對時變序列進(jìn)行頻率計算。謝平等[11]按照“分離”“合成”的思路對水文序列非一致性進(jìn)行診斷?？紤]時間因素對傳統(tǒng)頻率分布的影響，劉德地等[12]推導(dǎo)了非平穩(wěn)條件下洪水重現(xiàn)期的計算公式。顧西輝等[13]選取珠江流域28個水文站1951年至2010年歷年最大流量數(shù)據(jù)，在突變點分析的基礎(chǔ)上采用GEV、GPD、Gumbel、P-III、GLO、廣義對數(shù)正態(tài)分布(GNO)6種分布進(jìn)行擬合。研究發(fā)現(xiàn)，最優(yōu)分布以GEV居多。珠江流域內(nèi)西江、北江頻率設(shè)計值增大，且隨著重現(xiàn)期增加而增大，基于整體樣本的水利水電工程規(guī)劃設(shè)計，可能低估大洪水影響，導(dǎo)致防洪風(fēng)險。顧西輝等[14]又將GEV-CDN模型[4]應(yīng)用于珠江流域6個水文站近45年歷年最大日流量序列。研究發(fā)現(xiàn)，50年一遇以下洪水，平穩(wěn)和非平穩(wěn)模型成果差異不大；但50年一遇以上洪水，不同模型頻率設(shè)計值差異較大，存在潛在的防洪風(fēng)險。梁忠民等[15]針對非平穩(wěn)環(huán)境下工程水文設(shè)計值發(fā)生變化的問題，以等可靠度建立平穩(wěn)序列與非平穩(wěn)序列的聯(lián)系，提出了待建工程、已建工程推求和協(xié)調(diào)頻率設(shè)計值的思路。

考慮位置參數(shù)、尺度參數(shù)、形狀參數(shù)的廣義可變模型GAMLSS為變化環(huán)境水文頻率計算提供了直接可靠的流程和方法。顧西輝等[16-17]將GAMLSS模型應(yīng)用于全國839個氣象站極端降水發(fā)生率及新疆塔里木河洪水發(fā)生率的非平穩(wěn)性研究。高潔[18-19]將GAMLSS模型應(yīng)用于雅礱江流域。研究發(fā)現(xiàn)，全流域降水特性無顯著變化，但是流域中下游極端降水量增加，并伴隨著最長連續(xù)干期變長、連續(xù)濕期變短。流域中上游極端降水平穩(wěn)性較好，下游不同站點極端降水非平穩(wěn)性特征差異較大。小流域及支流洪水的平穩(wěn)性特征易受局地極端降水影響。雖然有建壩修庫等人類活動，但隨著支流入?yún)R，流域干流水文站洪峰流量回復(fù)平穩(wěn)狀態(tài)。

關(guān)于水文時間序列，基于平穩(wěn)假設(shè)和考慮非平穩(wěn)性，均開展了大量的研究工作。根據(jù)傳統(tǒng)規(guī)范要求，我國水電水利工程行業(yè)，水文設(shè)計選擇P-III分布作為頻率計算的基礎(chǔ)?？紤]極值理論和變化環(huán)境，傳統(tǒng)頻率計算成果與GEV分布、GAMLSS模型成果的差異尚缺乏定量分析。本文以雅礱江流域為研究對象，通過GEV統(tǒng)計參數(shù)和GAMLSS最優(yōu)模擬，量化分析GEV分布、GAMLSS模型成果與傳統(tǒng)P-III適線成果差異。雅礱江干流河道天然落差約3 870 m，水能資源豐富，是我國規(guī)劃建設(shè)的十四大水電基地之一。根據(jù)河流水電規(guī)劃，雅礱江流域水電技術(shù)可開發(fā)量約2 900萬kW，主要集中在流域中下游河段，分布著眾多水電水利工程。因此，充分認(rèn)識流域洪水時變特征、洪水頻率計算成果偏差，對于流域內(nèi)眾多涉水工程設(shè)計運行具有重要意義。

1 數(shù)據(jù)和方法

1.1 研究對象

雅礱江流域(96°52′E～102°48′E，26°32′N～33°58′N)是長江上游金沙江最大的一級支流，地處青藏高原東南部，位于金沙江和大渡河之間的狹長地帶。流域面積約13萬km2，谷嶺高差大[20]，地形地貌以高山峽谷為主，基于90 m空間分辨率的DEM數(shù)據(jù)，海拔高程位于966～5 888 m[20]，地勢西北高東南低。本研究收集了流域內(nèi)8個水文站：甘孜、雅江、瀘寧、道孚、朱巴、濯桑、莊房(二)、孫水關(guān)長系列年最大洪峰流量系列。

1.2 研究方法

1.2.1洪水系列趨勢檢驗和突變分析

針對洪水系列的單調(diào)變化和突變特征，分別進(jìn)行趨勢檢驗和突變點分析。

MK檢驗作為廣泛應(yīng)用于水文氣象時間序列趨勢分析的非參數(shù)化方法，得到了國際氣象組織的推薦[21]，被廣泛應(yīng)用。在突變點分析方面，分別采用Pettitt、AMOC和ECP三種檢驗方法。Pettitt法通過檢驗時間序列均值變化，確定序列躍變發(fā)生時間[22]。基于R語言平臺，changepoint程序包[23]的AMOC檢驗法是基于似然函數(shù)框架的突變點檢測，時間序列無需滿足正態(tài)分布假設(shè)。ECP程序包[24]是基于凝聚和分裂算法的聚類分析法，可用于時間序列的多突變點檢驗。

1.2.2洪水系列非平穩(wěn)性分析

Stasinopoulos和Rigby于2005年提出(半)參數(shù)回歸模型GAMLSS(基于位置參數(shù)、尺度參數(shù)和形狀參數(shù)的廣義可加模型)，該模型[25]可模擬統(tǒng)計參數(shù)隨時間的變化模式。本文為了簡化模型結(jié)構(gòu)，采用兩參數(shù)分布函數(shù)，基于對數(shù)正態(tài)(LOGNO)、Gamma(GA)、Gumbel(GU)和Weibull(WEI)分布，擬定位置參數(shù)、尺度參數(shù)與時間的關(guān)系，通過AIC準(zhǔn)則選擇最優(yōu)分布函數(shù)，計算水文序列的時變頻率設(shè)計值。

1.2.3洪水系列適線

我國傳統(tǒng)設(shè)計洪水頻率適線通常采用P-III型曲線，極值統(tǒng)計近年來也得到了廣泛關(guān)注[26]。GEV作為極值統(tǒng)計的一個重要分類，包含了Gumbel、Frechet和Weibull三種分布類型。

H(x;μ，σ，ξ)=exp{-[1+ξ(x-μ)/σ]-1/ξ}，[1+ξ(x-μ)]/σ>0

(1)

式中，μ為位置參數(shù)；σ為尺度參數(shù)；ξ為形狀參數(shù)。μ和σ反映洪水的量級和變幅。ξ決定了分布函數(shù)的尾部形態(tài)[7]。ξ>0的極值II型(Frechet)呈現(xiàn)厚尾分布，尾部較長且無上邊界。ξ<0極值的III型分布(Weibull)具有上端點μ-σ/ξ，呈現(xiàn)薄尾分布。當(dāng)統(tǒng)計變量x趨近于無窮，ξ=0的極值I型分布(Gumbel)為無上下邊界的薄尾分布。

表1 各水文站趨勢檢驗和突變點分析

2 成果分析

2.1 時間序列特性

根據(jù)MK趨勢檢驗，各水文站的年最大洪峰流量序列均無顯著變化趨勢(顯著度水平5%)。由于收集的洪峰流量系列長度僅在37年至65年，本研究檢測一個突變點，并根據(jù)變點位置將序列劃分為前后兩部分。經(jīng)Pettitt、AMOC(mean)和ECP突變檢驗，甘孜、濯桑、莊房(二)、孫水關(guān)4站通過顯著水平5%突變檢驗，雅江站滿足顯著水平10%突變檢驗。變化趨勢和突變時間點如表1所示。

2.2 平穩(wěn)性分析

對各水文站年最大洪峰流量系列進(jìn)行GAMLSS非平穩(wěn)性分析，在LOGNO、GA、GU和WEI四種分布中，以LOGNO分布占優(yōu)。雖然各水文站年最大洪峰流量的MK趨勢檢驗無顯著變化，但是與甘孜、濯桑、莊房、孫水關(guān)、雅江5站通過突變檢驗相協(xié)調(diào)。GAMLSS模型顯示，其中的甘孜、濯桑和孫水關(guān)3站的年最大洪峰流量呈現(xiàn)非平穩(wěn)特征，如圖1所示。

(1)甘孜站位置參數(shù)表征的年最大洪峰流量均值隨時間單調(diào)減小，MK檢驗呈不顯著減小趨勢，1986年以后均值突變減小。

(2)濯桑站尺度參數(shù)表征的不同重現(xiàn)期年最大流量離散程度逐漸減小，體現(xiàn)為常遇頻率設(shè)計值增加、稀遇頻率設(shè)計值減小，并逐漸穩(wěn)定。MK檢驗常遇頻率設(shè)計值呈不顯著增加。

(3)孫水關(guān)站位置參數(shù)表征的年最大洪峰流量均值在1980年～1990年出現(xiàn)拐點，先減小后增大。MK檢驗常遇頻率設(shè)計值呈不顯著增加，在1985年以后突變增加。

受限于觀測資料多為常遇頻率，MK趨勢檢驗和突變分析主要反映的是常遇頻率下洪水事件。GAMLSS成果可以體現(xiàn)不同重現(xiàn)期表征的常遇頻率和稀遇頻率設(shè)計值差異化的時變特征。年最大洪峰流量的非平穩(wěn)性分析(見表2)，基本映射出漸變性和突變性檢驗中尚不顯著的發(fā)展趨勢，具有一定的前瞻性。

時間序列趨勢檢驗無顯著變化，采用傳統(tǒng)的P-III適線和穩(wěn)態(tài)GEV分布擬合仍然合適。基于GAMLSS非平穩(wěn)性分析和考慮突變點的GEV分布，可反映變化環(huán)境下洪峰流量頻率設(shè)計值隨時間的變化。

表2 各水文站GAMLSS非平穩(wěn)性分析成果

2.3 頻率適線成果對比

P-III分布目估適線，通常微調(diào)Cv和Cs兩個參數(shù)，控制適線的傾斜程度和彎曲程度，使高水端上翹，傾向于頭部特大值，使稀遇頻率設(shè)計偏安全，得到比較保守的頻率設(shè)計值。GAMLSS兩參數(shù)分布由于參數(shù)較少，靈活性低，不能人為提高高水端的擬合優(yōu)度，可能出現(xiàn)低頻稀遇設(shè)計洪水洪峰流量小于實測頭部特大值，該方法調(diào)整裕度較少。GEV分布通過擬合位置參數(shù)、尺度參數(shù)和形狀參數(shù)計算頻率設(shè)計值。位置參數(shù)和尺度參數(shù)均與集水面積具有相關(guān)性，關(guān)于形狀參數(shù)的影響因素，學(xué)術(shù)界尚無定論。但是，形狀參數(shù)又對于GEV尾部形態(tài)具有關(guān)鍵作用，以正負(fù)為界，區(qū)別厚尾和薄尾特征，體現(xiàn)了分布極值的偏離程度和密集程度，對于研究極端事件具有重要意義。

2.3.1平穩(wěn)序列適線成果對比

年最大洪峰流量序列相對平穩(wěn)的雅江、瀘寧、朱巴、道孚、莊房5站經(jīng)三參數(shù)P-III目估適線、GEV分布擬合及兩參數(shù)GAMLSS模型優(yōu)選，除了道孚站在低頻稀遇洪水時GEV分布推算頻率設(shè)計值偏大，其他4站均為GEV分布和GAMLSS模型優(yōu)選成果較為接近，普遍小于P-III適線成果，且部分站點(如瀘寧、莊房)隨著重現(xiàn)期增大，稀遇洪水設(shè)計值的偏差增加。三種方法成果對比見圖2。

圖2 平穩(wěn)序列GAMLSS、P-III、GEV 3種方法計算成果比較

根據(jù)平穩(wěn)序列5個水文站的擬合成果，三種方法對于20年一遇以下洪水洪峰流量的擬合偏差低于10%，50年一遇至200年一遇洪水洪峰流量的擬合偏差≤15%，1 000年一遇以下洪水洪峰流量的擬合誤差基本不超過20%。對于5 000年、10 000年等稀遇特大洪水的擬合偏差在25%至30%之間(見圖2f)。

2.3.2非平穩(wěn)序列適線成果對比

通過GAMLSS模型推算年最大洪峰流量時間序列非平穩(wěn)的甘孜、濯桑和孫水關(guān)3站頻率設(shè)計值的時變特性。隨著重現(xiàn)期增大，頻率設(shè)計值變幅增大，體現(xiàn)了采用有限資料推算稀遇事件的外延性導(dǎo)致不確定性增加(見圖3a、3b、3c)，可通過頻率設(shè)計值時變范圍與設(shè)計值的比例來反映。

甘孜站1956年～2016年洪峰流量頻率設(shè)計值時變范圍與均值比例，對于不同重現(xiàn)期基本在14.8%～15.5%。濯桑站1966年～2016年洪峰流量頻率設(shè)計值隨著重新期增大，時變范圍比例從10年一遇洪水23%增至10 000年一遇洪水的57%。孫水關(guān)站時變范圍比例從10年一遇洪水45%降至10 000年一遇洪水的39%。

濯桑站和孫水關(guān)站由于集水面積較小，實測年最大洪峰流量量級相對較小(均值分別為287 m3/s和610 m3/s)；因此，頻率設(shè)計值時變范圍比例較大(23%～57%)。GAMLSS頻率設(shè)計值的時變特征是否具有尺度效應(yīng)，小流域非平穩(wěn)性普遍較大，有待更多站點驗證。

隨著重現(xiàn)期增大，甘孜站和孫水關(guān)站洪峰流量時變范圍與均值的比例基本穩(wěn)定。濯桑站洪峰流量時變范圍比例明顯增大，主要由于低頻稀遇洪水時變減小趨勢(見圖1b、3b)，體現(xiàn)為10年一遇以上洪水隨著重現(xiàn)期增大增幅減小，稀遇洪水致災(zāi)風(fēng)險隨時間降低。

P-III目估適線成果均在GAMSS測算區(qū)間內(nèi)。但是，GEV分布在低頻大值的估算上與P-III適線和GAMLSS模型成果(見圖3a和3c)產(chǎn)生偏差。GEV分布對于非平穩(wěn)系列的擬合較為敏感，需要考慮參數(shù)時變特性和尾部特點。

圖3 非平穩(wěn)序列GAMLSS、P-III、GEV三種方法計算成果比較

表3 各水文站突變前后及完整系列GEV參數(shù)

2.4 考慮突變的GEV變化分析

2.4.1GEV參數(shù)分析

Yang[7]發(fā)現(xiàn)中國極端降水GEV分布的位置參數(shù)、尺度參數(shù)與集水面積相關(guān)，長江流域極端降水的GEV形狀參數(shù)以負(fù)值居多。與該結(jié)論一致，流域內(nèi)8個水文站原始系列、突變前后系列的位置參數(shù)和尺度參數(shù)均與集水面積(見表3)呈現(xiàn)雙對數(shù)相關(guān)(見圖4)。本流域位于長江流域上游支流，年最大洪峰流量GEV分布也以形狀參數(shù)<0的薄尾分布為主。僅道孚站和孫水關(guān)站的形狀參數(shù)>0，呈現(xiàn)厚尾分布，由此造成兩站GEV擬合頻率設(shè)計值比P-III分布和GAMLSS上限值均偏大(見圖2c、3c)。GEV對于極大值的過分高估，是否由于形狀參數(shù)>0所致，值得研究。

圖4 GEV統(tǒng)計參數(shù)與集水面積雙對數(shù)線性相關(guān)

2.4.2突變前后GEV分析

選取形狀參數(shù)<0的甘孜、雅江、濯桑、莊房4個站點，基于突變點分析，對比完整系列、突變前、突變后3種系列的GEV位置參數(shù)、尺度參數(shù)和形狀參數(shù)，推算頻率設(shè)計值。

(1)甘孜站和濯桑站分別在1986年突變后和1989年突變前，系列的GEV形狀參數(shù)>0，由此導(dǎo)致GEV擬合值偏離較大，分析無效。

(2)雅江站在突變點1998年以后，各頻率設(shè)計洪水均有增加(見圖5a)。因此1998年以前設(shè)計的工程有必要復(fù)核現(xiàn)狀防洪安全。

(3)莊房站在2006年以前，各頻率設(shè)計洪水略大于現(xiàn)狀水平，2006年以前的設(shè)計仍可滿足現(xiàn)狀相應(yīng)的防洪標(biāo)準(zhǔn)(見圖5b)。

圖5 基于完整系列、突變前后系列推算設(shè)計洪水頻率值

3 結(jié) 論

基于MK趨勢檢驗，Pettitt、AMOC和ECP突變分析，研究雅礱江流域8個水文站年最大洪峰流量漸變和突變特征，基于GAMLSS模型分析流域洪水的平穩(wěn)性特征。按照平穩(wěn)序列和非平穩(wěn)序列，對各水文站年最大洪峰流量系列進(jìn)行P-III目估適線、GEV分布擬合和GAMLSS模型成果對比分析，探索GEV分布位置參數(shù)、尺度參數(shù)和形狀參數(shù)特點，結(jié)合突變點分析，研究流域防洪安全。研究發(fā)現(xiàn)：

(1)GAMLSS擬合以對數(shù)正態(tài)最優(yōu)分布居多。其中，甘孜站、濯桑站和孫水關(guān)站年最大洪峰流量時間序列非平穩(wěn)。GAMLSS模型成果與實測系列漸變、突變分析統(tǒng)籌考慮，體現(xiàn)了不同重現(xiàn)期設(shè)計值時變特征，揭示了實測點難以反映的低頻洪水特性。

(2)比較了3參數(shù)P-III目估適線、GEV分布和兩參數(shù)GAMLSS模型成果，針對平穩(wěn)序列5個站點，傾向于頭部特大值的P-III目估適線成果普遍大于以穿過點據(jù)中央為擬合原則的GEV分布和GAMLSS模擬成果。3種途徑推算20年一遇頻率洪水偏差<10%，200年一遇頻率洪水偏差<15%，1 000年一遇頻率洪水偏差基本不超過20%。

(3)非平穩(wěn)序列的3個站點，P-III適線成果基本在GAMLSS模型擬合時變范圍內(nèi)。GEV擬合非平穩(wěn)序列建議考慮參數(shù)的時變特性和尾部特點。

(4)雅礱江流域年最大洪峰流量序列GEV分布形狀參數(shù)以<0的薄尾分布為主。呈現(xiàn)厚尾分布的道孚站和孫水關(guān)站GEV頻率設(shè)計值遠(yuǎn)大于其他模型成果，存在極大值被高估的現(xiàn)象。雅江站1998年以前設(shè)計的工程有必要復(fù)核現(xiàn)狀防洪安全。