陳子燊，劉占明，黃 強

(中山大學地理學院水資源與環(huán)境系，廣州510275)

進入21世紀以來，全球氣候變化背景下的區(qū)域干旱現(xiàn)象非常引人注目.由于干旱的發(fā)生與發(fā)展過程具有漸進性和持續(xù)性、影響范圍廣等特征，如形成災害將對區(qū)域經(jīng)濟社會發(fā)展造成重大影響.從概率的角度分析并預測干旱發(fā)生的規(guī)律，可為區(qū)域水資源規(guī)劃、風險管理、防旱抗旱與制定區(qū)域社會與經(jīng)濟發(fā)展規(guī)劃或計劃等提供決策依據(jù).干旱類型通常可分為氣象干旱或氣候干旱、農(nóng)業(yè)干旱、水文干旱及社會經(jīng)濟干旱4類［1］.其中，水文干旱被認為是氣象干旱和農(nóng)業(yè)干旱的延續(xù)和發(fā)展，是最終最徹底的干旱［2］.

目前，國內(nèi)外應用游程理論和Copula理論與方法在量化干旱特征與變化規(guī)律的研究方面做了不少研究工作.其中，F(xiàn)leig等［3］采用3種方法:滑動平均方法(MA-procedure)、事件間隔方法(IT-method)和連續(xù)峰值算法(SPA)對徑流干旱特征值的提取做了綜合評價.周玉良等［4］以徑流量為水文干旱指標，通過徑流量距平百分率和徑流量累積頻率兩種方法識別水文干旱特征變量，并以適線法確定單個干旱特征變量的分布曲線，在此基礎上利用Copula函數(shù)構(gòu)建了干旱歷時與干旱烈度間的聯(lián)合分布.陸桂華等［5-10］分別利用游程理論提取降水資料或標準降水指數(shù)(SPI)的干旱特征值，應用Copula方法計算干旱歷時和干旱烈度的聯(lián)合分布概率，推算干旱極限事件重現(xiàn)水平.尹正杰等［11-14］分別利用截斷水平提取水文干旱事件，應用馬爾可夫平穩(wěn)概率和一階自回歸模式模擬隨機合成的徑流，研究干旱強度、干旱烈度、干旱歷時和強度及風險概率等統(tǒng)計特征.不過，由于干旱烈度至今尚沒有一個明確的指標劃分標準，而干旱強度指標是建立在概率分布基礎并可直接聯(lián)系旱情級別，建立干旱歷時和干旱強度之間的聯(lián)合分布將可直接分析干旱風險及給予概率意義上的預測.

1 水文干旱頻率分析方法

1.1 水文干旱指標

對水文干旱指標劃分的研究表明［15］，Seiler等［16］提出的標準化降水指標(SPI)優(yōu)于在我國廣泛應用的Z指數(shù).Nalbantis等［17］采用了與SPI相同的方法計算水文年多時間尺度的徑流干旱指數(shù)(Streamflow Drought Index，SDI).首先采用伽馬分布描述不同時間尺度的月徑流量的偏態(tài)分布，進而經(jīng)標準正態(tài)化得出SDI.由此計算獲得的SDI不僅可量化描述多時間尺度的水文干旱特征，而且適用于不同時間與空間尺度的干旱監(jiān)測和評價.具體計算方法如下:

假設某一時間尺度下的徑流量為x，則其滿足伽馬分布的概率密度函數(shù)f(x)為:

式中，α、β分別為形狀參數(shù)和尺度參數(shù).給定時間尺度的徑流量x的累積概率為

對Γ分布概率進行正態(tài)標準化處理求解得到SDI:

1.2 水文干旱特征值與重現(xiàn)期定義

設定SDI=－1為游程的截斷水平提取水文干旱特征值:干旱歷時D、干旱強度I和干旱烈度S(圖1).應用Copula函數(shù)構(gòu)建干旱歷時和干旱強度兩變量的聯(lián)合分布模式，推算其相應重現(xiàn)水平值，進而計算條件概率.

邊緣分布干旱歷時和干旱強度的重現(xiàn)期TD、TI公式分別為［10］:

式中，E(L)為干旱時間間隔的期望值.

聯(lián)合重現(xiàn)期To:

同現(xiàn)重現(xiàn)期Ta:

條件概率分布:

式中，P(D≥d|I≥i)為超值累積條件概率分布，可定義為遭遇概率或風險概率，其倒數(shù)即為條件(或遭遇)重現(xiàn)期.同理可計算P(I≥i|D≥d)的遭遇概率和條件重現(xiàn)期.

2 實例研究

2.1 研究背景與SDI指數(shù)統(tǒng)計分析

圖1 干旱事件特征描述Fig.1 Characterization of drought events

干旱歷來是廣東三大心腹大患(洪澇、干旱、風災)之一.近十幾年來全球氣候變化異常，降水變率增大，時空分布愈加不均，干旱趨于增強.其中，1998下半年至1999年廣東曾出現(xiàn)夏秋冬春連旱現(xiàn)象，其引發(fā)的重大災害成為社會普遍關注的問題.

珠江流域西江水系源于云南沾益縣馬雄山，東流至廣西梧州匯桂江后始稱西江，流入廣東經(jīng)珠江三角洲的磨刀門、雞啼門、虎跳門和崖門注入南海，平均年徑流量為2300×108m3.馬口水文站位于西、北江三角洲河網(wǎng)區(qū)的頂端，作為國家重點水文監(jiān)測站是西江進入河網(wǎng)區(qū)前必經(jīng)的第一關隘，掌握和預測馬口水文站的水情是做好西北江三角洲防洪防旱工作的關鍵.1990s以來，珠江河口三角洲枯季咸潮活動呈現(xiàn)出越來越頻繁、強度不斷增大、影響范圍越來越大和持續(xù)時間不斷增加等新特點［18］.受此影響，位于西江主要入海通道的磨刀門水道沿岸的中山、珠海和澳門等城市枯季用水受到極大影響.

采用1959年1月到2009年12月馬口水文站月徑流量計算了徑流干旱指數(shù).對SDI的三個時間尺度的統(tǒng)計分析結(jié)果表明，1960s是最近數(shù)十年間西江發(fā)生旱情最嚴重的時段，1970s-1980s旱情緩解，1990s以來干旱出現(xiàn)頻次又持續(xù)增大(表1).進一步從風險概率角度深入分析馬口站的水文干旱特征有助于有關部門對西江三角洲枯季水資源進行規(guī)劃管理，并可加深對咸潮發(fā)生與持續(xù)性的認識.

表1 不同時間尺度的水文干旱頻次的年代際變化Tab.1 Interdecadal variation of different time scales of hydrological drought occurrence frequencies

2.2 聯(lián)合分布計算與分析

對D和I樣本都使用了三參數(shù)的Weibull概率分布函數(shù)［19］:

式中，ξ、β、μ分別為形態(tài)參數(shù)、尺度參數(shù)和位置參數(shù).參數(shù)估計使用線性矩(L-矩)方法［20］.經(jīng)驗頻率分布Pi使用 Gringorten 公式計算［21］:Pi=(i－0.44)/(n+0.12).擬合結(jié)果采用均方根誤差(RMSE)、經(jīng)驗頻率和理論頻率擬合誤差平方和(Q)和概率點據(jù)相關系數(shù)(PPCC)檢驗擬合優(yōu)度.擬合的干旱歷時和干旱強度的PPCC分別達0.994和0.992，其Weibull概率分布見圖2.

圖2 干旱歷時(a)和干旱強度(b)Weibull概率分布Fig.2 Weibull distribution of drought duration(a)and drought intensity(b)

統(tǒng)計的干旱歷時和干旱強度之間的秩相關系數(shù)為0.617，說明二者之間存在高關聯(lián)性.選擇Archimedean 族4 種 Copula函數(shù):Gumbel-Hougaard Copula、Clayton Copula、Ali-Mikhail-Haq Copula和 Frank Copula，通過Genest-Rivest圖示方法［22］、均方根誤差(RMSE)準則法和AIC信息準則法等擬合優(yōu)度檢驗方法擇優(yōu)構(gòu)建干旱變量D和I的聯(lián)合概率分布模型，進而推算其重現(xiàn)水平和風險概率.

由相關性指標法計算的兩個變量的聯(lián)合概率分布Copula函數(shù)的參數(shù)θ等于2.611，利用AIC、RMSE檢驗其擬合優(yōu)度的結(jié)果見表2.

表2 4個Copula函數(shù)的參數(shù)及其擬合優(yōu)度指標對比Tab.2 Parameters of four Copula functions and evaluation indices of goodness of fit

根據(jù)擬合優(yōu)度評價指標，選擇兩樣本中AIC和RMSE最小、Genest-Rivest圖示方法Kt-Ke關系圖中點據(jù)和理論直線最接近45°對角線的Gumbel-Hougaard Copula作為聯(lián)合概率分布的連接函數(shù)(圖3).由于Gumbel-Hougaard Copula函數(shù)具有描述變量之間上尾正相關的性質(zhì)，適用于干旱事件的概率分析.擇優(yōu)構(gòu)建的干旱歷時和干旱強度聯(lián)合概率分布模式如下:

2.3 干旱歷時和干旱強度聯(lián)合分布的重現(xiàn)水平

由二變量Gumbel-Hougaard Copula聯(lián)合分布計算的干旱歷時和干旱強度的不同重現(xiàn)期可見，重現(xiàn)期為5～10年的干旱歷時約2～3個月，相應的干旱強度介于1.5～2.0之間，即可達到重旱級別(表3).重現(xiàn)期為20年情況下，干旱歷時約4個月，干旱強度即可達到特旱級別.50年、100年一遇的水文干旱歷時分別約為5個月、6個月.

表3 干旱歷時和干旱強度的不同重現(xiàn)水平值Tab.3 Different return level values of drought duration and drought intensity

由聯(lián)合分布推算的干旱歷時和干旱強度重現(xiàn)水平值分別略小于其相應的邊緣分布的重現(xiàn)水平值，同現(xiàn)分布條件下推算的干旱歷時和干旱強度重現(xiàn)水平值分別略大于其相應的邊緣分布的重現(xiàn)水平值，聯(lián)合分布和同現(xiàn)分布推算的重現(xiàn)水平值可視為單變量重現(xiàn)水平值的上下限，這樣有助于對水資源工程建設投入與風險管理進行綜合考慮.

圖3 4種Copula函數(shù)擬合對比Fig.3 Comparison of fitness of four theoretical and empirical Copulas

2.4 干旱歷時和干旱強度條件分布與遭遇概率

干旱歷時和干旱強度P(I≥i|D≥d)和P(D≥d|I≥i)兩個條件概率分布見圖4，干旱歷時D為2、3、4、5、6、7 月和干旱強度 I為 1.63、1.91、2.19、2.55、2.82、3.09(分別對應重現(xiàn)期 5、10、20、50、100、200 年)之間的組合遭遇概率計算結(jié)果見表4.

圖4 干旱歷時(a)和干旱強度(b)條件概率等值線分布Fig.4 Isograms distribution of conditional probabilities of drought duration(a)and drought intensity(b)

干旱歷時和干旱強度遭遇概率表明:1)對特定干旱強度或干旱歷時，隨干旱歷時的加長或干旱強度的加大，二者的遭遇概率明顯增大，表明出現(xiàn)干旱歷時或干旱級別的風險迅速增大.例如，20年一遇的干旱強度值為 2.2，其與干旱歷時為 2、3、…、7 個月的遭遇概率分別為 23.5%、45.6%、…、92.0%，條件重現(xiàn)期介于1～4年之間;2)二者相關性強:同時增大干旱強度和干旱歷時，二者大致仍保持在66% ～71%之間的高遭遇概率水平;3)隨著干旱歷時或干旱強度的加大，其與特定干旱特征值的遭遇概率隨之遞減.例如，干旱歷時從2個月增大至7個月，其與干旱強度1.63的遭遇概率從66.6%降至2.1%，條件重現(xiàn)期大約介于1.5～48.2年之間，這進一步佐證了高的干旱級別通常與長的干旱歷時相關聯(lián).綜上分析表明，對應于特定干旱歷時出現(xiàn)干旱級別(中旱、重旱、特旱)或特定干旱強度下出現(xiàn)干旱歷時的條件概率分布實際上給出了干旱特征概率意義上的診斷預測，此對于枯季西北江三角洲的水資源的調(diào)度管理與調(diào)淡水壓咸潮提供了有實質(zhì)意義的信息.

表4 干旱歷時和干旱強度遭遇概率Tab.4 Encounter probabilities of drought dura tion and drought intensity

3 結(jié)論

本文對西江下游馬口水文站61年的徑流干旱指數(shù)SDI提取干旱特征值，進而基于二變量Copula函數(shù)構(gòu)建水文干旱特征值——干旱歷時和干旱強度的聯(lián)合分布模式，獲得以下主要結(jié)論:

1)擇優(yōu)檢驗表明，采用Archimedean類的Gumbel-Hougaard Copula函數(shù)較好地構(gòu)建了干旱歷時和干旱強度的聯(lián)合分布模式.

2)5～10年一遇的水文干旱即屬重旱級別，干旱持續(xù)約2～3個月.20年一遇即達特旱級別，干旱持續(xù)約4個月.

3)聯(lián)合分布和同現(xiàn)分布推算的重現(xiàn)水平值可作為干旱歷時或干旱強度重現(xiàn)水平值的上下限，此有助于對枯季水資源工程建設與風險管理的綜合考慮.

4)干旱歷時和干旱強度條件概率分布表明，二者具有很強的相關性，其為特定干旱歷時與水文干旱級別(中旱、重旱、特旱)或特定干旱強度與干旱歷時之間的對應關系提供了概率意義上的干旱特征診斷與預測.

［1］Richard R，Heim Jr著.周躍武，馮建英譯.美國20世紀干旱指數(shù)評述.干旱氣象，2006，24(1):79-89.

［2］耿鴻江，沈必成.水文干旱的定義及其意義.干旱地區(qū)農(nóng)業(yè)研究，1992，10(4):91-94.

［3］Fleig AK，Tallaksen LM，Hisdal H et al.A global evaluation of stream flow drought characteristics.Hydrology and Earth System Sciences，2006，10:535-552.

［4］周玉良，袁瀟晨，金菊良等.基于Copula的區(qū)域水文干旱頻率分析.地理科學，2011，31(11):1383-1388.

［5］陸桂華，閆桂霞，吳志勇等.基于Copula函數(shù)的區(qū)域干旱分析方法.水科學進展，2010，21(2):188-193.

［6］許月萍，張慶慶，樓章華等.基于Copula方法的干旱歷時和烈度的聯(lián)合概率分析.天津大學學報，2010，43(10):928-932.

［7］馬明衛(wèi)，宋松柏.橢圓型Copulas函數(shù)在西安站干旱特征分析中的應用.水文，2010，30(4):36-42.

［8］張 雨，宋松柏.Copulas函數(shù)在多變量干旱聯(lián)合分布中的應用.灌溉排水學報，2010，29(3):64-68.

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

［10］Shiau J.Fitting drought duration and severity with Two-Dimensional Copulas.Water Resources Management，2006，20(5):795-815.

［11］尹正杰，黃 薇，陳 進.水庫徑流調(diào)節(jié)對水文干旱的影響分析.水文，2009，29(2):41-44.

［12］袁 超，宋松柏，荊 萍.極限水文干旱歷時概率分布解析法研究.西北農(nóng)林科技大學學報，2008，36(7):212-218.

［13］Salas JD，Chongjin F，Cancelliere A et al.Characterizing the severity and risk of drought in the Poudre River.Journal of Water Resources Planning and Management，2005，131(5):383-393.

［14］蔡明科.關中地區(qū)水文、氣象干旱特征對比研究.灌溉排水學報，2003，22(6):33-37.

［15］袁文平，周廣勝.標準化降水指標與Z指數(shù)在我國應用的對比分析.植物生態(tài)學報，2004，28(4):523-529.

［16］Seiler BA，Hayes M，Bressan L.Using the standardized precipitation index for flood risk monitoring.International Journal of Climatology，2002，22:1365-1376.

［17］Nalbantis I，Tsakiris G.Assessment of hydrological drought revisited.Water Resources Management，2009，23(5):881-897.

［18］胥加仕，羅承平.近年來珠江三角洲咸潮活動特點及重點研究領域探討.人民珠江，2005，2:21-23.

［19］Hosking JRM.L-moments:Analysis and estimation of distributions using linear combinations of order statistics.Royal Statistical Society，1990，52(1):105-124.

［20］Seguro JV，Lambert TW.Modern estimation of the parameters of the Weibull wind speed distribution for wind energy analysis.Journal of Wind Engineering and Industrial Aerodynamics，2000，85(1):75-84.

［21］郭生練，葉寧澤.論水文計算中的經(jīng)驗頻率公式.武漢水利電力學院學報，1992，25(2):38-45.

［22］Genest C，Rivest L.Statistical inference procedures for bivariate Archimedean copulas.Journal of American Statistical Association，1993，88:1034-1043.