熊 威，胡鐵松，張兵堂，程曉峰

(武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072)

0 引 言

受全球氣候變化和劇烈人類活動影響，我國旱澇交替事件呈現(xiàn)廣發(fā)、頻發(fā)的態(tài)勢，特別是近年來頻發(fā)的旱澇急轉(zhuǎn)事件，據(jù)不完全統(tǒng)計，近百年來長江流域旱澇急轉(zhuǎn)發(fā)生頻率達(dá)到3～4年一遇，淮河流域達(dá)到4年一遇[1,2]。因此辨識旱澇交替事件基本特征，改變傳統(tǒng)基于單一干旱或者洪澇災(zāi)害治理理論體系與思路，對水旱災(zāi)害治理具有重要價值，是我國災(zāi)害治理關(guān)注的熱點與難點問題。

旱澇交替事件是一類多變量隨機事件，旱澇交替過程的持續(xù)時間、發(fā)生強度以及旱澇交替程度三者之間存在一定的相關(guān)關(guān)系，需要多特征量才能完整描述旱澇交替過程。但是目前旱澇交替事件主要是側(cè)重于旱澇交替過程單一特征的描述，如發(fā)生強度，而沒有從一個完整過程的多屬性特征角度進行完整描述。如吳志偉[3,4]提出了長、短周期旱澇急轉(zhuǎn)指數(shù)重要描述了水旱交替事件強度隨時間演變特征；王勝等[5]研究了淮河流域水旱交替事件強度空間分布特征。

Copula函數(shù)理論是描述這種相關(guān)性的一種有效方法，并得到了廣泛的應(yīng)用[6-8]。本文以四湖流域為例，采用游程理論從降雨距平指標(biāo)中分離出旱澇交替事件的3個特征變量，進而應(yīng)用Frank Copula函數(shù)構(gòu)建兩變量聯(lián)合概率分布以及三變量聯(lián)合分布，并計算相應(yīng)的重現(xiàn)期，以期得到四湖流域水旱交替事件的持續(xù)時間、發(fā)生強度以及旱澇交替程度三者之間的相關(guān)關(guān)系，為應(yīng)對地區(qū)旱澇交替災(zāi)害防治工作提供一定的依據(jù)。

1 方 法

1.1 旱澇交替事件以及特征變量

旱澇交替事件作為一種隨機事件，為了全面描述其本質(zhì)屬性，故定義了3個特征變量：持續(xù)時間(歷時)、發(fā)生強度以及旱澇交替程度，如圖1所示。旱澇交替歷時D為一次旱澇交替事件的總歷時，即L1、L2、L3三者之和；旱澇交替程度S為歷時D內(nèi)降水距平指數(shù)(以下簡稱Pa)絕對值之和，即S1和S2之和，表示一次旱澇交替事件產(chǎn)生的影響；峰值I為歷時D內(nèi)Pa的

圖1 旱澇交替特征變量定義圖Fig.1 Definition of drought and flood alternating variables

最大絕對值，反映了一次旱澇交替事件中某一天所產(chǎn)生的最大影響。

確定了特征變量以后，下一步是選取適合的分布函數(shù)進行擬合。本文選取對數(shù)正態(tài)分布、Gamma分布和廣義Pareto分布來分別擬合旱澇交替歷時D、程度S和峰值I[9-11]，這3種分布函數(shù)計算公式分別如下：

(2)

(3)

式中：μ為對數(shù)正態(tài)分布的方差；σ為標(biāo)準(zhǔn)差；α為Gamma分布的位置參數(shù)；β為尺度參數(shù)；k為廣義Pareto分布函數(shù)的形狀參數(shù)；γ為尺度參數(shù)；x0為閾值。

1.2 Copula函數(shù)理論

Copula函數(shù)式定義域為[0,1]的多維聯(lián)合分布，它可以將多個隨機變量的邊際分布連接起來構(gòu)造聯(lián)合分布。其中，Archimedean Copula函數(shù)結(jié)構(gòu)簡單，計算簡便，實際應(yīng)用中較多。不同的Copula函數(shù)有不同的適用范圍，本文選擇Frank Copula來構(gòu)造聯(lián)合分布。二維Frank Copula函數(shù)的表達(dá)式為：

(4)

相應(yīng)的三維Frank Copula函數(shù)的表達(dá)式為：

(5)

其中，θ為Copula函數(shù)的參數(shù)，對于二維Copula函數(shù)的參數(shù)估計，本文采用普遍適用且相對簡便的非參數(shù)法[12]，即利用τ與θ的關(guān)系式計算θ。

對于三維Copula函數(shù)，非參數(shù)法不再適用，本文采用參數(shù)法中的IFM(Inference of Functions for Margins)方法確定三維Copula函數(shù)的參數(shù)[13]。

選定的 Copula 函數(shù)能否描述變量之間的相關(guān)性結(jié)構(gòu)，需要對 Copula 函數(shù)進行擬合優(yōu)度檢驗，本文采用Kolomogorov-Smirnov(K-S)檢驗方法[14]。

1.3 重現(xiàn)期

工程中常用重現(xiàn)期來描述水文事件的嚴(yán)重性。假定X1、X2和X3分別表示旱澇交替事件中具有相依關(guān)系的特征變量，u1、u2和u3分別表示旱澇交替歷時、程度和峰值的邊緣分布， 則單變量重現(xiàn)期計算公式[15]如下：

(6)

式中：TXi為單變量重現(xiàn)期；N為系列長度，a；n為發(fā)生旱澇交替事件的次數(shù)，i=1,2,3。

對于二維聯(lián)合分布的重現(xiàn)期，以X1和X2的組合重現(xiàn)期為例，有“X1>x1或X2>x2”和“X1>x1且X2>x2”兩種情況，對應(yīng)的二維變量聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期計算公式[16]分別為：

(8)

三維變量聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期的計算公式[17]分別為：

(9)

(10)

2 特征變量聯(lián)合分布模型

本文以四湖流域監(jiān)利站為例，以降水距平指數(shù)[18]作為診斷量，從該站1957-2013年的逐日降雨資料中挑選出52次旱澇交替事件。借助于游程理論[19]從Pa序列中分離出旱澇事件，挑選標(biāo)準(zhǔn)如下：若前期干旱歷時L1(洪澇歷時L3)在10 d以上，并在10 d內(nèi)經(jīng)過天數(shù)L2由旱轉(zhuǎn)為澇(澇轉(zhuǎn)為旱)，且洪澇歷時L3(干旱歷時L1)在10 d以上，則選為一次旱澇交替事件。

2.1 單變量邊緣分布

分別用對數(shù)正態(tài)分布來擬合旱澇交替歷時D，Gamma分布來擬合旱澇交替程度S，廣義極值分布來擬合旱澇交替程度I。對數(shù)正態(tài)分布采用極大似然法估計參數(shù)，其余均采用矩法估計參數(shù)。參數(shù)以及擬合結(jié)果見表1，結(jié)果表明，3種分布均通過了α=0.05顯著性檢驗(n=52，Dmax<0.196)，即對旱澇交替歷時、程度以及峰值的擬合效果良好。

表1 邊緣分布函數(shù)參數(shù)值Tab.1 Estimated parameters for marginal distribution function

2.2 Copula聯(lián)合分布函數(shù)及擬合檢驗

2.2.1 二維模型

二維Frank Copula函數(shù)采用非參數(shù)法進行參數(shù)估計，參數(shù)計算結(jié)果和擬合效果見表2。結(jié)果表明，3種擬合情況下的Frank Copula函數(shù)均通過α=0.05顯著性檢驗(n=52，Dmax<0.196)，即擬合效果良好。

表2 二維Frank Copula函數(shù)參數(shù)和K-S檢驗結(jié)果Tab.2 Parameters of two-dimension Frank Copulafunction and the K-S fitting results

2.2.2 三維模型

根據(jù)IFM方法計算得到三維Frank Copula 函數(shù)的參數(shù)θ值為3.072，并采用K-S檢驗進行擬合檢驗，其檢驗值為0.113，表明三維Frank Copula函數(shù)達(dá)到了α=0.05顯著性水平(n=52，Dmax<0.196)，即擬合效果良好。

3 旱澇交替特征值的頻率分析

3.1 旱澇交替事件概況

監(jiān)利站從1957年到2013年共發(fā)生52次旱澇交替事件，約1.1 a一遇，平均歷時為33.85 d，平均程度為49.31，平均峰值為5.18。其中，旱轉(zhuǎn)澇事件占絕大多數(shù)，占總數(shù)的62%，其余為澇轉(zhuǎn)旱事件。

從年內(nèi)分布來看，如表3所示，旱澇交替事件主要集中在夏季，占總事件數(shù)的46%，冬季發(fā)生的次數(shù)最少；夏季旱澇交度替事件平均歷時為30.13 d，平均程度為58.21，平均峰值為7.01，平均程度和平均峰值均大于其他季節(jié)。從年間分布來看，如圖2所示，旱澇交替事件2000年以前年代見發(fā)生次數(shù)大致相同，2000以后有增加趨勢；20世紀(jì)80年代以前旱澇交替事件分布較為均勻，80年代以后分布較為集中；旱澇交替歷時和程度在2000年以后變化越來越劇烈。

表3 旱澇交替事件年內(nèi)分布情況Tab.3 The probability of drought and flood alternating events

圖2 旱澇交替特征值變化趨勢Fig.2 The tendency of drought and flood alternating variables

3.2 二維變量聯(lián)合概率分布

根據(jù)式(1)計算得到二維聯(lián)合分布函數(shù)，分別將旱澇交替D-S、D-I、S-I的聯(lián)合不超越概率以及聯(lián)合超越概率繪于圖3～圖5中，從二維平面投影圖上可以查出三個特征變量兩兩組合的聯(lián)合不超越概率和聯(lián)合超越概率。由圖3可知，常見的旱澇交替事件歷時小于33 d，程度小于56，此類事件占所有事件的50%左右；而歷時小于39 d，程度小于74的旱澇交替事件約占所有事件的80%；歷時60 d以上，程度大于91的旱澇交替事件極少發(fā)生。

由圖4可知，常見的旱澇交替事件程度小于50，峰值小于5，此類事件占所有事件的50%左右；而程度小于72，峰值小于9的旱澇交替事件約占所有事件的80%；程度大于92，峰值大于12的旱澇交替事件極少發(fā)生。

由圖5可知，常見的旱澇交替事件歷時小于32 d，峰值小于7，此類事件占所有事件的50%左右；而歷時小于39 d，峰值小于10的旱澇交替事件約占所有事件的80%；歷時60 d以上，程度大于91的旱澇交替事件極少發(fā)生。

圖3 D-S聯(lián)合分布概率三維圖Fig.3 The joint probability of D-S

圖4 S-I聯(lián)合分布概率三維圖Fig.4 The joint probability of S-I

圖5 D-I聯(lián)合分布概率三維圖Fig.5 The joint probability of D-I

3.3 組合重現(xiàn)期

在相同單變量重現(xiàn)期下，兩個和三個特征變量同時發(fā)生的組合重現(xiàn)期的計算結(jié)果見表4?？梢钥闯?，單變量重現(xiàn)期介于兩種多維變量的重現(xiàn)期之間，故可根據(jù)組合重現(xiàn)期來估計單變量重現(xiàn)期的范圍。綜合對比4種組合重現(xiàn)期，隨著單變量重現(xiàn)期的增加，D-I的組合重現(xiàn)期增幅最大，其次是D-S，S-I的增幅最小，說明S與I同時發(fā)生的風(fēng)險最高，與前面的分析一致；D-I的組合重現(xiàn)期范圍最大，例如單變量重現(xiàn)期為20 a時，同現(xiàn)重現(xiàn)期達(dá)到了500 a以上，這與D與I呈負(fù)相關(guān)有關(guān)。由3.1可知，考慮比較容易發(fā)生的特征變量即單變量發(fā)生頻率達(dá)到80%時，旱澇交替事件的D小于33.57 d，S小于61.13，I小于7.44，對應(yīng)的組合重現(xiàn)期見表6，可以看到，監(jiān)利站大部分旱澇交替事件的聯(lián)合重現(xiàn)期介于2～5 a之間，同現(xiàn)重現(xiàn)期介于48～100 a之間。

3.4 典型年頻率分析

三個特征變量值越大，產(chǎn)生的影響越大，故將特征變量較大的10次旱澇交替事件列在表7中，可以看出，2011年的3個特征變量值均排在前列，旱澇交替產(chǎn)生的影響嚴(yán)重，這與實際相符：2011年春末至夏初，長江中下游地區(qū)發(fā)生一次旱澇急轉(zhuǎn)現(xiàn)象，給人民生命、財產(chǎn)帶來了重大損失[20]。故以2011年為旱澇交替典型年，其特征變量為典型特征變量，分別記為D0、S0和I0。對典型年進行頻率分析，結(jié)果見表8?？梢钥闯觯舭l(fā)生比2011年特征變量值更大的旱澇交替事件，D-S同時超過事件的發(fā)生概率與I-S同時超過事件的發(fā)生概率接近，分別為0.31%和0.35%，D-I同時超過事件的發(fā)生概率三者中最小，這間接反映了三變量之間的相關(guān)性：相關(guān)性越強，二維聯(lián)合分布概率越大。三變量值同時超過事件的發(fā)生概率很小，接近于0。

表4 單變量重現(xiàn)期以及相應(yīng)二維變量的組合重現(xiàn)期Tab.4 Calculated results of the return periods of two-dimension joint distribution

注：Ta為聯(lián)合重現(xiàn)期，To為同現(xiàn)重現(xiàn)期。

表5 單變量重現(xiàn)期以及相應(yīng)三維變量的組合重現(xiàn)期Tab.5 Calculated results of the return periods ofthree-dimension joint distribution

表6 給定特征變量的多維變量組合重現(xiàn)期Tab.6 The return periods of multi-dimension jointdistribution with given variables

表7 監(jiān)利站旱澇交替事件的特征變量Tab.7 Variables of drought and flood alternatingevents at Jianli Meteorological stations

表8 典型年旱澇交替事件頻率Tab.8 Probability of the typical drought and flood alternating event

注：D0、S0和I0為2011年時時澇交替事件的特征變量；D0=60，S0=119.48，I0=14.13。

4 結(jié) 語

本文基于游程理論從監(jiān)利站1957-2007年的降雨距平序列中挑選出了旱澇交替事件，定義了旱澇交替特征變量，并應(yīng)用二維以及三維Frank Copula函數(shù)構(gòu)建了旱澇交替特征變量的聯(lián)合分布模型，基于該聯(lián)合分布模型，分析了旱澇交替事件的頻率，結(jié)果表明：

(1)監(jiān)利站共發(fā)生52次旱澇交替事件，約1.1 a一遇，平均歷時為33.85 d，平均程度為49.31，平均峰值為5.18。旱轉(zhuǎn)澇事件占絕大多數(shù)，從年內(nèi)分布上看，主要集中在夏季，冬季發(fā)生的次數(shù)最少；從年間分布上看，旱澇交替事件呈現(xiàn)出歷時變長，程度變大，分布越來越密集的態(tài)勢。

(2)旱澇交替事件作為隨機水文事件，其特征變量間具有相依關(guān)系，用Frank Copula函數(shù)構(gòu)造的兩變量聯(lián)合分布以及三變量聯(lián)合分布擬合精度是滿意的。通過聯(lián)合分布計算得到監(jiān)利站旱澇交替事件多為歷時小于39 d、程度小于74，或者程度小于72、峰值小于9，或者歷時小于39 d、程度小于74的旱澇交替事件，這3種情況均占所有事件的80%左右。

(3)單變量重現(xiàn)期介于兩種多維變量的重現(xiàn)期之間，可根據(jù)組合重現(xiàn)期來估計單變量重現(xiàn)期的范圍；監(jiān)利站大部分旱澇交替事件的聯(lián)合重現(xiàn)期介于2～5 a之間，同現(xiàn)重現(xiàn)期介于48～100 a之間。

(4)2011發(fā)生的年旱澇交替事件歷時達(dá)到60 d、程度達(dá)到119.48、峰值達(dá)到14.13，歷史上極為罕見，若發(fā)生比其特征變量值更大的旱澇交替事件，D-S同時超過事件的發(fā)生概率為0.31%，I-S同時超過事件的發(fā)生概率為0.35%，D-I同時超過事件的發(fā)生概率為0.20%，三變量值同時超過事件的發(fā)生概率接近于0。

需要指出的是，本文僅從氣象干旱出發(fā)，沒有考慮農(nóng)業(yè)、水文以及社會等方面的綜合影響，今后應(yīng)進一步深入研究。

□

[1] 程 智, 徐 敏, 羅連升,等. 淮河流域旱澇急轉(zhuǎn)氣候特征研究[J]. 水文, 2012,32(1):73-79.

[2] 曾凡榮. 湖北水利志[M]. 中國水利水電出版社, 2000.

[3] 吳志偉, 李建平, 何金海,等. 大尺度大氣環(huán)流異常與長江中下游夏季長周期旱澇急轉(zhuǎn)[J]. 科學(xué)通報, 2006,51(14):1 717-1 724.

[4] 吳志偉. 長江中下游夏季風(fēng)降水“旱澇并存、旱澇急轉(zhuǎn)”現(xiàn)象的研究[D]. 南京：南京信息工程大學(xué), 2006.

[5] 王 勝, 田 紅, 徐 敏,等. 1961-2008年淮河流域主汛期極端降水事件分析[J]. 氣象科技, 2012,40(1):87-91.

[6] 閆寶偉, 郭生練, 肖 義,等. 基于兩變量聯(lián)合分布的干旱特征分析[J]. 干旱區(qū)研究, 2007,24(4):537-542.

[7] 張 翔, 冉啟香, 夏 軍,等. 基于Copula函數(shù)的水量水質(zhì)聯(lián)合分布函數(shù)[J]. 水利學(xué)報, 2011,42(4):483-489.

[8] 王 珍, 李久生, 栗巖峰. 基于三維Copula函數(shù)的滴灌硝態(tài)氮淋失風(fēng)險評估方法[J]. 農(nóng)業(yè)工程學(xué)報, 2013,29(19):79-87.

[9] Mayer P, Mayer H. Radar target detection-handbook of theory and practise[M].New York:Academic Press,1973:1-108.

[10] Emir Z, Atila S. A method of streamflow drought analysis[J]. Water Resources Research, 2010,23(23):156-168.

[11] Mannshardt-Shamseldin E C, Smith R L, Sain S R, et al. Downscaling extremes: a comparison of extreme value distributions in point-source and gridded precipitation data[J]. Annals of Applied Statistics, 2010,4(1):484-502.

[12] Genest C, Rivest L P. Statistical inference procedures for bivariate archimedean copulas[J]. Journal of the American Statistical Association, 2012,88(423):1 034-1 043.

[13] 楊益黨, 羅羨華. Copula函數(shù)的參數(shù)估計[J]. 新疆師范大學(xué)學(xué)報:自然科學(xué)版, 2007,26(2):15-18.

[14] Singh V P, Zhang L. Bivariate flood frequency analysis using the copula method[J]. Journal of Hydrologic Engineering, 2006,11(2):150-164.

[15] Shiau J T, Shen H W. Recurrence analysis of hydrologic droughts of differing severity[J]. Journal of Water Resources Planning & Management, 2014,127(1):30-40.

[16] Shiau J T. Return period of bivariate distributed hydrological events[J]. Stochastic Environmental Research & Risk Assessment, 2003,17:42-57.

[17] 李 計, 李 毅, 賀纏生. 基于Copula函數(shù)的黑河流域干旱頻率分析[J]. 西北農(nóng)林科技大學(xué)學(xué)報:自然科學(xué)版, 2013,(1):213-220.

[18] 鞠笑生, 楊賢為, 陳麗娟,等. 我國單站旱澇指標(biāo)確定和區(qū)域旱澇級別劃分的研究[J]. 應(yīng)用氣象學(xué)報, 1997,(1):26-33.

[19] Herbst P H, Bredenkamp D B, Barker H M G. A technique for the evaluation of drought from rainfall data[J]. Journal of Hydrology, 1966,4(66):264-272.

[20] 封國林, 楊涵洧, 張世軒,等. 2011年春末夏初長江中下游地區(qū)旱澇急轉(zhuǎn)成因初探[J]. 大氣科學(xué), 2012,36(5):1 009-1 026.