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

高 潔，楊 龍

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

0 引 言

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

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

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

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

1 數據和方法

1.1 研究對象

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

1.2 研究方法

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

針對洪水系列的單調變化和突變特征，分別進行趨勢檢驗和突變點分析。

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

1.2.2洪水系列非平穩性分析

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

1.2.3洪水系列適線

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

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

(1)

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

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

2 成果分析

2.1 時間序列特性

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

2.2 平穩性分析

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

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

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

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

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

時間序列趨勢檢驗無顯著變化，采用傳統的P-III適線和穩態GEV分布擬合仍然合適?；贕AMLSS非平穩性分析和考慮突變點的GEV分布，可反映變化環境下洪峰流量頻率設計值隨時間的變化。

表2 各水文站GAMLSS非平穩性分析成果

2.3 頻率適線成果對比

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

2.3.1平穩序列適線成果對比

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

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

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

2.3.2非平穩序列適線成果對比

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

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

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

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

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

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

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

2.4 考慮突變的GEV變化分析

2.4.1GEV參數分析

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

圖4 GEV統計參數與集水面積雙對數線性相關

2.4.2突變前后GEV分析

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

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

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

(3)莊房站在2006年以前，各頻率設計洪水略大于現狀水平，2006年以前的設計仍可滿足現狀相應的防洪標準(見圖5b)。

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

3 結 論

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

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

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

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

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