陳子燊，劉占明，黃 強(qiáng)

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

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

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

1 水文干旱頻率分析方法

1.1 水文干旱指標(biāo)

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

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

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

對(duì)Γ分布概率進(jìn)行正態(tài)標(biāo)準(zhǔn)化處理求解得到SDI:

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

設(shè)定SDI=－1為游程的截?cái)嗨教崛∷母珊堤卣髦?干旱歷時(shí)D、干旱強(qiáng)度I和干旱烈度S(圖1).應(yīng)用Copula函數(shù)構(gòu)建干旱歷時(shí)和干旱強(qiáng)度兩變量的聯(lián)合分布模式，推算其相應(yīng)重現(xiàn)水平值，進(jìn)而計(jì)算條件概率.

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

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

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

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

條件概率分布:

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

2 實(shí)例研究

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4 干旱歷時(shí)和干旱強(qiáng)度條件分布與遭遇概率

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

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

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

表4 干旱歷時(shí)和干旱強(qiáng)度遭遇概率Tab.4 Encounter probabilities of drought dura tion and drought intensity

3 結(jié)論

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

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

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

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

4)干旱歷時(shí)和干旱強(qiáng)度條件概率分布表明，二者具有很強(qiáng)的相關(guān)性，其為特定干旱歷時(shí)與水文干旱級(jí)別(中旱、重旱、特旱)或特定干旱強(qiáng)度與干旱歷時(shí)之間的對(duì)應(yīng)關(guān)系提供了概率意義上的干旱特征診斷與預(yù)測(cè).

［1］Richard R，Heim Jr著.周躍武，馮建英譯.美國(guó)20世紀(jì)干旱指數(shù)評(píng)述.干旱氣象，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ū)域水文干旱頻率分析.地理科學(xué)，2011，31(11):1383-1388.

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

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

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

［8］張 雨，宋松柏.Copulas函數(shù)在多變量干旱聯(lián)合分布中的應(yīng)用.灌溉排水學(xué)報(bào)，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］尹正杰，黃 薇，陳 進(jìn).水庫(kù)徑流調(diào)節(jié)對(duì)水文干旱的影響分析.水文，2009，29(2):41-44.

［12］袁 超，宋松柏，荊 萍.極限水文干旱歷時(shí)概率分布解析法研究.西北農(nóng)林科技大學(xué)學(xué)報(bào)，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］蔡明科.關(guān)中地區(qū)水文、氣象干旱特征對(duì)比研究.灌溉排水學(xué)報(bào)，2003，22(6):33-37.

［15］袁文平，周廣勝.標(biāo)準(zhǔn)化降水指標(biāo)與Z指數(shù)在我國(guó)應(yīng)用的對(duì)比分析.植物生態(tài)學(xué)報(bào)，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］胥加仕，羅承平.近年來(lái)珠江三角洲咸潮活動(dòng)特點(diǎn)及重點(diǎn)研究領(lǐng)域探討.人民珠江，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ì)算中的經(jīng)驗(yàn)頻率公式.武漢水利電力學(xué)院學(xué)報(bào)，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.