考慮非一致性條件的淮河干流洪水分析

洪雙玲， 石 朋， 瞿思敏， 馮 穎， 趙夢杰

(1.河海大學(xué) 水文水資源學(xué)院， 江蘇 南京 210098； 2.福建省水利水電勘測設(shè)計(jì)研究院，福建 福州 350000； 3.淮河水利委員會(huì)水文局(信息中心)， 安徽 蚌埠 233001)

1 研究背景

全球氣候變化使得降水時(shí)空分配發(fā)生改變，環(huán)境的變化和不確定性導(dǎo)致水文序列表現(xiàn)出非平穩(wěn)特征，因此基于一致性假設(shè)的洪水頻率分析很可能低估洪水風(fēng)險(xiǎn)[1-4]。目前非一致性水文頻率分析大多是基于還原或還現(xiàn)的途徑，只能反映過去或現(xiàn)狀條件下的水文過程，難以反映流域水文情勢的變化[5]。因此，研究非一致性條件下的水文頻率分析具有重要意義。非平穩(wěn)條件下，水文隨機(jī)變量的頻率分布會(huì)隨時(shí)間發(fā)生變化。Vogel等[6]基于兩參數(shù)對數(shù)正態(tài)分布構(gòu)建了分布函數(shù)參數(shù)與時(shí)間的線性關(guān)系，引入洪水放大因子概念，對比同一設(shè)計(jì)標(biāo)準(zhǔn)下，未來設(shè)計(jì)流量相比現(xiàn)在的變化情況，量化了氣候變化和人類活動(dòng)對極端事件分布影響在未來的變化情況。在非一致性情況下，水文要素邊際分布的概率參數(shù)隨時(shí)間發(fā)生變化，可以表示為協(xié)變量的函數(shù)，以模擬變化情況下水文要素的非平穩(wěn)特征[7-8]。

洪水過程是一個(gè)多要素水文過程，包含洪峰、時(shí)段洪量等多個(gè)特征變量，為考慮各水文變量之間的相關(guān)性，更全面地理解洪水過程，需要分析其多變量聯(lián)合分布特征[9]。Wahl等[10]發(fā)現(xiàn)美國沿海城市強(qiáng)降水和風(fēng)暴潮之間的相關(guān)關(guān)系表現(xiàn)出非平穩(wěn)性，證明了復(fù)合洪水風(fēng)險(xiǎn)的增加及其與天氣和氣候的聯(lián)系。水文要素之間的相關(guān)性往往是非線性的，其邊際分布也常是非正態(tài)的，在此情況下構(gòu)建多變量聯(lián)合分布具有極大的挑戰(zhàn)性[11]。Copula函數(shù)作為多元數(shù)據(jù)依賴結(jié)構(gòu)的有力解釋工具，在連接水文隨機(jī)變量復(fù)雜邊際分布、模擬變量間非線性關(guān)系方面具有優(yōu)勢[12-13]。馮穎等[14]基于考慮非一致性的GAMLSS(generalized additive models for location, scale and shape)模型，采用Copula函數(shù)構(gòu)建了淮河流域干、支流洪峰序列的聯(lián)合分布模型，研究了變化環(huán)境下干、支流洪水的遭遇規(guī)律。

淮河流域地處我國南北氣候過渡帶，夏季降水充沛，暴雨洪水事件頻發(fā)，給周邊地區(qū)帶來嚴(yán)重危害。近年來，受氣候變化影響，淮河流域的降水情勢發(fā)生改變，極端水文事件發(fā)生頻率增加，如何評估這種變化環(huán)境下的洪水風(fēng)險(xiǎn)，對于今后淮河流域重要水利工程的高效管理應(yīng)用具有重要的參考價(jià)值。

2 數(shù)據(jù)來源與研究方法

2.1 研究區(qū)域概況

淮河流域位于113°16′E～116°15′E，31°15′N～33°28′N之間，總面積為27×104km2，淮河以北屬暖溫帶半濕潤季風(fēng)氣候區(qū)，以南屬亞熱帶濕潤季風(fēng)氣候區(qū)，多年平均降水量約為920 mm，其分布狀況大致為由南向北遞減，山區(qū)多于平原，沿海大于內(nèi)陸?；春恿饔蛩导八摹⒂炅空军c(diǎn)分布見圖1。

圖1 淮河流域水系及水文、雨量站點(diǎn)分布

2.2 數(shù)據(jù)來源

研究所使用的數(shù)據(jù)包括90 m分辨率DEM(digital elevation model)數(shù)據(jù)集、流量數(shù)據(jù)和降水?dāng)?shù)據(jù)。流量數(shù)據(jù)來源于淮河流域內(nèi)干流4個(gè)水文站(表1)的日流量觀測數(shù)據(jù)集，降水?dāng)?shù)據(jù)來源于流域內(nèi)16個(gè)雨量站(表2)的日降水觀測數(shù)據(jù)集。洪水序列采用年最大值取樣法進(jìn)行取樣。

表1 淮河流域干流水文站點(diǎn)

表2 淮河流域雨量站

降水序列選擇控制流域內(nèi)16個(gè)雨量站的逐日降水資料，以泰森多邊形劃分每個(gè)控制流域內(nèi)各雨量站權(quán)重，依權(quán)重計(jì)算水文站點(diǎn)控制流域內(nèi)的面平均日降水量。

2.3 研究方法

2.3.1 非一致性檢驗(yàn) 變化環(huán)境下，水文序列出現(xiàn)很多非一致性問題，其非一致性檢驗(yàn)包括突變點(diǎn)檢驗(yàn)、趨勢檢驗(yàn)和周期檢驗(yàn)。由于滑動(dòng)平均差法能同時(shí)檢測出所有突變點(diǎn)的突變位置和突變強(qiáng)度[15]，故本文采用該方法與Pettitt法、CUSUM(cumulative sum)檢驗(yàn)法共同進(jìn)行突變點(diǎn)分析，采用Mann-Kendall法進(jìn)行趨勢分析。

2.3.2 邊際分布的確定 采用GAMLSS模型確定各洪水要素的頻率分布。GAMLSS模型是一種評估時(shí)間序列的半?yún)?shù)回歸模型，可以用來描述服從同一分布簇的隨機(jī)變量序列的統(tǒng)計(jì)參數(shù)與協(xié)變量之間的線性或非線性關(guān)系。本文只考慮模型參數(shù)隨單個(gè)協(xié)變量x變化的情況，對于t時(shí)刻，相應(yīng)變量服從分布函數(shù)UX(yt|θt)，其中g(shù)(·)定義為參數(shù)θt的單調(diào)變化函數(shù)。以雙參數(shù)單變量形式(包括均值u、方差σ)為例，其表達(dá)式如下：

g(ut)=a1xt+b1

(1)

g(σt)=a2xt+b2

(2)

模型參數(shù)采用極大似然法進(jìn)行估計(jì)，選擇面均降水P作為參數(shù)協(xié)變量，以AIC(Akaike information criterion)準(zhǔn)則和SBC(Schewarz Bayesian criterion)準(zhǔn)則對模型進(jìn)行擇優(yōu)。

2.3.3 Copula函數(shù) 根據(jù)Sklar定理，設(shè)X1，X2，…，Xd為具有邊際分布U1，U2，…，Ud的隨機(jī)變量，聯(lián)合分布為U(x1，x2,…，xd)，對于所有x∈R，R∈(-∞，∞)，存在一個(gè)d維CopulaC，描述多元聯(lián)合分布和相應(yīng)的單變量分布函數(shù)之間的關(guān)系式為：

U=(x1,x2,…,xd)

=P(X1≤x1,X2≤x2,…,Xd≤xd)

=C[U1(x1),U2(x2),…,Ud(xd)|θc]

=C(u1,u2,…,ud|θc)

(3)

式中：C(·)：[0,1]d→[0,1]為一個(gè)d維Copula函數(shù)；θc為Copula函數(shù)的參數(shù), 且ui～U(0,1)。

Copula函數(shù)種類眾多，水文領(lǐng)域中，對稱Archimedean Copula函數(shù)因其結(jié)構(gòu)簡單、構(gòu)造形式多樣及適應(yīng)性強(qiáng)等優(yōu)點(diǎn)被廣泛應(yīng)用，對于滿足交換性檢驗(yàn)的兩變量序列即可采用該Copula函數(shù)連接，Archimedean Copula函數(shù)分布簇中的Gumbel-Hougaard Copula、Clayton Copula和Frank Copula這3種二維函數(shù)在水文上應(yīng)用較多，故從中進(jìn)行優(yōu)選，其函數(shù)表達(dá)式見表3。

表3 水文領(lǐng)域中常用的3種Archimedean Copula函數(shù)

考慮非一致情況的邊際分布較為復(fù)雜，故Copula函數(shù)的參數(shù)估計(jì)采用半?yún)?shù)推斷法[16]，該方法用經(jīng)驗(yàn)頻率代替理論計(jì)算值，有效避免了理論值誤差對Copula參數(shù)估計(jì)的影響，表示為：

(4)

式中：UX(x)與UY(y)分別為兩變量X與Y的經(jīng)驗(yàn)累積分布函數(shù)，即有：

(5)

(6)

將經(jīng)驗(yàn)分布函數(shù)的累積概率值代入Copula對數(shù)似然函數(shù)中，使得對數(shù)似然函數(shù)值最大，即可得到Copula參數(shù)。根據(jù)AIC準(zhǔn)則和SBC準(zhǔn)則對Copula模型進(jìn)行擇優(yōu)，繪制Pickands依賴圖[17]，對模型擬合優(yōu)度進(jìn)行定性分析。

2.3.4 兩變量同現(xiàn)重現(xiàn)期 水文頻率分析的目的是確定某一洪水設(shè)計(jì)值下的設(shè)計(jì)標(biāo)準(zhǔn)，即重現(xiàn)期。設(shè)UX(x)和UY(y)分別為洪水要素X和洪水要素Y的分布函數(shù)。單要素重現(xiàn)期計(jì)算公式如下：

(7)

(8)

同現(xiàn)重現(xiàn)期表示洪水事件中兩要素均超過某一特定值，公式如下：

(9)

3 結(jié)果與分析

3.1 洪峰洪量序列非一致性檢測

3.1.1 突變點(diǎn)檢驗(yàn) 采用Pettitt法、CUSUM檢驗(yàn)法和滑動(dòng)平均差法3種突變點(diǎn)檢測方法對流域干流各站水文序列(洪峰Qm、15 d洪量W15)進(jìn)行突變點(diǎn)檢驗(yàn)，表4給出了突變檢驗(yàn)的結(jié)果。由表4可看出，不同的檢驗(yàn)方法可能得到不同的結(jié)果。對于同一站點(diǎn)，洪峰序列和洪量序列的突變年份大多相同，不同站點(diǎn)的突變年份一般有所不同。

表4 淮河流域干流水文站洪水序列突變點(diǎn)檢驗(yàn)結(jié)果

3.1.2 趨勢檢驗(yàn) 采用Mann-Kendall法對各水文站洪峰和洪量序列進(jìn)行趨勢分析，表5給出了各站序列變化趨勢。從表5中的結(jié)果來看，各水文站的洪峰序列和洪量序列均呈減小趨勢，其中息縣站洪峰和洪量減小趨勢較為顯著，潤河集站峰量減小趨勢較為不顯著。

表5 淮河流域干流水文站洪水序列趨勢分析結(jié)果

3.2 邊際分布

3.2.1 協(xié)變量選擇 氣候變化對水文的直接影響是降水條件的改變，降水又是洪水的主要成因，與洪水各要素的相關(guān)性較強(qiáng)，故采用水文站點(diǎn)控制流域內(nèi)的面均降水量作為洪峰Qm和最大15 d洪量W15邊際分布統(tǒng)計(jì)參數(shù)的協(xié)變量。以年最大值法選取年最大1 d面均降水量(P1)、年最大3 d面均降水量(P3)、年最大5 d面均降水量(P5)、年最大7 d面均降水量(P7)、年最大15 d面均降水量作為備選協(xié)變量(P15)，分別將其與洪峰和洪量進(jìn)行相關(guān)性分析。圖2為各面均降水量與洪峰、洪量的相關(guān)性熱圖，采用Kendall相關(guān)系數(shù)、Spearman相關(guān)系數(shù)、Pearson相關(guān)系數(shù)進(jìn)行度量。

根據(jù)圖2熱圖顯示，各水文站的洪峰和洪量與各面均降水量均存在較強(qiáng)的正相關(guān)性，選擇相關(guān)性最強(qiáng)的降水量序列作為峰量邊際分布時(shí)變參數(shù)的協(xié)變量。其中息縣站、王家壩站洪峰邊際分布選擇P15，淮濱站、潤河集站選擇P5；4個(gè)站點(diǎn)的最大15 d洪量邊際分布模型均選擇P15。

圖2 各面均降水量與各水文站洪峰、15 d洪量相關(guān)關(guān)系熱圖表6 各水文站洪峰序列最優(yōu)模型擬合結(jié)果

3.2.2 邊際分布擬合 GAMLSS模型包含多種分布，選擇水文上較為常用的4種雙參數(shù)概率分布： Gamma(GA)分布、Weibull (WEI) 分布、Lognormal (LOGNO) 分布以及Gumbel(GU)分布作為備選分布參數(shù)，每種分布考慮4種情況：(1)均值與方差均不變；(2)均值隨協(xié)變量變化，方差為常數(shù)；(3)均值為常數(shù)，方差隨協(xié)變量變化；(4)均值和方差均隨協(xié)變量變化。第1種為一致性情況，后3種為非一致性情況。

采用極大似然法進(jìn)行參數(shù)估計(jì)，根據(jù)AIC、SBC最小準(zhǔn)則選擇最優(yōu)模型，以worm圖和正態(tài)QQ圖定性分析模型的擬合優(yōu)度，各水文站洪峰和洪量的模型擬合結(jié)果分別見表6和表7。

由表6中AIC和SBC結(jié)果來看，洪峰邊際分布非一致性模型擬合效果均優(yōu)于一致性模型，其中，息縣站、淮濱站、王家壩站的最優(yōu)模型皆為均值隨降水協(xié)變量變化、方差不變，潤河集站的最優(yōu)模型為均值和方差皆隨降水協(xié)變量變化。

根據(jù)表7中結(jié)果，息縣站、淮濱站、王家壩站洪量邊際分布的最優(yōu)模型為均值和方差均隨降水協(xié)變量變化，只有潤河集站洪量邊際分布為均值隨降水協(xié)變量變化、方差不變。由趨勢分析結(jié)果可以看出，潤河集站峰量趨勢變化并不明顯，但是其峰量邊際分布的非一致性模型優(yōu)于一致性模型。

表7 各水文站洪量序列最優(yōu)模型擬合結(jié)果

峰量序列GAMLSS模型殘差定量評價(jià)指標(biāo)如表8所示。由表8可看出，各模型殘差均值基本為0，方差為1，且根據(jù)Filliben系數(shù)及KS檢驗(yàn)(Kolmogorov-Smirnovtest)表明，各模型均通過顯著性水平為0.05的假設(shè)，因此可以認(rèn)為所選的模型是合理的。各水文站一致性和非一致性模型擬合效果可通過worm圖及正態(tài)QQ圖直觀表示，如圖3～6所示。

表8 峰量序列GAMLSS模型殘差定量評價(jià)指標(biāo)

圖3、4的worm圖顯示，模型殘差點(diǎn)據(jù)基本落在95%置信區(qū)間內(nèi)；由圖5、6的正態(tài)QQ圖可以直觀看出，峰量邊際分布模型的殘差點(diǎn)據(jù)大致均勻地分布在45°直線兩側(cè)。這說明模型擬合效果較好，模型選擇合理。

圖3 水文站洪峰序列(Qm)邊際分布模型擬合殘差圖

圖4 水文站洪量序列(W15)邊際分布模型擬合殘差圖

圖5 水文站洪峰序列(Qm)邊際分布模型正態(tài)QQ圖

3.2.3 時(shí)變設(shè)計(jì)洪峰和洪量 根據(jù)擬合的邊際分布，分別計(jì)算各水文站50年一遇的設(shè)計(jì)洪峰流量和15 d設(shè)計(jì)洪量，結(jié)果如圖7、8所示。由圖7、8可以看出，非一致性條件下邊際分布統(tǒng)計(jì)參數(shù)逐年變化，其設(shè)計(jì)洪峰值也隨之發(fā)生變化；各水文站非一致性條件下設(shè)計(jì)洪峰、設(shè)計(jì)洪量與實(shí)測流量、實(shí)測洪峰變化趨勢一致。

圖6 水文站洪量序列(W15)邊際分布模型正態(tài)QQ圖

圖7 不同年份各水文站50年一遇時(shí)變設(shè)計(jì)洪峰流量

水文站模型AICSBC分布函數(shù)形式參 數(shù)μσ息縣 一致11451149GA25640.76非一致10701038WEIexp(6.40+0.0058cs(P15))2.65淮濱 一致10081012GA27010.739非一致941953WEIexp (6.70+0.007cs(P5))2.81王家壩 一致947951GA35530.723非一致885896WEIexp (6.74+0.006cs(P15))2.91潤河集 一致647654GA34790.654非一致620635LOGNOexp (6.60+0.009cs(P5))exp(-0.48-0.004cs(P5))

通過圖中數(shù)據(jù)對比可知，在大洪水年(1954、1968、2007年等)，各站在一致性模型下的設(shè)計(jì)洪峰和洪量均小于非一致性模型的計(jì)算結(jié)果，其余大多數(shù)年份，一致性模型下的設(shè)計(jì)洪峰和洪量大于非一致性模型的結(jié)果?；春痈闪鞲髡局?，除了潤河集站的最高設(shè)計(jì)洪峰為1982年，其余站點(diǎn)均于1968或2007年出現(xiàn)最高設(shè)計(jì)洪峰或次高設(shè)計(jì)洪峰，這是1968、2007年淮河流域罕見特大暴雨導(dǎo)致的結(jié)果。

3.3 聯(lián)合分布

水文領(lǐng)域中，常用Archimedean Copula 函數(shù)族擬合多變量間的相關(guān)關(guān)系，在選取Copula函數(shù)前需對所選數(shù)據(jù)進(jìn)行交換性檢驗(yàn)[18-19]，檢驗(yàn)方法采用Cramér-von Mises檢驗(yàn)法，檢驗(yàn)結(jié)果見表9。只有滿足交換性檢驗(yàn)，才可采用具有對稱結(jié)構(gòu)的Archimedean Copula函數(shù)族，否則需要選擇其他非對稱的Copula函數(shù)。

圖8 不同年份各水文站50年一遇時(shí)變15 d設(shè)計(jì)洪量

表9中結(jié)果顯示，各站峰量序列檢驗(yàn)P值均大于0.05，通過交換性檢驗(yàn)，可以采用具有對稱結(jié)構(gòu)的Archimedean Copula函數(shù)進(jìn)行擬合。根據(jù)峰量序列的邊際分布模型，采用Copula函數(shù)構(gòu)建峰量聯(lián)合分布。采用3種Copula函數(shù)擬合，得到峰量序列的Copula參數(shù)估計(jì)結(jié)果，如表10所示。

表9 各站點(diǎn)洪峰序列與洪量序列交換性檢驗(yàn)

由表10可見，息縣、淮濱和潤河集水文站峰量關(guān)系以Frank Copula函數(shù)擬合最優(yōu)，王家壩站以Gumbel Copula函數(shù)擬合最優(yōu)。Pickands依賴圖(圖9)顯示，息縣、淮濱、王家壩水文站的峰量聯(lián)合Copula函數(shù)擬合效果好，潤河集站擬合效果一般，這可能與潤河集站水文數(shù)據(jù)系列較短有關(guān)。

圖9 Copula函數(shù)擬合Pickands依賴圖

表10 峰量序列的Copula參數(shù)估計(jì)結(jié)果

3.4 洪水風(fēng)險(xiǎn)分析

基于Copula函數(shù)建立的峰量聯(lián)合分布，計(jì)算50年一遇的設(shè)計(jì)洪峰、設(shè)計(jì)洪量以及一致性條件和非一致性條件下峰量同現(xiàn)概率，結(jié)果見表11。

由表11可看出，一致性條件下，淮河干流發(fā)生50年一遇洪峰的同時(shí)遭遇大于50年一遇洪量的概率均小于0.02，從概率上可認(rèn)為流域內(nèi)不會(huì)同時(shí)遭遇50年一遇的洪峰及洪量。而在非一致性條件下，同現(xiàn)超過概率受降水影響，逐年變化。其中，王家壩站和潤河集站的同現(xiàn)超過概率均小于0.05；而息縣站和淮濱站的最大概率均超過0.1，且非一致性條件下概率大于一致性條件下概率的年份多為1954、1968、2007等大洪水年。以上統(tǒng)計(jì)結(jié)果表明，若在大洪水年采用一致性條件下的同現(xiàn)超過概率容易忽略峰高量大的洪水風(fēng)險(xiǎn)。

表11 50年一遇的峰量同現(xiàn)超過概率

4 討 論

淮河干流水文站水文極值基于一致性條件和非一致性條件的聯(lián)合頻率分析結(jié)果對比表明，以降水為協(xié)變量的水文極值序列設(shè)計(jì)值與實(shí)測系列變化趨勢一致，相同重現(xiàn)期下峰量同現(xiàn)超過概率在豐水年高于一致性情況，這與前人研究成果相符[20-22]。水文頻率分析是水文科學(xué)中的重要問題，而非一致性情況下多變量分析相對于一致性情況下單變量分析更為復(fù)雜。以時(shí)變邊際分布和Copula函數(shù)結(jié)合，對洪峰和15 d洪量進(jìn)行二元建模，計(jì)算非一致性條件下50年一遇的峰量同現(xiàn)超過概率，將洪水風(fēng)險(xiǎn)動(dòng)態(tài)化，結(jié)合未來氣候模式，可以更加準(zhǔn)確地預(yù)估未來洪水風(fēng)險(xiǎn)，為防洪工程規(guī)劃和決策提供參考。

基于GAMLSS 4種雙參數(shù)概率模型擬合各站洪峰和洪量序列時(shí)變邊際分布時(shí)，一致性情況下主要以Gamma分布擬合度最高，非一致性情況下主要以Weibull分布擬合度最高，說明在考慮非平穩(wěn)性情況下Weibull線型適用性更好。以與峰量序列相關(guān)性較強(qiáng)的降水為協(xié)變量的時(shí)變邊際分布模型，也反映了氣候變化下的水文響應(yīng)變化，具有一定的物理意義。時(shí)變邊際分布以兩參數(shù)線型為主，未考慮到三參數(shù)線型，如我國推薦使用的P-Ⅲ線型，后期工作需要進(jìn)一步擴(kuò)大研究范圍，將時(shí)變邊際分布與更多線型結(jié)合，進(jìn)行對比分析，得出最優(yōu)擬合線型；同時(shí)協(xié)變量需要引入水利工程因子、下墊面因子等，考慮人類活動(dòng)對水文機(jī)制的影響。

5 結(jié)論與展望

考慮到氣候變化對極端洪水事件的影響，經(jīng)過相關(guān)性分析，選擇降水為非一致性條件下洪水要素邊際分布模型的協(xié)變量，一定程度上反映了氣候變化對洪水要素頻率分布的非線性影響。

(1)以淮河干流水文站的年洪峰序列和年最大15 d洪量序列為例，運(yùn)用GAMLSS模型推求峰量序列的邊際分布，對比一致性條件和非一致性條件下的結(jié)果，可以發(fā)現(xiàn)，在大洪水年，基于以降水為協(xié)變量的非一致性洪峰和洪量邊際分布模型得到的設(shè)計(jì)值超出一致性條件下的設(shè)計(jì)值。

(2)運(yùn)用Copula函數(shù)構(gòu)建峰量聯(lián)合分布，分析50年一遇的峰量同現(xiàn)超過概率，結(jié)果表明，在非一致性條件下，同現(xiàn)超過概率逐年變化，且息縣站、淮濱站的最高概率均超過0.1，且非一致性條件下概率大于一致性條件下概率的年份多為1954、1968、1991、2007等大洪水年，說明在降水充沛的年份，洪水風(fēng)險(xiǎn)變高，峰量同現(xiàn)超過概率也相應(yīng)變高。