魯帆 江明 蔣云鐘 周毓彥 徐揚

摘要：海河流域天然河川徑流持續(xù)衰減，水文豐枯情勢顯著變化，亟需研究適用于非一致性水文序列的豐枯概率計算方法?；跇?biāo)準(zhǔn)化徑流指數(shù)、GAMLSS模型等方法，提出一種不同等級豐枯水事件期望發(fā)生次數(shù)和期望等待時間的計算方法，研究變化環(huán)境下海河流域天然河川徑流豐枯概率的演變規(guī)律。結(jié)果表明：① 徑流豐枯概率呈現(xiàn)出顯著的枯增豐減趨勢;② 同傳統(tǒng)的一致性分布等多類概率分布相比，以時間t為協(xié)變量的LOGNO分布擬合流域徑流系列的效果最優(yōu)，且基于該分布計算的期望發(fā)生次數(shù)更接近于歷史實際;③ 非一致性最優(yōu)模型不同情景條件下計算的流域極端枯水和極端豐水事件的期望等待時間分別為4.9～9.4 a、14.5～36.0 a，說明海河流域近期發(fā)生極端枯水的概率遠(yuǎn)大于極端豐水。

關(guān)鍵詞：河川徑流;豐枯概率;非一致性;GAMLSS模型;海河流域

中圖分類號：P333.1

文獻(xiàn)標(biāo)志碼：A

文章編號：1001-6791（2023）01-0012-09

收稿日期：2022-09-08;

網(wǎng)絡(luò)出版日期：2023-02-01

網(wǎng)絡(luò)出版地址：https：∥kns.cnki.net/kcms/detail∥32.1309.P.20230131.1539.006.html

基金項目：國家重點研發(fā)計劃資助項目（2021YFC3000200;2018YFC0406500）

作者簡介：魯帆（1981—），男，湖北天門人，正高級工程師，主要從事水文水資源研究。E-mail：lufan@iwhr.com

通信作者：江明，E-mail：1541427645@qq.com

河川徑流對氣候、下墊面等環(huán)境變化的響應(yīng)敏感^［1-2］，部分流域由于環(huán)境變化破壞了徑流序列的一致性，基于傳統(tǒng)一致性假設(shè)的水文頻率分析方法已不再適用^［3-4］，分析徑流豐枯變化趨勢與成因、深入研究環(huán)境變化對極值豐枯水事件時間間隔和發(fā)生概率的影響，對于流域水資源安全保障具有重要的實際意義。

已有研究常基于距平百分率、Z指數(shù)、標(biāo)準(zhǔn)化指數(shù)、概率矩陣等方法計算豐枯變化閾值和評價等級^［5-7］。傳統(tǒng)水文頻率分析則采用矩法、三點法、權(quán)函數(shù)、極大似然等方法估計分布參數(shù)，計算水文極值的發(fā)生概率和重現(xiàn)期。當(dāng)變量總體分布不便確定時，可采用核密度估計或基于徑流形成因子推求分布概率^［8］。許多研究考慮環(huán)境變化影響，引入時間、降水等協(xié)變量，構(gòu)建變參數(shù)概率分布模型分析水文頻率^［9-11］，其中GAMLSS（Generalized Additive Model for Location，Scale and Shape）模型研究較多^［12-13］。非一致性條件下水文系列服從的分布參數(shù)每年不再是常量，國內(nèi)外學(xué)者提出利用期望超過次數(shù)、期望等待時間、設(shè)計水平年限和等可靠度等方法推求變參數(shù)下指定設(shè)計標(biāo)準(zhǔn)（如百年一遇）的水文設(shè)計值^［14-16］，但上述方法的推演計算存在不確定性，應(yīng)考慮水文序列的變化趨勢和驅(qū)動機(jī)制，例如期望等待時間法在收斂前的預(yù)估時間越長，不確定性將越大^［17-18］ 。

近年來，海河流域徑流呈現(xiàn)出顯著的衰減趨勢，對比2001—2021年與1956—2000年2個時期，流域后一時期多年平均降水量比前一時期減幅為1.9%，天然河川徑流量減幅高達(dá)31.8%，已有研究認(rèn)為社會經(jīng)濟(jì)取用水增加、下墊面變化、地下水超采等人類活動是徑流減少的主要原因^［19-20］，洪水與徑流的設(shè)計值均發(fā)生了較大變化^［9］。

本文提出一套揭示變化環(huán)境下徑流豐枯概率演變規(guī)律的計算方法，分析非一致性條件下海河流域天然河川徑流豐枯變化及極值豐枯水事件發(fā)生次數(shù)、發(fā)生時間間隔的演變特征，并比較不同情景的計算結(jié)果，以期為未來水資源評價和旱澇風(fēng)險管理提供科學(xué)依據(jù)。

1 研究區(qū)概況

海河流域位于112°E—120°E、35°N—43°N之間，總面積約32萬km²，包括灤河水系、海河北系、海河南系和徒駭馬頰河水系4個水資源二級區(qū)，見圖1。流域西部為山西高原和太行山區(qū)，北部為內(nèi)蒙古高原和燕山山區(qū)，東部和東南部為廣闊平原，山地和高原約占流域總面積的60%。流域?qū)儆跍貛О霛駶?、半干旱大陸性季風(fēng)氣候區(qū)，年平均氣溫為1.5～14 ℃，氣溫自西北向東南遞增，年均氣溫差距較大，1961年來以每10 a超過0.30 ℃的速率顯著上升，升溫速率在中國僅次于青藏高原。1956—2021年多年平均年降水量為531.6 mm，年內(nèi)分配極不均勻，夏季降水量占全年的60%～80%。流域河川徑流變化劇烈，水土資源分布嚴(yán)重不匹配，本地水資源極度匱乏。多年來，通過開展生態(tài)清潔小流域建設(shè)、風(fēng)沙源治理、坡耕地治理、退耕還林等工程，流域水土流失綜合治理成效顯著，林草覆蓋率持續(xù)增加。遙感監(jiān)測結(jié)果表明：1980—2020年間，海河流域山區(qū)耕地面積減少4 360 km²，林地面積增加840 km²，城鎮(zhèn)居民用地增加5 450 km²。在多元因素共同作用下，山區(qū)植被生長狀況明顯改善，年均歸一化植被指數(shù)（NDVI）由2000年的0.36增加到2020年的0.43（土地利用數(shù)據(jù)來源于中國科學(xué)院地理科學(xué)與資源研究所資源環(huán)境科學(xué)與數(shù)據(jù)中心，https：∥www.resdc.cn;NDVI數(shù)據(jù)來自美國NASA的MODIS產(chǎn)品數(shù)據(jù)MOD13A3v006，https：∥lpdaac.usgs.gov/products/mod13a3v006/）。

2 研究方法

M-K趨勢檢驗和Pettitt突變檢驗法廣泛應(yīng)用于水文氣象序列的統(tǒng)計檢驗^［21］。本文首先基于上述方法分析海河流域降水及天然河川徑流的變化趨勢，利用標(biāo)準(zhǔn)化徑流指數(shù)（SRI）劃分徑流豐枯，然后針對徑流持續(xù)衰減的特征，通過GAMLSS模型確定徑流系列的最優(yōu)時變參數(shù)分布，得到最優(yōu)分布下不同等級豐枯水事件的時變概率，在此基礎(chǔ)上結(jié)合期望發(fā)生次數(shù)和期望等待時間法揭示非一致性條件下不同等級豐枯事件發(fā)生次數(shù)及極值豐枯水發(fā)生時間間隔的變化規(guī)律。

2.1 基于GAMLSS的豐枯概率計算

SRI是監(jiān)測和評估水文干旱的一種指標(biāo)，反映徑流豐枯的變化程度。該指數(shù)通過對某個時間段的徑流系列選取合適的概率分布，估計分布最優(yōu)參數(shù)，采用等概率轉(zhuǎn)換將徑流量對應(yīng)的累積概率進(jìn)行正態(tài)標(biāo)準(zhǔn)化，從而得到一定時間尺度范圍的SRI值。本文以SRI作為海河流域天然河川徑流量的豐枯表征指標(biāo)，其計算類似于表征氣象干旱指標(biāo)的標(biāo)準(zhǔn)化降水指數(shù)（SPI）^［6］，選用與SPI相同的等級劃分標(biāo)準(zhǔn)，將徑流豐枯劃分為7個等級，詳見表1。用Z_SRI=y（y=-2，-1.5，-1，1，1.5，2）分別表示發(fā)生極端枯水、嚴(yán)重枯水、中等枯水、中等豐水、較大豐水和極端豐水的臨界徑流閾值。

GAMLSS模型包含眾多分布函數(shù)，同時可靈活建立分布參數(shù)與時間、氣象等協(xié)變量之間的函數(shù)關(guān)系^［22-23］ 。假設(shè)水文樣本序列y_t （t=1，2，…，n）服從分布函數(shù)F_z（y_t|θ_t），其中θt=（θ_t1，θ_t2，θ_t3，…，θ_tp）是包含p個參數(shù)的向量，通常在建模中p≤4。用g_k表示參數(shù)向量θ_t與解釋變量X_k和隨機(jī)效應(yīng)項之間的單調(diào)連接函數(shù)，在實際建模過程中，通常忽略隨機(jī)效應(yīng)的影響，GAMLSS模型可表示為

式中：η_k和θ_t是長度為n的參數(shù)向量;β_k為長度為J_k的回歸參數(shù)向量;

X_k為n×J_k的解釋變量矩陣;h_jk為解釋變量的函數(shù)，模型原理介紹可參考文獻(xiàn)［23］。

采用AIC（Akaike Information Criterion）準(zhǔn)則優(yōu)選分布，值最小時對應(yīng)的模型最優(yōu)，同時借助殘差正態(tài)QQ圖和蠕蟲圖^［24］評價模型擬合優(yōu)度。基于最優(yōu)GAMLSS分布模型，結(jié)合不同豐枯等級對應(yīng)的臨界徑流閾值，可推導(dǎo)出非一致性條件下第t年不同等級豐枯水事件累積時變概率，見式（2）。

式中：Z_SRI為豐枯等級閾值;θ_t為第t年的分布參數(shù);p_t為徑流低于閾值Z_SRI的累積時變概率。

2.2 期望發(fā)生次數(shù)

期望發(fā)生次數(shù)定義為從某一年開始在T年內(nèi)某一等級豐枯水事件發(fā)生次數(shù)的期望值。記N（N=1，2，3，…，T）為在T年間出現(xiàn)某一等級以上豐枯水事件的可能次數(shù)。對于枯水事件，N=∑^T_t=1I（Z_t≤Z_SRI），其中I（·）為示性函數(shù)，Z_t為第t年的徑流量。一致性條件下，隨機(jī)變量N服從二項分布，記q₀為徑流低于某一等級枯水事件閾值的概率，則在T年內(nèi)對應(yīng)等級枯水事件的期望發(fā)生次數(shù)為Tq₀。非一致性條件下，各等級水文事件在每年的發(fā)生概率不再為常數(shù)，此時N不再服從二項分布，在T年內(nèi)某一等級以上枯水事件的期望發(fā)生次數(shù)為

式中：P（·）為發(fā)生括號內(nèi)豐枯事件的概率。

對于豐水事件，記q₁為一致性條件下徑流高于某一等級豐水事件閾值的發(fā)生概率，則在T年內(nèi)該等級以上豐水事件的期望發(fā)生次數(shù)為Tq₁，非一致性條件下對應(yīng)的期望發(fā)生次數(shù)為

2.3 期望等待時間

記T為首次出現(xiàn)某一等級豐枯事件的年份（T=1，2，3…，∞），期望等待時間E（T）定義為從某一年開始直到下一次某一等級豐枯事件出現(xiàn)的平均時間間隔。對于某一等級以上枯水事件，一致性條件下隨機(jī)變量T服從幾何分布，期望等待時間為1/q₀，非一致性條件下T的概率密度函數(shù)及對應(yīng)的期望等待時間分別為：

對于某一等級以上豐水事件，一致性條件下的期望等待時間為1/q₁，非一致性條件下對應(yīng)的期望等待時間為

3 結(jié)果與分析

3.1 徑流變化趨勢

本文采用的海河流域1956—2021年降水及徑流數(shù)據(jù)序列來源于第二次全國水資源調(diào)查評價和歷年的《中國水資源公報》，趨勢分析結(jié)果見圖2。結(jié)果表明，流域年降水量存在輕微的下降趨勢，M-K檢驗及Pettitt檢驗的p值分別為0.29和0.26，均未通過顯著性水平95%檢驗，可認(rèn)為該序列仍保持平穩(wěn)性，而天然河川徑流量下降趨勢明顯，M-K檢驗統(tǒng)計量為-3.58，通過了顯著性水平為99%的趨勢檢驗。Pettitt突變檢驗表明

，徑流量向下跳躍突變發(fā)生于1979年，通過了顯著性為95%的突變檢驗，1956—1979年多年平均天然河川徑流量為255.9億m³，1980—2021年則變?yōu)?58.9億m³，衰減幅度達(dá)38%。

3.2 概率計算模型優(yōu)選

選取GAMLSS模型中6種常用的兩參數(shù)分布為備選概率分布集，包括正態(tài)分布（Normal，NO）、對數(shù)正態(tài)分布（Log Normal，LOGNO）、耿貝爾分布（Gumbel，GU）、伽瑪分布（Gamma，GA）、威布爾分布（Weibull，WEI）和邏輯斯諦分布（Logistic，LO），分布參數(shù)設(shè)置為常變量和以時間t為協(xié)變量，對流域天然河川徑流量序列進(jìn)行頻率計算。為避免分布參數(shù)與時間t的回歸方程太過復(fù)雜，分布位置參數(shù)和尺度參數(shù)與時間t的多項式次數(shù)不超過2，分別設(shè)為常變量及與時間t成1次方、2次方的多項式函數(shù)共9種組合對序列進(jìn)行擬合，選取AIC值最小的概率分布和參數(shù)作為最優(yōu)模型。結(jié)果表明，LOGNO分布的位置參數(shù)μ（t）=1 310-1.3t+0.000 32t²，尺度參數(shù)為常數(shù)時的模型AIC值最小，位置參數(shù)反映整個水文序列均值情況，體現(xiàn)了海河流域2010年前徑流量持續(xù)減小、2010年后略有回升的趨勢，與實際情況吻合。

為進(jìn)一步驗證最優(yōu)非一致時變參數(shù)模型的擬合效果，防止模型過度擬合，圖3展示了最優(yōu)模型殘差的正態(tài)QQ圖和蠕蟲圖。QQ圖上的散點分布在1∶1直線附近，認(rèn)為最優(yōu)模型的殘差序列服從標(biāo)準(zhǔn)正態(tài)分布。蠕蟲圖中間的紅線是由圖中散點系列擬合的三次多項式曲線，所有散點都位于上下2條曲線之間的置信區(qū)間內(nèi)，從而判斷最優(yōu)模型的分布類型和參數(shù)選擇合理，模型擬合效果較好。

3.3 發(fā)生次數(shù)的比較

計算海河流域逐年SRI值（I_SR），同時根據(jù)前面的徑流量突變點分析結(jié)果，劃分1956—1979年（簡稱PⅠ時期）及1980—2021年（簡稱PⅡ時期）2個時期揭示徑流豐枯變化特征，I_SR年際變化及不同時期豐枯等級發(fā)生次數(shù)見圖4。結(jié)果表明：I_SR呈顯著的下降趨勢（M-K檢驗統(tǒng)計量-3.6），通過置信水平95%的顯著性檢驗，下降速率為0.34/（10 a），均值由PⅠ時期的0.84下降至PⅡ的-0.70。PⅠ時期豐水事件出現(xiàn)次數(shù)較多，而PⅡ時期豐水事件次數(shù)遠(yuǎn)低于枯水事件，枯水事件主要集中于1997—2010年間，甚至在2000—2001年連續(xù)出現(xiàn)極枯事件，河川徑流的豐枯概率呈現(xiàn)出顯著的枯增豐減趨勢。Ling等^［25］采用標(biāo)準(zhǔn)化降水蒸散指數(shù)（SPEI）研究海河流域氣象干旱變化特征，結(jié)論與上述SRI有同樣的變化規(guī)律，但SRI表征河川徑流干枯程度更為嚴(yán)重。1964年和2021年分別是PⅠ、PⅡ時期降水量最大的年份，前者降水量相比后者少4.6%，但徑流量卻高出1.7%，可見流域內(nèi)產(chǎn)流能力減弱，人類活動導(dǎo)致下墊面變化是河川徑流減少的重要影響因素，同時氣溫升高、人類活動取用水、地下水超采等均有一定影響^［19］ 。

利用優(yōu)選的LOGNO時變分布，計算不同等級豐枯事件在非一致性條件下的發(fā)生概率，采用期望發(fā)生次數(shù)法計算一致性和非一致性條件下海河流域歷史期豐枯事件的期望發(fā)生次數(shù)，并與實際發(fā)生次數(shù)比較。由表2可知，一致性假設(shè)情況下，不同等級豐水和枯水事件的發(fā)生次數(shù)在PⅠ時期分別被低估和高估，而在PⅡ時期則分別被高估和低估。與一致性假設(shè)相比，考慮非一致性計算豐枯事件的期望發(fā)生次數(shù)與歷史期的實際發(fā)生次數(shù)更接近，比較符合海河流域徑流演變特征。因此，非一致性模型描述海河流域歷史期豐枯特征要優(yōu)于一致性模型。

3.4 極值豐枯事件的期望等待時間

探求未來極值豐枯事件發(fā)生的期望等待時間，重點是給出未來t=1 a，2 a，3 a，…的分布參數(shù)。鑒于徑流變化存在一定的不確定性，本文立足于流域現(xiàn)狀，基于3.2節(jié)擬合的最優(yōu)時變LOGNO分布模型設(shè)置2種情景：一是假設(shè)未來年份維持非一致性模型的2021年分布（記為M₁分布），此情景代表流域徑流變化近期已接近穩(wěn)定，分布參數(shù)及逐年發(fā)生極值豐枯事件的概率為固定值;二是以時間t為協(xié)變量往后延伸計算（記為M₂分布），此情景下分布位置參數(shù)與時間t呈二次方，近期將呈現(xiàn)一定的增大趨勢，代表流域在近期降水回升和水源涵養(yǎng)加強(qiáng)背景下未來徑流將呈現(xiàn)一定的回升趨勢。選取歷史序列的下一年（2022年）作為起始年，計算不同分布下未來極端豐枯事件在第t年首次發(fā)生的概率及期望等待時間，如圖5所示。

圖5（a）和圖5（b）縱坐標(biāo)分別為極端豐水、極端枯水事件在未來第t年首次發(fā)生的概率。在M₁分布下極豐事件的縱坐標(biāo)隨著t增大而減小，M₂分布下極豐事件的縱坐標(biāo)則先增加后減少，而極枯事件縱坐標(biāo)均隨著t增大而減小，但兩者衰減速率有差異。圖5（c）和圖5（d）縱坐標(biāo)分別為預(yù)估極端豐水、極端枯水事件的期望等待時間，M₁、M₂分布下極端豐水事件的期望等待時間分別為36.0 a、14.5 a，而在過去65 a中平均9.2 a出現(xiàn)一次極端豐水;在M₁、M₂分布下極端枯水事件的期望等待時間分別為9.4 a、4.9 a，而在過去65 a中平均21.6 a出現(xiàn)一次極端枯水，極枯事件的平均間隔時間將大幅縮短。本文中不同等級豐枯水事件期望等待時間的計算均能在未來15～30 a內(nèi)收斂，模型預(yù)測的不確定性相對較小。綜合上述分析，預(yù)估海河流域近期發(fā)生枯水的頻率較高，發(fā)生豐水的頻率相對較小，預(yù)估極枯、極豐事件下一次發(fā)生的期望等待時間分別為4.9～9.4 a和14.5～36.0 a。目前已有文獻(xiàn)普遍認(rèn)為下墊面與人類活動影響是海河流域徑流減少的主要原因^［19-20］，由于這種影響帶有一定的不可逆性，徑流豐枯很難在短期內(nèi)恢復(fù)到20世紀(jì)80年代以前的水平，因此，本文關(guān)于期望等待時間的計算結(jié)果與已有文獻(xiàn)的認(rèn)識相符。此外，由于不同氣候變化排放情景和模式的預(yù)估結(jié)果間可能存在較大差異，本文在分析徑流豐枯概率變化時沒有考慮未來氣候變化情景。

4 結(jié)論

本文針對海河流域徑流量變化劇烈、水量銳減等問題，提出了一套研究變化環(huán)境下徑流豐枯特征的計算方法，基于1956—2021年流域降水及天然河川徑流量數(shù)據(jù)，研究了天然河川徑流豐枯概率及未來發(fā)生時間間隔的演變規(guī)律，主要結(jié)論如下：

（1） 基于GAMLSS的時變矩模型中，以時間t為協(xié)變量的LOGNO分布為擬合徑流系列的最優(yōu)分布，且位置參數(shù)表征水文系列均值存在顯著的非一致性，與時間t成二次多項式，尺度參數(shù)為常數(shù)。

（2） 流域標(biāo)準(zhǔn)化徑流指數(shù)持續(xù)減小，徑流豐枯概率呈現(xiàn)出顯著的枯增豐減趨勢。

（3） 非一致性條件下不同豐枯水等級在PⅠ與PⅡ 2個時期的期望發(fā)生次數(shù)更接近于實際發(fā)生情況，對于豐枯水事件擬合結(jié)果非一致性情況優(yōu)于一致性假設(shè)情況。

（4） 預(yù)估極枯、極豐水文事件下一次發(fā)生的期望等待時間分別為4.9～9.4 a、14.5～36.0 a，發(fā)生枯水的頻率相對較高，存在一定的不確定性。

參考文獻(xiàn)：

［1］王國慶，張建云，管曉祥，等.中國主要江河徑流變化成因定量分析［J］.水科學(xué)進(jìn)展，2020，31（3）：313-323.（WANG G Q，ZHANG J Y，GUAN X X，et al.Quantifying attribution of runoff change for major rivers in China［J］.Advances in Water Science，2020，31（3）：313-323.（in Chinese））

［2］周祖昊，劉佳嘉，嚴(yán)子奇，等.黃河流域天然河川徑流量演變歸因分析［J］.水科學(xué)進(jìn)展，2022，33（1）：27-37.（ZHOU Z H，LIU J J，YAN Z Q，et al.Attribution analysis of the natural runoff evolution in the Yellow River basin［J］.Advances in Water Science，2022，33（1）：27-37.（in Chinese））

［3］SUN P，WEN Q Z，ZHANG Q，et al.Nonstationarity-based evaluation of flood frequency and flood risk in the Huai River basin，China［J］.Journal of Hydrology，2018，567：393-404.

［4］梁忠民，胡義明，王軍.非一致性水文頻率分析的研究進(jìn)展［J］.水科學(xué)進(jìn)展，2011，22（6）：864-871.（LIANG Z M，HU Y M，WANG J.Advances in hydrological frequency analysis of non-stationary time series［J］.Advances in Water Science，2011，22（6）：864-871.（in Chinese））

［5］白君瑞，徐宗學(xué)，班春廣，等.基于Z指數(shù)的雅魯藏布江流域徑流豐枯變化及其特征分析［J］.北京師范大學(xué)學(xué)報（自然科學(xué)版），2019，55（6）：715-723.（BAI J R，XU Z X，BAN C G，et al.Runoff variation and characteristics in Yarlung Zangbo River by Z-index［J］.Journal of Beijing Normal University（Natural Science），2019，55（6）：715-723.（in Chinese））

［6］氣象干旱等級：GB/T 20481—2017［S］.北京：中國標(biāo)準(zhǔn)出版社，2017.（Grades of meteorological drought：GB/T 20481—2017［S］.Beijing：Standards Press of China，2017.（in Chinese））

［7］涂新軍，謝育廷，吳海鷗，等.基于概率矩陣的干旱等級變化評估及應(yīng)用［J］.水科學(xué)進(jìn)展，2021，32（4）：520-533.（TU X J，XIE Y T，WU H O，et al.An evaluation method for drought level changes based on probability matrix and its application［J］.Advances in Water Science，2021，32（4）：520-533.（in Chinese））

［8］熊立華，鄺韻琪，于坤霞，等.年徑流頻率分析的一次二階矩法及其應(yīng)用［J］.水科學(xué)進(jìn)展，2017，28（3）：390-397.（XIONG L H，KUANG Y Q，YU K X，et al.Derived moment approach for analyzing the annual runoff frequency and its application［J］.Advances in Water Science，2017，28（3）：390-397.（in Chinese））

［9］SONG X Y，LU F，WANG H，et al.Penalized maximum likelihood estimators for the nonstationary Pearson type 3 distribution［J］.Journal of Hydrology，2018，567：579-589.

［10］LU F，SONG X Y，XIAO W H，et al.Detecting the impact of climate and reservoirs on extreme floods using nonstationary frequency models［J］.Stochastic Environmental Research and Risk Assessment，2020，34（1）：169-182.

［11］XIONG B，XIONG L H，GUO S L，et al.Nonstationary frequency analysis of censored data：a case study of the floods in the Yangtze River from 1470 to 2017［J］.Water Resources Research，2020，56（8）：e2020WR027112.

［12］顧西輝，張強(qiáng)，孔冬冬，等.新疆塔里木河流域洪水發(fā)生率非平穩(wěn)性特征及氣候變化影響研究［J］.自然資源學(xué)報，2016，31（9）：1499-1513.（GU X H，ZHANG Q，KONG D D，et al.Nonstationarity of flooding processes in the Tarim River basin and climate-related impacts［J］.Journal of Natural Resources，2016，31（9）：1499-1513.（in Chinese））

［13］XU Y R，LU F，ZHOU Y Y，et al.Dryness-wetness encounter probabilities analysis for lake ecological water replenishment considering non-stationarity effects［J］.Frontiers in Environmental Science，2022，10：806794.

［14］YAN L，XIONG L H，GUO S L，et al.Comparison of four nonstationary hydrologic design methods for changing environment［J］.Journal of Hydrology，2017，551：132-150.

［15］王銀堂，李伶杰，胡慶芳，等.考慮局部趨勢的非一致性水文頻率分析方法［J］.水科學(xué)進(jìn)展，2017，28（3）：406-414.（WANG Y T，LI L J，HU Q F，et al.Nonstationary hydrologic frequency analysis method considering local trends［J］.Advances in Water Science，2017，28（3）：406-414.（in Chinese））

［16］梁忠民，胡義明，王軍，等.基于等可靠度法的變化環(huán)境下工程水文設(shè)計值估計方法［J］.水科學(xué)進(jìn)展，2017，28（3）：398-405.（LIANG Z M，HU Y M，WANG J，et al.Estimation of design flood using equivalent reliability method under changing environment［J］.Advances in Water Science，2017，28（3）：398-405.（in Chinese））

［17］DU T，XIONG L H，XU C Y，et al.Return period and risk analysis of nonstationary low-flow series under climate change［J］.Journal of Hydrology，2015，527：234-250.

［18］READ L K，VOGEL R M.Reliability，return periods，and risk under nonstationarity［J］.Water Resources Research，2015，51（8）：6381-6398.

［19］鮑振鑫，張建云，嚴(yán)小林，等.基于四元驅(qū)動的海河流域河川徑流變化歸因定量識別［J］.水科學(xué)進(jìn)展，2021，32（2）：171-181.（BAO Z X，ZHANG J Y，YAN X L，et al.Quantitative assessment of the attribution of runoff change caused by four factors in the Haihe River basin［J］.Advances in Water Science，2021，32（2）：171-181.（in Chinese））

［20］YANG L，ZHAO G J，TIAN P，et al.Runoff changes in the major river basins of China and their responses to potential driving forces［J］.Journal of Hydrology，2022，607：127536.

［21］WANG M H，JIANG S H，REN L L，et al.Nonstationary flood and low flow frequency analysis in the upper reaches of Huaihe River basin，China，using climatic variables and reservoir index as covariates［J］.Journal of Hydrology，2022，612：128266.

［22］JEHANZAIB M，ALI SHAH S，YOO J，et al.Investigating the impacts of climate change and human activities on hydrological drought using non-stationary approaches［J］.Journal of Hydrology，2020，588：125052.

［23］RIGBY R A，STASINOPOULOS D M.Generalized additive models for location，scale and shape （with discussion）［J］.Journal of the Royal Statistical Society：Series C （Applied Statistics），2005，54（3）：507-554.

［24］魯帆，肖偉華，戴雁宇，等.黃河干流年徑流量非一致性頻率計算［J］.水力發(fā)電學(xué)報，2020，39（12）：76-84.（LU F，XIAO W H，DAI Y Y，et al.Study on non-stationary frequency calculations for the annual runoff of Yellow River mainstream［J］.Journal of Hydroelectric Engineering，2020，39（12）：76-84.（in Chinese））

［25］LING M H，GUO X M，SHI X L，et al.Temporal and spatial evolution of drought in Haihe River basin from 1960 to 2020［J］.Ecological Indicators，2022，138：108809.

Evolution law of wet and dry probability of natural river runoff in Haihe River basin under changing environment

The study is financially supported by the National Key R&D Program of China（No.2021YFC3000200;No.2018YFC0406500）.

LU Fan¹，JIANG Ming¹，JIANG Yunzhong¹，ZHOU Yuyan¹，XU Yang^1，2

（1. Department of Water Resources，China Institute of Water Resources and Hydropower Research，Beijing 100038，China;2. School of Resources and Earth Science，China University of Mining and Technology，Xuzhou 221116，China）

Abstract：The natural river flow in the Haihe River basin has been continuously reducing in recent years，which leads to significant changes particularly in the wet and dry conditions of hydrology.As such，research on the calculation method of wet and dry probability catering for nonstationary hydrological series is needed.Based on the standardized runoff index and GAMLSS model，a new method computing expected number of occurrences （ENO） and expected waiting time （EWT） was proposed according to different levels of wet and dry hydrological events.Then，we further investigate on the evolution law of probability of wet and dry of natural river flow in Haihe River basin under changing environment.Our main finding are outlined as follows：① The probability of wet and dry of surface runoff showed a significant trend of increasing low flow and decreasing high flow.②? Compared with other probability distributions including the traditional stationary hypothesis，the LOGNO distribution with time（t） as the covariate has the best performance when fitting the surface runoff series.The ENO of wet and dry hydrological events in the historical period calculated based on this distribution is closer to the actual situation.③ The EWT of extremely dry and wet hydrological events in the Haihe River basin calculated under different scenarios of the nonstationary optimal model were 4.9—9.4 a and 14.5—36.0 a respectively，indicating that the considerably higher likelihood of extreme dry event in oncoming future.

Key words：river runoff;wet and dry probability;nonstationary;GAMLSS model;Haihe River basin